pysawsim/test/constant_rate.py

   1 # Copyright (C) 2008-2010  W. Trevor King <wking@drexel.edu>
   2 #
   3 # This program is free software: you can redistribute it and/or modify
   4 # it under the terms of the GNU General Public License as published by
   5 # the Free Software Foundation, either version 3 of the License, or
   6 # (at your option) any later version.
   7 #
   8 # This program is distributed in the hope that it will be useful,
   9 # but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  11 # GNU General Public License for more details.
  12 #
  13 # You should have received a copy of the GNU General Public License
  14 # along with this program.  If not, see <http://www.gnu.org/licenses/>.
  15 #
  16 # The author may be contacted at <wking@drexel.edu> on the Internet, or
  17 # write to Trevor King, Drexel University, Physics Dept., 3141 Chestnut St.,
  18 # Philadelphia PA 19104, USA.
  19
  20 """Test a Hookian chain with constant unfolding rate.
  21
  22 With the constant velocity experiment and a Hookian domain, the
  23 unfolding force is proportional to time, so we expect exponential
  24 population decay, even with multiple unfolding domains.
  25
  26 For example, with a spring constant
  27
  28 .. math:: k = df/dx
  29
  30 and a pulling velocity
  31
  32 .. math:: v = dx/dt
  33
  34 we have a loading rate
  35
  36 .. math:: df/dt = df/dx dx/dt = kv \;,
  37
  38 so
  39
  40 .. math:: f = kvt  + f_0 \;.
  41
  42 With an unfolding rate constant
  43
  44 .. math:: K = 1 \text{frac/s}
  45
  46 The population follows
  47
  48 .. math::
  49   dp/dt = Kp = 1 \text{s^{-1}}
  50   p(t) = exp(-tK) = exp(-(f-f_0)K/kv) = p(f) \;.
  51
  52 Therefore, a histogram of unfolding vs. force :math:`p(f)` normalized
  53 to :math:`p(0)=1` should follow
  54
  55 .. math::
  56   -w N dp/dF = w N pK/kv = w N K/kv exp(-FK/kv)
  57
  58 where :math:`w` is the binwidth and :math:`N` is the number of tested
  59 domains.
  60
  61 >>> test()
  62 """
  63
  64 from numpy import arange, exp, log
  65
  66 from ..histogram import Histogram
  67 from ..fit import HistogramModelFitter
  68 from ..manager import get_manager
  69 from ..sawsim import SawsimRunner
  70 from ..sawsim_histogram import sawsim_histogram
  71
  72
  73 def probability_distribution(x, params):
  74     """Exponential decay decay probability distribution.
  75
  76     .. math:: 1/\\tau \cdot exp(-x/\\tau)
  77     """
  78     p = params  # convenient alias
  79     p[0] = abs(p[0])  # cannot normalize negative tau.
  80     return (1.0/p[0]) * exp(-x/p[0])
  81
  82
  83 class ExponentialModelFitter (HistogramModelFitter):
  84     """Exponential decay model fitter.
  85     """
  86     def model(self, params):
  87         """An exponential decay model.
  88
  89         Notes
  90         -----
  91         .. math:: y \\propto \cdot e^{-t/\tau}
  92         """
  93         self._model_data.counts = (
  94             self.info['binwidth']*self.info['N']*probability_distribution(
  95                 self._model_data.bin_centers, params))
  96         return self._model_data
  97
  98     def guess_initial_params(self):
  99         return [self._data.mean]
 100
 101     def guess_scale(self, params):
 102         return [self._data.mean]
 103
 104     def model_string(self):
 105         return 'p(x) = A exp(-t/tau)'
 106
 107     def param_string(self, params):
 108         pstrings = []
 109         for name,param in zip(['tau'], params):
 110             pstrings.append('%s=%g' % (name, param))
 111         return ', '.join(pstrings)
 112
 113
 114
 115 def constant_rate(sawsim_runner, num_domains=1, unfolding_rate=1,
 116                   spring_constant=1, velocity=1, N=100):
 117     loading_rate = float(spring_constant * velocity)
 118     tau = loading_rate / unfolding_rate
 119     w = 0.1 * tau  # calculate bin width (in force)
 120     A = w*num_domains*N / tau
 121     theory = Histogram()
 122     # A exp(-x/tau) = 0.001
 123     #        -x/tau = log(0.001/A)
 124     #             x = -log(0.001/A)*tau = log(A/0.001)*tau = log(1e3*A)*tau
 125     theory.bin_edges = arange(start=0, stop=log(1e3*A)*tau, step=w)
 126     theory.bin_centers = theory.bin_edges[:-1] + w/2
 127     theory.counts = w*num_domains*N*probability_distribution(
 128         theory.bin_centers, [tau])
 129     theory.analyze()
 130
 131     max_time_step = tau/10.0
 132     max_force_step = loading_rate * max_time_step
 133     param_string = (
 134         '-d %(max_time_step)g -F %(max_force_step)g -v %(velocity)g '
 135         '-s cantilever,hooke,%(spring_constant)g -N1 '
 136         '-s folded,null -N%(num_domains)d -s unfolded,null '
 137         '-k folded,unfolded,const,%(unfolding_rate)g -q folded '
 138         ) % locals()
 139     sim = sawsim_histogram(sawsim_runner, param_string, N=N,
 140                            bin_edges=theory.bin_edges)
 141
 142     e = ExponentialModelFitter(sim)
 143     params = e.fit()
 144     sim_tau = abs(params[0])
 145     assert (sim_tau - tau)/tau, 'simulation tau = %g != %g = tau' % (
 146         sim_tau, tau)
 147     return sim.residual(theory)
 148
 149
 150 def test(num_domains=[1,50], unfolding_rate=[1,2], spring_constant=[1,2],
 151          velocity=[1,2], threshold=0.2):
 152     m = get_manager()()
 153     sr = SawsimRunner(manager=m)
 154     try:
 155         for nd in num_domains:
 156             for ur in unfolding_rate:
 157                 for sc in spring_constant:
 158                     for v in velocity:
 159                         residual = constant_rate(
 160                             sawsim_runner=sr, num_domains=nd,
 161                             unfolding_rate=ur, spring_constant=sc, velocity=v)
 162                         assert residual < threshold, residual
 163     finally:
 164         m.teardown()