Standardize max step calculation and fit-parameter requirements in pysawsim.test.

[sawsim.git] / pysawsim / test / constant_rate.py

diff --git a/pysawsim/test/constant_rate.py b/pysawsim/test/constant_rate.py
old mode 100755 (executable)
new mode 100644 (file)
index 97f755a..c71bdf5
--- a/pysawsim/test/constant_rate.py
+++ b/pysawsim/test/constant_rate.py
@@ -23,7 +23,7 @@ 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
+Analytically, with a spring constant
 
 .. math:: k = df/dx
 
@@ -33,21 +33,17 @@ and a pulling velocity
 
 we have a loading rate
 
-.. math:: df/dt = df/dx dx/dt = kv \;,
+.. math:: df/dt = df/dx dx/dt = kv \\;,
 
 so
 
-.. math:: f = kvt  + f_0 \;.
+.. math:: f = kvt  + f_0 \\;.
 
-With an unfolding rate constant
-
-.. math:: K = 1 \text{frac/s}
-
-The population follows
+With an unfolding rate constant :math:`K`, the population follows
 
 .. math::
-  dp/dt = Kp = 1 \text{s^{-1}}
-  p(t) = exp(-tK) = exp(-(f-f_0)K/kv) = p(f) \;.
+  dp/dt = Kp
+  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
@@ -57,6 +53,8 @@ to :math:`p(0)=1` should follow
 
 where :math:`w` is the binwidth and :math:`N` is the number of tested
 domains.
+
+>>> test()
 """
 
 from numpy import arange, exp, log
@@ -71,7 +69,7 @@ from ..sawsim_histogram import sawsim_histogram
 def probability_distribution(x, params):
     """Exponential decay decay probability distribution.
 
-    .. math:: 1/\\tau \cdot exp(-x/\\tau)
+    .. math:: 1/\\tau \\cdot exp(-x/\\tau)
     """
     p = params  # convenient alias
     p[0] = abs(p[0])  # cannot normalize negative tau.
@@ -86,7 +84,7 @@ class ExponentialModelFitter (HistogramModelFitter):
 
         Notes
         -----
-        .. math:: y \\propto \cdot e^{-t/\tau}
+        .. math:: y \\propto \\cdot e^{-t/\\tau}
         """
         self._model_data.counts = (
             self.info['binwidth']*self.info['N']*probability_distribution(
@@ -114,7 +112,7 @@ def constant_rate(sawsim_runner, num_domains=1, unfolding_rate=1,
                   spring_constant=1, velocity=1, N=100):
     loading_rate = float(spring_constant * velocity)
     tau = loading_rate / unfolding_rate
-    w = 0.1 * tau  # calculate bin width (in force)
+    w = 0.2 * tau  # calculate bin width (in force)
     A = w*num_domains*N / tau
     theory = Histogram()
     # A exp(-x/tau) = 0.001
@@ -126,8 +124,8 @@ def constant_rate(sawsim_runner, num_domains=1, unfolding_rate=1,
         theory.bin_centers, [tau])
     theory.analyze()
 
-    max_time_step = tau/10.0
-    max_force_step = loading_rate * max_time_step
+    max_force_step = w/10.0
+    max_time_step = max_force_step / loading_rate
     param_string = (
         '-d %(max_time_step)g -F %(max_force_step)g -v %(velocity)g '
         '-s cantilever,hooke,%(spring_constant)g -N1 '
@@ -140,8 +138,9 @@ def constant_rate(sawsim_runner, num_domains=1, unfolding_rate=1,
     e = ExponentialModelFitter(sim)
     params = e.fit()
     sim_tau = abs(params[0])
-    assert (sim_tau - tau)/tau, 'simulation tau = %g != %g = tau' % (
-        sim_tau, tau)
+    for s,t,n in [(sim_tau, tau, 'tau')]:
+        assert abs(s - t)/w < 3, (
+            'simulation %s = %g != %g = %s (bin width = %g)' % (n,s,t,n,w))
     return sim.residual(theory)