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