In [2]:
%pylab inline
import sys
sys.path.append('FOLDER') ## Put the folder in which contract.py sits
import contract
from scipy.integrate import odeint
simu_length=20.
simu_points=2000.
ngp=256
In [3]:
##LOAD THE DENSITY DATA HERE
from scipy.io import loadmat
f=loadmat('Fixed_Profiles.mat')
pp=f['Profile']
##LOAD VELOCITY DATA HERE
from scipy.io import loadmat
f=loadmat('Velocity_timecourse.mat')
xt=f['Position_timecourse'].flatten()
xt=[xx[0] for xx in xt]
vt=f['Velocity_timecourse']
vt=[xx.flatten() for xx in vt[0]]
In [5]:
## Defines how the time evolution is done.
def compute(a,b,r0,w_init=None):
l=sqrt(a/b)
s0=-1./sqrt(a*b)
#if no w_init is given take it from movie
if w_init==None:
ww=[i for i in xrange(206) if not isnan(pp[0,i])]
w_init=max(ww)-min(ww)
w_init=w_init*8.
#set parameters
contract.l=l
contract.s0=s0
contract.rho0=1.
contract.rho_init=r0*ones(ngp)
contract.w_init=w_init
t=linspace(0,simu_length,simu_points)
return odeint(contract.dt_packed,contract.init_packed(),t)
In [6]:
def normalization():
##obtain width of first experimental frame
ww=[i for i in xrange(206) if not isnan(pp[0,i])]
w_init=max(ww)-min(ww)
w_init=w_init*8.
##obtain experimental mean level in the first frame
ppc=[i for i in xrange(206) if not isnan(pp[0][i])]
pm=mean(ppc)
exp_mean0=mean([pp[0,i] for i in ppc])*w_init
##obtain theory mean level in the frame cooresponding to the first experimental frame
th_mean0=mean(o[0][1:ngp+1])*o[0][0]
return th_mean0/exp_mean0
In [7]:
##This defines the plotting routine
##frames are the frames which you wish to show
##offset omit the first offset frames (this is to get rid of the initial lag phase)
def and_plot_velocity(frames=(0,10,15,25,30),framerate=3.,offset=0):
##obtain time offset between simulation and experiment
simu_frame_offset=offset/simu_length*simu_points
r=normalization()
colors={0:'r',1:'g',2:'b'}
c=0
for j in frames:
#pick simulation frame to be shown
#(Here the total sim length is div minutes with 2000 time points stored)
k=(j)/framerate/simu_length*simu_points+simu_frame_offset
#center data to show
#ppc=[i for i in xrange(206) if not isnan(pp[j][i])]
#pm=mean(ppc)
#plot
#plot((arange(206)-pm)*8,pp[j],'-',color=colors[c])
if k>=0:
v=contract.get_v(o[k][0],o[k][1:])
plot(linspace(0,o[k][0]/2.,ngp),v,'--',lw=3,color=colors[c])
plot(linspace(0,-o[k][0]/2.,ngp),-v,'--',lw=3,color=colors[c])
#xlim(-650, 650)
#ylim(.0055, .0185)
#pick new color for the next curve
c=(c+1)%3
In [8]:
##This defines the plotting routine
##frames are the frames which you wish to show
##offset omit the first offset frames (this is to get rid of the initial lag phase)
def and_plot(frames=(0,10,15,25,30),framerate=3.,offset=0):
##obtain time offset between simulation and experiment
simu_frame_offset=offset/simu_length*simu_points
r=normalization()
colors={0:'r',1:'g',2:'b'}
c=0
for j in frames:
#pick simulation frame to be shown
#(Here the total sim length is div minutes with 2000 time points stored)
k=(j)/framerate/simu_length*simu_points+simu_frame_offset
#center data to show
ppc=[i for i in xrange(206) if not isnan(pp[j][i])]
pm=mean(ppc)
#plot
plot((arange(206)-pm)*8,pp[j],'-',color=colors[c])
if k>=0:
plot(linspace(0,o[k][0]/2.,ngp),o[k][1:]/r,'--',lw=3,color=colors[c])
plot(linspace(0,-o[k][0]/2.,ngp),o[k][1:]/r,'--',lw=3,color=colors[c])
xlim(-650, 650)
# ylim(.0055, .0185)
#pick new color for the next curve
c=(c+1)%3
In [10]:
a,b,r0=array([ 2.26413429e+00, 7.09505347e-05, 3.26194452e-01]) #
## these are the values coming from a leastsq procedure
a=a/3.
b=b/3.
b=b/4.
o=compute(a,b,r0,w_init=647*2)
and_plot(offset=-0.8,frames=(3,6,9,12,15))
xlabel('Channel Position ($\mu$m)', size=16)
ylabel('Relative Density', size=16)
Out[10]:
In [11]:
def fit_fun(args):
offset=-0.4
framerate=3.
simu_frame_offset=offset/simu_length*simu_points
from scipy.interpolate import interp1d
o=compute(*args,w_init=2*694)
oo=[]
r=normalization()
for j in range(0,15):
k=(j)/framerate/simu_length*simu_points+simu_frame_offset
if k>=0:
#wrap data in function
ppc=[i for i in xrange(206) if not isnan(pp[j][i])]
pm=mean(ppc)
p=nan_to_num(pp[j])
x=(arange(206)-pm)*8
data=interp1d((arange(206)-pm)*8,p)
edat=data(x)
#wrap theory result
xt=linspace(0,o[k][0]/2.,ngp)
yt=o[k][1:ngp+1]/r
theory=interp1d(xt,yt)
tdat=[]
for xx in x:
try:
tdat.append(theory(abs(xx)))
except:
tdat.append(0.)
tdat=array(tdat)
oo.append(tdat-edat)
return array(oo).flatten()
oo=fit_fun((a,b,r0))
In [12]:
from scipy.optimize import leastsq,fmin
best_fit=leastsq(fit_fun,(a,b,r0))
print best_fit
In [13]:
a,b,r0=best_fit[0]
o=compute(a,b,r0,w_init=2*647)
and_plot(offset=-0.8,frames=(3,6,9,12,15))
xlabel('Channel Position ($\mu$m)', size=16)
ylabel('Relative Density', size=16)
In [14]:
offset=-.4+.5
framerate=3.
simu_frame_offset=offset/simu_length*simu_points
r=normalization()
colors={0:'r',1:'b',2:'g'}
c=0
for j in (5,10):
#pick simulation frame to be shown
#(Here the total sim length is div minutes with 2000 time points stored)
k=(j)/framerate/simu_length*simu_points+simu_frame_offset
plot(xt[j]-mean(xt[j]),(vt[j]-mean(vt[j])),'-', lw=3.,color=colors[c])
if k>=0:
v=contract.get_v(o[k][0],o[k][1:])
plot(linspace(0,o[k][0]/2.,ngp),v,'--',lw=3,color=colors[c],label='t=%i s'%(j*20))
plot(linspace(0,-o[k][0]/2.,ngp),-v,'--',lw=3,color=colors[c])
#pick new color for the next curve
c=(c+1)%3
legend()
grid()
ylim(-120,120)
xlabel('Channel Position ($\mu$m)', size=16)
ylabel('Velocity ($\mu$m / min)', size=16)
In [15]:
def fit_fun2(args):
a,b,r0,offset,w_init=args
framerate=3.
simu_frame_offset=offset/simu_length*simu_points
from scipy.interpolate import interp1d
o=compute(a,b,r0,w_init=2*694)
oo=[]
r=normalization()
for j in range(0,30):
k=(j)/framerate/simu_length*simu_points+simu_frame_offset
if k>=0:
#wrap data in function
ppc=[i for i in xrange(206) if not isnan(pp[j][i])]
pm=mean(ppc)
p=nan_to_num(pp[j])
x=(arange(206)-pm)*8
data=interp1d((arange(206)-pm)*8,p)
edat=data(x)
#wrap theory result
xt=linspace(0,o[k][0]/2.,ngp)
yt=o[k][1:ngp+1]/r
theory=interp1d(xt,yt)
tdat=[]
for xx in x:
try:
tdat.append(theory(abs(xx)))
except:
tdat.append(0.)
tdat=array(tdat)
oo.append(tdat-edat)
return array(oo).flatten()
#oo=fit_fun2((a,b,r0,-0,2.*450))
In [16]:
from scipy.optimize import leastsq,fmin
a,b,r0=array([ 9.00036043e-01, 3.73472775e-06, 3.03218838e-01]) ### these are the values coming from a leastsq procedure
best_fit2=leastsq(fit_fun2,(a,b,r0,-.4,2*694),maxfev=10000000,xtol=1.49012e-12,ftol=1.49012e-12,epsfcn=1.0e-3)
print best_fit2
In [17]:
a,b,r0,offset,w_init=best_fit2[0]
o=compute(a,b,r0,w_init=w_init)
#offset=-.4
and_plot(offset=offset,frames=(3,9,15,21,27))
xlabel('Channel Position ($\mu$m)', size=16)
ylabel('Relative Density', size=16)
In [42]:
def make_movie_velocity(frame_max=30,folder='/Users/peterfoster/Desktop/Profile_Fitting_Simulations/velocity_movie/',framerate=3.,offset=-0.4):
simu_frame_offset=offset/simu_length*simu_points
r=normalization()
for j in xrange(frame_max):
#pick simulation frame to be shown
k=(j)/framerate/simu_length*simu_points+simu_frame_offset
#center data to show
ppc=[i for i in xrange(206) if not isnan(pp[j][i])]
pm=mean(ppc)
#plot
if k>=0:
v=contract.get_v(o[k][0],o[k][1:])
plot(linspace(0,o[k][0]/2.,ngp),o[k][1:]/r,lw=4,color='b')
plot(linspace(0,-o[k][0]/2.,ngp),o[k][1:]/r,lw=4,color='b')
ylim(0,0.02)
xlim(-700, 700)
twinx()
plot(linspace(0,o[k][0]/2.,ngp),v,'--',lw=3,color='r')
plot(linspace(0,-o[k][0]/2.,ngp),-v,'--',lw=3,color='r')
#pick new color for the next curve
grid()
ylim(-120,120)
xlim(-700, 700)
#xlabel('Channel Position ($\mu$m)', size=16)
#ylabel('Velocity ($\mu$m / min)', size=16)
fname=folder+'fig_v_%03i.png'%j
savefig(fname)
clf()
make_movie_velocity()
In [29]:
def make_movie_density(frame_max=30,folder='/Users/peterfoster/Desktop/Profile_Fitting_Simulations/data_overlap_movie/',framerate=3.,offset=-0.4):
simu_frame_offset=offset/simu_length*simu_points
r=normalization()
for j in xrange(frame_max):
#pick simulation frame to be shown
k=(j)/framerate/simu_length*simu_points+simu_frame_offset
#center data to show
ppc=[i for i in xrange(206) if not isnan(pp[j][i])]
pm=mean(ppc)
#plot
plot((arange(206)-pm)*8,pp[j],'-',lw=3,color='b')
twinx()
if k>=0:
plot(linspace(0,o[k][0]/2.,ngp),o[k][1:]/r,'--',lw=4,color='r')
plot(linspace(0,-o[k][0]/2.,ngp),o[k][1:]/r,'--',lw=4,color='r')
xlim(-650, 650)
ylim(.0055, .02)
xlabel('Channel Position ($\mu$m)', size=16)
ylabel('Relative Density', size=16)
fname=folder+'fig_rho_%03i.png'%j
savefig(fname)
clf()
make_movie_density()