From 780d75e14fb1a05686db24bb0599d214a938dd21 Mon Sep 17 00:00:00 2001 From: "W. Trevor King" Date: Tue, 9 Sep 2008 13:14:07 -0400 Subject: [PATCH] Added testing scripts, comparing sawsim with theory TODO, adjust bell & kramers testing to read param files like const_rate does. --- testing/README | 22 ++++++ testing/bell_rate/fit_bell | 67 +++++++++++++++++ testing/bell_rate/run_test.sh | 103 ++++++++++++++++++++++++++ testing/common/fit.py | 89 ++++++++++++++++++++++ testing/common/fit.pyc | Bin 0 -> 4361 bytes testing/common/fitHist.pyc | Bin 0 -> 4362 bytes testing/const_rate/const_rate.sh | 101 +++++++++++++++++++++++++ testing/const_rate/fit_exponential | 35 +++++++++ testing/const_rate/params | 5 ++ testing/const_rate/run_test.sh | 13 ++++ testing/kramers_rate/fit_kramers | 46 ++++++++++++ testing/kramers_rate/run_test.sh | 115 +++++++++++++++++++++++++++++ testing/test_cluster.sh | 11 +++ 13 files changed, 607 insertions(+) create mode 100644 testing/README create mode 100755 testing/bell_rate/fit_bell create mode 100755 testing/bell_rate/run_test.sh create mode 100755 testing/common/fit.py create mode 100644 testing/common/fit.pyc create mode 100644 testing/common/fitHist.pyc create mode 100755 testing/const_rate/const_rate.sh create mode 100755 testing/const_rate/fit_exponential create mode 100644 testing/const_rate/params create mode 100755 testing/const_rate/run_test.sh create mode 100755 testing/kramers_rate/fit_kramers create mode 100755 testing/kramers_rate/run_test.sh create mode 100755 testing/test_cluster.sh diff --git a/testing/README b/testing/README new file mode 100644 index 0000000..5e63141 --- /dev/null +++ b/testing/README @@ -0,0 +1,22 @@ +A bunch of validation checks. Run automatically before allowing +merges/commits to the main branch. The smaller, faster subset is run +before allowing commits/merges to any branch. + +Since this software is being developed on a cluster, many of the tests +take advantage of the more powerful environment. You should still be +able to run the tests if you're on a standard machine though, it will +just take more time :p. To allow for the differences in invocation, +the most tests are broken out into two scripts, a X.cl to be run on a +cluster and an X.sh to be run on a standard machine. + +Each test directory contains the script run_test.sh, which gives +automated testing scripts a common file name to look for. Exporting +the global variable ISACLUSTER allows more efficient execution on a +cluster. This can be accomplised (in Bash) with the command + + $ export ISACLUSTER=1 + +Note that it is also necessary to `touch' the testing directory after +you add or remove a test to update bin/run-test.sh. + +The common directory stores code shared among the tests. diff --git a/testing/bell_rate/fit_bell b/testing/bell_rate/fit_bell new file mode 100755 index 0000000..c6433c3 --- /dev/null +++ b/testing/bell_rate/fit_bell @@ -0,0 +1,67 @@ +#!/usr/bin/python +# +# fit data coming in on stdin to bell-theory histogram +# output: +# +# a:\t +# b:\t +# N:\t +# bin:\t +# +# usage: fit_bell +# for example: +# $ cat datafile | fit_bell + +import sys +sys.path.append('../common/') +import fit +from scipy import exp + +class Model (fit.Model) : + def printModel(self) : + print "bin*N*a*exp(b*x + a/b(1-exp(b*x)))" + def printParams(self, xopt) : + print "a:\t%g" % xopt[0] + print "b:\t%g" % xopt[1] + def model(self, x, params) : + a = params[0] + b = params[1] + return self.binwidth*self.N*a*exp(b*x + a/b*(1.-exp(b*x))) + def guessInitialParams(self) : + print "guessing initial params" + init = self.data[0] + # find the coords (Xmax,Ymax) of the peak of the histogram + Yarr = list(zip(*self.data)[1]) + Xarr = list(zip(*self.data)[0]) + Ymax = max(Yarr) + Xmax = Xarr[Yarr.index(Ymax)] + # two guessing methods for B: + if Ymax < 1.3*init[1] : + # in the first, we have a high unforced rate, and + # the histogram is more-or-less exponential decay. + print "guessing with method 1" + Aguess = init[1]/self.N/self.binwidth + i = 0 + s = 0 + while s < self.N/2.0 : + s += Yarr[i] + i += 1 + print "half dead by", self.data[i] + Aguess = init[1]/self.N/self.binwidth + Bguess = 1.0/self.data[i][0] + else : + # in the second, we have a low unforced rate, and the histogram has + # a well defined 'bump' at some likely unfolding force + # h(f) = binwidth N a e^{b x + a/b[1-e^(b f)]} + # dh/df = h [b - a e^(b f)] + # so dh/df = 0 at f=fmax when + # b = a e^(b fmax) + # hmm, trancendental equation for b + print "Xmax = ", Xmax + Bguess = 1.0/Xmax + Aguess = init[1]/self.N/self.binwidth + print "guessing ", (Aguess, Bguess) + return (Aguess, Bguess) + +if __name__ == "__main__" : + fit.run(2, Model) diff --git a/testing/bell_rate/run_test.sh b/testing/bell_rate/run_test.sh new file mode 100755 index 0000000..931eeba --- /dev/null +++ b/testing/bell_rate/run_test.sh @@ -0,0 +1,103 @@ +#!/bin/bash +# +# Use a single hooke domain with a bell unfolding rate, so unfolding force +# is proportional to time, and we expect a peaked death histogram. +# +# e.g. +# spring constant k=df/dx=1 \__ +# velocity v=dx/dt=1 / `--> df/dt = kv = 1 --> f = t (+ f0) +# unfolding rate K= a e^(bf) (frac/s) +# let f0 = 0 +# let P(f) be the population existing at force f (normed so P(0)=0). +# the probability of dying (normed # deaths) between f and f+df is given by +# dD = K P(f) * dt = K P(f) df/kv since dt = df/kv +# this is what whe plot in the normalized histogram with dD in a bin of width +# df. The continuous analog is the function dD/df = K/kv P(f). +# +# dP/dt = - KP +# dP/P = -K dt = -a e^{bf} df/kv = -z e^{bf} df where z ≡ a/kv +# ln(P) = -z/b e^{bf) + c +# P = C e^{-z/b e^{bf}} +# Applying our f0=0 boundary condition we see that C = P0 e^{z/b}. +# Plugging P(f) into our formula for dD/df yields +# dD/df = K/kv P(f) = z e^{bf} P0 e^{z/b} e^{-z/b e^{bf}} +# dD/df = P0 z e^{bf + z/b(1-e^{bf})} +# +# The histogram scaling isn't dD/df, though, it's dD/d(bin) +# dD/d(bin) = dD/df df/d(bin) = binwidth P0 z e^{bf + z/b(1-e^{bf})} +# +# test with lots of runs of v-dx=1 +# sawsim -v1 -mhooke -a1 -kbell -K1,1 -T7.24296e22 -F1 +# ^- v=1 ^- k=1 ^-a=1 ^-1/kB +# (kB = 1.3806503 10-23 J/K, we use T=1/kB to get the bell K=e^{-f}) +# +# 70 runs of that takes ~1 second: +# time xtimes 200 sawsim -v1 -mhooke -a1 -kconst -K1 -F1 > /dev/null +# +# usage: bell_rate.sh + +# since we're on the cluster, we'll go crazy and get really good statistics ;) +N=10000 + +SAWSIM=../../bin/sawsim +DATA=bell_rate.d +HIST=bell_rate.hist +GPFILE=bell_rate.gp +PNGFILE=bell_rate.png + +# set up the theory +# dD/d(bin) = binwidth P0 z e^{bf + z/b(1-e^{bf})} +# where z ≡ a/kv, and b = Δx/kBT +KB=1.3806503e-23 # J/K +STYLE="protein" +if [ "$STYLE" == "ones" ] +then + df=1e-1 # histogram bin width + DX=1 + T=7.24296e22 # 1/kB + A=1 + K=1 + V=1 +else + df=1e-12 + DX=1e-9 + T=300 + A=1e-3 + K=0.05 + V=1e-6 +fi +B=`echo "print ($DX/($KB*$T))" | python` +Z=`echo "print (float($A)/($K*$V))" | python` +THEORY_FN="$df * $N * $Z * exp($B*x + $Z/$B *(1 - exp($B*x)))" + +# run the experiments +> $DATA +if [ "$ISACLUSTER" -eq 1 ] +then + qcmd "xtimes $N $SAWSIM -v$V -mhooke -a$K -kbell -K$A,$DX -T$T -F1 | grep -v '^#' >> $DATA" +else + xtimes $N $SAWSIM -v$V -mhooke -a$K -kbell -K$A,$DX -T$T -F1 | grep -v '^#' >> $DATA +fi + +# histogram the data +cat $DATA | stem_leaf -v-1 -b$df | sed 's/://' > $HIST + +# fit the data +FIT=`cat $HIST | ./fit_bell` +#echo "$FIT" +FITZ=`echo "$FIT" | sed -n 's/a:[[:space:]]*//p'` +FITB=`echo "$FIT" | sed -n 's/b:[[:space:]]*//p'` +echo -e "a:\t$FITZ\t(expected $Z)" +echo -e "b:\t$FITB\t(expected $B)" + +# generate gnuplot script +echo "set terminal png" > $GPFILE +echo "set output '$PNGFILE'" >> $GPFILE +echo "set logscale y" >> $GPFILE +echo "set xlabel 'Time'" >> $GPFILE +echo "set ylabel 'Deaths'" >> $GPFILE +echo "theory(x) = $THEORY_FN" >> $GPFILE +echo "fit(x) = $df * $N * $FITZ * exp($FITB*x + $FITZ/$FITB*(1-exp($FITB*x)))" >> $GPFILE +echo "plot theory(x), '$HIST' with points, fit(x)" >> $GPFILE + +gnuplot $GPFILE diff --git a/testing/common/fit.py b/testing/common/fit.py new file mode 100755 index 0000000..1f0aa12 --- /dev/null +++ b/testing/common/fit.py @@ -0,0 +1,89 @@ +# -*- coding: utf-8 -*- + +"""Fit histogram data coming in on stdin to a user-defined model. You +should write a python script sublcassing fit-hist.Model.""" + +import scipy.optimize, scipy +import sys + +kB=1.3806503e-23 # m²·kg/s²·K + +class Model (object) : + + """For minimizing fit error: returns the normalized rms difference + between a dataset and a function y(X). In order to use this + class, you must sublcass it and overwrite Model.model(x,params). + It is recommended that you also overwrite + Model.guess_initial_params(), and then not pass initializing + params to Model.model(). This helps ensure good automation for + unsupervised testing. Developing good heuristics for + Model.guessInitialParams() is usually the hardest part writing a + theoretical histogram test. + + For bonus points, you can override Model.printModel() and + Model.printParams(), to give your users an easy to understand idea + of what's going on. It's better to be verbose and clear. + Remember, this is the testing branch. It needs to be solid. Save + your wonderful, sneaky hacks for the main code ;). + + """ + def __init__(self, numParams, data) : + """Initialize the Model instance with a set of data, for + example: data = [(0,4),(1,3),(2,3),(4,0)]""" + self.data = data # an iterable filled with (x,y) pairs + # N (# unfolding events) = the sum of the y values (# hits per bin) + self.N = sum(zip(*data)[1]) # zip(* takes data to [[x0,x1,..],[y..]] + # assumes constant binwidth + self.binwidth = self.data[1][0] - self.data[0][0] + self.numParams = numParams + def nrms(self, params) : + """Returns the normalized root mean squagred between the data + and a model for a given set of parameters""" + s = 0.0 + for x,y in self.data : + yResidual = y-self.model(x, params) + s += yResidual**2 + s /= float(len(self.data)) + return s + def fit(self, x0=None) : + """Minimize error given an initial guess vector x0, the heart of + the class.""" + if x0 == None : + x0 = self.guessInitialParams() + xopt = scipy.optimize.fmin_powell(self.nrms, x0) + return xopt + def printModel(self) : + "OVERRIDE: Give the user some idea of what model we're fitting with" + print "model string not defined" + def printParams(self) : + "OVERRIDE: Give the user nicely named format for a parameter vector" + for i,p in zip(range(len(params)), params) : + print "param_%d\t=\t%g\n" % (i,p) + def model(self, x, params) : + "OVERRIDE: Implement your model" + assert len(params) == self.numParams + raise Exception, "model method not defined" + return 1.0 + def guessInitialParams(self) : + """OVERRIDE: Provide some intelligent heuristic to pick a + starting point for minimizing your model""" + return [1.0] * self.numParams + +def readData(fileObj) : + """Read histogram data from a whitespace-seperated, two-column + ASCII text-file into a list of tuples""" + data = [] + for line in fileObj : + if line[0] == "#" : continue + data.append([float(x) for x in line.split()]) + return data + +def run(numParams, modelClass) : + import sys + data = readData(sys.stdin) + mod = modelClass(numParams, data) + popt = mod.fit() + mod.printModel() + mod.printParams(popt) + print "N:\t%g" % mod.N + print "bin:\t%g" % mod.binwidth diff --git a/testing/common/fit.pyc b/testing/common/fit.pyc new file mode 100644 index 0000000000000000000000000000000000000000..396d9d25504bf2b1931414ce0d9ad0fff0a2b2ea GIT binary patch literal 4361 zcma)9ZFAek5k8Q5F%`LS?9^j7NqaMCVydNRDoz_$88?k3*Vcz94%Dd|r5+9h97$Ms zp>YR^^4R?l=VO0Ue^mPi`aFA($hec4Qo`{7-0kf?yZh`i|M%s!-T(aSQLLIz8{dDz zWe+eUN*!Vh)N!B=TQX^#Pul9RtrFZWsidvUvKj>{SyDeJ^;)TK4_8#O%$vhim8{6* z689f;S1J4@hCioHPg3j3aTVn{iCm=PGEa-KP77TY+PVZ2SL#UD)>NB`8Ks3ubY3PV z>uddOS+ChiS!anptx{*OZhGb>+!`6EPeo3Fki5s`5p(D z?{Vm;WK{rJQ{h4P5+x7YfTOH5tWWdwMWe7bRaI8ox-zb=3aj13=%TFhC`(@ehbp%^ zNk=17nIbl88lNNMPK_ym1mU#C>8MC>H>!)+r9kGa`?n21w+Dc$#8gBV*g|q@ML?VZ zubw_DYn|724xxsxD5^X$l?U77D8A^Sa9vgdSap*6gjn<9a#ZW1{W8fCVe z3z5}c(6~0%4q?7aqipEKbvJr)DC{m^xYJW8hv1Mkf|3aEi)qTDPBil`FiTHNHnrLm zwyumGmnHD1U71IM$Ed7CMvd1~Q=O!i+8XQN0tnbOCnhVW^oQ)6n7Tp`rLooLI~Ua7 z^O`*tEV>&6p|*9DWoP2LNd(A&DfD#`u@oE0@vyLjm(w`P7E*$G^w(qogX5?yYOAMZ zTDZnVaRfiYBUPH1#`tN4*P;^?fX8IF*HK=bE7v1lWb3vCe<@LOda zeGcJX%W|Yo;mBJSU{Lq6U?}k(ap!zA9T^RUj!Fa-_r{rtDjyZ0$xVJ_Dx6HW;Ggl_ zV6Bg;sE8*bPZ!1{cw0iOwq=$kczh6@_z)L$PfMyds|n>);ScY}ME z)l0t?Dc=UaEm8IuLlX|jRnbYDfaoQtvC*fgn;~UDzz`%>Zkza$;6s$1jC?Goy+$NFe7*DJ2KBv zTAZedo45{ciaKw?gBwEfUP!hMxxS0b5AF(P{r+vHcUjy9v-ZKels1$E9u9wrxn0Ac zS}N$YZni!QVoF=VMSd0dBTP^bUbWTib3wVKcApN^-&^X{vU&yJJ@$ZC5PscKH(tS8 z;v|57?KA9renmaO6C4f$=$n>;NX?be0EaYwxIKz@Iaqt*lLoppcALOEHMda5a41WZu!y|@F%m&&nLz(zsg+H)wC`1>CSuE*7BFb$HRChf zS@H)|_a@LnJ`fYWz+`rpJrH8oGp0{B#3ekx`}93YUnI3ON94tp8nx8*hmW4Wub%vW zh8JmUP%#VSQo=MtrEaqABEgzqv5b9rIQ%H--0OTaUfcJXL5;=fs48OW{ebUDT5qDm09h!T+qkutFc6;ikvSJdg~2Eg{RR`@1G|BY&le=mbXAh&d*%#U z#6W!yfEWf;)biwc80?Yu(KkuMgqD!hu;5}a%mh6Q+o8tAqxM#R&=`JUsSjMTx1-@- zTErAPo(F6xg_h>puuO^BBLdo=3sDu{S+X?>2$X zayYl}aaEqMWAPzgIHX=Wrok8O6^s8gjgOo5h~$g12?G$^(f~dL5W>r&{Q+V;Be@eBf}x2*BNgK>TsBp@o*S2sf9MecRBk7hxa)A ziUU!T6XtNrOcK z4v4cIAxR_7o?jMjM)?F+W5W{28jC>Gues@(5H!*kd@VtG_4$4E4=I}5`ULx$urxR@ zG*Pn9k+K}>QB~$T0w>`&mz_qj*|Y||bL30{?mjIy<1(xBLYVu@gPpxSG{SIztmpUUM~K-KIFmr{;(I@_R@c4#An5zCL;;q@3Z`pVY``I>pZlc{coR>#1k= zX_h)Ed?AZlm(YbDbJ*Y@fqjcJM#9`3zBXRMvkI2(!n{YAFVx3XwGzCeR)f3q4yKJ2 zJ6tuwkSKKwr=WWjvhxiHyBA3TfXRv|JxsbPl@lIA$?R|*aJEM0fGYxck@s(rE#Xbh z*JZAnaXEyvTZ@5l=J$oP-zmgj}{B-h`6?Rdw+VzC4oPinsB%9;^kIF-rXoe~s6| zO;VWq=S4ZZj^}cURQTij;_~w78TtcB88I>3!YZ%Y9j@oCe` Qt+cM9JzZ&Czt`#f7iD_5RsaA1 literal 0 HcmV?d00001 diff --git a/testing/common/fitHist.pyc b/testing/common/fitHist.pyc new file mode 100644 index 0000000000000000000000000000000000000000..91c493848c7d2f3621d3ec8c49ffdf30827a2282 GIT binary patch literal 4362 zcma)9+j1Mn5uF8i5oj`&C5I)&b~06tBQOCW&~jphF2%7#MdE`QWtoa%Ix5t1cL1!p z7iede1XWTVEI;N0@-z8>{7$|g=kzQISSeQpLgG{Et#|~CT(@pRtavGRMJ*vS&agfEUE95dacw~M=L5>=FQQn zN>*fYh5HY?s}%kc!yi(oC#iMixQcR}L@v^CnWx2Ar-d#HZC!$iD|Mu6YpU(UjMBm+ zIxiEG^|k)0tk>+Mtg}R)RjD&rH$8U~?u@H+>a?wovpBMrB1fs)=3)IWW%t2obEkjj z?%+4S9^4r;pB=OP@g27e&16lCiF{fZ@8hyhFgREQJ246zG*D3BWQiyMi|+sy%oi+Z zzQ+OPdmK6{SrtImRCw6ELdnB6;3z8%>(e}a(I~70$Gh+%MK{%~(Iw}&}jp`zHDUdnu{&fq`?E~N{F%{7Twve1! z5fEp^q%WXbo@XrbwWwn?z28 zMwu<=LS(fUG_H-cLzwTwfl-L0M+3cCv!?(`JOAvk1>pdw#GWR00Q>RsmaPI{UJLirmhe~X>9ey z&IR@Ny=G4Yi|!UdsBIl(*}1rG5&?2x3VoeKEX780JS;5XVuu%11?La+4pM3MbPo z_-8yfSnJ~|D&mRA(}gh!-j)!nZJDJB9v?=hKEy@cvy$qK>a3@2VWN|BJ&EEIaf(RI zBQQsdOxK?wCIy2SDW)#22TZyU!^@Te6%CQsE#=z0{fF?~&DVjvSmK!aN8q}f%^haE zyTQH7>ZM8*DvG1DpjFl{ye^byR@zWMWM=f?-gp&gJ_Bn2|dU z9GT}hEzZ)!O6Y4iHc)?SsaMPD6~OoS17bn> zbxYlN1&@iN;DGC&V(0U#>M5S!a2P@Vv=p3&|6RPqi^d}vjtYN(NK0@_ZV>(Uc^!=_ zP{0L|GDL}`gb>}rppJYAqPs9>5@DzJBXT>+4RXPb+JpW8y3j3zK|r@2RN^W3Ga1jN z=;SO+j_`;OD=s+?jZG79h%6apW#q(*nJKzl>gD~2v$r|$q;Q?XTNvg(3+WmOhyrEf z7H;fS4D`)Pu%_07&EQ&a(|d@{=k)-W{VfK>H=OiMOO;nuK!n>Y@vj0^1S){gbXlb< z0{tr2!b@~g3&_)aDVtw{wNcoNFKG>C0K6e{E%XoOh#fGU4agA&%uS{&N0;F9p;V~; zL6-rGq*ZD*5W5b_!U!7f8Xk(g@A0t!XyddA({k!WAxE4>sJp{yd1kUK z@EBj`>?VdcLFE$UzyM#D@xBpU4Q>Q6{d{o=np7H4K6;GcJrv>gJ|6w?FAu|TfA8V0 zendv5ne2E#ez_6)B}I^Wo4`9Ww@}M)D9e;EiM;_a5wGj`JMfu7jm7DxEMpVC zgWHgcCFN2Xrwji5fbU6KZ=%BhSty;mxV2X>5T5suIhTcn!6*>@5)ed@qGX{fWjWNNs?2o+PQvdlJB?zqZ4LV8$e9G(eO7MAWme~fF!vXS5BB%b z=*`^rD9u1<1w0$+3{3-h&DDr?o8lm!njap>Zz5?rgowFGCK-PINr8-}Q~a)wXTz^| zo_T7YW~q~^7qYT-iCg#!4qF@~s&8?|Fqr$o7se}iR>90Y*!CFnh2prXR)TlbYH)Yn zzO=F80auMMBubsYBWNFm>3sLW)w*@U}LFFTwEoUPF_Ac_Dkmx5c56Wxeb5C8?KX|&mKhvuO!9;{=Ac(PS)d@ zgsUJqMZH^yTfvVlhG$cr7xp>yII!A@?C>^bZ^B7{s=D|SzC4m3insB%9;^i$7^QfJ zKgVk!CaKE(i;5g>;<=n6<^ANoxV$|62F-!wi df/dt = kv --> f = kvt (+ f0) +# unfolding rate K=1 frac/s +# +# Population as a function of force/time follows +# p(t) = exp(-tK) = exp(-FK/kv) = p(F) +# so histogram of unfolding vs. force p(F) normalized to p(0)=1 should follow +# -binwidth N0 dp/dF = binwidth N0 pK/kv = binwidth N0 K/kv exp(-FK/kv) +# +# test with lots of runs of +# sawsim -v1 -mhooke -a1 -kconst -K1 -F1 +# which is not a problem, because 200 runs of that takes ~1 second: +# time xtimes 200 sawsim -v1 -mhooke -a1 -kconst -K1 -F1 > /dev/null +# +# usage: const_rate.sh + +if [ $# -ne 3 ] +then + echo "usage: const_rate.sh " + exit 1 +fi + +V=$1 +K=$2 +RATE=$3 +DEBUG=0 + +# since we're on the cluster, we'll go crazy and get really good statistics ;) +N=10000 +df=`python -c "print 0.1/$RATE"` + +SAWSIM=../../bin/sawsim +DATA=const_rate_${V}_${K}_$RATE.d +HIST=const_rate_${V}_${K}_$RATE.hist +COMPFILE=const_rate_${V}_${K}_$RATE.compare +GPFILE=const_rate_${V}_${K}_$RATE.gp +PNGFILE=const_rate_${V}_${K}_$RATE.png + +# set up the theory +# our predicted histogram (see notes above) is +# hist(F) = binwidth N0 K/kv exp(-FK/kv) = A exp(F/B) +theoryA=`python -c "print $df*$N*$RATE/float($K*$V)"` +theoryB=`python -c "print float(-$K*$V)/$RATE"` +THEORY_FN="$theoryA * exp(x / ($theoryB))" + + +# run the experiments +> $DATA +if [ "$ISACLUSTER" -eq 1 ] +then + if [ $DEBUG -eq 1 ]; then + echo "qcmd \"xtimes $N $SAWSIM -v$V -mhooke -a$K -kconst -K$RATE -F1 | grep -v '^#' >> $DATA\"" + fi + qcmd "xtimes $N $SAWSIM -v$V -mhooke -a$K -kconst -K$RATE -F1 | grep -v '^#' >> $DATA" +else + if [ $DEBUG -eq 1 ]; then + echo "xtimes $N $SAWSIM -v$V -mhooke -a$K -kconst -K$RATE -F1 | grep -v '^#' >> $DATA" + fi + xtimes $N $SAWSIM -v$V -mhooke -a$K -kconst -K$RATE -F1 | grep -v '^#' >> $DATA +fi + +# histogram the data +if [ $DEBUG -eq 1 ]; then + echo "cat $DATA | stem_leaf -v-1 -b$df | sed 's/://' > $HIST" +fi +cat $DATA | stem_leaf -v-1 -b$df | sed 's/://' > $HIST + +# fit the data to an exponential decay +if [ $DEBUG -eq 1 ]; then + echo "cat $HIST | ./fit_exponential" +fi +FIT=`cat $HIST | ./fit_exponential` || exit 1 +B=`echo "$FIT" | sed -n 's/time constant:[[:space:]]*//p'` +A=`echo "$FIT" | sed -n 's/initial value:[[:space:]]*//p'` +echo -e "force scale:\t$B\t(expected $theoryB)" > $COMPFILE +echo -e "initial value:\t$A\t(expected $theoryA)" >> $COMPFILE +cat $COMPFILE + +# generate gnuplot script +if [ $DEBUG -eq 1 ]; then + echo "generate Gnuplot script" +fi +echo "set terminal png" > $GPFILE +echo "set output '$PNGFILE'" >> $GPFILE +echo "set logscale y" >> $GPFILE +echo "set xlabel 'Force'" >> $GPFILE +echo "set ylabel 'Deaths'" >> $GPFILE +echo "theory(x) = $THEORY_FN" >> $GPFILE +echo "fit(x) = $A * exp(x / $B)" >> $GPFILE +echo "plot theory(x), '$HIST' with points, fit(x)" >> $GPFILE + +if [ $DEBUG -eq 1 ]; then + echo "Generate PNG" +fi +gnuplot $GPFILE diff --git a/testing/const_rate/fit_exponential b/testing/const_rate/fit_exponential new file mode 100755 index 0000000..7506812 --- /dev/null +++ b/testing/const_rate/fit_exponential @@ -0,0 +1,35 @@ +#!/usr/bin/python +# +# fit data coming in on stdin to an exponential decay +# output: +# +# initial value:\t +# time constant:\t +# +# usage: fit_exponential +# for example: +# $ cat datafile | fit_exponential + +import sys +sys.path.append('../common/') +import fit +from scipy import exp + +class Model (fit.Model) : + def printModel(self) : + print "a*exp(x/b)" + def printParams(self, params) : + print "initial value:\t%g" % params[0] + print "time constant:\t%g" % params[1] + def model(self, x, params) : + a = params[0] + b = params[1] + return a*exp(x/b) + def guessInitialParams(self) : + init = self.data[0] + fin = self.data[-1] + slope = (fin[1]-init[1])/(fin[0]-init[0]) + return (init[1],slope) + +if __name__ == "__main__" : + fit.run(2, Model) diff --git a/testing/const_rate/params b/testing/const_rate/params new file mode 100644 index 0000000..150a6f0 --- /dev/null +++ b/testing/const_rate/params @@ -0,0 +1,5 @@ +# velocity spring constant unfolding rate +1 1 1 +2 1 1 +1 2 1 +1 1 2 diff --git a/testing/const_rate/run_test.sh b/testing/const_rate/run_test.sh new file mode 100755 index 0000000..2407efc --- /dev/null +++ b/testing/const_rate/run_test.sh @@ -0,0 +1,13 @@ +#!/bin/bash +# +# run const_rate.sh with a variety of parameters +# +# usage: run_test.sh + +while read LINE +do + echo ./const_rate.sh $LINE + ./const_rate.sh $LINE +done < <(cat params | grep -v '^#') + +exit 0 diff --git a/testing/kramers_rate/fit_kramers b/testing/kramers_rate/fit_kramers new file mode 100755 index 0000000..faa9ee4 --- /dev/null +++ b/testing/kramers_rate/fit_kramers @@ -0,0 +1,46 @@ +#!/usr/bin/python +# -*- coding: utf-8 -*- +# +# fit data coming in on stdin to Kramers' k(F) over a cusp. +# k ≅ ω²/2πγ √(πβE_b) e^(-βE_b) +# E_b(F) = ½ω²Δx² - FΔx - (½+1/ω²)F² +# fitting gamma, omega, T, and Δx +# output: +# +# gamma:\t +# omega:\t +# temp:\t +# delta x:\t +# +# usage: fit_kramers +# for example: +# $ cat datafile | fit_kramers + +import sys +sys.path.append('../common/') +import fit +from scipy import pi, sqrt, exp + +kB=1.3806503e-23 # m²·kg/s²·K + +class Model (fit.Model) : + def printModel(self) : + print "k ≅ ω²/2πγ √(πβE_b) e^(-βE_b) with E_b(F) = ½ω²Δx² - FΔx - (½+1/ω²)F²" + def printParams(self, xopt) : + print "gamma:\t%g" % xopt[0] + print "omega:\t%g" % xopt[1] + print "temp:\t%g" % xopt[2] + print "delta x:\t%g" % xopt[3] + def model(self, x, params) : + F = x + gamma = params[0] + omega = params[1] + T = params[2] + dX = params[3] + betaEb = (0.5*(omega**2)*(dX**2) - F*dX - (0.5 + 1.0/omega**2)*F**2)/(kB*T) + return (omega**2/(2*pi*gamma) * sqrt(pi*betaEb) * exp(-betaEb)) + def guessInitialParams(self) : + return [1, 1, 1.0/kB, 1.0] + +if __name__ == "__main__" : + fit.run(4, Model) diff --git a/testing/kramers_rate/run_test.sh b/testing/kramers_rate/run_test.sh new file mode 100755 index 0000000..cb52a72 --- /dev/null +++ b/testing/kramers_rate/run_test.sh @@ -0,0 +1,115 @@ +#!/bin/bash +# -*- coding: utf-8 -*- +# +# Test the Kramers integral k(F) implementation on the cusp potential. +# +# The test is escape from a cusp potential from Kramers' 1940 paper. +# He solves the cusp potential in the text immediately after Eq. 17 & Fig. 2 +# We can translate Kramers' notation to Hänggi's more modern notation with +# Kramers Hängii Meaning Why I think so +# 2πω ω Curvature at minimum Hänggi's footnote to Eq. 1.5 +# η γ Viscosity/damping Compare H. Eq.1.5 with K. Eq. 25 +# Q E_b Energy barrier height K. Fig. 1 (p 291) +# T k_B·T≡1/β Absolute temperature K: “k_B = 1” (text after Eq. 5) +# r k Reaction rate seems more reasonable than k ;) +# With these translations, Kramers' Eq. is +# k ≅ ω²/2πγ √(πβE_b) e^(-βE_b) +# for E(x) = ½ω²(|x|-Δx)². +# The barrier is at x_b=0 with minimums at x_ac=±Δx. +# The barrier height E_b = E(0) = ½ω²Δx² +# +# As a force is applied to tip the potential, +# E(F,x) = ½ω²(|x|-Δx)²-Fx. +# Focusing only on the area left of the barrier (x<0), we have +# E_L(F,x) = ½ω²(-x-Δx)²-Fx = ½ω²(x+Δx)²-Fx +# E_L'(F,x) = ω²(x+Δx)-F +# E_L"(F,x) = ω², +# just like the F=0 case. Kramers' solution still holds then, but with a +# different barrier height. The minimum moves to +# x_a = F/ω² - Δx. +# so +# Eb(F) = E(F,0) - E(F,x_a) = E(0,0) - E(F,x_a) +# = ½ω²Δx² - ½ω²[(F/ω² - Δx)+Δx]²-F(F/ω² - Δx) +# = ½ω²Δx² - ½F² - F²/ω² - FΔx +# = ½ω²Δx² - FΔx - (½+1/ω²)F² +# ----+----- ^- reduction due to decreasing barrier distance +# ^- Bell model +# +# Usage: kramers_rate.sh + +# Since we're on the cluster, we'll go crazy and get really good statistics ;) +FMAX=1 +K_UTILS=../../bin/k_model_utils +DATA=kramers_rate.d +HIST=kramers_rate.hist +GPFILE=kramers_rate.gp +PNGFILE=kramers_rate.png + +# set up the physics +# E(x) = ½ω²(|x|-Δx)² +OMEGA=3 +GAMMA=100 +KB=1.3806503e-23 # m²·kg/s²·K +T=`python -c "print (1 / ($KB))"` # pick units where k_B = 1 +DX=1 +EB=`python -c "print (0.5 * (($OMEGA)*($DX))**2)"` + +# set up the theory +E_FN="0.5*($OMEGA)**2 * (abs(x)-($DX))**2" +BETA=`python -c "print (1.0 / (($KB)*($T)))"` +# k ≅ ω²/2πγ √(πβE_b) e^(-βE_b) +BETAEB=`python -c "print (($BETA)*($EB))"` +THEORY_K="BETAEB(F) = ($BETA)*(($EB) - (0.5+1/($OMEGA)**2)*F**2 - F*($DX)) +K(F) = ($OMEGA)**2/(2*pi*($GAMMA)) * sqrt(pi*BETAEB(F)) * exp(-BETAEB(F))" + +# Get the diffusion coefficient from the Einstein-Smoluchowski relation +# D = µ_p·k_B·T = k_B·T/γ +# where µ_p=1/γ is the particle mobility +# See, for example, +# http://en.wikipedia.org/wiki/Einstein_relation_(kinetic_theory) +D=`python -c "print (1/(($BETA)*($GAMMA)))"` +# pick safe integration bounds +XMIN=`python -c "print (-2*($DX))"` +XMAX="$DX" +# set up the cubic-spline approximation to a sharp cusp… get close anyway… +SPLINE_PTS=`python -c "from sys import stdout +dx = (($XMAX)-($XMIN))/100.0 +stdout.write('{') +for i in range(100) : + x = ($XMIN)+i*dx + stdout.write('(%g,%g),' % (x, ($E_FN))) +stdout.write('(%g,%g)}\n' % (($XMAX), ($E_FN)))"` + +# Run the simulation +echo "$K_UTILS -kkramers_integ -K \"$D,$SPLINE_PTS,$XMIN,$XMAX\" -x$XMIN -X$XMAX -T$T -F$FMAX > $DATA" +$K_UTILS -kkramers_integ -K "$D,$SPLINE_PTS,$XMIN,$XMAX" -x$XMIN -X$XMAX -T$T -F$FMAX > $DATA + +# fit the data +FIT=`cat $DATA | ./fit_kramers` +FITGAMMA=`echo "$FIT" | sed -n 's/gamma:[[:space:]]*//p'` +FITOMEGA=`echo "$FIT" | sed -n 's/omega:[[:space:]]*//p'` +FITT=`echo "$FIT" | sed -n 's/temp:[[:space:]]*//p'` +FITDX=`echo "$FIT" | sed -n 's/delta x:[[:space:]]*//p'` +echo -e "gamma:\t$FITGAMMA\t(expected $GAMMA)" +echo -e "omega:\t$FITOMEGA\t(expected $OMEGA)" +echo -e "temp:\t$FITT\t(expected $T)" +echo -e "delta x:\t$FITDX\t(expected $DX)" + +FITEB=`python -c "print (0.5 * (($FITOMEGA)*($FITDX))**2)"` +FITE_FN="0.5*($FITOMEGA)**2 * (abs(x)-($FITDX))**2" +FITBETA=`python -c "print (1.0 / (($KB)*($FITT)))"` +FITBETAEB=`python -c "print (($FITBETA)*($FITEB))"` +FIT_K="FITBETAEB(F) = ($FITBETA)*(($FITEB) - (0.5+1/($FITOMEGA)**2)*F**2 - F*($FITDX)) +fitK(F) = ($FITOMEGA)**2/(2*pi*($FITGAMMA)) * sqrt(pi*FITBETAEB(F)) * exp(-FITBETAEB(F))" + +# generate gnuplot script +echo "set terminal png" > $GPFILE +echo "set output '$PNGFILE'" >> $GPFILE +echo "set logscale y" >> $GPFILE +echo "set xlabel 'Time'" >> $GPFILE +echo "set ylabel 'Deaths'" >> $GPFILE +echo "$THEORY_K" >> $GPFILE +echo "$FIT_K" >> $GPFILE +echo "plot K(x), fitK(x), '$DATA' with points" >> $GPFILE + +gnuplot $GPFILE diff --git a/testing/test_cluster.sh b/testing/test_cluster.sh new file mode 100755 index 0000000..d76d77f --- /dev/null +++ b/testing/test_cluster.sh @@ -0,0 +1,11 @@ +#!/bin/bash +# +# Check to make sure the PBS queue is working before going +# whole hog on the simulations. + +if [ "$ISACLUSTER" -eq 1 ] +then + qcmd printenv +fi + +exit 0 -- 2.26.2