more fixes

py/mains/merging_time.py

+#!/usr/bin/python
+
+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')
+from uns_simu import *
+from ccod import * 
+import argparse
+
+
+# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+# commandLine, parse the command line 
+def commandLine():
+    dbname=None
+    ncores=None
+    pngfile=None
+    txtfile=None
+    cod1=None
+    cod2=None
+    
+    threshold=10000
+     # help
+    parser = argparse.ArgumentParser(description="Display merging time",
+                                    formatter_class=argparse.ArgumentDefaultsHelpFormatter)
+    # options
+    parser.add_argument('simname', help='Simulation name')
+    parser.add_argument('--cod1', help='first cod file',default=cod1)
+    parser.add_argument('--cod2', help='second cod file',default=cod2)
+    parser.add_argument('--pngfile', help='png filename, if None interactive plot',default=pngfile)
+    parser.add_argument('--dbname',help='UNS database file name', default=dbname)
+    parser.add_argument('--verbose',help='verbose mode', default=False)
+
+     # parse
+    args = parser.parse_args()
+    
+    if args.pngfile=="None":
+        args.pngfile=None
+    # start main funciton
+    process(args)
+
+
+# -----------------------------------------------------
+# process, is the core function
+def process(args):
+    cod = CCod(simname=args.simname,verbose_debug=args.verbose,dbname=args.dbname)
+    cod.computeMergingTime(halo_1=args.cod1,halo_2=args.cod2,pngfile=args.pngfile,txtfile=None)
+
+    
+
+# -----------------------------------------------------
+# main program
+if __name__ == '__main__':
+    commandLine()

py/mains/uns_select_byid.py

@@ -85,13 +85,14 @@ def getComponentByIds(snap,massratio,idgal):
     m1    = float(m1)         # m1 mass of gal1
     m2    = m1 / m1om2        # m2 mass of gal2
     mtot  = m1 + m2           # mtot total mass
-    
+    print ("m1om2=%f m1=%f m2=%f mtot=%f"%(m1om2,m1,m2,mtot),file=sys.stderr)
     #compute #bodies per galaxy
     nbody    = np.zeros(2)             # array to store #bodies ofr each particles
     nbody[0] = m1 * snap.nbody / mtot  # nbody gal1
     nbody[1] = m2 * snap.nbody / mtot  # nbody gal2
     nbody=nbody.astype(int)            # convert to int
-
+    
+    print("nbody :",nbody)
 
     offset = snap.ids.min() # first ID value for the component
     if (idgal==1):

py/modules/simulations/ccod.py

@@ -3,6 +3,7 @@ from __future__ import print_function
 from uns_simu import *
 from csnapshot import *
 from cfalcon import *
+from ctree import * 
 
 from multiprocessing import Process, Lock,Pool
 import multiprocessing
@@ -11,19 +12,29 @@ import time
 import os
 import signal
 
+import matplotlib.pyplot as plt
+
+
 #
 # ----
 #
 class CCod:
     "compute Center Of Density on UNS snapshots"
     __sql3=None
-    __simname=None
+    simname=None
     __r=None
     __slist=None # snap list
     __verbose=False
     __dbname=None
     __cod_file_base=None
-
+    __ctree_threshold=10000
+    __fastcod=True
+    __halo_part=np.zeros(10)
+    __is_multiple_halo=False
+    __ok_halo_N=False
+    __halo_comp=None
+    __halo_N=None
+    __COD_DIR_NAME="codFAST"
     #
     # ----
     #
@@ -45,29 +56,45 @@ class CCod:
         if self.__vdebug:
             self.__sql3.printInfo(simname)
         self.__slist = self.__sql3.getSnapshotList()
-
+        self.__cod_file_base=self.__r["dir"]+"/ANALYSIS/"+self.__COD_DIR_NAME
+        
     #
     # ----
     #        
     # compute COD according "select"ion component
-    def compute(self,select,ncores=None,cod_file=None):
+    def compute(self,select,ncores=None,cod_file=None, fastcod=True,threshold=10000):
         """
         compute Center of Density on list of list of components
 
-        select: list of components, example :
+        select: string
+                list of components, example :
                 simple list : select=gas,disk
                 two list    : select=gas,disk:stars  (note ':' to separate lists)
 
-        ncores=#cores used for multiprocessing, if None then all availables cores will be used
+        ncores : integer
+                #cores used for multiprocessing, if None then all availables cores will be used
 
-        cod_file=cod's file name, if None it will be automatically build to be stored in sim's directory
+        cod_file : string
+                  cod's file name, if None it will be automatically build to be stored in sim's directory
+
+        fastcod  : boolean (default true)
+                    True : compute fas_cod using octree 
+                    rue : compute fas_cod using octree 
+                    
+        threshold : number of max dens particles estimated by ctree depht keept
+                    if > 0 use this value as number of particles
+                    if < 0 use as percentage of particles
+                    default (10000 particles)
         """
 
         if self.__r is None:
             print ("Simulation",self.simname," does not exist")
             return False
+
+        self.__ctree_threshold=threshold
+        self.__fastcod=fastcod
         if cod_file is None:
-            self.__cod_file_base=self.__r["dir"]+"/ANALYSIS/codp"
+            self.__cod_file_base=self.__r["dir"]+"/ANALYSIS/"+self.__COD_DIR_NAME
             # create directory
             if not os.path.isdir(self.__cod_file_base):
                 try:
@@ -92,11 +119,14 @@ class CCod:
              
         if self.__vdebug:
             print("Best select :",self.__best_select)
+        
 
         # get half of list of snapshot to check if all components exist
+
         half_snap=self.__slist[len(self.__slist)/2]
         if self.__vdebug:
             print("Half snapshot =",half_snap)
+        #sys.exit()
         test_snap=CSnapshot(half_snap, self.__best_select)
         ok=test_snap.nextFrame()
 
@@ -104,6 +134,9 @@ class CCod:
         self.__new_select=""
         for colon_s in select.split(":"):
             ok,data=test_snap.getData(colon_s,"pos")
+            if not ok:
+                ok,comp,ii=self.__checkSelectCompN(colon_s)
+                print(ok,comp,ii)
             if ok:
                 if len(self.__new_select)==0:
                         sep=""
@@ -113,8 +146,10 @@ class CCod:
         if self.__vdebug:
             print("new select :",self.__new_select)
         test_snap.close()
+        self.__is_multiple_halo = self.__getHaloParticlesNumber() #
         self.__startProcess(select,ncores,cod_file)
-    
+        self.computeMergingTime()
+
     #
     # ----
     #
@@ -144,6 +179,7 @@ class CCod:
         
         # lock to control access to file
         lock=Lock()
+
         # loop on all #cores and start a process
         for p in range(ncores):
             p = Process(target=self.__codProcess, args=(q,n,lock,))  # create
@@ -187,7 +223,7 @@ class CCod:
                 #my_snap=queue_list.get()  # do not use this, could block if nothing to get
                 #time.sleep(0.01)
                 uns_snap=CSnapshot(my_snap,self.__best_select)
-                ok=uns_snap.nextFrame("mxv")
+                ok=uns_snap.nextFrame("mxvI")
                 okt,time=uns_snap.getData("time")
 
                 # checkCodTime
@@ -201,18 +237,40 @@ class CCod:
                     if not exist_time:
                         if self.__vdebug:
                             print (">> computing density on: ", select)
+
+                        # check select if it is halo_N string
+                        #
+                        self.__ok_halo_N,self.__halo_comp,self.__halo_N=self.__checkSelectCompN(select)
+                        
+                        
+                        if self.__ok_halo_N: # then it's an halo_X like
+                            select=self.__halo_comp # rewrite halo_N to halo
+                            
                         ok,pos=uns_snap.getData(select,"pos")
+                        
                         cxv=np.zeros(6)
                         if ok and pos.size/3 > 32: # more than 32 particles at that time                   
                             ok,mass=uns_snap.getData(select,"mass")
                             ok,vel =uns_snap.getData(select,"vel")
 
-                            c=CFalcon() # new falcon object
-                            ok,rho,hsml=c.getDensity(pos,mass) # compute density
-
-                            # compute cod
-
-                            cxv=uns_snap.center(pos,vel,rho*mass)
+                            
+                            if self.__ok_halo_N: # WE MUST RE ORDER PARTICLES
+                                print (">> select =",select)
+                                ok,id=uns_snap.getData(select,"id")
+                                print ("ok, id.size=",ok,id.size)
+                                ok,mass,pos,vel=self.__extractHalo(id,mass,pos,vel)
+                                print("after extracting halo mass.size:", mass.size)
+                                #!!!!sys.exit()
+                            if not self.__fastcod:
+                                #c=CFalcon() # new falcon object
+                                ok,rho,hsml=c.getDensity(pos,mass) # compute density
+                                
+                                # compute cod
+                            
+                                cxv=uns_snap.center(pos,vel,rho*mass)
+                            else : # fastcod
+                                ctree=CTree(pos,None,mass)  # instiate a ctree object
+                                cxv=ctree.fastCod3(self.__ctree_threshold)
                             print ("TCXV : ", time,cxv)
 
                         else: # no data
@@ -242,6 +300,152 @@ class CCod:
 
         #print "Core [",n,"] got snap : ",queue_list.get()
 
+    #
+    # ----
+    #
+    def __extractHalo(self,id,mass,pos,vel):
+        print("id =",id.size)
+        id_sort=np.argsort(id) # sort according to ID order
+        print("id_sort=",id_sort.size)
+        offset=0
+        myidsort=np.zeros(1)
+        if self.__halo_N==1:   # first halo
+            mylast=self.__halo_part[0]
+            myidsort = id_sort[0:mylast]
+        else :                 # second halo
+            mylast=self.__halo_part[1]
+            offset = self.__halo_part[0]+mylast
+            myidsort = id_sort[self.__halo_part[0]:offset]
+        print("myidsort : ",myidsort,mylast,self.__halo_part)
+        pos = np.reshape(pos,(-1,3)) # reshape pos to nbody,3
+        pos = np.reshape(pos[myidsort],-1) # select ids
+        vel = np.reshape(vel,(-1,3)) # reshape vel to nbody,3
+        vel = np.reshape(vel[myidsort],-1) # select ids
+        
+        return True,mass[myidsort],pos,vel
+    
+    #
+    # ----
+    #
+    # check in select string component is indexed like halo_1 halo_2 etc..
+    def __checkSelectCompN(self,ss,max=10):
+        xx=ss.find("_")
+        if (xx!=-1):
+            comp=ss[0:xx]
+            #print "comp=",comp
+            for i in range(max):
+                if (ss==comp+"_%d"%(i)):
+                    return True,comp,i-1
+        return False,None,None
+
+    #
+    # ----
+    #
+    # try to read model_param.txt file in SIMULATION directory and fill up an array with #particles per galaxy/halo
+    def __getHaloParticlesNumber(self,model_param="model_param.txt"):
+
+        ok=False
+        infile=self.__r['dir']+"/"+model_param
+        if os.path.exists(infile):
+            print("Model file exist :",infile)
+            fh = open(infile)
+            for line in fh:
+                line=line.strip()
+                if line.find('_halo_nbody')>=0 :
+                    ok=True
+                    sline=line.split()
+                    ind=int(((sline[0].split("_"))[0].split("l"))[1])
+                    print (sline,ind,type(self.__halo_part))
+                    
+                    self.__halo_part[ind-1] = int(sline[1])
+
+        self.__halo_part=self.__halo_part[np.where(self.__halo_part>0)]
+        print (self.__halo_part,self.__halo_part.size)
+        return ok
+
+    
+
+    #
+    # ----
+    # Compute mering time between two halos galaxy
+    def computeMergingTime(self,halo_1=None,halo_2=None,simname=None,outdir="./",txtfile="merging_time.txt",pngfile="merging_time.png",dmax=1.0):
+        if halo_1 is None: # we are in a simulation
+            if self.__cod_file_base is None: # no cod computation ?
+                return 0
+            else :
+                c_halo_1=self.__cod_file_base+"/"+self.simname+".halo_1.cod"
+                c_halo_2=self.__cod_file_base+"/"+self.simname+".halo_2.cod"
+                outdir=self.__cod_file_base
+                simname=self.simname
+        else:
+            c_halo_1=halo_1
+            c_halo_2=halo_2
+            simname="None"
+        if os.path.isfile(c_halo_1) and os.path.isfile(c_halo_2):
+            print(c_halo_1,c_halo_2)
+            # load files in numpy arrays
+            data1=np.loadtxt(c_halo_1,dtype=np.float32)
+            data2=np.loadtxt(c_halo_2,dtype=np.float32)
+            # sort according time
+            time1_sort=data1[:,0].argsort()
+            time2_sort=data2[:,0].argsort()
+            n = min(time1_sort.size,time2_sort.size) # get minimum in case of #lines
+            time1_sort = time1_sort[0:n]
+            time2_sort = time2_sort[0:n]
+            # compute distance
+            x2=(data1[time1_sort][:,1]-data2[time2_sort][:,1])**2
+            y2=(data1[time1_sort][:,2]-data2[time2_sort][:,2])**2
+            z2=(data1[time1_sort][:,3]-data2[time2_sort][:,3])**2
+            distance=np.sqrt(x2+y2+z2)
+            # select times only when distance > dmax
+            a=(distance>=dmax)
+            x=np.where(a)
+            if x[0].size==0:
+                print("\n\n\n NO MERGING TIME FOUND \n\n\n")
+                return -1
+            id_merge = x[0][-1]+1 # id matching to merging time
+            merging_time=data1[id_merge][0]
+            print("Merging time =",merging_time,"\nmerging distance =",
+                  distance[id_merge],"\nmax distance from merging=",
+                  distance[data1[time1_sort][:,0]>merging_time].max())
+
+            # plot
+            # first plot: plot distance between 2 halos on whole simulation
+            plt.subplot(211)
+            plt.plot(data1[time1_sort][:,0],distance[time1_sort],'b.')
+            plt.axvline(merging_time,color='g',dashes=(10,3),label="merging")
+            plt.title('Distance between 2 Halos : '+simname)
+            plt.xlabel('Gyears')
+            plt.ylabel('Kpc')
+            ax = plt.axis() # get axis coordinates
+            plt.text(merging_time+0.01,(ax[3]-ax[2])/2,"%.3f"%(merging_time)) # print merging time
+            #plt.show()
+            # second plot: plot distance between 2 halos from merging time
+            plt.subplot(212)
+            #plt.title('Distance between 2 Halos')
+            plt.xlabel('Gyears')
+            plt.ylabel('Kpc')
+            select_t=data1[time1_sort][:,0]>merging_time
+            plt.plot(data1[select_t][:,0],distance[select_t],'b.')
+            plt.axvline(merging_time,color='g',dashes=(10,3))
+            ax = plt.axis() # get axis coordinates
+            plt.text(merging_time+0.01,(ax[3]-ax[2])/2,"%.3f"%(merging_time)) # print merging time
+            plt.tight_layout()
+            #plt.show()
+
+            if txtfile is not None:
+                m_txt=outdir+"/"+txtfile
+                print("Saving txt ",m_txt)
+                f=open(m_txt,"w")
+                f.write(str(merging_time)+"\n")
+                f.close()
+
+            if pngfile is not None:
+                img_png=outdir+"/"+pngfile
+                print("Saving plot :",img_png)
+                plt.savefig(img_png)
+            else:
+                plt.show()
 
     #
     # ----

py/modules/simulations/ctree.py

@@ -12,15 +12,17 @@ class CTree:
     __pos=None
     __mass=None
     __vel=None
+    __time=0.0
 # -----------------------------------------------------
 #
-    def __init__(self,pos,vel=None,mass=None,fcells=0.9,rsize=0.4):
+    def __init__(self,pos,vel=None,mass=None,fcells=0.9,rsize=0.4,time=0.0):
         self.__pos=pos
         self.__mass=mass
         self.__vel=vel
+        self.__time=float(time)
         if (type(pos[0])==np.float32):
             self.__tree=ut.CTreeF(pos,mass,fcells,rsize)
-            
+
         else: # assume 64 bits (double)
             self.__tree=ut.CTreeD(pos,mass,fcells,rsize)
 
@@ -142,7 +144,7 @@ class CTree:
 
 # -----------------------------------------------------
 #  
-    def fastCod3(self,threshold=10000):
+    def fastCod3(self,threshold=10000,outfile=None):
         """
         compute cod using octree
 
@@ -191,12 +193,16 @@ class CTree:
              v=v[shuffle_idx_threshold, ]
              v=np.reshape(v,(-1,))
         cxv=CSnapshot(None).center(p,v,m*rho)
-        
-        unso=CunsOut("/home/jcl/x_x","gadget2")
-        
-        unso.setArrayF("gas","pos",p)
-        unso.setArrayF("gas","mass",m)
-        unso.setArrayF("gas","rho",rho)
-        unso.setArrayF("gas","hsml",hsml)
-        unso.save()
-        return cxv[0:3]
+        if v is None: # no velocities
+            cxv[3]=0.0
+            cxv[4]=0.0
+            cxv[5]=0.0
+        if outfile is not None:
+            unso=CunsOut(outfile,"nemo")
+            unso.setValueF("time",self.__time)
+            unso.setArrayF("all","pos",p)
+            unso.setArrayF("all","mass",m)
+            unso.setArrayF("all","rho",rho)
+            unso.setArrayF("all","hsml",hsml)
+            unso.save()
+        return cxv

py/modules/uns_simu.py

@@ -86,7 +86,7 @@ class UnsSimu:
         if name is None:
             name=self.name
         r = self.getInfo(name)
-        if (r and r['type']=="Gadget"):
+        if (r and ( r['type']=="Gadget" or r['type']=="Gadget3") ):
             if ( os.path.exists(r['dir'] )):
                 base = r['dir']+"/"+r['base']
                 print ("base : ", base)

py/test/test_analysis.py

+#!/usr/bin/python
+
+# coding: utf-8
+
+## Example of howto get gravity using falcon library directly from Python 
+
+# In[1]:
+
+#from py_unstools import *         # import py_unstools package
+from py_unsio import *            # import py_unsio package
+import numpy as np                # arrays are treated as numpy arrays
+import time
+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')
+from uns_simu import *
+from cuns_analysis import *
+
+
+
+# instantiate CUnsAnalysis object
+a=CUnsAnalysis()
+
+t0=time.clock()
+
+simname="mdf001%50" 
+components="stars" 
+times="all" 
+verbose=False
+uns = CunsIn(simname,components,times,verbose)
+bits=""         # select properties, "" means all
+ok=uns.nextFrame(bits)   # load data from disk
+
+# Load positions and masses
+comp="stars" 
+
+ok,pos = uns.getArrayF(comp,"pos")
+ok,vel = uns.getArrayF(comp,"vel")
+ok,mass = uns.getArrayF(comp,"mass")
+
+print "Nbody = ",pos.size/3
+
+# softening
+eps=0.05
+
+# get density
+ok,rho,hsml=a.falcon.getDensity(pos,mass,eps)
+t2=time.clock()
+print ok,rho,hsml
+
+print "POS =", pos 
+com = a.snapshot.center(pos,vel,rho*mass,center=False)
+
+print "POS =", pos
+print "VEL =", vel
+
+
+print "Coordinates center :", com
+

py/mains/test_cod.py → py/test/test_cod.py

@@ -13,12 +13,16 @@ import argparse
 def commandLine():
     dbname=None
     ncores=None
+    fastcod=True
+    threshold=10000
      # help
     parser = argparse.ArgumentParser(description="UNS test COD python",
                                     formatter_class=argparse.ArgumentDefaultsHelpFormatter)
     # options
     parser.add_argument('simname', help='Simulation name')
     parser.add_argument('select', help='Selected component')
+    parser.add_argument('--fastcod', help='compute density by selecting particles from octree',default=fastcod)
+    parser.add_argument('--threshold', help='number of particles used for fastcod (<0 percentage)',default=threshold,type=int)
     parser.add_argument('--ncores', help='Use ncores, None means all',default=ncores,type=int)
     parser.add_argument('--dbname',help='UNS database file name', default=dbname)
     parser.add_argument('--verbose',help='verbose mode', default=False)
@@ -34,7 +38,7 @@ def commandLine():
 # process, is the core function
 def process(args):
     cod = CCod(simname=args.simname,verbose_debug=args.verbose,dbname=args.dbname)
-    cod.compute(args.select,args.ncores)
+    cod.compute(select=args.select,ncores=args.ncores,fastcod=args.fastcod,threshold=args.threshold)
     
 
 # -----------------------------------------------------

py/test/test_ctree.py

+#!/usr/bin/python
+
+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')
+from uns_simu import *
+from ctree import * 
+from csnapshot import *
+import argparse
+
+
+# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+# commandLine, parse the command line 
+def commandLine():
+    dbname=None
+    ncores=None
+    threshold=10000
+     # help
+    parser = argparse.ArgumentParser(description="UNS test COD python",
+                                    formatter_class=argparse.ArgumentDefaultsHelpFormatter)
+    # options
+    parser.add_argument('inputfile', help='Input file')
+    parser.add_argument('outfile', help='output file')
+    parser.add_argument('select', help='Selected component')
+    parser.add_argument('--threshold', help='number of particles used for cod (<0 percentage)',default=threshold,type=int)
+    parser.add_argument('--ncores', help='Use ncores, None means all',default=ncores,type=int)
+    parser.add_argument('--dbname',help='UNS database file name', default=dbname)
+    parser.add_argument('--verbose',help='verbose mode', default=False)
+
+     # parse
+    args = parser.parse_args()
+    if args.outfile=="None":
+        args.outfile=None