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 Analytically, 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 :math:`K`, the population follows
  43
  44 .. math::
  45   dp/dt = Kp
  46   p(t) = exp(-tK) = exp(-(f-f_0)K/kv) = p(f) \\;.
  47
  48 Therefore, a histogram of unfolding vs. force :math:`p(f)` normalized
  49 to :math:`p(0)=1` should follow
  50
  51 .. math::
  52   -w N dp/dF = w N pK/kv = w N K/kv exp(-FK/kv)
  53
  54 where :math:`w` is the binwidth and :math:`N` is the number of tested
  55 domains.
  56
  57 >>> test()
  58 """
  59
  60 from numpy import arange, exp, log
  61
  62 from ..histogram import Histogram
  63 from ..fit import HistogramModelFitter
  64 from ..manager import get_manager
  65 from ..sawsim import SawsimRunner
  66 from ..sawsim_histogram import sawsim_histogram
  67
  68
  69 def probability_distribution(x, params):
  70     """Exponential decay decay probability distribution.
  71
  72     .. math:: 1/\\tau \\cdot exp(-x/\\tau)
  73     """
  74     p = params  # convenient alias
  75     p[0] = abs(p[0])  # cannot normalize negative tau.
  76     return (1.0/p[0]) * exp(-x/p[0])
  77
  78
  79 class ExponentialModelFitter (HistogramModelFitter):
  80     """Exponential decay model fitter.
  81     """
  82     def model(self, params):
  83         """An exponential decay model.
  84
  85         Notes
  86         -----
  87         .. math:: y \\propto \\cdot e^{-t/\\tau}
  88         """
  89         self._model_data.counts = (
  90             self.info['binwidth']*self.info['N']*probability_distribution(
  91                 self._model_data.bin_centers, params))
  92         return self._model_data
  93
  94     def guess_initial_params(self):
  95         return [self._data.mean]
  96
  97     def guess_scale(self, params):
  98         return [self._data.mean]
  99
 100     def model_string(self):
 101         return 'p(x) = A exp(-t/tau)'
 102
 103     def param_string(self, params):
 104         pstrings = []
 105         for name,param in zip(['tau'], params):
 106             pstrings.append('%s=%g' % (name, param))
 107         return ', '.join(pstrings)
 108
 109
 110
 111 def constant_rate(sawsim_runner, num_domains=1, unfolding_rate=1,
 112                   spring_constant=1, velocity=1, N=100):
 113     loading_rate = float(spring_constant * velocity)
 114     tau = loading_rate / unfolding_rate
 115     w = 0.1 * tau  # calculate bin width (in force)
 116     A = w*num_domains*N / tau
 117     theory = Histogram()
 118     # A exp(-x/tau) = 0.001
 119     #        -x/tau = log(0.001/A)
 120     #             x = -log(0.001/A)*tau = log(A/0.001)*tau = log(1e3*A)*tau
 121     theory.bin_edges = arange(start=0, stop=log(1e3*A)*tau, step=w)
 122     theory.bin_centers = theory.bin_edges[:-1] + w/2
 123     theory.counts = w*num_domains*N*probability_distribution(
 124         theory.bin_centers, [tau])
 125     theory.analyze()
 126
 127     max_time_step = tau/10.0
 128     max_force_step = loading_rate * max_time_step
 129     param_string = (
 130         '-d %(max_time_step)g -F %(max_force_step)g -v %(velocity)g '
 131         '-s cantilever,hooke,%(spring_constant)g -N1 '
 132         '-s folded,null -N%(num_domains)d -s unfolded,null '
 133         '-k folded,unfolded,const,%(unfolding_rate)g -q folded '
 134         ) % locals()
 135     sim = sawsim_histogram(sawsim_runner, param_string, N=N,
 136                            bin_edges=theory.bin_edges)
 137
 138     e = ExponentialModelFitter(sim)
 139     params = e.fit()
 140     sim_tau = abs(params[0])
 141     assert (sim_tau - tau)/tau, 'simulation tau = %g != %g = tau' % (
 142         sim_tau, tau)
 143     return sim.residual(theory)
 144
 145
 146 def test(num_domains=[1,50], unfolding_rate=[1,2], spring_constant=[1,2],
 147          velocity=[1,2], threshold=0.2):
 148     m = get_manager()()
 149     sr = SawsimRunner(manager=m)
 150     try:
 151         for nd in num_domains:
 152             for ur in unfolding_rate:
 153                 for sc in spring_constant:
 154                     for v in velocity:
 155                         residual = constant_rate(
 156                             sawsim_runner=sr, num_domains=nd,
 157                             unfolding_rate=ur, spring_constant=sc, velocity=v)
 158                         assert residual < threshold, residual
 159     finally:
 160         m.teardown()