XPCS Pipeline

"This notebook corresponds to version {{ version }} of the pipeline tool: https://github.com/NSLS-II/pipelines"

NSLS2 data retrieval imports



In [3]:

    
from databroker import DataBroker as db, get_images, get_table, get_events
from filestore.api import register_handler, deregister_handler
from filestore.retrieve import _h_registry, _HANDLER_CACHE



In [4]:

    
hdr = db[{{ uid }}]









    



---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-4-27d71802241f> in <module>()
----> 1 hdr = db[{{ uid }}]

NameError: name 'uid' is not defined



In [5]:

    
cd /home/yuzhang/chx-pipelines/Develops/









    



/home/yuzhang/chx-pipelines/Develops



In [6]:

    
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm
import time
from ipywidgets import interact



In [ ]:

    
cd /XF11ID/analysis/Analysis_Pipelines/Develop/



In [7]:

    
%run develop.py









    



/home/yuzhang/.conda/envs/user_analysis/lib/python3.4/site-packages/IPython/html.py:14: ShimWarning: The `IPython.html` package has been deprecated. You should import from `notebook` instead. `IPython.html.widgets` has moved to `ipywidgets`.
  "`IPython.html.widgets` has moved to `ipywidgets`.", ShimWarning)



In [8]:

    
%matplotlib notebook
#%matplotlib inline

Users put

uid for Bluesky Scan
filename for acquiring data directly by EigerSofteare



In [9]:

    
BlueScan = True
DirectAcq = False
detector = 'eiger_4M_cam_img_image_lightfield'  #for 4M

Users put uid here



In [10]:

    
if BlueScan:
    uid = '54614d43'
    #uid = '95782687'
    uid = '95782687'
    uid= 'ff9f20c0'
    uid='71720966'
    uid='1663d34a'
    uid = 'f505e052-3baa-47d4-bdc4-61c2eb1bcc7a'  #sid= 551, 1%PEG, 
    uid='ee6975a1-9161'   #1% wt PEG
    
else:
    uid = '/XF11ID/data/2015/11/23/d01ab510-3cf3-4719-bee3_795_master.h5'

Get data



In [16]:

    
if BlueScan:
    hdr = db[uid]
    ev, = get_events(  hdr, [detector] )
    imgs = ev['data'][detector]
else:    
    imgs =  Images(uid)
print (imgs)
Nimg=len(imgs)









    



hdf5 path = /XF11ID/data/2015/11/23/5f7b93e2-71ef-41b0-9c2c_138_master.h5
<Frames>
Length: 10000 frames
Frame Shape: 2167 x 2070
Pixel Datatype: uint16

Get data path



In [17]:

    
if BlueScan:
    from datetime import datetime
    dt = datetime.fromtimestamp(hdr['start'].time)
    path ='/XF11ID/analysis' + '/%s/%s/%s/' % (dt.year, dt.month, dt.day)
else:
    path ='/XF11ID/analysis/2015/11/23/' 
path









    Out[17]:





'/XF11ID/analysis/2015/11/23/'

Note: experiment information

The physical size of the pixels
Wavelegth of the X-rays - (units in Angstroms)
Detector to sample distance
Exposure time - (units in seconds)
acqusition period - (units in seconds)
dead time - (units in seconds)
time per frame = (exposure time + dead_time or acqusition period) - (units in seconds)



In [18]:

    
imgs.md









    Out[18]:





{'beam_center_x': 852.0,
 'beam_center_y': 1830.0,
 'count_time': 0.0049899998,
 'detector_distance': 4.8400002,
 'frame_time': 0.0049999999,
 'framerate': 200.00000447034847,
 'incident_wavelength': 1.3776,
 'pixel_mask': array([[1, 1, 1, ..., 1, 1, 0],
        [1, 1, 1, ..., 1, 1, 1],
        [1, 1, 1, ..., 1, 1, 1],
        ..., 
        [1, 1, 1, ..., 1, 1, 1],
        [1, 1, 1, ..., 1, 1, 1],
        [1, 1, 1, ..., 1, 1, 1]], dtype=uint32),
 'x_pixel_size': 7.5000004e-05,
 'y_pixel_size': 7.5000004e-05}



In [19]:

    
# The physical size of the pixels
dpix = imgs.md['x_pixel_size'] * 1000.  
lambda_ = imgs.md['incident_wavelength']    # wavelegth of the X-rays in Angstroms
Ldet = 4812.        # detector to sample distance (mm)

exposuretime= imgs.md['count_time']
acquisition_period = imgs.md['frame_time']

# deadtime= 0   # 60e-6 
# timeperframe = exposuretime + deadtime
timeperframe = acquisition_period  

timeperframe, exposuretime









    Out[19]:





(0.0049999999, 0.0049899998)

load a mask

load the mask saved in the mask_pipeline



In [20]:

    
mask = np.load( path +  str(uid)+ "_mask.npy")

Reverse the mask in y-direction due to the coordination difference between python and Eiger software



In [21]:

    
maskr = mask[::-1,:]

Plot the mask



In [22]:

    
fig, ax = plt.subplots()
im=ax.imshow(maskr, origin='lower' ,vmin=0,vmax=1,cmap='viridis')
fig.colorbar(im)
plt.show()

Interactive way to browse through images.

Note : Provide the number of images that you want to browse



In [23]:

    
def view_image(i):    
    fig, ax = plt.subplots()
    ax.imshow(imgs[i]*mask, interpolation='nearest', cmap='viridis',
                  origin='lower', norm= LogNorm(vmin=0.001, vmax=1e1) )
    ax.set_title("Browse the Image Stack")
    plt.show()



In [24]:

    
#interact(view_image, i=(0, Nimg-1))

Movie



In [25]:

    
def view_image(sleeps=1, ims=0, ime = 1):    
    fig, ax = plt.subplots()  
    for i in range( ims, ime  ):
        im=ax.imshow(imgs[i]*mask,  interpolation='nearest', cmap='viridis',
                  origin='lower', norm= LogNorm( vmin=0.01, vmax=10 ) )
        ax.set_title("images_%s"%i)
        
        time.sleep( sleeps )
        plt.draw()
    #fig.colorbar(im)
        
#view_image(.2, 0, 2)

hey, let's see if any images are bad!

load the image intensity (kymograph) saved in the mask_pipeline



In [26]:

    
kymo_sum = np.load( path +  str(uid)+"_kymo_sum.npy" )



In [27]:

    
fig, axes = plt.subplots(  )
axes.plot( kymo_sum, '-go'  ) 
ax.set_ylabel('Intensity')
ax.set_xlabel('Frame')
ax.set_title('Kymograph_sum') 
plt.show()

Get the Averaged Image Data

load the average intensity saved in the mask_pipeline



In [28]:

    
avg_img = np.load( path + str(uid)+"_avg_img.npy" )
avg_imgm =  avg_img * mask

Reverse the image in y-direction due to the coordination difference between python and Eiger software



In [29]:

    
avg_imgr  = avg_img[::-1,:] 
avg_imgmr  = avg_imgm[::-1,:]

Plot the averged image with the mask



In [30]:

    
fig, ax = plt.subplots()
im = ax.imshow(avg_imgmr, cmap='viridis',origin='lower',
               norm= LogNorm(vmin=0.001, vmax=1e1))
ax.set_title("Masked Averaged Image_Reversed")
fig.colorbar(im)
plt.show()

Import all the required packages for Data Analysis

scikit-xray - data analysis tools for X-ray science
- https://github.com/scikit-xray/scikit-xray
xray-vision - plotting helper functions for X-ray science
- https://github.com/Nikea/xray-vision

Get the approximate center and see the statistic to make sure



In [31]:

    
imgs.md['beam_center_x'], imgs.md['beam_center_y']









    Out[31]:





(852.0, 1830.0)



In [32]:

    
#center = (imgs.md['beam_center_x'], imgs.md['beam_center_y'])
center = [ 2167 - 336, 849]  #for not reversed
center = [ 336, 849]  #for reversed
center = [ 2167- 1830, 846]
 
center









    Out[32]:





[337, 846]

check the center with a ring

to be done-->> a better beam center finding algorithm



In [33]:

    
fig, ax = plt.subplots()
im = ax.imshow(avg_imgr, cmap='viridis',origin='lower', norm= LogNorm(vmin=0.001, vmax=1e1))

radius = 54
circle=plt.Circle( [center[1], center[0]], radius, color='b', alpha=1.0, lw=2, edgecolor='r',fill=False)
plt.gcf().gca().add_artist(circle)


ax.set_title("Masked Averaged Image_Reversed")
fig.colorbar(im)


rwidth = 100 
x1,x2 = [center[1] - rwidth, center[1] + rwidth]
y1,y2 = [center[0] - rwidth, center[0] + rwidth]
ax.set_xlim( [x1,x2])
ax.set_ylim( [y1,y2])

plt.show()









    














    











    



/home/yuzhang/.conda/envs/user_analysis/lib/python3.4/site-packages/matplotlib/patches.py:107: UserWarning: Setting the 'color' property will overridethe edgecolor or facecolor properties. 
  warnings.warn("Setting the 'color' property will override"

Circular Average : compute the radial integartion from the center of the speckle pattern



In [37]:

    
bin_centers, ring_averages= circular_average(avg_imgr,  center, pixel_size=(dpix, dpix), mask= maskr)

#  convert to q (reciprocal space)
two_theta = utils.radius_to_twotheta(Ldet, bin_centers)
q_val = utils.twotheta_to_q(two_theta, lambda_)



In [38]:

    
fig,axes = plt.subplots(figsize=(8, 6))
axes.semilogy(q_val, ring_averages, '-o')
axes.set_title('Circular Average')
axes.set_ylabel('Ring Avearge')
axes.set_xlabel('Q ('r'$\AA^{-1}$ )')
axes.set_xlim(0.001, 0.02)
axes.set_ylim(0.001, 10.0)
plt.show()



In [39]:

    
fig,axes = plt.subplots(figsize=(8, 6))
axes.semilogy(bin_centers/dpix, ring_averages, '-o')
axes.set_title('Circular Average')
axes.set_ylabel('Ring Avearge')
axes.set_xlabel('Bin Centers, (pixel)')
axes.set_xlim(30,  250)
axes.set_ylim(0.001, 10.0)
plt.show()

Create label array (Q rings)



In [44]:

    
inner_radius = 58 # radius of the first ring
width = 2       # width of each ring
spacing =  (166 - 58)/9 - 2    # spacing between rings
num_rings = 6   # number of rings

#  find the edges of the required rings
edges = roi.ring_edges(inner_radius, width, spacing, num_rings)
edges









    Out[44]:





array([[  58.,   60.],
       [  70.,   72.],
       [  82.,   84.],
       [  94.,   96.],
       [ 106.,  108.],
       [ 118.,  120.]])



In [45]:

    
two_theta = utils.radius_to_twotheta(Ldet, edges*dpix)
q_ring_val = utils.twotheta_to_q(two_theta, lambda_)

q_ring_center = np.average(q_ring_val, axis=1)
q_ring_center









    Out[45]:





array([ 0.00419415,  0.0050472 ,  0.00590025,  0.00675329,  0.00760634,
        0.00845939])



In [46]:

    
rings = roi.rings(edges, center, avg_imgmr.shape)

ring_mask = rings*maskr

Extract the labeled array



In [47]:

    
ring_mask









    Out[47]:





array([[0, 0, 0, ..., 0, 0, 0],
       [0, 0, 0, ..., 0, 0, 0],
       [0, 0, 0, ..., 0, 0, 0],
       ..., 
       [0, 0, 0, ..., 0, 0, 0],
       [0, 0, 0, ..., 0, 0, 0],
       [0, 0, 0, ..., 0, 0, 0]])



In [48]:

    
qind, pixelist = roi.extract_label_indices(   ring_mask  )
noqs = len( np.unique(qind) )
nopr = np.bincount(qind, minlength=(noqs+1))[1:]



In [49]:

    
qind









    Out[49]:





array([6, 6, 6, ..., 6, 6, 6])



In [50]:

    
nopr









    Out[50]:





array([ 128,  312,  459,  583,  848, 1068])

Number of pixels in each q ring

check center



In [51]:

    
pixel = roi.roi_pixel_values(avg_imgmr, ring_mask, [2] )
fig,ax=plt.subplots()
ax.plot( pixel[0][0] ,'bo', ls='-' )









    














    











    Out[51]:





[<matplotlib.lines.Line2D at 0x7fe28c307d68>]



In [52]:

    
fig, axes = plt.subplots( figsize=(8, 6))
#axes.semilogy(q_val, ring_averages, '-o')
axes.plot(q_val, ring_averages, '-o')
axes.set_title('Circular Average with the Q ring values')
axes.set_ylabel('Ring Avearge')
axes.set_xlabel('Bin Centers 'r'$\AA^{-1}$')
axes.set_xlim(0.00, 0.02)
axes.set_ylim(0, 6)
for i in range(num_rings):
    axes.axvline(q_ring_center[i])
plt.show()



In [53]:

    
plt.close('all')



In [55]:

    
# plot the figure
fig, axes = plt.subplots(figsize=(8,8))
axes.set_title("Labeled Array on Averaged Data")
im,im_label = show_label_array_on_image(axes, avg_imgmr, ring_mask, imshow_cmap='viridis',
                        cmap='Paired',
                         vmin=0.01, vmax=5,  origin="lower")
#rwidth = 200 
#x1,x2 = [center[1] - rwidth, center[1] + rwidth]
#y1,y2 = [center[0] - rwidth, center[0] + rwidth]
#axes.set_xlim( [x1,x2])
#axes.set_ylim( [y1,y2])

#fig.colorbar(im)
rwidth = 200 

x1,x2 = [center[1] - rwidth, center[1] + rwidth]
y1,y2 = [center[0] - rwidth, center[0] + rwidth]
axes.set_xlim( [x1,x2])
axes.set_ylim( [y1,y2])
fig.colorbar(im_label)
plt.show()

Kymograph(waterfall plot) of the max-intensity ring



In [56]:

    
imgs_ =imgs
imgsr = Reverse_Coordinate(imgs_, mask)



In [57]:

    
%run two_time.py

users put the number of ring with max intensity here



In [58]:

    
max_inten_ring =2



In [59]:

    
#kymo = roi.kymograph(imgsr, ring_mask, num = max_inten_ring)



In [60]:

    
t0 = time.time()
data_pixel =   Get_Pixel_Array( imgsr, pixelist).get_data()
run_time(t0)









    



Total time: 4.07 min



In [61]:

    
pixelist_qi =  np.where( qind == 2)[0]         
data_pixel_qi = data_pixel[:,pixelist_qi]



In [62]:

    
fig, ax = plt.subplots(figsize=(8,6))
ax.set_ylabel('Pixel')
ax.set_xlabel('Frame')
ax.set_title('Kymograph')

im = ax.imshow(data_pixel_qi.T, cmap='viridis', vmax=5.0)
fig.colorbar( im   )
ax.set_aspect(30.)
plt.show()

Mean intensities for each ring



In [63]:

    
mean_inten = get_mean_intensity( data_pixel, qind)



In [64]:

    
times = np.arange(  mean_inten[1].shape[0]   )  #*timeperframe  # get the time for each frame

fig, ax = plt.subplots(figsize=(8, 6))
ax.set_title("Mean intensity of each ring")
for i in range(num_rings):
    ax.plot(times, mean_inten[i+1], '--o', label="Ring "+str(i+1))
    ax.set_xlabel("Frame")
    ax.set_ylabel("Mean Intensity")
ax.legend() 
plt.show()

One time Correlation

Note : Enter the number of levels and number of buffers for Muliti tau one time correlation number of buffers has to be even. More details in https://github.com/scikit-xray/scikit-xray/blob/master/skxray/core/correlation.py



In [65]:

    
nopr









    Out[65]:





array([ 128,  312,  459,  583,  848, 1068])



In [66]:

    
edges









    Out[66]:





array([[  58.,   60.],
       [  70.,   72.],
       [  82.,   84.],
       [  94.,   96.],
       [ 106.,  108.],
       [ 118.,  120.]])



In [67]:

    
qind









    Out[67]:





array([6, 6, 6, ..., 6, 6, 6])



In [ ]:



In [ ]:



In [ ]:



In [68]:

    
ring_mask = np.array( ring_mask, dtype=int)



In [69]:

    
good_start = 0
good_end = 7600
imgs_ =imgs[good_start: good_end-1]
imgsr2 = Reverse_Coordinate(imgs_, mask)



In [70]:

    
num_lev = 15
num_buf = 8

t0 = time.time()
g2, lag_steps = corr.multi_tau_auto_corr(num_lev, num_buf,ring_mask, imgsr2)
run_time(t0)









    



Total time: 1.61 min



In [71]:

    
noframes = good_end - good_start
num_buf = #500
num_lev = int(np.log( noframes/(num_buf-1))/np.log(2) +1) +1

t0 = time.time()
g2, lag_steps = corr.multi_tau_auto_corr(num_lev, num_buf,ring_mask, imgsr2)
run_time(t0)









    



  File "<ipython-input-71-d48ecbb64399>", line 2
    num_buf = #500
                  ^
SyntaxError: invalid syntax



In [72]:

    
lags = lag_steps*timeperframe



In [73]:

    
tg2 = np.hstack( [ lags.reshape( len(lags),1), g2] )
#np.save( path + 'g2_%s-%s--%s'%(uid,good_start, good_end), tg2)
np.savetxt( path + 'g2_%s-%s--%s.txt'%(uid,good_start, good_end), tg2)

Plot the one time correlation functions



In [74]:

    
sx = int(round(np.sqrt(num_rings)) )
if num_rings%sx == 0: 
    sy = int(num_rings/sx)
else:
    sy=int(num_rings/sx+1)
#fig = plt.figure(figsize=(14, 10))
fig = plt.figure()
plt.title('uid= %s'%uid,fontsize=20, y =1.08)  
plt.axis('off')
for sn in range(num_rings):
    ax = fig.add_subplot(sx,sy,sn+1 )
    ax.set_ylabel("g2") 
    ax.set_title(" Q= " + '%.5f  '%(q_ring_center[sn]) + r'$\AA^{-1}$')
    y=g2[:, sn]
    ax.semilogx(lags, y, '-o', markersize=6) 
    #ax.set_ylim([min(y)*.95, max(y[1:])*1.05 ])
    ax.set_ylim([    min(y) , max(y[1:]) ])
    ax.set_xlim([   min(lags)+0*1e-6, max(lags)])
plt.show()
fig.tight_layout()



In [76]:

    
FILENAME = 'ee6975a1-9161'
outDir =  '/home/yuzhang/XPCS_Anlysis/Results/'
noframes = 7600

g2_yg = cpopen( 'sid_%s-g2_%sframes'%(FILENAME,noframes),outDir)

sx = int(round(np.sqrt(num_rings)) )
if num_rings%sx == 0: 
    sy = int(num_rings/sx)
else:
    sy=int(num_rings/sx+1)
#fig = plt.figure(figsize=(14, 10))
fig = plt.figure()
plt.title('uid= %s'%uid,fontsize=20, y =1.08)  
plt.axis('off')
for sn in range(num_rings):
    ax = fig.add_subplot(sx,sy,sn+1 )
    ax.set_ylabel("g2") 
    ax.set_title(" Q= " + '%.5f  '%(q_ring_center[sn]) + r'$\AA^{-1}$')
    y=g2[:, sn]
    ax.semilogx(lags, y, 'ro', markersize=6) 
    
    y2=g2_yg[:, sn]
    ax.semilogx(lags, y2, '--b', markersize=6) 
    
    #ax.set_ylim([min(y)*.95, max(y[1:])*1.05 ])
    ax.set_ylim([    min(y) , max(y[1:]) ])
    ax.set_xlim([   min(lags)+0*1e-6, max(lags)])
plt.show()
fig.tight_layout()

Fit g2



In [77]:

    
from lmfit import  Model
mod = Model(corr.auto_corr_scat_factor)



In [78]:

    
rate = []

sx = int( round (np.sqrt(num_rings)) )
if num_rings%sx==0:
    sy = int(num_rings/sx)
else:
    sy = int(num_rings/sx+1)
    
#fig = plt.figure(figsize=(14, 10))
fig = plt.figure()
plt.title('uid= %s'%uid, fontsize=20, y =1.02)  
for sn in range(num_rings):
    ax = fig.add_subplot(sx, sy, sn+1 )
    y=g2[1:, sn]
    result1 = mod.fit(y, lags=lags[1:], beta=.1,
                      relaxation_rate =.5, baseline=1.0)
    rate.append(result1.best_values['relaxation_rate'])
    
    ax.semilogx(lags[1:], y, 'ro')
    ax.semilogx(lags[1:], result1.best_fit, '-b')
    ax.set_title(" Q= " + '%.5f  '%(q_ring_center[sn]) + r'$\AA^{-1}$')  
    ax.set_ylim([min(y)*.95, max(y[1:]) *1.05])
    txts = r'$\gamma$' + r'$ = %.3f$'%(rate[sn]) +  r'$ s^{-1}$'
    ax.text(x =0.015, y=.55, s=txts, fontsize=14, transform=ax.transAxes)              
 
fig.tight_layout()

Plot the relaxation rates vs (q_ring_center)**2



In [79]:

    
fig, ax=plt.subplots()
ax.plot(q_ring_center**2, rate, 'ro', ls='--')
ax.set_ylabel('Relaxation rate 'r'$\gamma$'"($s^{-1}$)")
ax.set_xlabel("$q^2$"r'($\AA^{-2}$)')
plt.show()

Fitted the Diffusion Coefficinet D0



In [80]:

    
D0 = np.polyfit(q_ring_center**2, rate, 1)
gmfit = np.poly1d(D0)
print ('The fitted diffusion coefficient D0 is:  %.2E   A^2S-1'%D0[0])









    



The fitted diffusion coefficient D0 is:  7.99E+03   A^2S-1



In [81]:

    
fig,ax = plt.subplots()
ax.plot(q_ring_center**2, rate, 'ro', ls='')
ax.plot(q_ring_center**2,  gmfit(q_ring_center**2),  ls='-')
ax.set_ylabel('Relaxation rate 'r'$\gamma$'"($s^{-1}$)")
ax.set_xlabel("$q^2$"r'($\AA^{-2}$)')
plt.show()

Two_time Correlation



In [82]:

    
%run two_time.py



In [83]:

    
good_start= 0
good_end = 10000
imgs_ =imgs[good_start: good_end]
imgsr = Reverse_Coordinate(imgs_, mask)



In [84]:

    
g12b = auto_two_Array( imgsr, ring_mask, data_pixel = data_pixel )









    



Total time: 0.34 min



In [ ]:



In [85]:

    
g12_num = 0  #0: the firs ring
data = g12b[:,:,g12_num]
fig, ax = plt.subplots()
im=ax.imshow( data, origin='lower' , cmap='viridis', 
             norm= LogNorm( vmin= 1, vmax= 1.6 ), 
        extent=[0, data.shape[0]*timeperframe, 0, data.shape[0]*timeperframe ] )
ax.set_title('0-%s frames--Qth= %s'%(Nimg,g12_num))
ax.set_xlabel( r'$t_1$ $(s)$', fontsize = 18)
ax.set_ylabel( r'$t_2$ $(s)$', fontsize = 18)
fig.colorbar(im)
plt.show()

find the bad frames and mask it in two-time correlation



In [ ]:

    
#np.where( g12b[:,:,0] == g12b[:,:,0].min() )



In [86]:

    
g12b_mask = make_g12_mask(g12b,[7629])
g12bm = masked_g12(g12b,[7629])

get one time correlation



In [87]:

    
g2b = get_one_time_from_two_time(g12bm, timeperframe= timeperframe)



In [88]:

    
sx = int(round(np.sqrt(num_rings)) )
if num_rings%sx == 0: 
    sy = int(num_rings/sx)
else:
    sy=int(num_rings/sx+1)
#fig = plt.figure(figsize=(14, 10))
fig = plt.figure()
plt.title('uid= %s'%uid,fontsize=20, y =1.08)  
plt.axis('off')
for sn in range(num_rings):
    ax = fig.add_subplot(sx,sy,sn+1 )
    ax.set_ylabel("g2") 
    ax.set_title(" Q= " + '%.5f  '%(q_ring_center[sn]) + r'$\AA^{-1}$')     
    y=g2b[:, sn]
    ax.semilogx( np.arange(len(y))*timeperframe, y, '-o', markersize=6) 
    #ax.semilogx(lags, y, '-o', markersize=6) 
    #ax.set_ylim([min(y)*.95, max(y[1:])*1.05 ])
    ax.set_ylim( [1.0, max(y[1:])*1.05 ] )
    
plt.show()
fig.tight_layout()

plot g2 by multi-tau and g2 from two-time



In [89]:

    
sx = int(round(np.sqrt(num_rings)) )
if num_rings%sx == 0: 
    sy = int(num_rings/sx)
else:
    sy=int(num_rings/sx+1)
#fig = plt.figure(figsize=(14, 10))
fig = plt.figure()
plt.title('uid= %s'%uid,fontsize=20, y =1.06)  
plt.axis('off')
for sn in range(num_rings):
    ax = fig.add_subplot(sx,sy,sn+1 )
    ax.set_ylabel("g2") 
    ax.set_title(" Q= " + '%.5f  '%(q_ring_center[sn]) + r'$\AA^{-1}$')     
    y=g2b[:, sn]
    ax.semilogx( np.arange(len(y))*timeperframe, y, '--r', markersize=6) 
    
    y2=g2[:, sn]
    ax.semilogx(lags, y2, 'o', markersize=6) 
    
    #ax.semilogx(lags, y, '-o', markersize=6) 
    #ax.set_ylim([min(y)*.95, max(y[1:])*1.05 ])
    ax.set_ylim( [1.0, max(y2[1:])*1.05 ] )
plt.show()
fig.tight_layout()

work one bad images in multi-tau code



In [90]:

    
good_start= 0
good_end = 500
imgs_ =imgs[good_start: good_end]
imgsr = Reverse_Coordinate(imgs_, mask)
noframes = good_end - good_start

num_buf = 8
num_lev = int(np.log( noframes/(num_buf-1))/np.log(2) +1) +1


t0 = time.time()
g2, lag_steps = corr.multi_tau_auto_corr(num_lev, num_buf,ring_mask, imgsr)
run_time(t0)

lags = lag_steps*timeperframe









    



Total time: 0.11 min



In [94]:

    
#g2_t, lag_steps_t =autocor_one_time( num_buf,  ring_mask, imgsr,       num_lev=None, bad_images=None  )



In [393]:

    
pwd









    Out[393]:





'/home/yuzhang/chx-pipelines/Develops'



In [585]:

    
%run two_time.py



In [586]:

    
good_start= 0
good_end = 100
imgs_ =imgs[good_start: good_end]
imgsr = Reverse_Coordinate(imgs_, mask)
noframes = good_end - good_start

num_buf = 4



In [587]:

    
g2_t2, lag_steps_t =autocor_one_time( num_buf,  ring_mask, imgsr, 
                                    num_lev=None, bad_images=[ 94,95,96,97,98,99,100]  )









    



The lev number is 7
Doing g2 caculation of 100 frames---
##########bad image: 94 here!
bad image: 95 here!
bad image: 96 here!
bad image: 97 here!
bad image: 98 here!
bad image: 99 here!
Total time: 0.03 min



In [588]:

    
good_start= 0
good_end = 95
imgs_ =imgs[good_start: good_end]
imgsr = Reverse_Coordinate(imgs_, mask)
noframes = good_end - good_start

num_buf = 4



In [589]:

    
g2_t3, lag_steps_t3 =autocor_one_time( num_buf,  ring_mask, imgsr, 
                                    num_lev=None, bad_images=[  ]  )









    



The lev number is 6
Doing g2 caculation of 95 frames---
#####Total time: 0.03 min



In [590]:

    
g2_t2[:,0]









    Out[590]:





array([ 1.78143577,  1.02329077,  1.0171744 ,  1.03811612,  1.02307037,
        1.02151458,  1.02397441,  1.0215932 ,  1.02466236,  1.02632637,
        1.02589129,  1.02554867])



In [591]:

    
g2_t3[:,0]









    Out[591]:





array([ 1.78238823,  1.02458214,  1.01637848,  1.03652324,  1.02307037,
        1.02151458,  1.02397441,  1.0215932 ,  1.02466236,  1.02632637,
        1.02589129,  1.02554867])



In [602]:

    
%run two_time.py



In [603]:

    
good_start= 0
good_end = 8
imgs_ =imgs[good_start: good_end]
imgsr = Reverse_Coordinate(imgs_, mask)
noframes = good_end - good_start

num_buf = 4



In [604]:

    
g2_t2, lag_steps_t =autocor_one_time( num_buf,  ring_mask, imgsr, 
                                    num_lev=None, bad_images=[ 1, 10, 20, 45  ]  )









    



The lev number is 3
Doing g2 caculation of 8 frames---
#bad image: 1 here!
{1: [1, 2, 1, 1], 2: [1, 1], 3: [0, 0]}
Total time: 0.00 min



In [595]:

    
g2_t2[:,0]









    Out[595]:





array([ 1.78711109,  1.0612573 ,  1.04040582,  1.00682279,  1.01833726])



In [575]:

    
lag_steps_t









    Out[575]:





array([  0.,   1.,   2.,   3.,   4.,   6.,   8.,  12.,  16.,  24.,  32.])



In [ ]:



In [566]:

    
g2_t[:,0]









    Out[566]:





array([ 1.76746229,  1.07993428,  1.0665884 ,  1.07064173,  1.05660431,
        1.05903657,  1.05855856,  1.05027663,  1.04623448,  1.04321904,
        1.03755408,  1.03741158,  1.03787578,  1.03441453,  1.03325287,
        1.03169963,  1.03257515,  1.03067179,  1.03077536,  1.03145597,
        1.03129382,  1.03313767,  1.03096145,  1.02980621,  1.02984819,
        1.0261275 ,  1.02480121,  1.02520176,  1.02300024,  1.02342068,
        1.02289539])



In [154]:

    
dly, dict_dly = delays( num_lev =3, num_buf=4, time=1 )



In [418]:

    
dly









    Out[418]:





array([  0.,   1.,   2.,   3.,   4.,   6.,   8.,  12.,  16.])



In [493]:

    
dict_dly









    Out[493]:





{1: array([ 3.,  3.,  3.,  3.]), 2: array([ 4.,  6.]), 3: array([  8.,  12.])}



In [495]:

    
{ key: [0]* len(  dict_dly[key] ) for key in list(dict_dly.keys())  }









    Out[495]:





{1: [0, 0, 0, 0], 2: [0, 0], 3: [0, 0]}



In [ ]:



In [420]:

    
lev_leng = [  len(  dict_dly[i] ) for i in list(dict_dly.keys())   ]



In [421]:

    
lev_leng









    Out[421]:





[4, 2, 2]



In [422]:

    
mm = dict_dly



In [424]:

    
mm[1][:] = 3



In [425]:

    
mm









    Out[425]:





{1: array([ 3.,  3.,  3.,  3.]), 2: array([ 4.,  6.]), 3: array([  8.,  12.])}



In [ ]:



In [ ]:



In [ ]:



In [323]:

    
g2_t2[:,0]









    Out[323]:





array([ 1.98867232,  1.00308826,  1.02341422,  0.95997843,  1.00161315,
        1.0143698 ])



In [ ]:



In [ ]:



In [ ]:



In [310]:

    
x = (np.ravel(imgsr[0]))[pixelist]









    



---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-310-d0ad802d5630> in <module>()
----> 1 x = (np.ravel(imgsr[0]))[pixelist]

/home/yuzhang/chx-pipelines/Develops/two_time.py in __getitem__(self, key)
    408         if self.mask is not None:
    409             img =self.indexable[key] * self.mask
--> 410         else:
    411             img = self.indexable[key]
    412 

NameError: name 'len' is not defined



In [177]:

    
x=np.linspace(0,10,11)



In [178]:

    
y1=   np.ma.zeros(x.shape)
y1.mask = True



In [184]:

    
y2=   np.ma.zeros(x.shape)
y2.mask = True



In [180]:

    
x









    Out[180]:





array([  0.,   1.,   2.,   3.,   4.,   5.,   6.,   7.,   8.,   9.,  10.])



In [181]:

    
y1









    Out[181]:





masked_array(data = [-- -- -- -- -- -- -- -- -- -- --],
             mask = [ True  True  True  True  True  True  True  True  True  True  True],
       fill_value = 1e+20)



In [182]:

    
(x +y1)/2.









    Out[182]:





masked_array(data = [-- -- -- -- -- -- -- -- -- -- --],
             mask = [ True  True  True  True  True  True  True  True  True  True  True],
       fill_value = 1e+20)



In [185]:

    
((y2 +y1)/2. ).data









    Out[185]:





array([ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.])



In [ ]:



In [345]:

    
(y1 == np.ma.zeros(x.shape)).all()









    Out[345]:





masked



In [355]:

    
((x+y1).data ==0).all()









    Out[355]:





True



In [356]:

    
(x.data ==0)









    Out[356]:





False



In [360]:

    
if ((x+y1).data ==0):
    print ('y')









    



---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-360-4b8667f1b361> in <module>()
----> 1 if ((x+y1).data ==0):
      2     print ('y')

ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()



In [ ]:



In [ ]:



In [285]:

    
x+y1









    Out[285]:





masked_array(data = [-- -- -- -- -- -- -- -- -- -- --],
             mask = [ True  True  True  True  True  True  True  True  True  True  True],
       fill_value = 1e+20)



In [286]:

    
x









    Out[286]:





array([  0.,   1.,   2.,   3.,   4.,   5.,   6.,   7.,   8.,   9.,  10.])



In [297]:

    
np.bincount(qind,   weights= x )









    Out[297]:





array([   0.,  144.,  254.,  111.,  103.,  111.,   66.])



In [300]:

    
np.bincount(qind,   weights= y1 + 100*x )









    Out[300]:





array([ 0.,  0.,  0.,  0.,  0.,  0.,  0.])



In [ ]:



In [ ]:



In [ ]:



In [295]:

    
np.bincount(np.int_(x),   weights= x )









    Out[295]:





array([   0.,  493.,  204.,   51.,   36.,    5.])



In [296]:

    
np.bincount(np.int_(x),   weights= x+ y1 )









    Out[296]:





array([ 0.,  0.,  0.,  0.,  0.,  0.])



In [ ]:



In [ ]:



In [250]:

    
y1.shape









    Out[250]:





(3398,)



In [266]:

    
yb1









    Out[266]:





array([  0.,   1.,   2.,   3.,   4.,   5.,   6.,   7.,   8.,   9.,  10.])



In [267]:

    
yb2









    Out[267]:





array([  0.,   1.,   2.,   3.,   4.,   5.,   6.,   7.,   8.,   9.,  10.])



In [ ]:



In [ ]:



In [ ]:



In [235]:

    
%run two_time.py



In [ ]:



In [ ]:



In [ ]:

a test for comparison of one and two-time



In [211]:

    
good_start= 0
good_end = 7600
imgs_ =imgs[good_start: good_end]
imgsr = Reverse_Coordinate(imgs_, mask)
noframes = good_end - good_start

num_buf = 8
num_lev = int(np.log( noframes/(num_buf-1))/np.log(2) +1) +1
  


t0 = time.time()
g2, lag_steps = corr.multi_tau_auto_corr(num_lev, num_buf,ring_mask, imgsr)
run_time(t0)

lags = lag_steps*timeperframe









    



Total time: 1.59 min



In [205]:

    
#g12b_norm = auto_two_Array( imgsr, ring_mask, data_pixel = None )









    



Total time: 0.20 min



In [212]:

    
g12b_norm, g12b_not_norm, norms = auto_two_Array_g1_norm( imgsr,
                            ring_mask, data_pixel = None )









    



Total time: 1.76 min



In [213]:

    
g12_num = 0  #0: the firs ring
data = g12b_test[:,:,g12_num]
fig, ax = plt.subplots()
im=ax.imshow( data, origin='lower' , cmap='viridis', 
             norm= LogNorm( vmin= 1, vmax= 1.6 ), )
        #extent=[0, data.shape[0]*timeperframe, 0, data.shape[0]*timeperframe ] )
ax.set_title('0-%s frames--Qth= %s'%(Nimg,g12_num))
ax.set_xlabel( r'$t_1$ $(s)$', fontsize = 18)
ax.set_ylabel( r'$t_2$ $(s)$', fontsize = 18)
fig.colorbar(im)
plt.show()



In [ ]:



In [214]:

    
g2b_norm = get_one_time_from_two_time(g12b_norm, timeperframe= timeperframe)



In [215]:

    
g2b_not_norm = get_one_time_from_two_time(g12b_not_norm,
                             norms=norms, nopr=nopr, timeperframe= timeperframe)



In [226]:

    
#num_rings =6


sx = int(round(np.sqrt(num_rings)) )
if num_rings%sx == 0: 
    sy = int(num_rings/sx)
else:
    sy=int(num_rings/sx+1)
#fig = plt.figure(figsize=(14, 10))
fig = plt.figure()
plt.title('uid= %s'%uid,fontsize=20, y =1.06)  
plt.axis('off')
for sn in range(num_rings):
    ax = fig.add_subplot(sx,sy,sn+1 )
    ax.set_ylabel("g2") 
    ax.set_title(" Q= " + '%.5f  '%(q_ring_center[sn]) + r'$\AA^{-1}$')    
    
    y=g2b_norm[:, sn]
    ax.semilogx( np.arange(len(y))*timeperframe, y, '-bs', markersize=3) 
    
    y3=g2b_not_norm[:, sn]
    ax.semilogx( np.arange(len(y))*timeperframe, y3, '-r', lw=4, markersize=3)     
    
    y2=g2[:, sn]
    ax.semilogx(lags, y2, '-ko', markersize= 3) 
    
    #ax.semilogx(lags, y, '-o', markersize=6) 
    #ax.set_ylim([min(y)*.95, max(y[1:])*1.05 ])
    ax.set_ylim( [1.0, max(y2[1:])*1.05 ] )
plt.show()
fig.tight_layout()



In [ ]:

get fit of g2 from two-time



In [354]:

    
rate_g12 = []


sx = int(round(np.sqrt(num_rings)) )
if num_rings%sx == 0: 
    sy = int(num_rings/sx)
else:
    sy=int(num_rings/sx+1)
#fig = plt.figure(figsize=(14, 10))
fig = plt.figure()
plt.title('uid= %s'%uid,fontsize=20, y =1.06)  
plt.axis('off')
for sn in range(num_rings):
    ax = fig.add_subplot(sx,sy,sn+1 )
    ax.set_ylabel("g2") 
    ax.set_title(" Q= " + '%.5f  '%(q_ring_center[sn]) + r'$\AA^{-1}$')     
    y=g2b[:, sn]
    ax.semilogx( np.arange(len(y))*timeperframe, y, 'bo', markersize=6) 
    
    
    y2= y[1:7600]
    result2 = mod.fit(y2, lags= np.arange(len(y2))*timeperframe, beta=.1,
                      relaxation_rate =.5, baseline=1.0)
    
    rate_g12.append(result2.best_values['relaxation_rate'])
    ax.semilogx( np.arange(len(y2))*timeperframe, result2.best_fit, '-r')
    
 
    ax.set_title(" Q= " + '%.5f  '%(q_ring_center[sn]) + r'$\AA^{-1}$')  
    ax.set_ylim([min(y)*.95, max(y[1:]) *1.05])
    txts = r'$\gamma$' + r'$ = %.3f$'%(rate_g12[sn]) +  r'$ s^{-1}$'
    ax.text(x =0.015, y=.55, s=txts, fontsize=14, transform=ax.transAxes)    
    
    ax.set_ylim( [1.0, max(y[1:])*1.01 ] )
 
fig.tight_layout()



In [377]:

    
D0_g12 = np.polyfit(q_ring_center**2, rate_g12, 1)
gmfit_g12 = np.poly1d(D0_g12)
print ('The fitted diffusion coefficient D0_g12 is:  %.2E   A^2S-1'%D0_g12[0])









    



The fitted diffusion coefficient D0_g12 is:  1.30E+04   A^2S-1



In [378]:

    
fig,ax = plt.subplots()
ax.plot(q_ring_center**2, rate_g12, 'ro', ls='')
ax.plot(q_ring_center**2,  gmfit_g12(q_ring_center**2),  ls='-')
ax.set_ylabel('Relaxation rate 'r'$\gamma$'"($s^{-1}$)")
ax.set_xlabel("$q^2$"r'($\AA^{-2}$)')
plt.show()



In [ ]:



In [ ]:

get one time correlation @different age



In [173]:

    
g12_num =0

g2b = show_g12_age_cuts( g12bm, g12_num=g12_num, slice_num =3, slice_width= 500, 
                slice_start=4000, slice_end= 20000-4000,
    timeperframe= timeperframe,vmin= 1.0, vmax =  1.22) #vmax=g12b[:,:,g12_num].max() )









    



the cut age centers are: [  4000.  10000.  16000.]

get taus histgram



In [230]:

    
taus = show_g12_tau_cuts(g12b_norm, g12_num=0,  slice_num = 5, slice_width= 1.0, 
    slice_start=3, slice_end= 5000-1,draw_scale_tau =1.0, vmin=1.01,vmax=1.55 )









    



the cut tau centers are: [  3.00000000e+00   1.25200000e+03   2.50100000e+03   3.75000000e+03
   4.99900000e+03]



In [232]:

    
his = his_tau(taus, hisbin = 10, plot=True)



In [ ]: