Skip to content

GitLab

Commit 9155402a authored by LAMBERT Jean-charles's avatar LAMBERT Jean-charles
Browse files
Options

more stuffs for analysis

parent ea8e7b6c
Expand all Hide whitespace changes
Inline Side-by-side
Showing with 307 additions and 176 deletions
lib/projects/nemodep/c2dplot.cc
View file @ 9155402a
......@@ -165,7 +165,7 @@ template <class T> void C2dplot<T>::computeImage(const int xaxis, const int yaxi
// look for indexes in the range axis
findIndexes(xaxis,yaxis);
std::cerr << "XYZ particles => Nb indexes found = " << indexes.size() << "\n";
//std::cerr << "XYZ particles => Nb indexes found = " << indexes.size() << "\n";
xmin=ymin=min(range[xaxis][0],range[yaxis][0]);
xmax=ymax=max(range[xaxis][1],range[yaxis][1]);
......@@ -177,8 +177,8 @@ template <class T> void C2dplot<T>::computeImage(const int xaxis, const int yaxi
// launch computation on threads
startWorkers(nbody,pos,xaxis,yaxis,zmin,zmax);
#endif
std::cerr << "Work done cpu time : "<< timing.cpu() << "\n";
std::cerr << "Work done elapsed : "<< timing.elapsed() << "\n";
//std::cerr << "Work done cpu time : "<< timing.cpu() << "\n";
//std::cerr << "Work done elapsed : "<< timing.elapsed() << "\n";
// transfert function
float tr[6] = { xmin, (xmax-xmin)/(float)dimx,
......@@ -359,7 +359,7 @@ template <class T> void C2dplot<T>::startWorkers(const int nbody, T * data, cons
}
}
std::cerr << "[-->] zmin="<<zmin<<" zmax="<<zmax<<" nbody="<<nbody<< "\n";
//std::cerr << "[-->] zmin="<<zmin<<" zmax="<<zmax<<" nbody="<<nbody<< "\n";
}
template void C2dplot<float >::startWorkers(const int nbody, float * data, const int xaxis, const int yaxis, float& zmin, float& zmax);
template void C2dplot<double>::startWorkers(const int nbody, double * data, const int xaxis, const int yaxis, float& zmin, float& zmax);
......@@ -393,8 +393,9 @@ template <class T> void C2dplot<T>::worker(const int ithread, const int offset,
float hsml_size;
if (hsml) {
hsml_size=ceil(((hsml[ii])/(xmax-xmin))*(dimx-1));
if (hsml_size>700)
std::cerr << "hsml_size="<<hsml_size << " hsml="<< hsml[ii] << "\n";
if (hsml_size>700) {
//std::cerr << "hsml_size="<<hsml_size << " hsml="<< hsml[ii] << "\n";
}
}
assert(x<dimx);
......@@ -428,7 +429,7 @@ template <class T> void C2dplot<T>::worker(const int ithread, const int offset,
zmax = std::max(zmax,tab[ithread][x*dimx+y]);
}
std::cerr << "Before gaussian Zmin = "<<zmin<< " zmax = "<<zmax << "\n";
//std::cerr << "Before gaussian Zmin = "<<zmin<< " zmax = "<<zmax << "\n";
// loop to store cells from the array which have a weight
std::map<int, int> HSML;
pvec.clear();
......@@ -456,7 +457,7 @@ template <class T> void C2dplot<T>::worker(const int ithread, const int offset,
// sort vector of particles
//std::sort(pvec.begin(),pvec.end(),CPartProp::mySort);
std::cerr << "HSML size =" << HSML.size() << "\n";
//std::cerr << "HSML size =" << HSML.size() << "\n";
int cpt=0;
for (std::map<int, int>::iterator ii=HSML.begin(); ii !=HSML.end(); ii++) {
......
lib/projects/nemodep/crectify.cc
View file @ 9155402a
......@@ -38,6 +38,7 @@ CRectify::CRectify(bool _verbose)
{
verbose=_verbose;
only_eigen_values=false;
w_sum_ok=false;
init();
}
......@@ -62,7 +63,6 @@ bool CRectify::rectify(const int _nbody,const float _time,
const bool _use_rho, const bool _rectify,
std::string _cod_file,std::string _rect_file, const float _radius,const float _dmin, const float _dmax)
{
bool ok = true;
nbody= _nbody;
pos = _pos;
vel = _vel;
......@@ -81,7 +81,7 @@ bool CRectify::rectify(const int _nbody,const float _time,
process();
return ok;
return w_sum_ok;
}
// PUBLIC CDataIndex
......@@ -266,7 +266,11 @@ void CRectify::findMoment()
}
}
}
DIVMS(w_qpole, w_qpole, w_sum);
if (w_sum>0) {
DIVMS(w_qpole, w_qpole, w_sum);
w_sum_ok=true;
}
}
//------------------------------------------------------------------------------
// computeVectors
......
lib/projects/nemodep/crectify.h
View file @ 9155402a
......@@ -106,10 +106,18 @@ public:
for (int i=0; i<6; i++) {
_rect_array[ii++] = codf[i];
}
for (int i=0; i<3; i++) {
_rect_array[ii++] = frame[i][0];
_rect_array[ii++] = frame[i][1];
_rect_array[ii++] = frame[i][2];
if (ok) {
for (int i=0; i<3; i++) {
_rect_array[ii++] = frame[i][0];
_rect_array[ii++] = frame[i][1];
_rect_array[ii++] = frame[i][2];
}
} else {
for (int i=0; i<3; i++) {
_rect_array[ii++] = 0.0;
_rect_array[ii++] = 0.0;
_rect_array[ii++] = 0.0;
}
}
return ok;
......@@ -140,6 +148,7 @@ private:
std::string cod_file, rect_file;
bool rect, use_rho;
int nbody;
bool w_sum_ok; // w_sum is positif
std::vector<CDataIndex> rho_di;
void init();
void processRho();
......
py/mains/mdf_mt.py
View file @ 9155402a
......@@ -21,6 +21,11 @@ def commandLine():
# parse
args = parser.parse_args()
process(args)
# -----------------------------------------------------
# process, is the core function
def process(args):
data=np.genfromtxt(args.infile,dtype=None,names=['strvar','fltvar','fltvar','fltvar'])
found=False
for i in data:
......@@ -35,6 +40,8 @@ def commandLine():
print ("none", 0., 0., 0.)
#embed()
# -----------------------------------------------------
# main program
commandLine() # parse command line
if __name__ == '__main__':
commandLine()
py/mains/merging_time.py
View file @ 9155402a
#!/usr/bin/python
from __future__ import print_function
import sys
#sys.path.append('/home/jcl/works/GIT/uns_projects/py/modules/')
#sys.path.append('/home/jcl/works/GIT/uns_projects/py/modules/simulations')
#sys.path=['/home/jcl/works/GIT/uns_projects/py/modules/','/home/jcl/works/GIT/uns_projects/py/modules/simulations']+sys.path
from uns_simu import *
from simulations.ccod import *
import argparse
......@@ -18,24 +17,28 @@ def commandLine():
cod1=None
cod2=None
dmax=1.0
seger=False
threshold=10000
# help
# help
parser = argparse.ArgumentParser(description="Display merging time",
formatter_class=argparse.ArgumentDefaultsHelpFormatter)
# options
parser.add_argument('simname', help='Simulation name')
parser.add_argument('--dmax', help='max distance of separation (kpc)',default=dmax,type=float)
parser.add_argument('--cod1', help='first cod file',default=cod1)
parser.add_argument('--cod2', help='second cod file',default=cod2)
parser.add_argument('--seger', help='use sergey files',default=seger)
parser.add_argument('--seger', help='use sergey files',dest="seger",action='store_true',default=False)
#parser.add_argument('--no-seger', help='use base simulation files',dest="seger",action='store_false',default=True)
parser.add_argument('--pngfile', help='png filename, if None interactive plot',default=pngfile)
parser.add_argument('--txtfile', help='merging time filename, if None not saved',default=txtfile)
parser.add_argument('--plot', help='Enable ploting',dest='plot',action='store_true',default=False)
#parser.add_argument('--no-plot', help='Disable ploting',dest='plot',action='store_false',default=True)
parser.add_argument('--dbname',help='UNS database file name', default=dbname)
parser.add_argument('--verbose',help='verbose mode', default=False)
parser.add_argument('--verbose',help='verbose mode on', dest="verbose", action="store_true", default=False)
#parser.add_argument('--no-verbose',help='verbose mode off', dest="verbose", action="store_false", default=True)
# parse
args = parser.parse_args()
# parse
args = parser.parse_args()
if args.pngfile=="None":
args.pngfile=None
......@@ -48,7 +51,10 @@ def commandLine():
def process(args):
cod = CCod(simname=args.simname,verbose_debug=args.verbose,dbname=args.dbname)
try:
cod.computeMergingTime(halo_1=args.cod1,halo_2=args.cod2,pngfile=args.pngfile,txtfile=None,dmax=args.dmax,seger=args.seger)
mt=cod.computeMergingTime(halo_1=args.cod1,halo_2=args.cod2,pngfile=args.pngfile,txtfile=args.txtfile,dmax=args.dmax,seger=args.seger,plot=args.plot)
if mt<0.0:
mt=0.0
print("%f\n"%(mt))
except KeyboardInterrupt:
sys.exit()
......
py/mains/script_analysis.txt
View file @ 9155402a
......@@ -7,9 +7,9 @@ from simulations.ccom import *
from simulations.crectify import *
from simulations.cmovie import *
run_cod=0
run_com=0
run_rectify=0
run_cod=1
run_com=1
run_rectify=1
run_movie=1
##
## COD
......
py/modules/simulations/c2dplot.py
View file @ 9155402a
......@@ -72,7 +72,6 @@ class C2dplot:
# try to load positions
ok,pos=uns_snap.getData(select,"pos")
print("pos=",pos)
if ok:
# load time
ok,time=uns_snap.getData("time")
......@@ -87,7 +86,6 @@ class C2dplot:
nrange[:,0]=-float(rrange)
nrange[:,1]= float(rrange)
nrange=np.reshape(nrange,nrange.size) #reshape in 1D
print (nrange,nrange.size)
legend=""
comp_prop=select
......@@ -136,9 +134,9 @@ class C2dplot:
hsml=np.array([],dtype='f')
print(pos.size,outdev,no,weight)
print(outdev,no,pos,nrange,title,comp_prop,filename,time.item(),xy,xz,zy,\
True,weight,psort,hsml,itf,cb,legend,cmap)
#print(pos.size,outdev,no,weight)
#print(outdev,no,pos,nrange,title,comp_prop,filename,time.item(),xy,xz,zy,\
# True,weight,psort,hsml,itf,cb,legend,cmap)
self.__c.compute_image(outdev,no,pos,nrange,title,comp_prop,filename,time.item(),xy,xz,zy,\
True,weight,psort,hsml,itf,cb,legend,cmap)
......
py/modules/simulations/ccod.py
View file @ 9155402a
This diff is collapsed.
py/modules/simulations/ccom.py
View file @ 9155402a
......@@ -84,17 +84,17 @@ class CCom:
print("Unable to create directory [%s]\n"%(self.__analysis.com_dir))
self.__analysis.lock[self.__analysis.lock_id].release()
sys.exit()
else:
#
if not os.path.isfile(self.__analysis.com_file):
f=open(self.__analysis.com_file,"w")
try:
f.write("%s\n"%(self.__analysis.sim_info['name']))
f.close()
except:
print("Unable to create file [%s], aborting"%(self.__analysis.com_file))
self.__analysis.lock[self.__analysis.lock_id].release()
sys.exit()
#
if not os.path.isfile(self.__analysis.com_file):
f=open(self.__analysis.com_file,"w")
try:
f.write("%s\n"%(self.__analysis.sim_info['name']))
f.close()
except:
print("Unable to create file [%s], aborting"%(self.__analysis.com_file))
self.__analysis.lock[self.__analysis.lock_id].release()
sys.exit()
# release process
self.__analysis.lock[self.__analysis.lock_id].release()
......@@ -145,8 +145,8 @@ class CCom:
if self.__vdebug:
print("Time [%f] does not exist!!\n"%(time))
for select in data.com_select.split(":"):
if self.__vdebug:
print("COM select[%s]\n"%(select))
#if self.__vdebug:
# print("COM select[%s]\n"%(select))
icxv[cpt][0] = cpt
ok,mass=uns_snap.getData(select,"mass")
if mass.size==0 : # no mass
......
py/modules/simulations/cmovie.py
View file @ 9155402a
......@@ -4,6 +4,7 @@ from uns_simu import *
from csnapshot import *
from c2dplot import *
import py_unstools # rectify swig
import subprocess
from multiprocessing import Lock
import time
......@@ -60,16 +61,71 @@ class CMovie:
ok,time=uns_snap.getData("time")
print("Core [%d] time <%f>"%(data.core_id,time))
element=0
# loop on all existing component
for comp,prop,percen,newradius,extdir,fdir in self.__comp_data:
t_file=fdir+"/.time_completed"
radius=self.__newradius_0[element] # in case of newradius==0 (computed in INIT phase)
element=element+1
if not self.__checkTimeExist(time,t_file): # time not computed yet
if newradius==1: # we must compute newradius for this frame
if percen < 0 : # we use percen*(-1) as range
radius=percen*-1
else:
status,radius=uns_snap.getMaxRadius(select=comp,percen=percen)
if self.__vdebug:
print("Computed radius [%f]\n"%(radius))
radius = radius*1.1 # add 10% more radius
filename=uns_snap.getFileName()
basefile=os.path.basename(filename)
no=int(basefile.split("_")[1]) # compute no
c=C2dplot() # object
if prop=="none":
c.draw(uns_snap=uns_snap,select=comp,outdev=fdir+"/frame",no=no,prop=prop,cb=0,rrange=30.,times=time.item())
if prop=="age":
pixel=40
else:
pixel=20
myframe=fdir+"/frame.%05d.gif"%(no)
print("MY FRAME = ",myframe)
c=C2dplot(pixel=pixel) # object
if prop=="none" or prop=="rho":
c.draw(uns_snap=uns_snap,select=comp,outdev=fdir+"/frame",title=basefile,no=no,prop=prop,cb=0,rrange=radius,times=time.item())
else :
c.draw(uns_snap=uns_snap,select=comp,outdev=fdir+"/frame",title=basefile,no=no,prop=prop,cb=1,rrange=radius,times=time.item(),psort=1,itf=0)
if not os.path.isfile(myframe): # no frame created
subprocess.call(["convert","-size","850x680","xc:black",myframe]) # create a black frame
self.__saveNewTime(time,t_file) # save new time in ".time_completed" file
else:
if self.__vdebug:
print("MOVIE: skip time[%.3f] from <%s>\n"%(time,t_file))
#
# __checkTimeExist
#
def __saveNewTime(self,time,t_file):
"""
write in t_file a new time completed
*IN*
time : current time
t_file : filename
"""
# lock process
self.__analysis.lock[self.__analysis.lock_id].acquire()
if self.__vdebug:
print("Saving a new movie time[%.3f] in <%s>\n"%(time,t_file))
try:
f=open(t_file,"a")
f.write("%e\n"%(time))
f.close()
except:
print("\n\nUnable to save time[%.3f] in <%s>, aborting....\n"%(time,t_file))
sys.exit()
# release process
self.__analysis.lock[self.__analysis.lock_id].release()
#
# __checkTimeExist
#
......@@ -89,16 +145,20 @@ class CMovie:
try:
f=open(t_file,"r")
while True:
atime=float(f.readline())
if (atime-0.001)<time and (atime+0.001)>time:
f.close()
return True
atime=map(np.float,(f.readline().split()))
if (len(atime)>0):
if (atime[0]-0.001)<time and (atime[0]+0.001)>time:
f.close()
return True
else:
f.close()
return False
except EOFError:
pass
else:
if self.__vdebug:
print("FILE <%s> does not exist...\n"%(t_file))
return False
......@@ -154,8 +214,8 @@ class CMovie:
print("CMovie#Warning : component <%s> from select <%s> does not exist...\n"\
%(comp,colon_s))
else: # comp exist
fdir=self.__movie_dir+"/work/rsnap_"+newradius+"_percen_"+\
percen+"_prop_"+prop+"/"+comp
fdir=self.__movie_dir+"/w/r_"+newradius+"_p_"+\
percen+"_pp_"+prop+"/"+comp
self.__comp_data.append([comp,prop,float(percen),int(newradius),extdir,fdir])
......@@ -175,12 +235,13 @@ class CMovie:
s=CSnapshot(data.slist[-1],comp)
ok=s.nextFrame("xv")
if ok:
status,radius=s.getMaxRadius(select=comp,percen=99.)
# compute radius for at percen particles for the last snapshot
status,radius=s.getMaxRadius(select=comp,percen=percen)
else:
print ("Cannot compute newradius_0 for [%s] on file "%(comp,data.slist[-1]))
s.close() # close test snapshot
print("Comp [%s] newradius_0 = <%f>"%(comp,radius))
self.__newradius_0.append([radius]) # store newradius_0
self.__newradius_0.append(radius) # store newradius_0
......
py/modules/simulations/crectify.py
View file @ 9155402a
......@@ -102,7 +102,7 @@ class CRectify:
rcut=float(cod_R[1])
codf=cod_R[0]
else:
rcut=100.
rcut=50.
codf=cod_R[0]
# check components exist?
......@@ -170,15 +170,19 @@ class CRectify:
ok1,pos = uns_snap.getData(select,"pos" )
ok2,vel = uns_snap.getData(select,"vel" )
ok3,mass = uns_snap.getData(select,"mass")
print("ok1 ok2 ok3 ",ok1,ok2,ok3)
#print("ok1 ok2 ok3 ",ok1,ok2,ok3)
if pos.size>1 and vel.size>1 and mass.size>1 :
if self.__vdebug:
print("Rectifing time (%f) select <%s> cod[%s]\n"%(time.item(),select,cod_file))
ok,eigen_v=c.computeEigenVectors(pos,vel,mass,mass,16,time.item(),False,False,str(cod_file),rcut)
print("Time [%f] eigen:\n"%(time),eigen_v)
data.lock[data.lock_id].acquire()
self.__saveEigenVectors(pre_rect_file,eigen_v)
data.lock[data.lock_id].release()
try:
ok,eigen_v=c.computeEigenVectors(pos,vel,mass,mass,16,time.item(),False,False,str(cod_file),rcut)
print("Time [%f] eigen:\n"%(time),eigen_v)
data.lock[data.lock_id].acquire()
self.__saveEigenVectors(pre_rect_file,eigen_v)
data.lock[data.lock_id].release()
except:
print("UNABLE TO COMPUTE EIGENVECTORS for time (%f) select <%s> cod[%s]\n"%(time.item(),select,cod_file))
......
py/modules/simulations/ctree.py
View file @ 9155402a
......@@ -54,7 +54,7 @@ class CTree:
else:
threshold = min(threshold,self.__tree.getNbody())
print("Threshold = ",threshold)
#print("Threshold = ",threshold)
ok,levels=self.getLevels()
#idx=np.argsort(levels)
......@@ -63,15 +63,15 @@ class CTree:
p=np.reshape(self.__pos,(-1,3))
p=p[idx[0:threshold],]
print(p.shape,p.size)
#print(p.shape,p.size)
m=self.__mass[idx[0:threshold]]
p=np.reshape(p,(-1,))
print(p.shape,p.size)
print (p[:,])
#print(p.shape,p.size)
#print (p[:,])
c=CFalcon() # new falcon object
ok,rho,hsml=c.getDensity(p,m) # compute density
print (ok)
#print (ok)
v=None
if self.__vel is not None:
......@@ -108,7 +108,7 @@ class CTree:
else:
threshold = min(threshold,self.__tree.getNbody())
print("Threshold = ",threshold)
#print("Threshold = ",threshold)
ok,radius=self.__tree.get_closest_distance_to_mesh(self.__tree.getNbody())
#idx=np.argsort(radius)
......@@ -117,15 +117,15 @@ class CTree:
p=np.reshape(self.__pos,(-1,3))
p=p[idx[0:threshold],]
print(p.shape,p.size)
#print(p.shape,p.size)
m=self.__mass[idx[0:threshold]]
p=np.reshape(p,(-1,))
print(p.shape,p.size)
print (p[:,])
#print(p.shape,p.size)
#print (p[:,])
c=CFalcon() # new falcon object
ok,rho,hsml=c.getDensity(p,m) # compute density
print (ok)
#print (ok)
v=None
if self.__vel is not None:
......@@ -163,7 +163,7 @@ class CTree:
else:
threshold = min(threshold,self.__tree.getNbody())
print("Threshold = ",threshold)
#print("Threshold = ",threshold)
ok,levels=self.getLevels() # get octree levels for each particles
......@@ -179,15 +179,15 @@ class CTree:
p=np.reshape(self.__pos,(-1,3)) # convert flat array to [nbody,3]
p=p[shuffle_idx_threshold,] # get randomized thresholded indexes for pos
print(p.shape,p.size)
#print(p.shape,p.size)
m=self.__mass[shuffle_idx_threshold] # get randomized thresholded indexes for mass
p=np.reshape(p,(-1,)) # flatten pos array
print(p.shape,p.size)
print (p[:,])
#print(p.shape,p.size)
#print (p[:,])
c=CFalcon() # new falcon object
ok,rho,hsml=c.getDensity(p,m) # compute density
print (ok)
#print (ok)
v=None
if self.__vel is not None:
......
py/modules/uns_simu.py
View file @ 9155402a
......@@ -19,7 +19,7 @@ class UnsSimu:
if (dbname is not None ) :
self.dbname = dbname
if (not os.path.isfile(dbname)) :
print('sqlite3 db file ['+dbname+'] does not exist...\n')
print('sqlite3 db file ['+dbname+'] does not exist...\n',file=sys.stderr)
self.__status = 0;
else :
self.__conn = sqlite3.connect(dbname)
......@@ -31,7 +31,7 @@ class UnsSimu:
c = self.__conn.cursor()
c.execute("select name from eps") # try a request to check if sdb3 is valid
except sqlite3.Error:
print ("File [",self.dbname,"] is not a valid sqlite3 DB")
print ("File [",self.dbname,"] is not a valid sqlite3 DB",file=sys.stderr)
sys.exit()
......@@ -50,7 +50,7 @@ class UnsSimu:
try:
gp=open(inputfile,"r")
except IOError:
print ("no file ["+inputfile+"], will use default :"+self.dbname)
print ("no file ["+inputfile+"], will use default :"+self.dbname,file=sys.stderr)
return False,self.dbname
for line in gp:
......@@ -61,7 +61,7 @@ class UnsSimu:
if ('dbname' in gparam and gparam['dbname']!="") :
if (not os.path.isfile(gparam['dbname'])) :
print ("FILE [",gparam['dbname'],"] does not exist !!")
print ("FILE [",gparam['dbname'],"] does not exist !!",file=sys.stderr)
return False, self.dbname
else :
return True,gparam['dbname']
......@@ -75,7 +75,7 @@ class UnsSimu:
r = self.getInfo(name)
if (r):
for keys in r.keys():
print ( keys," : ",r[keys]) #for row in cursor:
print ( keys," : ",r[keys],file=sys.stderr) #for row in cursor:
def getInfo(self,name=None):
if name is None:
......@@ -98,7 +98,7 @@ class UnsSimu:
if (r and ( r['type']=="Gadget" or r['type']=="Gadget3") ):
if ( os.path.exists(r['dir'] )):
base = r['dir']+"/"+r['base']
print ("base : ", base)
print ("base : ", base,file=sys.stderr)
snap_list = sorted(glob.glob(base+'_?'))
snap_list = snap_list+sorted(glob.glob(base+'_??'))
snap_list = snap_list+sorted(glob.glob(base+'_???'))
......
Markdown is supported
0%