--- /dev/null
+# Copyright (C) 2008-2010 W. Trevor King <wking@drexel.edu>
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program. If not, see <http://www.gnu.org/licenses/>.
+#
+# The author may be contacted at <wking@drexel.edu> on the Internet, or
+# write to Trevor King, Drexel University, Physics Dept., 3141 Chestnut St.,
+# Philadelphia PA 19104, USA.
+
+"""Test a Hookian chain with constant unfolding rate.
+
+With the constant velocity experiment and a Hookian domain, the
+unfolding force is proportional to time, so we expect exponential
+population decay, even with multiple unfolding domains.
+
+For example, with a spring constant
+
+.. math:: k = df/dx
+
+and a pulling velocity
+
+.. math:: v = dx/dt
+
+we have a loading rate
+
+.. math:: df/dt = df/dx dx/dt = kv \;,
+
+so
+
+.. math:: f = kvt + f_0 \;.
+
+With an unfolding rate constant
+
+.. math:: K = 1 \text{frac/s}
+
+The population follows
+
+.. math::
+ dp/dt = Kp = 1 \text{s^{-1}}
+ p(t) = exp(-tK) = exp(-(f-f_0)K/kv) = p(f) \;.
+
+Therefore, a histogram of unfolding vs. force :math:`p(f)` normalized
+to :math:`p(0)=1` should follow
+
+.. math::
+ -w N dp/dF = w N pK/kv = w N K/kv exp(-FK/kv)
+
+where :math:`w` is the binwidth and :math:`N` is the number of tested
+domains.
+"""
+
+from numpy import arange, exp, log
+
+from ..histogram import Histogram
+from ..fit import HistogramModelFitter
+from ..manager import get_manager
+from ..sawsim import SawsimRunner
+from ..sawsim_histogram import sawsim_histogram
+
+
+def probability_distribution(x, params):
+ """Exponential decay decay probability distribution.
+
+ .. math:: 1/\\tau \cdot exp(-x/\\tau)
+ """
+ p = params # convenient alias
+ p[0] = abs(p[0]) # cannot normalize negative tau.
+ return (1.0/p[0]) * exp(-x/p[0])
+
+
+class ExponentialModelFitter (HistogramModelFitter):
+ """Exponential decay model fitter.
+ """
+ def model(self, params):
+ """An exponential decay model.
+
+ Notes
+ -----
+ .. math:: y \\propto \cdot e^{-t/\tau}
+ """
+ self._model_data.counts = (
+ self.info['binwidth']*self.info['N']*probability_distribution(
+ self._model_data.bin_centers, params))
+ return self._model_data
+
+ def guess_initial_params(self):
+ return [self._data.mean]
+
+ def guess_scale(self, params):
+ return [self._data.mean]
+
+ def model_string(self):
+ return 'p(x) = A exp(-t/tau)'
+
+ def param_string(self, params):
+ pstrings = []
+ for name,param in zip(['tau'], params):
+ pstrings.append('%s=%g' % (name, param))
+ return ', '.join(pstrings)
+
+
+
+def constant_rate(sawsim_runner, num_domains=1, unfolding_rate=1,
+ spring_constant=1, velocity=1, N=400):
+ loading_rate = spring_constant * velocity
+ tau = loading_rate / unfolding_rate
+ w = 0.1 * tau # calculate bin width (in force)
+ A = w*num_domains*N / tau
+ theory = Histogram()
+ # A exp(-x/tau) = 0.001
+ # -x/tau = log(0.001/A)
+ # x = -log(0.001/A)*tau = log(A/0.001)*tau = log(1e3*A)*tau
+ theory.bin_edges = arange(start=0, stop=log(1e3*A)*tau, step=w)
+ theory.bin_centers = theory.bin_edges[:-1] + w/2
+ theory.counts = w*num_domains*N*probability_distribution(
+ theory.bin_centers, [tau])
+ theory.analyze()
+
+ param_string = (
+ '-v %(velocity)s -s cantilever,hooke,%(spring_constant)s -N1 '
+ '-s folded,null -N%(num_domains)s -s unfolded,null '
+ '-k folded,unfolded,const,%(unfolding_rate)s -q folded '
+ ) % locals()
+ sim = sawsim_histogram(sawsim_runner, param_string, N=N,
+ bin_edges=theory.bin_edges)
+
+ e = ExponentialModelFitter(sim)
+ params = e.fit()
+ sim_tau = abs(params[0])
+ assert (sim_tau - tau)/tau, 'simulation tau = %g != %g = tau' % (
+ sim_tau, tau)
+ return sim.residual(theory)
+
+
+def test(num_domains=[1,50], unfolding_rate=[1,2], spring_constant=[1,2],
+ velocity=[1,2], threshold=0.3):
+ m = get_manager()()
+ sr = SawsimRunner(manager=m)
+ try:
+ for nd in num_domains:
+ for ur in unfolding_rate:
+ for sc in spring_constant:
+ for v in velocity:
+ residual = constant_rate(
+ sawsim_runner=sr, num_domains=nd,
+ unfolding_rate=ur, spring_constant=sc, velocity=v)
+ assert residual < threshold, residual
+ finally:
+ m.teardown()
+++ /dev/null
-all : const_rate.png
-
-view : all
- qiv const_rate.png
-
-clean :
- rm -f const_rate_*.gp const_rate_*.hist const_rate_*.d const_rate_*.png const_rate_*.compare
-
-const_rate.png : run_test.sh
- ./$<
+++ /dev/null
-* Version 0.8 PASSED
-
-==> const_rate_1_1_1_1.compare <==
-force scale: -1.01187 (expected -1.0)
-initial value: 991.59 (expected 1000.0)
-
-==> const_rate_1_1_1_2.compare <==
-force scale: -0.504875 (expected -0.5)
-initial value: 995.942 (expected 1000.0)
-
-==> const_rate_1_1_2_1.compare <==
-force scale: -1.95167 (expected -2.0)
-initial value: 510.585 (expected 500.0)
-
-==> const_rate_1_2_1_1.compare <==
-force scale: -2.02266 (expected -2.0)
-initial value: 498.409 (expected 500.0)
-
-==> const_rate_50_1_1_1.compare <==
-force scale: -1.01996 (expected -1.0)
-initial value: 49270 (expected 50000.0)
-
-==> const_rate_50_1_1_2.compare <==
-force scale: -0.507475 (expected -0.5)
-initial value: 49484.4 (expected 50000.0)
-
-==> const_rate_50_1_2_1.compare <==
-force scale: -2.02617 (expected -2.0)
-initial value: 24796.2 (expected 25000.0)
-
-==> const_rate_50_2_1_1.compare <==
-force scale: -2.03676 (expected -2.0)
-initial value: 24674.6 (expected 25000.0)
-
-* Version 0.9 PASSED
-
-==> const_rate_1_1_1_1.compare <==
-force scale: -1.0122 (expected -1.0)
-initial value: 996.249 (expected 1000.0)
-
-==> const_rate_1_1_1_2.compare <==
-force scale: -0.498754 (expected -0.5)
-initial value: 1006.67 (expected 1000.0)
-
-==> const_rate_1_1_2_1.compare <==
-force scale: -1.97399 (expected -2.0)
-initial value: 504.787 (expected 500.0)
-
-==> const_rate_1_2_1_1.compare <==
-force scale: -2.01052 (expected -2.0)
-initial value: 497.477 (expected 500.0)
-
-==> const_rate_50_1_1_1.compare <==
-force scale: -1.00178 (expected -1.0)
-initial value: 49942.3 (expected 50000.0)
-
-==> const_rate_50_1_1_2.compare <==
-force scale: -0.498233 (expected -0.5)
-initial value: 50155.7 (expected 50000.0)
-
-==> const_rate_50_1_2_1.compare <==
-force scale: -2.00883 (expected -2.0)
-initial value: 24920.6 (expected 25000.0)
-
-==> const_rate_50_2_1_1.compare <==
-force scale: -1.99658 (expected -2.0)
-initial value: 25031.5 (expected 25000.0)
+++ /dev/null
-#!/bin/bash
-#
-# Use a single hooke domain with a constant unfolding rate, so
-# unfolding force is proportional to time, and connect multiple
-# identical constant unfolding domains. 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)
-#
-# usage: const_rate.sh <num_domains> <velocity> <spring_constant> <unfolding_rate>
-
-if [ $# -ne 4 ]
-then
- echo "usage: $0 <num_domains> <velocity> <spring_constant> <unfolding_rate>"
- exit 1
-fi
-
-NDOM=$1
-V=$2
-K=$3
-RATE=$4
-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_${NDOM}_${V}_${K}_${RATE}.d
-HIST=const_rate_${NDOM}_${V}_${K}_${RATE}.hist
-COMPFILE=const_rate_${NDOM}_${V}_${K}_${RATE}.compare
-GPFILE=const_rate_${NDOM}_${V}_${K}_${RATE}.gp
-PNGFILE=const_rate_${NDOM}_${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*$NDOM*$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
- # Sawsim >= 0.6
- qcmd xtimes $N $SAWSIM -v$V -s cantilever,hooke,$K -N1 \
- -s folded,null -N$NDOM -s unfolded,null \
- -k folded,unfolded,const,$RATE -q folded \
- | grep -v '^#' | cut -f1 >> $DATA
-else
- # Sawsim >= 0.6
- xtimes $N $SAWSIM -v$V -s cantilever,hooke,$K -N1 \
- -s folded,null -N$NDOM -s unfolded,null \
- -k folded,unfolded,const,$RATE -q folded \
- | grep -v '^#' | cut -f1 >> $DATA
-fi
-
-# histogram the data
-if [ $DEBUG -eq 1 ]; then
- echo "cat $DATA | stem_leaf -v-1 -b$df -c | sed 's/://' > $HIST"
-fi
-cat $DATA | stem_leaf -v-1 -b$df -c | 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
+++ /dev/null
-#!/usr/bin/python
-#
-# fit data coming in on stdin to an exponential decay
-# output:
-# <convergence information>
-# initial value:\t<some number>
-# time constant:\t<another number>
-#
-# 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)
+++ /dev/null
-# num domains velocity spring constant unfolding rate
-1 1 1 1
-1 2 1 1
-1 1 2 1
-1 1 1 2
-50 1 1 1
-50 2 1 1
-50 1 2 1
-50 1 1 2
+++ /dev/null
-#!/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
projects / sawsim.git / commitdiff