From dc9801a90b0d0f9b383a4018db3d6ce31d65ba7e Mon Sep 17 00:00:00 2001 From: "W. Trevor King" Date: Fri, 5 Nov 2010 13:03:14 -0400 Subject: [PATCH] Move constant unfolding rate test into pysawsim/test/constant_rate.py. --- pysawsim/test/constant_rate.py | 159 +++++++++++++++++++++++++++++ testing/const_rate/Makefile | 10 -- testing/const_rate/NOTES | 67 ------------ testing/const_rate/const_rate.sh | 101 ------------------ testing/const_rate/fit_exponential | 35 ------- testing/const_rate/params | 9 -- testing/const_rate/run_test.sh | 13 --- 7 files changed, 159 insertions(+), 235 deletions(-) create mode 100755 pysawsim/test/constant_rate.py delete mode 100644 testing/const_rate/Makefile delete mode 100644 testing/const_rate/NOTES delete mode 100755 testing/const_rate/const_rate.sh delete mode 100755 testing/const_rate/fit_exponential delete mode 100644 testing/const_rate/params delete mode 100755 testing/const_rate/run_test.sh diff --git a/pysawsim/test/constant_rate.py b/pysawsim/test/constant_rate.py new file mode 100755 index 0000000..bccf3f1 --- /dev/null +++ b/pysawsim/test/constant_rate.py @@ -0,0 +1,159 @@ +# Copyright (C) 2008-2010 W. Trevor King +# +# 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 . +# +# The author may be contacted at 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() diff --git a/testing/const_rate/Makefile b/testing/const_rate/Makefile deleted file mode 100644 index b638aa8..0000000 --- a/testing/const_rate/Makefile +++ /dev/null @@ -1,10 +0,0 @@ -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 - ./$< diff --git a/testing/const_rate/NOTES b/testing/const_rate/NOTES deleted file mode 100644 index 7a066b5..0000000 --- a/testing/const_rate/NOTES +++ /dev/null @@ -1,67 +0,0 @@ -* 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) diff --git a/testing/const_rate/const_rate.sh b/testing/const_rate/const_rate.sh deleted file mode 100755 index baad49c..0000000 --- a/testing/const_rate/const_rate.sh +++ /dev/null @@ -1,101 +0,0 @@ -#!/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 - -if [ $# -ne 4 ] -then - echo "usage: $0 " - 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 diff --git a/testing/const_rate/fit_exponential b/testing/const_rate/fit_exponential deleted file mode 100755 index 7506812..0000000 --- a/testing/const_rate/fit_exponential +++ /dev/null @@ -1,35 +0,0 @@ -#!/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 deleted file mode 100644 index 0a589c3..0000000 --- a/testing/const_rate/params +++ /dev/null @@ -1,9 +0,0 @@ -# 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 diff --git a/testing/const_rate/run_test.sh b/testing/const_rate/run_test.sh deleted file mode 100755 index 2407efc..0000000 --- a/testing/const_rate/run_test.sh +++ /dev/null @@ -1,13 +0,0 @@ -#!/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 -- 2.26.2