From 2eb3c5536c4c9a70d0f0dacc1ad4e952d500202a Mon Sep 17 00:00:00 2001 From: "W. Trevor King" Date: Sat, 1 May 2010 15:42:08 -0400 Subject: [PATCH] Converted order-dep from Gnuplot to Asymptote. This lead to some changes in the Asymptote build system: * New pyfit approach: site_init.link_pyfit() over old Alias('pyfit'). The Alias approach made path extraction ugly and brittle. The new approach mirrors link_wtk_graph(). * Adjust graphFile() and pyfit to allow single column (i.e. implicit x, explicit y) datafiles. * Adjust pyfit to handle single-parameter fitting. I also increased the horizontal spacing between the v-dep-sd histogram insets. --- tex/site_cons/site_init.py | 4 ++ tex/src/figures/asy/wtk_graph.asy | 14 +++-- tex/src/figures/order-dep/SConscript | 54 ++++++++++++---- tex/src/figures/order-dep/fig.gp | 80 ----------------------- tex/src/figures/order-dep/order-dep.asy | 84 +++++++++++++++++++++++++ tex/src/figures/order-dep/order_dep.py | 32 ++++++++++ tex/src/figures/script/SConscript | 1 - tex/src/figures/script/pyfit | 5 ++ tex/src/figures/v-dep/SConscript | 11 ++-- tex/src/figures/v-dep/v-dep-sd.asy | 36 +++++------ tex/src/sawsim/discussion.tex | 6 +- 11 files changed, 200 insertions(+), 127 deletions(-) delete mode 100644 tex/src/figures/order-dep/fig.gp create mode 100644 tex/src/figures/order-dep/order-dep.asy create mode 100644 tex/src/figures/order-dep/order_dep.py diff --git a/tex/site_cons/site_init.py b/tex/site_cons/site_init.py index 43377b1..c1d5f21 100644 --- a/tex/site_cons/site_init.py +++ b/tex/site_cons/site_init.py @@ -65,3 +65,7 @@ def recursive_glob(env, *args, **kwargs): 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')) diff --git a/tex/src/figures/asy/wtk_graph.asy b/tex/src/figures/asy/wtk_graph.asy index 072d292..ca7bec1 100644 --- a/tex/src/figures/asy/wtk_graph.asy +++ b/tex/src/figures/asy/wtk_graph.asy @@ -41,11 +41,15 @@ void graphFile(picture pic=currentpicture, string file="datafile", 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); } diff --git a/tex/src/figures/order-dep/SConscript b/tex/src/figures/order-dep/SConscript index 3e1e590..60a34f3 100644 --- a/tex/src/figures/order-dep/SConscript +++ b/tex/src/figures/order-dep/SConscript @@ -1,19 +1,47 @@ +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') diff --git a/tex/src/figures/order-dep/fig.gp b/tex/src/figures/order-dep/fig.gp deleted file mode 100644 index ec867d5..0000000 --- a/tex/src/figures/order-dep/fig.gp +++ /dev/null @@ -1,80 +0,0 @@ -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 diff --git a/tex/src/figures/order-dep/order-dep.asy b/tex/src/figures/order-dep/order-dep.asy new file mode 100644 index 0000000..a4958ec --- /dev/null +++ b/tex/src/figures/order-dep/order-dep.asy @@ -0,0 +1,84 @@ +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); diff --git a/tex/src/figures/order-dep/order_dep.py b/tex/src/figures/order-dep/order_dep.py new file mode 100644 index 0000000..3e740a7 --- /dev/null +++ b/tex/src/figures/order-dep/order_dep.py @@ -0,0 +1,32 @@ +#!/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) diff --git a/tex/src/figures/script/SConscript b/tex/src/figures/script/SConscript index 41d7526..48ea899 100644 --- a/tex/src/figures/script/SConscript +++ b/tex/src/figures/script/SConscript @@ -2,7 +2,6 @@ Import('env') pyfit = env.File('pyfit') -env.Alias('pyfit', pyfit) # Pass back the modified environment. Return('env') diff --git a/tex/src/figures/script/pyfit b/tex/src/figures/script/pyfit index e8d6ece..f0d6cf7 100755 --- a/tex/src/figures/script/pyfit +++ b/tex/src/figures/script/pyfit @@ -81,6 +81,8 @@ def read_data(datalines, xcol=0, ycol=1, wcol=-1): 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) @@ -93,6 +95,8 @@ def fit(data, fn, initial_ps, _locals={}, verbose=False): 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 @@ -158,6 +162,7 @@ ignored. 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): diff --git a/tex/src/figures/v-dep/SConscript b/tex/src/figures/v-dep/SConscript index d24a713..78c5c15 100644 --- a/tex/src/figures/v-dep/SConscript +++ b/tex/src/figures/v-dep/SConscript @@ -1,6 +1,6 @@ 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'] @@ -12,10 +12,7 @@ data_dirs = [env.Dir(d) for d in subdirs(env, '*')] 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 = [] @@ -24,7 +21,7 @@ for f in data_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, @@ -34,7 +31,7 @@ for f in data_files: 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, diff --git a/tex/src/figures/v-dep/v-dep-sd.asy b/tex/src/figures/v-dep/v-dep-sd.asy index 2555798..5bf72ef 100644 --- a/tex/src/figures/v-dep/v-dep-sd.asy +++ b/tex/src/figures/v-dep/v-dep-sd.asy @@ -3,10 +3,10 @@ import wtk_graph; //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; @@ -15,29 +15,29 @@ real fn_logxliny(real x, real[] params) { 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, @@ -49,8 +49,8 @@ picture hist_picture(string datafile, /* ^-- 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; } @@ -60,11 +60,11 @@ real xmax = 900; 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); diff --git a/tex/src/sawsim/discussion.tex b/tex/src/sawsim/discussion.tex index 232dd21..df48689 100644 --- a/tex/src/sawsim/discussion.tex +++ b/tex/src/sawsim/discussion.tex @@ -156,7 +156,7 @@ unfolding order and polymer length. \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 @@ -171,8 +171,8 @@ unfolding order and polymer length. 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} -- 2.26.2