<jupyter_start><jupyter_text>Mentions getting the ISF from single particle tracking data:https://journals.aps.org/prl/pdf/10.1103/PhysRevLett.121.078001<jupyter_code>isf = np.zeros((5,1000),dtype=np.complex128)
isf_std = np.zeros_like(isf)
isf_stderr = np.zeros_like(isf)
num_particles = np.zeros((1000))
q = 2*np.sqrt(2)
count=0

qs = np.array([1.98,3.02,3.96,5.0,6.04])

for j in range(5):
    q = qs[j]*np.sqrt(np.pi)
    
    qx = np.arange(0, q, 0.1)
    qy = np.sqrt(q**2 - qx**2)
    
    for i in range(1000):
        count=0
        temp = []
        for particle in imxA.keys():
            if not np.isnan(imxA[particle].values[i]):
                temp2 = []
                for qi in range(len(qx)):
                    temp2.append(np.exp(-1j*(qx[qi]*imxA[particle].values[i]+qy[qi]*imyA[particle].values[i])))
                    temp2.append(np.exp(-1j*(qx[qi]*imxA[particle].values[i]-qy[qi]*imyA[particle].values[i])))
                    temp2.append(np.exp(1j*(qx[qi]*imxA[particle].values[i]+qy[qi]*imyA[particle].values[i])))
                    temp2.append(np.exp(1j*(qx[qi]*imxA[particle].values[i]-qy[qi]*imyA[particle].values[i])))
                temp.append(np.mean(temp2))
                if j==0:
                    num_particles[i] = num_particles[i]+1
       
        for particle in imxA1.keys():
            if not np.isnan(imxA1[particle].values[i]):
                temp2 = []
                for qi in range(len(qx)):
                    temp2.append(np.exp(-1j*(qx[qi]*imxA1[particle].values[i]+qy[qi]*imyA1[particle].values[i])))
                    temp2.append(np.exp(-1j*(qx[qi]*imxA1[particle].values[i]-qy[qi]*imyA1[particle].values[i])))
                    temp2.append(np.exp(1j*(qx[qi]*imxA1[particle].values[i]+qy[qi]*imyA1[particle].values[i])))
                    temp2.append(np.exp(1j*(qx[qi]*imxA1[particle].values[i]-qy[qi]*imyA1[particle].values[i])))
                temp.append(np.mean(temp2))
                if j==0:
                    num_particles[i] = num_particles[i]+1
                    
        for particle in imxB.keys():
            if not np.isnan(imxB[particle].values[i]):
                temp2 = []
                for qi in range(len(qx)):
                    temp2.append(np.exp(-1j*(qx[qi]*imxB[particle].values[i]+qy[qi]*imyB[particle].values[i])))
                    temp2.append(np.exp(-1j*(qx[qi]*imxB[particle].values[i]-qy[qi]*imyB[particle].values[i])))
                    temp2.append(np.exp(1j*(qx[qi]*imxB[particle].values[i]+qy[qi]*imyB[particle].values[i])))
                    temp2.append(np.exp(1j*(qx[qi]*imxB[particle].values[i]-qy[qi]*imyB[particle].values[i])))
                temp.append(np.mean(temp2))
                if j==0:
                    num_particles[i] = num_particles[i]+1
                    
        for particle in imxB1.keys():
            if not np.isnan(imxB1[particle].values[i]):
                temp2 = []
                for qi in range(len(qx)):
                    temp2.append(np.exp(-1j*(qx[qi]*imxB1[particle].values[i]+qy[qi]*imyB1[particle].values[i])))
                    temp2.append(np.exp(-1j*(qx[qi]*imxB1[particle].values[i]-qy[qi]*imyB1[particle].values[i])))
                    temp2.append(np.exp(1j*(qx[qi]*imxB1[particle].values[i]+qy[qi]*imyB1[particle].values[i])))
                    temp2.append(np.exp(1j*(qx[qi]*imxB1[particle].values[i]-qy[qi]*imyB1[particle].values[i])))
                temp.append(np.mean(temp2))
                if j==0:
                    num_particles[i] = num_particles[i]+1
          
        isf[j,i] = np.mean(temp)
        isf_std[j,i] = np.std(temp)
        isf_stderr[j,i] = np.std(temp)/np.sqrt(len(temp))
        #isf[i] = isf[i]/count
num_particles
with open(r"isf_m04_at5qs_new.pkl",'wb') as output_file:
    pickle.dump(isf, output_file)
with open(r"isf_stderr_m04_at5qs_new.pkl",'wb') as output_file:
    pickle.dump(isf_stderr, output_file)
with open(r"isf_std_m04_at5qs_new.pkl",'wb') as output_file:
    pickle.dump(isf_std, output_file)
with open(r"for_isf_m04_trackpytimes.pkl",'wb') as output_file:
    pickle.dump(imxA.index, output_file)
############################
# QUICKER CALCULATION      #
############################

isf = np.zeros((5,1000),dtype=np.complex128)
isf_std = np.zeros_like(isf)
isf_stderr = np.zeros_like(isf)
q = 2*np.sqrt(2)
count=0

qs = np.array([1.98,3.02,3.96,5.0,6.04])

for j in range(5):
    q = qs[j]*(5./3) #*np.sqrt(np.pi)
    
    qx = np.arange(0, q, 0.6)
    qy = np.sqrt(q**2 - qx**2)
    
    for i in range(90):
        count=0
        temp = []
        for particle in imxA.keys():
            if not np.isnan(imxA[particle].values[i]):
                temp2 = []
                for qi in range(len(qx)):
                    temp2.append(np.exp(-1j*(qx[qi]*imxA[particle].values[i]+qy[qi]*imyA[particle].values[i])))
                    temp2.append(np.exp(-1j*(qx[qi]*imxA[particle].values[i]-qy[qi]*imyA[particle].values[i])))
                    temp2.append(np.exp(1j*(qx[qi]*imxA[particle].values[i]+qy[qi]*imyA[particle].values[i])))
                    temp2.append(np.exp(1j*(qx[qi]*imxA[particle].values[i]-qy[qi]*imyA[particle].values[i])))
                temp.append(np.mean(temp2))
       
        for particle in imxA1.keys():
            if not np.isnan(imxA1[particle].values[i]):
                temp2 = []
                for qi in range(len(qx)):
                    temp2.append(np.exp(-1j*(qx[qi]*imxA1[particle].values[i]+qy[qi]*imyA1[particle].values[i])))
                    temp2.append(np.exp(-1j*(qx[qi]*imxA1[particle].values[i]-qy[qi]*imyA1[particle].values[i])))
                    temp2.append(np.exp(1j*(qx[qi]*imxA1[particle].values[i]+qy[qi]*imyA1[particle].values[i])))
                    temp2.append(np.exp(1j*(qx[qi]*imxA1[particle].values[i]-qy[qi]*imyA1[particle].values[i])))
                temp.append(np.mean(temp2))
                    
        for particle in imxB.keys():
            if not np.isnan(imxB[particle].values[i]):
                temp2 = []
                for qi in range(len(qx)):
                    temp2.append(np.exp(-1j*(qx[qi]*imxB[particle].values[i]+qy[qi]*imyB[particle].values[i])))
                    temp2.append(np.exp(-1j*(qx[qi]*imxB[particle].values[i]-qy[qi]*imyB[particle].values[i])))
                    temp2.append(np.exp(1j*(qx[qi]*imxB[particle].values[i]+qy[qi]*imyB[particle].values[i])))
                    temp2.append(np.exp(1j*(qx[qi]*imxB[particle].values[i]-qy[qi]*imyB[particle].values[i])))
                temp.append(np.mean(temp2))
                    
        for particle in imxB1.keys():
            if not np.isnan(imxB1[particle].values[i]):
                temp2 = []
                for qi in range(len(qx)):
                    temp2.append(np.exp(-1j*(qx[qi]*imxB1[particle].values[i]+qy[qi]*imyB1[particle].values[i])))
                    temp2.append(np.exp(-1j*(qx[qi]*imxB1[particle].values[i]-qy[qi]*imyB1[particle].values[i])))
                    temp2.append(np.exp(1j*(qx[qi]*imxB1[particle].values[i]+qy[qi]*imyB1[particle].values[i])))
                    temp2.append(np.exp(1j*(qx[qi]*imxB1[particle].values[i]-qy[qi]*imyB1[particle].values[i])))
                temp.append(np.mean(temp2))
          
        isf[j,i] = np.mean(temp)
        isf_std[j,i] = np.std(temp)
        isf_stderr[j,i] = np.std(temp)/np.sqrt(len(temp))
        #isf[i] = isf[i]/count
cmap = plt.get_cmap("viridis") 


fig_width = 8
markersize = 6 #was 1
fig, ax = plt.subplots(nrows=1, figsize=(fig_width,fig_width/1.618))
data = 0.25 * (fit04_06_A.isf_data + fit04_06_A1.isf_data + fit04_06_B.isf_data + fit04_06_B1.isf_data)
theory = 0.25 * (fit04_06_A.theory + fit04_06_A1.theory + fit04_06_B.theory + fit04_06_B1.theory)
times = fit04_06_A.lagtime
q = fit04_06_A.q    
xlabel_str = "Lag time (s)"
ylabel_str = "ISF"

#q_index_to_plot = np.arange(12,100,7)
q_index_to_plot = np.array([21,32,42,53,64])
clrs = np.linspace(0,1,len(q_index_to_plot))
for i,qv in enumerate([21,32,42,53,64]):
    plt_color = cmap(clrs[i])
    qv = int(qv)

    ax.semilogx(times, data[:,qv], 'o', ms=markersize, color=plt_color, label="Data for q index %i, q=%.4f" % (qv, q[qv]))
    ax.semilogx(times, theory[:,qv], '-', color=plt_color, lw=2, alpha=0.8, label="Fit for q_index %i" % qv)
    ax.semilogx(times, theory[:,qv], color='k', linestyle=(0, (1, 1)), lw=1)
    
    #ax.semilogx(imxA.index, (isf[i,:]*np.conj(isf[i,:]))**1, 's', ms=4, color=plt_color, mec='k', mew=0.2)
    #ax.errorbar(imxA.index, (isf[i,:]*np.conj(isf[i,:]))**1, yerr=isf_stderr[i,:], fmt='none', color=plt_color, lw=1)
    
    ax.semilogx(imxA.index, np.real(isf[i,:]), 's', ms=4, color=plt_color, mec='k', mew=0.2)
    ax.errorbar(imxA.index, np.real(isf[i,:]), yerr=isf_stderr[i,:], fmt='none', color=plt_color, lw=1)

plt.legend()

#plt.semilogx(imxA.index, isf, 'ro')
#plt.errorbar(imxA.index, isf, yerr=isf_stderr, fmt='none', color='r')<jupyter_output><empty_output><jupyter_text>Getting van Hove distribution<jupyter_code>from scipy.optimize import curve_fit
# Gaussian function
def gaus_func(x, a, x0, sigma):
    return a*np.exp(-(x-x0)**2/(2*sigma**2))

def gaus_exp_func(x, a, x0, sigma, b, lam):
    return a*np.exp(-(x-x0)**2/(2*sigma**2)) + b*np.exp(-1*abs(x-x0)/lam)
bins = np.arange(-2,2,0.015)

posX = tA.set_index(['frame','particle'])['x'].unstack()
posY = tA.set_index(['frame','particle'])['y'].unstack()

lagframes = 10


vhX = tp.vanhove(posX, lagframes, ensemble=True, mpp=0.13, bins=bins)
vhY = tp.vanhove(posY, lagframes, ensemble=True, mpp=0.13, bins=bins)

prob10 = 0.5*(vhX+vhY)

popt10, pcovR = curve_fit(gaus_func, bins[:-1], prob10)
popt10E, pcovRE = curve_fit(gaus_exp_func, bins[:-1], prob10, p0=[2,0,0.1,1,0.1])
popt10
plt.figure()
#plt.semilogy(vhX.index, vhX, 'ro')
#plt.semilogy(vhX.index, vhY, 'bo')
#plt.semilogy(vhX.index, 0.5*(vhX+vhY), 'go')
plt.semilogy(bins[:-1], 0.5*(vhX+vhY), 'k.')

plt.plot(bins[:-1], gaus_func(bins[:-1], popt10[0], popt10[1], popt10[2]), '-', lw=3, alpha=0.6, c='tab:purple')
plt.plot(bins[:-1], gaus_exp_func(bins[:-1], popt10E[0], popt10E[1], popt10E[2], popt10E[3], popt10E[4]), '-', lw=3, alpha=0.6, c='tab:pink')
bins = np.arange(-2,2,0.015)
cmap = plt.get_cmap("viridis")
lag_times_in_frames = np.array([2,6,24,120])
clrs = np.linspace(0,1,len(lag_times_in_frames))
probs = np.zeros((len(lag_times_in_frames), len(bins)-1))
gauss_fit = np.zeros_like(probs)
gauss_width = np.zeros_like(lag_times_in_frames)
plt.figure()
for i,lagframes in enumerate(lag_times_in_frames):
    print(i)
    
    posXA = tA.set_index(['frame','particle'])['x'].unstack()
    posYA = tA.set_index(['frame','particle'])['y'].unstack()
    vhXA = tp.vanhove(posXA, int(lagframes), ensemble=True, mpp=0.13, bins=bins)
    vhYA = tp.vanhove(posYA, int(lagframes), ensemble=True, mpp=0.13, bins=bins)
    
    posXA1 = tA1.set_index(['frame','particle'])['x'].unstack()
    posYA1 = tA1.set_index(['frame','particle'])['y'].unstack()
    vhXA1 = tp.vanhove(posXA1, int(lagframes), ensemble=True, mpp=0.13, bins=bins)
    vhYA1 = tp.vanhove(posYA1, int(lagframes), ensemble=True, mpp=0.13, bins=bins)
    

    posXB = tB.set_index(['frame','particle'])['x'].unstack()
    posYB = tB.set_index(['frame','particle'])['y'].unstack()
    vhXB = tp.vanhove(posXB, int(lagframes), ensemble=True, mpp=0.13, bins=bins)
    vhYB = tp.vanhove(posYB, int(lagframes), ensemble=True, mpp=0.13, bins=bins)
    
    posXB1 = tB1.set_index(['frame','particle'])['x'].unstack()
    posYB1 = tB1.set_index(['frame','particle'])['y'].unstack()
    vhXB1 = tp.vanhove(posXB1, int(lagframes), ensemble=True, mpp=0.13, bins=bins)
    vhYB1 = tp.vanhove(posYB1, int(lagframes), ensemble=True, mpp=0.13, bins=bins)

    probs[i] = (1./8)*(vhXA.values+vhYA.values+vhXA1.values+vhYA1.values+vhXB.values+vhYB.values+vhXB1.values+vhYB1.values)
        
    popt, pcovR = curve_fit(gaus_func, bins[:-1], probs[i])
    gauss_fit[i] = gaus_func(bins[:-1], popt[0], popt[1], popt[2])
    gauss_width[i] = popt[2]
    
    
    #popt10, pcovR = curve_fit(gaus_func, bins[:-1], prob10)
    #popt10E, pcovRE = curve_fit(gaus_exp_func, bins[:-1], prob10, p0=[2,0,0.1,1,0.1])
    
    #plt.semilogy(bins[:-1], 0.5*(vhX+vhY), 'o', ms=4, color=cmap(clrs[i])
    
for i,lagframes in enumerate(lag_times_in_frames):
    plt_color = cmap(clrs[i])
    plt.semilogy(bins[:-1], probs[i], 'o', ms=4, color=plt_color)
    plt.plot(bins[:-1], gauss_fit[i], '-', lw=2, color=plt_color)
for i,lagframes in enumerate(lag_times_in_frames):
    plt_color = cmap(clrs[i])
    plt.semilogy(bins[:-1], probs[i], 'o', ms=8, mec=plt_color, mfc='w', mew=2)
    plt.plot(bins[:-1], gauss_fit[i], '-', lw=4, color=plt_color)
plt.figure()
for i,lagframes in enumerate(lag_times_in_frames):
    plt_color = cmap(clrs[i])
    plt.semilogy(bins[:-1], probs[i], 'o', ms=8, mec=plt_color, mfc='w', mew=2)
    plt.plot(bins[:-1], gauss_fit[i], '-', lw=4, color=plt_color)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.20577367205542726, 'fraction_numerical': 0.04056966897613549, 'mean_word_length': 3.2110210696920585, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 9, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '2', 'dir': 'jupyter-structured-clean-dedup', 'id': '264073', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 5902, 'n_tokens_neox': 5395, 'n_words': 634, 'path': 'examples/Combine Trajectories m04 to get ISF.ipynb', 'repo_name': 'rmcgorty/trackpy', 'size': '12990'}