[v2][_commonspeckle] what is the difference between chx_correlationp and chx_correlationp2

[ambarb]

From @yugangzhang 's google doc on pyCHX, chx_correlationp is the parallel computation of g2 using compressed data format.

Looks like chx_correlationp2 was added to debug g2.

Can we delete this module from .v2._commonspeckle

OR

Do we need to integrate parts of this module into the chx_correlationp or as tests for chx_correlationp

@yugangzhang or @afluerasu if this is quick for you to answer, this will save tons of time.

bash shell diff on two files

[abarbour@box64-3 ~/Repos/pyCHX/pyCHX/v2/_commonspeckle ✓ ] $ diff chx_correlationp.py chx_correlationp2.py 
4c4,6
< This module is for parallel computation of time correlation
---
> This module is for parallel computation of time correlation 
> Feb 20, 2018
> The chx_correlationp2 is for dedug g2
9c11
< from pyCHX.v2._commonspeckle.chx_libs import tqdm #common #TODO why import from chx module?
---
> from pyCHX.v2._commonspeckle.chx_libs import tqdm #common #TODO why not from chx module??
160,165c162,163
<                 else:                    
<                     S = norm.shape
<                     if len(S)>1:
<                         fra_pix[ pxlist] = v[w]/ norm[i,pxlist]   #-1.0     
<                     else:    
<                         fra_pix[ pxlist] = v[w]/ norm[pxlist]   #-1.0 
---
>                 else: 
>                     fra_pix[ pxlist] = v[w]/ norm[pxlist]   #-1.0  
170,174c168,169
<                     S = norm.shape
<                     if len(S)>1:               
<                         fra_pix[ pxlist] = v[w]/ imgsum[i]/  norm[i,pxlist] 
<                     else:    
<                         fra_pix[ pxlist] = v[w]/ imgsum[i]/  norm[pxlist] 
---
>                     fra_pix[ pxlist] = v[w]/ imgsum[i]/  norm[pxlist]            
>         
179c174,175
<         # get the current image time        
---
>         # get the current image time
>         
181c177,178
<         s.current_img_time += 1        
---
>         s.current_img_time += 1
>         
257,265c254,256
<         S = norm.shape
<         if len(S)>1:
<             norms = [ norm[ :, np.in1d(  pixelist, 
<                 extract_label_indices( np.array(ring_mask==i, dtype = np.int64))[1])]
<                     for i in np.unique( ring_mask )[1:] ] 
<         else:
<             norms = [ norm[ np.in1d(  pixelist, 
<                 extract_label_indices( np.array(ring_mask==i, dtype = np.int64))[1])]
<                     for i in np.unique( ring_mask )[1:] ] 
---
>         norms = [ norm[ np.in1d(  pixelist, 
>             extract_label_indices( np.array(ring_mask==i, dtype = np.int64))[1])]
>                 for i in np.unique( ring_mask )[1:] ] 
395,400c386,387
<                 else:          
<                     S = norm.shape
<                     if len(S)>1:
<                         fra_pix[ pxlist] = v[w]/ norm[i,pxlist]   #-1.0     
<                     else:    
<                         fra_pix[ pxlist] = v[w]/ norm[pxlist]   #-1.0 
---
>                 else:                     
>                     fra_pix[ pxlist] = v[w]/ norm[pxlist]   #-1.0 
405,409c392
<                     S = norm.shape
<                     if len(S)>1:               
<                         fra_pix[ pxlist] = v[w]/ imgsum[i]/  norm[i,pxlist] 
<                     else:    
<                         fra_pix[ pxlist] = v[w]/ imgsum[i]/  norm[pxlist]           
---
>                     fra_pix[ pxlist] = v[w]/ imgsum[i]/  norm[pxlist]           
529,538c512,517
<         S = norm.shape
<         if len(S)>1:
<             norms = [ norm[ :, np.in1d(  pixelist, 
<                 extract_label_indices( np.array(ring_mask==i, dtype = np.int64))[1])]
<                     for i in np.unique( ring_mask )[1:] ] 
<         else:
<             norms = [ norm[ np.in1d(  pixelist, 
<                 extract_label_indices( np.array(ring_mask==i, dtype = np.int64))[1])]
<                     for i in np.unique( ring_mask )[1:] ] 
<     inputs = range( len(ring_masks) ) 
---
>         norms = [ norm[ np.in1d(  pixelist, 
>             extract_label_indices( np.array(ring_mask==i, dtype = np.int64))[1])]
>                 for i in np.unique( ring_mask )[1:] ] 
> 
>     inputs = range( len(ring_masks) )    
>     
571,574c550,555
<         g2_err = np.zeros_like(g2)         
<         #g2_G = np.zeros((  int( (num_lev + 1) * num_buf / 2),  len(pixelist)) )  
<         #g2_P = np.zeros_like(  g2_G )
<         #g2_F = np.zeros_like(  g2_G )         
---
>         g2_err = np.zeros_like(g2)  
>         
>         g2_G = np.zeros((  int( (num_lev + 1) * num_buf / 2),  len(pixelist)) )  
>         g2_P = np.zeros_like(  g2_G )
>         g2_F = np.zeros_like(  g2_G ) 
>         
576a558
>     nopr_ =   np.lib.pad(  np.cumsum(nopr), [1], mode = 'constant',constant_values=(0))[:-1] 
612a595,598
>                 
>             g2_G[:, nopr_[i]: nopr_[i+1] ] =  s_Gall_qi 
>             g2_P[:, nopr_[i]: nopr_[i+1]] =   s_Pall_qi
>             g2_F[:, nopr_[i]: nopr_[i+1]] =   s_Fall_qi                        
618c604
<         return  g2[:Gmax,:],lag_steps_err[:Gmax], g2_err[:Gmax,:]/np.sqrt(nopr)
---
>         return  g2[:Gmax,:],lag_steps_err[:Gmax], g2_err[:Gmax,:]/np.sqrt(nopr), g2_G, g2_P, g2_F
624,656d609
< def cal_GPF( FD, ring_mask, bad_frame_list=None, 
<             good_start=0, num_buf = 8, num_lev = None, imgsum=None, norm=None,
<            cal_error=True ):
<     '''calculation G,P,D by using a multi-tau algorithm
<        for a compressed file with parallel calculation
<        if return_g2_details: return g2 with g2_denomitor, g2_past, g2_future
<     '''
<     FD.beg = max(FD.beg, good_start)
<     noframes = FD.end - FD.beg +1   # number of frames, not "no frames"
<     for i in range(FD.beg, FD.end):
<         pass_FD(FD,i)    
<     if num_lev is None:
<         num_lev = int(np.log( noframes/(num_buf-1))/np.log(2) +1) +1
<     print ('In this g2 calculation, the buf and lev number are: %s--%s--'%(num_buf,num_lev))
<     if  bad_frame_list is not None:
<         if len(bad_frame_list)!=0:
<             print ('%s Bad frames involved and will be discarded!'%len(bad_frame_list) )            
<             noframes -=  len(np.where(np.in1d( bad_frame_list, 
<                                               range(good_start, FD.end)))[0])   
<     print ('%s frames will be processed...'%(noframes-1))   
<     if np.min(ring_mask)==0:
<         qstart=1
<     else:
<         qstart=0
<     ring_masks = [   np.array(ring_mask==i, dtype = np.int64) 
<               for i in np.unique( ring_mask )[qstart:] ]  
<     qind, pixelist = roi.extract_label_indices(  ring_mask  )   
<     noqs = len(np.unique(qind))    
<     nopr = np.bincount(qind, minlength=(noqs+1))[qstart:]
<     if norm is not None:
<         norms = [ norm[ np.in1d(  pixelist, 
<             extract_label_indices( np.array(ring_mask==i, dtype = np.int64))[1])]
<                 for i in np.unique( ring_mask )[qstart:] ] 
658,691d610
<     inputs = range( len(ring_masks) ) 
<     pool =  Pool(processes= len(inputs) )
<     internal_state = None       
<     print( 'Starting assign the tasks...')    
<     results = {}    
<     if norm is not None: 
<         for i in  tqdm( inputs ): 
<             results[i] =  apply_async( pool, lazy_one_timep, ( FD, num_lev, num_buf, ring_masks[i],
<                         internal_state,  bad_frame_list, imgsum,
<                                     norms[i], cal_error  ) ) 
<     else:
<         #print ('for norm is None')    
<         for i in tqdm ( inputs ): 
<             results[i] = apply_async( pool, lazy_one_timep, ( FD, num_lev, num_buf, ring_masks[i], 
<                         internal_state,   bad_frame_list,imgsum, None, cal_error 
<                                      ) )             
<     pool.close()    
<     print( 'Starting running the tasks...')    
<     res =   [ results[k].get() for k in   tqdm( list(sorted(results.keys())) )   ]  
<             
<     #lag_steps  = res[0][1]         
<     g2_G = np.zeros((  int( (num_lev + 1) * num_buf / 2),  len(pixelist)) )  
<     g2_P = np.zeros_like(  g2_G )
<     g2_F = np.zeros_like(  g2_G )         
<     Gmax = 0
<     lag_steps_err  = res[0][1]  
<     #print('Here')
<     for i in inputs:
<         g2_G[:,qind==1+i] = res[i][2]#[:len_lag]   
<         g2_P[:,qind==1+i] = res[i][3]#[:len_lag]  
<         g2_F[:,qind==1+i] = res[i][4]#[:len_lag]            
<     del results
<     del res
<     return g2_G, g2_P, g2_F
694,703c613,639
< def get_g2_from_ROI_GPF( G,P,F,roi_mask):
<     '''YG. 2018.10.26. Get g2 from G, P, F by giving bins (roi_mask)
<     Input:
<         G: <I(t) * I(t+tau)>t
<         P: <  I(t)  >t
<         F: <  I(t+tau) >t
<         roi_mask: the roi mask
<         
<     Output:
<        g2 and g2_err
---
> 
> 
> 
> 
> 
> 
> 
> 
> 
> 
> 
> 
> 
> 
> 
> 
> 
> 
> 
> 
> 
> 
> 
> 
> 
> 
> 
705d640
<     '''
707,746d641
<     qind, pixelist = roi.extract_label_indices(  roi_mask  )   
<     noqs = len(np.unique(qind))    
<     g2 = np.zeros( [G.shape[0], noqs])
<     g2_err = np.zeros( [G.shape[0], noqs])  
<     for i in range(1,1+noqs):
<         ## G[0].shape is the same as roi_mask shape
<         if len( G.shape ) >2:
<             s_Gall_qi = G[:,roi_mask==i]  
<             s_Pall_qi = P[:,roi_mask==i]  
<             s_Fall_qi = F[:,roi_mask==i]  
<         ## G[0].shape is the same length as pixelist
<         else:
<             s_Gall_qi = G[:,qind==i]  
<             s_Pall_qi = P[:,qind==i]  
<             s_Fall_qi = F[:,qind==i]            
<             
<         #print( s_Gall_qi.shape,s_Pall_qi.shape,s_Fall_qi.shape )
<         avgGi = (np.average( s_Gall_qi, axis=1)) 
<         devGi = (np.std( s_Gall_qi, axis=1))
<         avgPi = (np.average( s_Pall_qi, axis=1)) 
<         devPi = (np.std( s_Pall_qi, axis=1))
<         avgFi = (np.average( s_Fall_qi, axis=1)) 
<         devFi = (np.std( s_Fall_qi, axis=1))
<         if len(np.where(avgPi == 0)[0]) != 0:
<             g_max1 = np.where(avgPi == 0)[0][0]
<         else:
<             g_max1 = avgPi.shape[0]    
<         if len(np.where(avgFi == 0)[0]) != 0:
<             g_max2 = np.where(avgFi == 0)[0][0]
<         else:
<             g_max2 = avgFi.shape[0]    
<         g_max = min( g_max1, g_max2)             
<         #print()
<         g2[:g_max,i-1] =   avgGi[:g_max]/( avgPi[:g_max] * avgFi[:g_max] )           
<         g2_err[:g_max,i-1] = np.sqrt( 
<             ( 1/ ( avgFi[:g_max] * avgPi[:g_max] ))**2 * devGi[:g_max] ** 2 +
<             ( avgGi[:g_max]/ ( avgFi[:g_max]**2 * avgPi[:g_max] ))**2 * devFi[:g_max] ** 2 +
<             ( avgGi[:g_max]/ ( avgFi[:g_max] * avgPi[:g_max]**2 ))**2 * devPi[:g_max] ** 2 )
<             
<     return g2, g2_err       
748,749d642
< 
<    

[ambarb]

@afluerasu will check to see if we can just delete this or if anything is of use for tests or better code.

for now, let's ignore p2 and plan to delete.

[afluerasu]

Looks like the normalization is different between the two. (normalizing or not by the total intensity in each frame - imgsum). Both schemes are valid and useful sometimes and we should definitely think about keeping a flexible way of doing different normalizations of the correlation functions is needed. Keeping, however, correlation and correlation2 as two almost identical copies of the code is certainly not the best implementation ... :) I'll look at exactly what we are doing and try to propose a flexible solutions..

else:
< S = norm.shape < if len(S)>1: < fra_pix[ pxlist] = v[w]/ norm[i,pxlist] #-1.0
< else:
< fra_pix[ pxlist] = v[w]/ norm[pxlist] #-1.0

            else: 
                fra_pix[ pxlist] = v[w]/ norm[pxlist]   #-1.0  

170,174c168,169 < S = norm.shape < if len(S)>1:
< fra_pix[ pxlist] = v[w]/ imgsum[i]/ norm[i,pxlist] < else:
< fra_pix[ pxlist] = v[w]/ imgsum[i]/ norm[pxlist]

                fra_pix[ pxlist] = v[w]/ imgsum[i]/  norm[pxlist]