def link_wtk_graph(env):
return env.Command('wtk_graph.asy', '../asy/wtk_graph.asy',
SCons.Script.Copy('$TARGET', '$SOURCE'))
+
+def link_pyfit(env):
+ return env.Command('pyfit', '../script/pyfit',
+ SCons.Script.Copy('$TARGET', '$SOURCE'))
markroutine markroutine=marknodes,
string t="Title",
bool dots=false) {
- file fin=input(file).line();
- real[][] a=fin.dimension(0,0);
- a=transpose(a);
- real[] x=a[xcol];
- real[] y=a[ycol];
+ file fin = input(file).line();
+ real[][] a = fin.dimension(0,0);
+ a = transpose(a);
+ real[] x;
+ real[] y = a[ycol];
+ if (xcol >= 0)
+ x = a[xcol];
+ else
+ x = sequence(y.length);
graphData(pic=pic, x=x, y=y, xscale=xscale, yscale=yscale, dx=dx, dy=dy, p=p,
mpath=mpath, markroutine=markroutine, t=t, dots=dots);
}
+import os.path
+
+from site_cons.site_init import link_wtk_graph, link_pyfit
+
+
+FIGURES = ['order-dep']
+
# Get the passed in environment.
Import('env')
-#order_data = Glob('data/order.avg-*')
-#hist_data = Glob('data/hist*.hist')
-#fig = env.Command('fig.pdf',
-# ['fig.gp'] + order_data + hist_data,
-# 'gnuplot fig.gp 2>&1 | grep -A10 "^Fitted Kwlcs:" > fit.params',
-# chdir=True)
-envb = env.Clone()
-envb['GNUPLOTCOM'] = \
- ('cd ${TARGET.dir} && $GNUPLOT $GNPLOTLAGS ${SOURCE.file} 2>&1 |'
- 'grep -A10 "^Fitted Kwlcs:" > fit.params')
-fig = envb.Gnuplot('fig.gp')
-env.SideEffect('fit.log', fig)
-env.SideEffect('fit.params', fig)
+data_dir = env.Dir('data')
+data_files = [f for f in env.Glob(os.path.join(str(data_dir), 'order.avg-*'))
+ if 'fit' not in str(f)]
+order_dep_py = File('order_dep.py')
+
+pyfit = link_pyfit(env)
+
+fit_files = []
+for f in data_files:
+ N = int(str(f).rsplit('-', 1)[-1])
+ fit = env.Command(
+ str(f)+'.fit',
+ [f, pyfit, order_dep_py],
+ "python %s -m order_dep:atheory -f 'atheory(i=x, N=%d, Kwlc=A)' -p 1 -x -1 -y 0 -v $SOURCE > $TARGET"
+ % (pyfit[0].get_abspath(), N))
+ fit_dat = env.Command(
+ str(fit[0])+'.dat',
+ fit,
+ "sed -n 's/^[A-Z]: -//p' $SOURCE > $TARGET")
+ # Note the minus --^ which ensures a Kwlc solution >= 0.
+ fit_files.append(fit_dat)
+
+order_dep_data = data_files + fit_files
+
+wtk_graph = link_wtk_graph(env)
+
+for fig in FIGURES:
+ asyfile = '%s.asy' % fig # static .asy file
+ pyfig = fig.replace('-', '_')
+ data = '%s_data' % (pyfig)
+ asydata = []
+ if data in globals(): # generated data dependencies
+ asydata = globals()[data]
+ env.Asymptote([asyfile, wtk_graph] + asydata)
# Pass back the modified environment.
Return('env')
+++ /dev/null
-set terminal pdf enhanced font 'Times,7' size 9cm, 6cm
-set output 'fig.pdf'
-unset key
-set xtics nomirror
-set ytics nomirror
-# border 3 = bottom & left
-set border 3
-set lmargin 9
-
-set multiplot
-
-set size 1.0,0.6
-set origin 0.0,0.0
-set tmargin 0
-set style data points
-set samples 500
-set xlabel 'Unfolding peak index {/Times-Italic i=N_u}'
-set ylabel 'Force (pN)'
-g = 0.5772156649
-kB = 1.3806504e-23
-T = 300
-p = 3.7e-10
-v = 1e-6
-Kc = 0.05
-Kwlc = 0.18
-dx = 0.225e-9
-ko = 5e-5
-f = kB*T/dx
-k(i,N) = (N-i)*ko
-kappa(i,N) = 1.0/(1.0/Kc+i/Kwlc)
-kappaK(i,N,Kwlcs) = 1.0/(1.0/Kc+i/Kwlcs)
-b(k,K) = k*f/(K*v)
-u(k,K) = (-log(b(k,K)) - g)*f
-U(i,N) = (i<=N-1) ? u(k(i,N), kappa(i,N))*1e12 : 1/0
-UK(i,N,Kwlcs) = (i<=N-1) ? u(k(i,N), kappaK(i,N,Kwlcs))*1e12 : 1/0
-
-fit U(x,4) 'data/order.avg-4' using 0:($1*1e12) via Kwlc
-Kwlc4 = Kwlc
-fit U(x,8) 'data/order.avg-8' using 0:($1*1e12) via Kwlc
-Kwlc8 = Kwlc
-fit U(x,12) 'data/order.avg-12' using 0:($1*1e12) via Kwlc
-Kwlc12 = Kwlc
-fit U(x,16) 'data/order.avg-16' using 0:($1*1e12) via Kwlc
-Kwlc16 = Kwlc
-print ''
-print 'Fitted Kwlcs:'
-print 'Kwlc4', ' ', Kwlc4
-print 'Kwlc8', ' ', Kwlc8
-print 'Kwlc12', ' ', Kwlc12
-print 'Kwlc16', ' ', Kwlc16
-
-plot 'data/order.avg-4' using ($1*1e12), UK(x,4,Kwlc4), \
- 'data/order.avg-8' using ($1*1e12), UK(x,8,Kwlc8), \
- 'data/order.avg-12' using ($1*1e12), UK(x,12,Kwlc12), \
- 'data/order.avg-16' using ($1*1e12), UK(x,16,Kwlc16)
-
-set size 0.33,0.4
-set tmargin -1
-set yrange [1e-4:1e-1]
-set xrange [100:380]
-set logscale y
-set xtics 150,100
-set ytics format "10^{%L}"
-set xlabel 'Force (pN)'
-set ylabel 'Frequency'
-set boxwidth 20
-set style data boxes
-set style fill solid 0.3
-#set style fill pattern 4
-
-set origin 0.0,0.6
-plot 'data/hist3i.hist' using ($1*1e12+10):($2/3200)
-
-set origin 0.33,0.6
-plot 'data/hist3ii.hist' using ($1*1e12+10):($2/3200)
-
-set origin 0.67,0.6
-plot 'data/hist3iii.hist' using ($1*1e12+10):($2/3200)
-
-unset multiplot
--- /dev/null
+import wtk_graph;
+
+picture pic;
+size(pic, 15cm,5cm,IgnoreAspect);
+
+scale(pic, Linear, Linear);
+real xscale=1;
+real fscale=1e12;
+
+real g = 0.5772156649;
+real kB = 1.3806504e-23;
+real T = 300;
+real p = 3.7e-10;
+real v = 1e-6;
+real Kc = 0.05;
+real Kwlc = 0.18;
+real dx = 0.225e-9;
+real ko = 5e-5;
+real f = kB*T/dx;
+
+real kappa(real i, real Kc=Kc, real Kwlc=Kwlc) {
+ return 1.0/(1.0/Kc + i/Kwlc);
+}
+
+real theory(real i, int N, real Kc=Kc, real Kwlc=Kwlc, real ko=ko, real dx=dx,
+ real v=v, real T=T, real kB=kB, real g=g) {
+ real f = kB*T/dx;
+ return f*(log((kappa(i, Kc, Kwlc)*v)/((N-i)*ko*f)) - g);
+}
+
+// todo: colorpen(N,i)
+graphFile(pic, "data/order.avg-4", xcol=-1, ycol=0,
+ xscale=xscale, yscale=fscale, p=psoft, dots=true);
+graphFile(pic, "data/order.avg-8", xcol=-1, ycol=0,
+ xscale=xscale, yscale=fscale, p=pmed, dots=true);
+graphFile(pic, "data/order.avg-12", xcol=-1, ycol=0,
+ xscale=xscale, yscale=fscale, p=phard, dots=true);
+graphFile(pic, "data/order.avg-16", xcol=-1, ycol=0,
+ xscale=xscale, yscale=fscale, p=black, dots=true);
+real fn(real x, real[] params) {return theory(i=x, N=4, Kwlc=params[0]);}
+fitFile(pic, "data/order.avg-4.fit.dat", f=fn, xmin=0, xmax=3,
+ xscale=xscale, yscale=fscale, p=psoft);
+real fn(real x, real[] params) {return theory(i=x, N=8, Kwlc=params[0]);}
+fitFile(pic, "data/order.avg-8.fit.dat", f=fn, xmin=0, xmax=7,
+ xscale=xscale, yscale=fscale, p=pmed);
+real fn(real x, real[] params) {return theory(i=x, N=12, Kwlc=params[0]);}
+fitFile(pic, "data/order.avg-12.fit.dat", f=fn, xmin=0, xmax=11,
+ xscale=xscale, yscale=fscale, p=phard);
+real fn(real x, real[] params) {return theory(i=x, N=16, Kwlc=params[0]);}
+fitFile(pic, "data/order.avg-16.fit.dat", f=fn, xmin=0, xmax=15,
+ xscale=xscale, yscale=fscale, p=black);
+
+xaxis(pic, sLabel("Unfolding peak index $i=N_u$"), BottomTop, LeftTicks);
+yaxis(pic, sLabel("Force (pN)"), LeftRight, RightTicks);
+
+picture hist_picture(string datafile,
+ real xmin=-infinity, real xmax=infinity,
+ real ymin=-infinity, real ymax=infinity) {
+ picture pic;
+ size(pic, 4cm, 3cm, IgnoreAspect);
+ scale(pic, Linear, Log);
+ histFile(pic, datafile, bin_scale=fscale, low=ymin);
+ /* ^-- yscale b/c bins are in Force */
+ xlimits(pic, xmin, xmax);
+ ylimits(pic, ymin, ymax);
+ xaxis(pic, sLabel("Force (pN)"), BottomTop, LeftTicks(N=3), above=true);
+ yaxis(pic, sLabel("Frequency"), LeftRight, RightTicks, above=true);
+ return pic;
+}
+
+/* ensure consistent ranges across all histograms */
+real xmin = 100;
+real xmax = 380;
+real ymin = 1e-4;
+real ymax = 1;
+
+add(pic.fit(), (0.0), S);
+add(hist_picture("data/hist3i.hist", xmin, xmax, ymin, ymax).fit(),
+ 5cm*W, N);
+add(hist_picture("data/hist3ii.hist", xmin, xmax, ymin, ymax).fit(),
+ (0,0), N);
+add(hist_picture("data/hist3iii.hist", xmin, xmax, ymin, ymax).fit(),
+ 5cm*E, N);
+label(sLabel("Pulling speed dependence"), point(N)+(0,3cm), N);
--- /dev/null
+#!/usr/bin/python
+
+from math import log
+
+
+g = 0.5772156649
+kB = 1.3806504e-23
+T = 300
+p = 3.7e-10
+v = 1e-6
+Kc = 0.05
+Kwlc = 0.18
+dx = 0.225e-9
+ko = 5e-5
+f = kB*T/dx
+
+def kappa(i, Kc=Kc, Kwlc=Kwlc):
+ return 1.0/(1.0/Kc + i/Kwlc)
+
+def theory(i, N, Kc=Kc, Kwlc=Kwlc, ko=ko, dx=dx, v=v, T=T, kB=kB, g=g):
+ f = kB*T/dx;
+ return f*(log((kappa(i, Kc, Kwlc)*v)/((N-i)*ko*f)) - g)
+
+def atheory(*args, **kwargs):
+ """Call theory with abs() of all arguments, since the fitting
+ routine may attempt negative values.
+ """
+ for i,x in enumerate(args):
+ args[i] = abs(x)
+ for k,v in kwargs.items():
+ kwargs[k] = abs(v)
+ return theory(*args, **kwargs)
Import('env')
pyfit = env.File('pyfit')
-env.Alias('pyfit', pyfit)
# Pass back the modified environment.
Return('env')
ys.append(row[ycol])
if wcol >= 0:
ws.append(row[wcol])
+ if xcol < 0:
+ xs = range(len(ys))
if wcol < 0:
ws = [1]*len(xs)
return (xs, ys, ws)
for x,y,sqw in zip(xs,ys,sqws)]
p,cov,info,mesg,ier = leastsq(residuals, initial_ps,
full_output=True, maxfev=10000)
+ if not hasattr(p, '__len__'):
+ p = [p] # with only one param, leastsq doesn't return a list
if verbose == True:
print "Fitted params:",p
print "Covariance mx:",cov
assert len(initial_ps) < 26, ('More parameters than letters! (%d)'
% len(params))
params = fit(data, fn, initial_ps, mod_includes, verbose=options.verbose)
+
assert len(params) == len(initial_ps)
for i,p in enumerate(params):
import os.path
-from site_cons.site_init import subdirs, link_wtk_graph
+from site_cons.site_init import subdirs, link_wtk_graph, link_pyfit
FIGURES = ['v-dep', 'v-dep-sd']
data_files = [env.File(os.path.join(str(d), 'v_dep'))
for d in data_dirs]
-pyfit = env.Alias('pyfit')
-# Extract the first source Node required to build the 'pyfit' Alias
-# (i.e. the File('pyfit') node).
-pyfit = pyfit[0].get_binfo().bsources[0]
+pyfit = link_pyfit(env)
fit_files = []
sd_fit_files = []
str(f)+'.fit',
[f, pyfit],
"python %s -m math:log -f 'A*log(x,10)+B' -v $SOURCE > $TARGET"
- % str(pyfit))
+ % pyfit[0].get_abspath())
fit_dat = env.Command(
str(fit[0])+'.dat',
fit,
str(f)+'-sd.fit',
[f, pyfit],
"python %s -m math:log -f 'A*log(x,10)+B' -y 2 -v $SOURCE > $TARGET"
- % str(pyfit))
+ % pyfit[0].get_abspath())
sd_fit_dat = env.Command(
str(sd_fit[0])+'.dat',
sd_fit,
//size(15cm, 10cm, IgnoreAspect);
-picture sd;
-size(sd, 13cm, 5cm, IgnoreAspect);
+picture pic;
+size(pic, 13cm, 5cm, IgnoreAspect);
-scale(sd, Log, Linear);
+scale(pic, Log, Linear);
real xscale=1e9;
real yscale=1e12;
return params[0] * log10(x) + params[1];
}
-graphFile(sd, "v_dep-5e-5_0.1e-9/v_dep", ycol=2,
+graphFile(pic, "v_dep-5e-5_0.1e-9/v_dep", ycol=2,
xscale=xscale, yscale=yscale, p=psoft,
t="$k_{u0}=5\cdot10^{-5}\mbox{ s$^{-1}$}$, $\Delta x_u=0.100\mbox{ nm}$",
dots=true);
-graphFile(sd, "v_dep-1e-5_0.225e-9/v_dep", ycol=2,
+graphFile(pic, "v_dep-1e-5_0.225e-9/v_dep", ycol=2,
xscale=xscale, yscale=yscale, p=pmed,
t="$k_{u0}=1\cdot10^{-5}\mbox{ s$^{-1}$}$, $\Delta x_u=0.225\mbox{ nm}$",
dots=true);
-graphFile(sd, "v_dep-5e-5_0.225e-9/v_dep", ycol=2,
+graphFile(pic, "v_dep-5e-5_0.225e-9/v_dep", ycol=2,
xscale=xscale, yscale=yscale, p=phard,
t="$k_{u0}=5\cdot10^{-5}\mbox{ s$^{-1}$}$, $\Delta x_u=0.225\mbox{ nm}$",
dots=true);
-fitFile(sd, "v_dep-5e-5_0.1e-9/v_dep-sd.fit.dat", f=fn_logxliny,
+fitFile(pic, "v_dep-5e-5_0.1e-9/v_dep-sd.fit.dat", f=fn_logxliny,
xmin=1e-9, xmax=1e-5, xscale=xscale, yscale=yscale, p=psoft);
-fitFile(sd, "v_dep-1e-5_0.225e-9/v_dep-sd.fit.dat", f=fn_logxliny,
+fitFile(pic, "v_dep-1e-5_0.225e-9/v_dep-sd.fit.dat", f=fn_logxliny,
xmin=1e-9, xmax=1e-5, xscale=xscale, yscale=yscale, p=pmed);
-fitFile(sd, "v_dep-5e-5_0.225e-9/v_dep-sd.fit.dat", f=fn_logxliny,
+fitFile(pic, "v_dep-5e-5_0.225e-9/v_dep-sd.fit.dat", f=fn_logxliny,
xmin=1e-9, xmax=1e-5, xscale=xscale, yscale=yscale, p=phard);
-xequals(sd, x=1e-6*xscale, p=dashed);
+xequals(pic, x=1e-6*xscale, p=dashed);
-xaxis(sd, sLabel("Pulling speed (nm/s)"), BottomTop, LeftTicks);
-yaxis(sd, sLabel("Unfolding force (pN)"), LeftRight, RightTicks);
-add(sd, legend(sd), point(sd, E), 20E, UnFill);
+xaxis(pic, sLabel("Pulling speed (nm/s)"), BottomTop, LeftTicks);
+yaxis(pic, sLabel("Unfolding force (pN)"), LeftRight, RightTicks);
+add(pic, legend(pic), point(pic, E), 20E, UnFill);
picture hist_picture(string datafile,
real xmin=-infinity, real xmax=infinity,
/* ^-- yscale b/c bins are in Force */
xlimits(pic, xmin, xmax);
ylimits(pic, ymin, ymax);
- xaxis(pic, sLabel("Force (pN)"), BottomTop, LeftTicks);
- yaxis(pic, sLabel("Frequency"), LeftRight, RightTicks);
+ xaxis(pic, sLabel("Force (pN)"), BottomTop, LeftTicks(N=2), above=true);
+ yaxis(pic, sLabel("Frequency"), LeftRight, RightTicks, above=true);
return pic;
}
real ymin = 1e-4;
real ymax = 1;
-add(sd.fit(), (0,0), S);
+add(pic.fit(), (0,0), S);
add(hist_picture("fig4i-5e-5_0.1e-9.hist", xmin, xmax, ymin, ymax).fit(),
- 4cm*W, N);
+ 5cm*W, N);
add(hist_picture("fig4i-1e-5_0.225e-9.hist", xmin, xmax, ymin, ymax).fit(),
(0,0), N);
add(hist_picture("fig4i-5e-5_0.225e-9.hist", xmin, xmax, ymin, ymax).fit(),
- 4cm*E, N);
+ 5cm*E, N);
label(sLabel("Pulling speed width dependence"), point(N)+(0,3cm), N);
\begin{figure}
\begin{center}
- \includegraphics{figures/order-dep/fig}
+ \asyfig{figures/order-dep/order-dep}
\caption{The dependence of the unfolding force on the temporal
unfolding order for four polymers with $4$, $8$, $12$, and $16$
identical protein domains. Each point in the figure is the
to right, for the polymer with eight protein domains. The
parameters used for generating the data were the same as those
used for \cref{fig:sawsim:sim-sawtooth}, except for the number of
- domains. The histograms in the insets were normalized in the same
- way as in \cref{fig:sawsim:sim-hist}.\label{fig:sawsim:order-dep}}
+ domains. The histogram insets were normalized in the same way as
+ in \cref{fig:sawsim:sim-hist}.\label{fig:sawsim:order-dep}}
\end{center}
\end{figure}
projects / thesis.git / commitdiff