From: W. Trevor King Date: Tue, 9 Sep 2008 17:14:07 +0000 (-0400) Subject: Added testing scripts, comparing sawsim with theory X-Git-Tag: v0.5~5 X-Git-Url: http://git.tremily.us/?a=commitdiff_plain;h=780d75e14fb1a05686db24bb0599d214a938dd21;p=sawsim.git Added testing scripts, comparing sawsim with theory TODO, adjust bell & kramers testing to read param files like const_rate does. --- 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 0000000..396d9d2 Binary files /dev/null and b/testing/common/fit.pyc differ diff --git a/testing/common/fitHist.pyc b/testing/common/fitHist.pyc new file mode 100644 index 0000000..91c4938 Binary files /dev/null and b/testing/common/fitHist.pyc differ diff --git a/testing/const_rate/const_rate.sh b/testing/const_rate/const_rate.sh new file mode 100755 index 0000000..5b59829 --- /dev/null +++ b/testing/const_rate/const_rate.sh @@ -0,0 +1,101 @@ +#!/bin/bash +# +# Use a single hooke domain with a constant unfolding rate, so unfolding force +# is proportional to time, and we expect an exponential population decay. +# +# e.g. +# spring constant k=df/dx \__ +# velocity v=dx/dt / `--> 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