myprojects-hg-reborn

Fork

#! /usr/bin/env python

import scipy as s
import numpy as n
import pylab as p

time_red=p.load("time_red.dat")  #time I used when comparing the snapshots


delta_t=time_red[1]-time_red[0]


time_coll=time_red[1:len(time_red)]

print "the snapshot time series is, ", time_red


N_mon=5000. #initial number of monomers

rho=0.01

box_vol=N_mon/rho


by=2

#NB: it all starts from the second configuration!!!! I need two snapshots to take a view of what is
#going on.

my_selection=s.arange(1,len(time_red))*by

k_ij=s.zeros(len(my_selection))

number_coll=s.zeros(len(my_selection))

n_i=s.zeros(len(my_selection))

n_j=s.zeros(len(my_selection))



sel_ker=s.arange(2)
sel_ker[0]=2
sel_ker[1]=3


sel_ker=sel_ker*1.  #to get a floating point array

#I will use sel_ker defined above to select the [1,1] collisions

print "sel_ker is, ", sel_ker

print "my_selection is, ", my_selection


delta_sel=my_selection[1]-my_selection[0]

for i in xrange(len(my_selection)):

    cluster_name="cluster_dist%05d"%(my_selection[i]-delta_sel)

    #print "cluster_name is, ", cluster_name
    
    size_dist=p.load(cluster_name)




    cluster_name="recorded_collisions%05d"%my_selection[i]
    extension=".dat"
    cluster_name=cluster_name+extension
    recorded_collisions=p.load(cluster_name)




    cluster_name="colliding_sizes%05d"%my_selection[i]
    extension=".dat"
    cluster_name=cluster_name+extension
    colliding_sizes=p.load(cluster_name)

    cluster_name="kernel_elements%05d"%my_selection[i]
    extension=".dat"
    cluster_name=cluster_name+extension

    kernel_elements=p.load(cluster_name)

    if (set(sel_ker)<=set(size_dist)):

        n_i[i]=len(s.where(size_dist==sel_ker[0])[0])
        n_j[i]=len(s.where(size_dist==sel_ker[1])[0])


        for j in xrange(len(kernel_elements)):
            if (colliding_sizes[j,:]==sel_ker).all():
                k_ij[i]=kernel_elements[j]
                number_coll[i]=recorded_collisions[j]

    else:

        k_ij[i]=-2.


print "the number of collisions are, ", number_coll



print "k_ij is, ", k_ij   


#Now I get rid of the negative k_ij elements (which corresponds to times where either n_i or n_j where 0)
sel_pos=s.where(k_ij>=0.)

print "len(sel_pos[0]) is, ", len(sel_pos[0])


k_ij=k_ij[sel_pos]

#and then I get rid of the corresponding times


time_coll=time_coll[sel_pos]

#and then do the same for n_i and n_j

n_i=n_i[sel_pos]
n_j=n_j[sel_pos]

#and finally for the number of collisions

number_coll=number_coll[sel_pos]



print "len(n_i) is, ", len(n_i)



k_ij_mean=k_ij.mean()

print "the mean value of k_ij is, ", k_ij_mean
    


print "s.shape(k_ij) is, ", s.shape(k_ij)


print "s.shape(time_red) is, ", s.shape(time_red)


print "n_i is, ", n_i

print "n_j is, ", n_j


fig = p.figure()
axes = fig.gca()


axes.plot(time_coll,k_ij, "bo", label="k_ij")
p.xlabel('Time')
p.ylabel('Istantaneous Value k_ij')
#p.title("Evolution Mean-free path")
p.grid(True)
cluster_name="kernel_element_vs_time.pdf"
axes.legend()
p.savefig(cluster_name)
 
p.clf()




print "len(time_coll) is, ", len(time_coll)
print "len(number_coll) is, ", len(number_coll)





fig = p.figure()
axes = fig.gca()


axes.plot(time_coll,number_coll, "bo", label="i-j collisions")
p.xlabel('Time')
p.ylabel('Number of collisions')
#p.title("Evolution Mean-free path")
p.grid(True)
cluster_name="i-j_collisions_vs_time.pdf"
axes.legend()
p.savefig(cluster_name)
 
p.clf()





fig = p.figure()
axes = fig.gca()


axes.plot(time_coll,n_i, "bo", label="n_i")

axes.plot(time_coll,n_j, "k^", label="n_j")


p.xlabel('Time')
p.ylabel('Number of i-mers and j-mers')
#p.title("Evolution Mean-free path")
p.grid(True)
cluster_name="number_i-mers_and_j-mers_vs_time.pdf"
axes.legend()
p.savefig(cluster_name)
 
p.clf()





fig = p.figure()
axes = fig.gca()


axes.plot(time_coll,n_i/box_vol, "bo", label="n_i [concentration]")

axes.plot(time_coll,n_j/box_vol, "k^", label="n_j [concentration]")


p.xlabel('Time')
p.ylabel('Concentration of i-mers and j-mers')
#p.title("Evolution Mean-free path")
p.grid(True)
cluster_name="concentration_i-mers_and_j-mers_vs_time.pdf"
axes.legend()
p.savefig(cluster_name)
 
p.clf()


p.save("n_i_concentration.dat", n_i/box_vol)


p.save("n_j_concentration.dat", n_j/box_vol)


#now a new mean of k11:

#time_11=time_red[1:len(time_red)]


time_late=time_coll[s.where(time_coll>120.)]

k_ij_late=k_ij[s.where(time_coll>120.)]


print "the mean value of the kernel element at late times is, ", k_ij_late.mean()

time_early=time_coll[s.where(time_coll<=120.)]

k_ij_early=k_ij[s.where(time_coll<=120.)]


print "the mean value of the kernel element at early times is, ", k_ij_early.mean()



def residuals_eval(fun1, fun2, par_inf, par_sup, n_steps):
    par_seq=s.linspace(par_inf, par_sup, n_steps) #interval where I am varying my parameter

    residual_vec=s.zeros(n_steps)

    #fun1 is taken to be n_{ij}=N_{ij}/box_vol and fun2 is n_in_j, the proportionality constant
    #(to be determined) is the kernel

    for i in xrange(n_steps):

        residual_vec[i]=s.sum((fun1-par_seq[i]*fun2)**2.)

    par_opt=par_seq[s.where(residual_vec==min(residual_vec))]

    return par_opt, residual_vec




k_ij_inf=0.2

k_ij_sup=12.

n_steps=120

#delta_t=4.

coll_dens=number_coll/(box_vol*delta_t)

print "coll_dens is, ", coll_dens

n_in_j=n_i/box_vol*n_j/box_vol

print "n_in_j is, ", n_in_j

fitted_kernel=residuals_eval(coll_dens, n_in_j, k_ij_inf, k_ij_sup, n_steps)

print "fitted_kernel is, ", fitted_kernel[0]





fig = p.figure()
axes = fig.gca()


axes.plot(time_coll,coll_dens, "bo", label="n_ij")

axes.plot(time_coll,fitted_kernel[0]*n_in_j, "k^", label="k_{ij}n_in_j")


p.xlabel('Time')
p.ylabel('i-mers and j-mers collisions per unit volume and time')
#p.title("Evolution Mean-free path")
p.grid(True)
cluster_name="i-mers_and_j-mers_collisions_and_fit_vs_time.pdf"
axes.legend()
p.savefig(cluster_name)
 
p.clf()

p.save("residuals.dat",fitted_kernel[1])


fig = p.figure()
axes = fig.gca()


axes.plot(s.linspace(k_ij_inf,k_ij_sup, n_steps),fitted_kernel[1], "bo", label="residuals")

p.xlabel('Value of k_{ij}')
p.ylabel('Squared sum of the residuals')
#p.title("Evolution Mean-free path")
p.grid(True)
cluster_name="residuals_vs_k_ij.pdf"
axes.legend()
p.savefig(cluster_name)
 
p.clf()

#Now I do as above but I specify a time-span

time_inf=0.
time_sup=100.

time_span=s.where((time_coll>time_inf) & (time_coll<time_sup))


fitted_kernel=residuals_eval(coll_dens[time_span], n_in_j[time_span], k_ij_inf, k_ij_sup, n_steps)

print "fitted_kernel is, ", fitted_kernel[0]




fig = p.figure()
axes = fig.gca()


axes.plot(time_coll[time_span],coll_dens[time_span], "bo", label="n_ij")

axes.plot(time_coll[time_span],fitted_kernel[0]*n_in_j[time_span], "k^", label="k_{ij}n_in_j")


p.xlabel('Time')
p.ylabel('i-mers and j-mers collisions per unit volume and time')
#p.title("Evolution Mean-free path")
p.grid(True)
cluster_name="i-mers_and_j-mers_collisions_and_fit_vs_time_time_span.pdf"
axes.legend()
p.savefig(cluster_name)
 
p.clf()



fig = p.figure()
axes = fig.gca()


axes.plot(s.linspace(k_ij_inf,k_ij_sup, n_steps),fitted_kernel[1], "bo", label="residuals")

p.xlabel('Value of k_{ij}')
p.ylabel('Squared sum of the residuals')
#p.title("Evolution Mean-free path")
p.grid(True)
cluster_name="residuals_vs_k_ij_time_span.pdf"
axes.legend()
p.savefig(cluster_name)
 
p.clf()

p.save("residuals_span.dat", fitted_kernel[1])

print "So far so good"