new 2d plot in devlopment

py/mains/uns_c2dplot.py → py/mains/uns_2dpgplot.py

@@ -7,7 +7,7 @@ import sys
 import argparse,textwrap
 sys.path=['/home/jcl/works/GIT/uns_projects/py/modules/','/home/jcl/works/GIT/uns_projects/py/modules/simulations']+sys.path
 
-from simulations.c2dplot import *
+from simulations.c2dpgplot import *
 from simulations.csnapshot import *
 
 # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
@@ -15,7 +15,7 @@ from simulations.csnapshot import *
 def commandLine():
 
     # help
-    parser = argparse.ArgumentParser(description="Compute 2D image using uns_projects c2dplot engine",
+    parser = argparse.ArgumentParser(description="Compute 2D image using uns_projects c2dpgplot engine",
                                     formatter_class=argparse.ArgumentDefaultsHelpFormatter)
 
     parser.add_argument('snapshot', help="uns input snapshot",default=None)
@@ -42,7 +42,7 @@ def process(args):
         uns=CSnapshot(args.snapshot,args.component,verbose_debug=args.verbose)
         ok=uns.nextFrame("")
         if ok:
-            c=C2dplot()
+            c=C2dpgplot()
             c.draw(uns_snap=uns,select=args.component,outdev=args.dev,rrange=args.range,\
                    prop=args.prop,psort=args.psort,pfname=args.pfname,cb=args.cb,cmap=args.cmap)
         else:

py/mains/uns_select_byid.py

-#!/usr/bin/python
+#!/usr/bin/env python
 
 from __future__ import print_function
 import sys

py/modules/simulations/c2dplot.py → py/modules/simulations/c2dpgplot.py

@@ -13,7 +13,7 @@ import os
 #
 # ----
 #
-class C2dplot:
+class C2dpgplot:
     """
     Draw 2d image using uns_project C2dplot engine
     """

py/modules/simulations/ccod.py

@@ -61,6 +61,10 @@ class CCod:
             self.__slist = self.__sql3.getSnapshotList()
             if self.__r is not None:
                 self.__cod_file_base=self.__r["dir"]+"/ANALYSIS/"+self.__COD_DIR_NAME
+            else:
+                message="From CCOD UNS simulation [%s] does not belong to unsio Database"%(self.simname)
+                raise Exception(message)
+            self.__is_multiple_halo = self.__getHaloParticlesNumber() #
         else :
             print("COD analysis is activated",file=sys.stderr)
             self.__r=analysis.sim_info
@@ -72,7 +76,454 @@ class CCod:
                 else:
                     print("SINGLE HALO present\n",file=sys.stderr)
             self.__smartAnalysis(analysis)
+
+    #
+    # isMultipleHalo
+    #                
+    def isMultipleHalo(self):
+        """
+        return if simulation has multiple halo
+        """
+        return self.__is_multiple_halo
+    #
+    # ----
+    #        
+    # compute COD according "select"ion component
+    def compute(self,select,ncores=None,cod_file=None, fastcod=True,threshold=10000):
+        """
+        compute Center of Density on list of list of components
+
+        select: string
+                list of components, example :
+                simple list : select=gas,disk
+                two list    : select=gas,disk:stars  (note ':' to separate lists)
+
+        ncores : integer
+                #cores used for multiprocessing, if None then all availables cores will be used
+
+        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",file=sys.stderr)
+            return False
+
+        self.__ctree_threshold=threshold
+        self.__fastcod=fastcod
+        if cod_file is None:
+            self.__cod_file_base=self.__r["dir"]+"/ANALYSIS/"+self.__COD_DIR_NAME
+            # create directory
+            if not os.path.isdir(self.__cod_file_base):
+                try:
+                    os.makedirs( self.__cod_file_base)
+                except:
+                    print ("Cannot create directory <%s>, aborting"%(self.__cod_file_base),file=sys.stderr)
+                    sys.exit()
+        if self.__vdebug:
+            print ("COD FILE BASE =",self.__cod_file_base,file=sys.stderr)
+        # Figure out best selected components
+        self.__best_select=""
+        for colon_s in select.split(":"):
+            for coma_s in colon_s.split(","):
+                if self.__vdebug:
+                    print (">> :",coma_s,file=sys.stderr)
+                if self.__best_select.find(coma_s)==-1: # substring not found
+                    if len(self.__best_select)==0:
+                        sep=""
+                    else:
+                        sep=","
+                    self.__best_select += sep+ coma_s
+             
+        if self.__vdebug:
+            print("Best select :",self.__best_select,file=sys.stderr)
+        
+
+        # 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,file=sys.stderr)
+        #sys.exit()
+        test_snap=CSnapshot(half_snap, self.__best_select)
+        ok=test_snap.nextFrame()
+
+        # rebuild select string with existing components
+        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,file=sys.stderr)
+            if ok:
+                if len(self.__new_select)==0:
+                        sep=""
+                else:
+                        sep=":"
+                self.__new_select += sep+ colon_s
+        if self.__vdebug:
+            print("new select :",self.__new_select,file=sys.stderr)
+        test_snap.close()
+        self.__is_multiple_halo = self.__getHaloParticlesNumber() #
+        self.__startProcess(select,ncores,cod_file)
+        self.computeMergingTime()
+
+    #
+    # ----
+    #
+    # launch parallel process
+    def __startProcess(self,select,ncores=None,code_file=None):
+        """
+        Start COD computation in parallel on nores
+        """
+        # compute cores
+        if ncores is None:
+            ncores=multiprocessing.cpu_count()
+        if self.__vdebug:
+            print ("#cores used :",ncores,file=sys.stderr)
+        
+        # Feed QUEUE with list of snapshots
+        q = multiprocessing.Queue()
+        for snap in self.__slist: # loop on list os snsphots
+            #print "SNAP put :",snap
+            q.put(snap) # out a snapshot in the list
+
+        # start process jobs
+        processes=[]  # list of processes
+        n=0
+        
+        # manage signint signal
+        #original_sigint_handler = signal.signal(signal.SIGINT, signal.SIG_IGN) # IGNOROE before process creation
+        
+        # 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
+            p.start()  # start
+            processes.append(p) # append list of process, used for joining
+            n += 1
+
+        # wait all processes to complete
+        try:
+            for p in processes:
+                print ("waiting.. ",p,file=sys.stderr)
+                p.join()
+        except KeyboardInterrupt: # allow to interupt all workers with  CTRL+C
+            for p in processes:
+                print ("Terminating.. ",p,file=sys.stderr)
+                p.terminate()
+                p.join()
+        while not q.empty():
+            q.get(block=False)
+    #
+    # ----
+    #
+    def __codProcess(self,queue_list,n,lock):
+        """
+        Compute COD on core
+        """
+
+        my_snap=self.__findLastSnapshotWithNoCod(queue_list,n,lock) 
+        stop=False
+        cpt=0
+        first=True
+        while (not stop and my_snap is not None): # and not queue_list.empty()  ) :
+            try:
+                print ("Core [",n,"] waiting...",file=sys.stderr)
+                if first:         # first time my_snap populated with self.__findLastSnapshotWithNoCod() 
+                    first=False
+                else:
+                    my_snap=queue_list.get(True,0.05) # we must use this line
+                                                      # True means bloc during 0.05 sec,
+                                                      # if nothing, then get raise Queue.Empty exception
+                #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("mxvI")
+                okt,time=uns_snap.getData("time")
+
+                # checkCodTime
+                for select in self.__new_select.split(":"):
+                    mycod_file=self.__cod_file_base+"/"+self.simname+"."+select+".cod"
+                    lock.acquire()
+                    #exist_time=self.__isTimeInCod(mycod_file,time)
+                    exist_time,tcxv=self.getCodFromFile(mycod_file,time,n)
+                    lock.release()
+                    #print ("MYCOD =",mycod_file,time,exist_time)
+                    if not exist_time:
+                        if self.__vdebug:
+                            print (">> computing density on: ", select,file=sys.stderr)
+
+                        # 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
+
+                        if self.__ok_halo_N and not  self.__is_multiple_halo:# no multiple Halo, no need to compute
+                            continue    
+
+                        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")
+
+                            
+                            if self.__ok_halo_N and self.__is_multiple_halo : # WE MUST RE ORDER PARTICLES
+                                print (">> select =",select,file=sys.stderr)
+                                ok,id=uns_snap.getData(select,"id")
+                                print ("ok, id.size=",ok,id.size,file=sys.stderr)
+                                ok,mass,pos,vel=self.__extractHalo(id,mass,pos,vel)
+                                print("after extracting halo mass.size:", mass.size,file=sys.stderr)
+                                #!!!!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)  # instantiate a ctree object
+                                cxv=ctree.fastCod3(self.__ctree_threshold)
+                            print ("TCXV : ", time,cxv,file=sys.stderr)
+
+                        else: # no data
+                            pass 
+                        lock.acquire()
+                        out=open(mycod_file,"a+")
+                        #cxv=np.insert(cxv,0,time)
+                        #out.write(str(cxv).replace("[","").replace("]","")+"\n")
+                        out.write("%e %e %e %e %e %e %e\n"%(time,cxv[0],cxv[1],cxv[2],cxv[3],cxv[4],cxv[5]))
+                        out.close()
+                        lock.release()
+                    
+                cpt+=1
+                print ("Core [",n,"] got snap : ",my_snap,cpt,file=sys.stderr)
+            except Queue.Empty:
+                stop = True # mo more data
+                if self.__verbose:
+                    print ("Queue.empty execption trapped...",file=sys.stderr)
+            #except (KeyboardInterrupt, SystemExit):
+            #    print ("Exiting...")
+            #    #lock.release()
+            #    sys.exit()
+            #    terminate()
+            #    stop = True # 
+
+        print ("Core [",n,"] DONE !",cpt,file=sys.stderr)
+
+        #print "Core [",n,"] got snap : ",queue_list.get()
+
+    #
+    # ----
+    #
+    def __extractHalo(self,id,mass,pos,vel):
+        id_sort=np.argsort(id) # sort according to ID order
+        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):
+        """
+        check string 'ss' given as parameter is indexed like halo_X (X 
+        """
+        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,file=sys.stderr)
+            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),file=sys.stderr)
+                    
+                    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,file=sys.stderr)
+        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,seger=False, plot=True):
         
+        merging_time=-1
+
+        if halo_1 is None: # we are in a simulation
+            if self.__cod_file_base is None: # no cod computation ?
+                return 0
+            else :
+                if not seger:
+                    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:  # seger
+                    c_halo_1=self.__cod_file_base+"/../seger.analysis/g1g2_center/g1_halo.cod0.01.txt"
+                    c_halo_2=self.__cod_file_base+"/../seger.analysis/g1g2_center/g2_halo.cod0.01.txt"
+                    outdir=None
+                    simname=self.simname+"(seger)"
+                
+        else:
+            c_halo_1=halo_1
+            c_halo_2=halo_2
+            simname="None"
+        if self.__vdebug:
+            print("c_halo_1 [%s]\nc_halo_2 [%s]\n"%(c_halo_1,c_halo_2),file=sys.stderr)
+        if os.path.isfile(c_halo_1) and os.path.isfile(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",file=sys.stderr)
+                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(),file=sys.stderr)
+
+            if txtfile is not None and outdir is not None :
+                m_txt=outdir+"/"+txtfile
+                print("Saving txt ",m_txt,file=sys.stderr)
+                f=open(m_txt,"w")
+                f.write(str(merging_time)+"\n")
+                f.close()
+
+            if plot:
+                # 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[time1_sort[select_t]][:,0],distance[time1_sort[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()
+
+
+                if pngfile is not None and outdir is not None:
+                    img_png=outdir+"/"+pngfile
+                    print("Saving plot :",img_png,file=sys.stderr)
+                    plt.savefig(img_png)
+                else:
+                    plt.show()
+        return merging_time
+
+    #
+    #
+    #
+    def getMergingTime(self, txtfile="merging_time.txt"):
+        """
+        Get merging time from file
+        """
+        ok=False
+        mtime=-1
+        if self.__cod_file_base is None: # no cod computation ?
+            return ok,-1
+
+        mtf=self.__cod_file_base+"/"+txtfile
+        if os.path.isfile(mtf):
+            f=open(mtf,"r")
+            try:
+                mtime=float(f.readline())
+                ok=True
+            except EOFError:
+                print("Unable to read file <%s>\n"%(mtf),file=sys.stderr)
+                pass
+        print("File <%s> does not exist\n"%(mtf),file=sys.stderr)
+        return ok,mtime
+    
+    #
+    # isMultipleHalo
+    #                
+    def isMultipleHalo(self):
+        """
+        return if simulation has multiple halo
+        """
+        return self.__is_multiple_halo
     #
     # ----
     #        
@@ -504,6 +955,21 @@ class CCod:
         print("File <%s> does not exist\n"%(mtf),file=sys.stderr)
         return ok,mtime
     
+
+        
+    #
+    # ----
+    #
+    def getCodFromComp(self, component, time,n=-1):
+        """
+        check if current time exist in cod file build with componen
+        """
+        cod_file=self.__cod_file_base+"/"+self.simname+"."+component+".cod"
+        if self.__vdebug:
+            print("getCodFromComp file [%s]"%(cod_file),file=sys.stderr)
+        
+        return self.getCodFromFile(cod_file,time,n)
+        
         
     #
     # ----