hooke/plugin/polymer_fit.py

   1 # -*- coding: utf-8 -*-
   2 #
   3 # Copyright (C) 2008-2010 Alberto Gomez-Casado
   4 #                         Fabrizio Benedetti
   5 #                         Massimo Sandal <devicerandom@gmail.com>
   6 #                         W. Trevor King <wking@drexel.edu>
   7 #
   8 # This file is part of Hooke.
   9 #
  10 # Hooke is free software: you can redistribute it and/or modify it
  11 # under the terms of the GNU Lesser General Public License as
  12 # published by the Free Software Foundation, either version 3 of the
  13 # License, or (at your option) any later version.
  14 #
  15 # Hooke is distributed in the hope that it will be useful, but WITHOUT
  16 # ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  17 # or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General
  18 # Public License for more details.
  19 #
  20 # You should have received a copy of the GNU Lesser General Public
  21 # License along with Hooke.  If not, see
  22 # <http://www.gnu.org/licenses/>.
  23
  24 """The ``polymer_fit`` module proviews :class:`PolymerFitPlugin` and
  25 serveral associated :class:`hooke.command.Command`\s for Fitting
  26 velocity-clamp data to various polymer models (WLC, FJC, etc.).
  27 """
  28
  29 import copy
  30 import Queue
  31 import StringIO
  32 import sys
  33
  34 import numpy
  35 from scipy.optimize import newton, bisect
  36
  37 from ..command import Command, Argument, Failure
  38 from ..config import Setting
  39 from ..curve import Data
  40 from ..plugin import Plugin
  41 from ..util.fit import PoorFit, ModelFitter
  42 from ..util.callback import is_iterable
  43 from .curve import CurveArgument
  44 from .vclamp import scale
  45
  46
  47 kB = 1.3806503e-23  # Joules/Kelvin
  48
  49
  50 def coth(z):
  51     """Hyperbolic cotangent.
  52
  53     Examples
  54     --------
  55     >>> coth(1.19967874)  # doctest: +ELLIPSIS
  56     1.199678...
  57
  58     Notes
  59     -----
  60     From `MathWorld`_
  61
  62     .. math::
  63       \coth(z) \equiv \frac{\exp(z) + \exp(-z)}{\exp(z) - \exp(-z)}
  64          = \frac{1}{\tanh(z)}
  65
  66     .. _MathWorld:
  67       http://mathworld.wolfram.com/HyperbolicCotangent.html
  68     """
  69     return 1.0/numpy.tanh(z)
  70
  71 def arccoth(z):
  72     """Inverse hyperbolic cotangent.
  73
  74     Examples
  75     --------
  76     >>> arccoth(1.19967874)  # doctest: +ELLIPSIS
  77     1.199678...
  78     >>> arccoth(coth(numpy.pi))  # doctest: +ELLIPSIS
  79     3.1415...
  80
  81     Notes
  82     -----
  83     Inverting from the :func:`definition <coth>`.
  84
  85     .. math::
  86       z \equiv \atanh\left( \frac{1}{\coth(z)} \right)
  87     """
  88     return numpy.arctanh(1.0/z)
  89
  90 def Langevin(z):
  91     """Langevin function.
  92
  93     Examples
  94     --------
  95     >>> Langevin(numpy.pi)  # doctest: +ELLIPSIS
  96     0.685...
  97
  98     Notes
  99     -----
 100     From `Wikipedia`_ or Hatfield [#]_
 101
 102     .. math::
 103       L(z) \equiv \coth(z) - \frac{1}{z}
 104
 105     .. _Wikipedia:
 106       http://en.wikipedia.org/wiki/Brillouin_and_Langevin_functions#Langevin_Function
 107
 108     .. [#]: J.W. Hatfield and S.R. Quake
 109       "Dynamic Properties of an Extended Polymer in Solution."
 110       Phys. Rev. Lett. 1999.
 111       doi: `10.1103/PhysRevLett.82.3548 <http://dx.doi.org/10.1103/PhysRevLett.82.3548>`
 112     """
 113     return coth(z) - 1.0/z
 114
 115 def inverse_Langevin(z, extreme=1.0 - 1e-8):
 116     """Inverse Langevin function.
 117
 118     Parameters
 119     ----------
 120     z : float or array_like
 121         object whose inverse Langevin will be returned
 122     extreme : float
 123         value (close to one) for which we assume the inverse is +/-
 124         infinity.  This avoids problems with Newton-Raphson
 125         convergence in regions with low slope.
 126
 127     Examples
 128     --------
 129     >>> inverse_Langevin(Langevin(numpy.pi))  # doctest: +ELLIPSIS
 130     3.14159...
 131     >>> inverse_Langevin(Langevin(numpy.array([1,2,3])))  # doctest: +ELLIPSIS
 132     array([ 1.,  2.,  3.])
 133
 134     Notes
 135     -----
 136     We approximate the inverse Langevin via Newton's method
 137     (:func:`scipy.optimize.newton`).  For the starting point, we take
 138     the first three terms from the Taylor series (from `Wikipedia`_).
 139
 140     .. math::
 141       L^{-1}(z) = 3z + \frac{9}{5}z^3 + \frac{297}{175}z^5 + \dots
 142
 143     .. _Wikipedia:
 144       http://en.wikipedia.org/wiki/Brillouin_and_Langevin_functions#Langevin_Function
 145     """
 146     if is_iterable(z):
 147         ret = numpy.ndarray(shape=z.shape, dtype=z.dtype)
 148         for i,z_i in enumerate(z):
 149             ret[i] = inverse_Langevin(z_i)
 150         return ret
 151     if z > extreme:
 152         return numpy.inf
 153     elif z < -extreme:
 154         return -numpy.inf
 155     # Catch stdout since sometimes the newton function prints exit
 156     # messages to stdout, e.g. "Tolerance of %s reached".
 157     stdout = sys.stdout
 158     tmp_stdout = StringIO.StringIO()
 159     sys.stdout = tmp_stdout
 160     ret = newton(func=lambda x: Langevin(x)-z,
 161                   x0=3*z + 9./5.*z**3 + 297./175.*z**5,)
 162     sys.stdout = stdout
 163     output = tmp_stdout.getvalue()
 164     return ret
 165
 166 def FJC_fn(x_data, T, L, a):
 167     """The freely jointed chain model.
 168
 169     Parameters
 170     ----------
 171     x_data : array_like
 172         x values for which the WLC tension should be calculated.
 173     T : float
 174         temperature in Kelvin.
 175     L : float
 176         contour length in meters.
 177     a : float
 178         Kuhn length in meters.
 179
 180     Returns
 181     -------
 182     f_data : array
 183         `F(x_data)`.
 184
 185     Examples
 186     --------
 187     >>> FJC_fn(numpy.array([1e-9, 5e-9, 10e-9], dtype=numpy.float),
 188     ...        T=300, L=15e-9, a=2.5e-10) # doctest: +ELLIPSIS
 189     array([  3.322...-12,   1.78...e-11,   4.889...e-11])
 190
 191     Notes
 192     -----
 193     The freely jointed chain model is
 194
 195     .. math::
 196       F(x) = \frac{k_B T}{a} L^{-1}\left( \frac{x}{L} \right)
 197
 198     where :math:`L^{-1}` is the :func:`inverse_Langevin`, :math:`a` is
 199     the Kuhn length, and :math:`L` is the contour length [#]_.  For
 200     the inverse formulation (:math:`x(F)`), see Ray and Akhremitchev [#]_.
 201
 202
 203     .. [#]: J.W. Hatfield and S.R. Quake
 204       "Dynamic Properties of an Extended Polymer in Solution."
 205       Phys. Rev. Lett. 1999.
 206       doi: `10.1103/PhysRevLett.82.3548 <http://dx.doi.org/10.1103/PhysRevLett.82.3548>`
 207
 208     .. [#] C. Ray and B.B. Akhremitchev.
 209       `"Fitting Force Curves by the Extended Freely Jointed Chain Model" <http://www.chem.duke.edu/~boris/research/force_spectroscopy/fit_efjc.pdf>`
 210     """
 211     return kB*T / a * inverse_Langevin(x_data/L)
 212
 213 class FJC (ModelFitter):
 214     """Fit the data to a freely jointed chain.
 215
 216     Examples
 217     --------
 218     Generate some example data
 219
 220     >>> T = 300  # Kelvin
 221     >>> L = 35e-9  # meters
 222     >>> a = 2.5e-10  # meters
 223     >>> x_data = numpy.linspace(10e-9, 30e-9, num=20)
 224     >>> d_data = FJC_fn(x_data, T=T, L=L, a=a)
 225
 226     Fit the example data with a two-parameter fit (`L` and `a`).
 227
 228     >>> info = {
 229     ...     'temperature (K)': T,
 230     ...     'x data (m)': x_data,
 231     ...     }
 232     >>> model = FJC(d_data, info=info, rescale=True)
 233     >>> outqueue = Queue.Queue()
 234     >>> Lp,a = model.fit(outqueue=outqueue)
 235     >>> fit_info = outqueue.get(block=False)
 236     >>> model.L(Lp)  # doctest: +ELLIPSIS
 237     3.500...e-08
 238     >>> a  # doctest: +ELLIPSIS
 239     2.499...e-10
 240
 241     Fit the example data with a one-parameter fit (`L`).  We introduce
 242     some error in our fixed Kuhn length for fun.
 243
 244     >>> info['Kuhn length (m)'] = 2*a
 245     >>> model = FJC(d_data, info=info, rescale=True)
 246     >>> Lp = model.fit(outqueue=outqueue)
 247     >>> fit_info = outqueue.get(block=False)
 248     >>> model.L(Lp)  # doctest: +ELLIPSIS
 249     3.199...e-08
 250     """
 251     def Lp(self, L):
 252         """To avoid invalid values of `L`, we reparametrize with `Lp`.
 253
 254         Parameters
 255         ----------
 256         L : float
 257             contour length in meters
 258
 259         Returns
 260         -------
 261         Lp : float
 262             rescaled version of `L`.
 263
 264         Examples
 265         --------
 266         >>> x_data = numpy.linspace(1e-9, 20e-9, num=10)
 267         >>> model = FJC(data=x_data, info={'x data (m)':x_data})
 268         >>> model.Lp(20e-9)  # doctest: +ELLIPSIS
 269         -inf
 270         >>> model.Lp(19e-9)  # doctest: +ELLIPSIS
 271         nan
 272         >>> model.Lp(21e-9)  # doctest: +ELLIPSIS
 273         -2.99...
 274         >>> model.Lp(100e-9)  # doctest: +ELLIPSIS
 275         1.386...
 276
 277         Notes
 278         -----
 279         The rescaling is designed to limit `L` to values strictly
 280         greater than the maximum `x` data value, so we use
 281
 282         .. math::
 283             L = [\exp(L_p))+1] * x_\text{max}
 284
 285         which has the inverse
 286
 287         .. math::
 288             L_p = \ln(L/x_\text{max}-1)
 289
 290         This will obviously effect the interpretation of the fit's covariance matrix.
 291         """
 292         x_max = self.info['x data (m)'].max()
 293         return numpy.log(L/x_max - 1)
 294
 295     def L(self, Lp):
 296         """Inverse of :meth:`Lp`.
 297
 298         Parameters
 299         ----------
 300         Lp : float
 301             rescaled version of `L`
 302
 303         Returns
 304         -------
 305         L : float
 306             contour length in meters
 307
 308         Examples
 309         --------
 310         >>> x_data = numpy.linspace(1e-9, 20e-9, num=10)
 311         >>> model = FJC(data=x_data, info={'x data (m)':x_data})
 312         >>> model.L(model.Lp(21e-9))  # doctest: +ELLIPSIS
 313         2.100...e-08
 314         >>> model.L(model.Lp(100e-9))  # doctest: +ELLIPSIS
 315         9.999...e-08
 316         """
 317         x_max = self.info['x data (m)'].max()
 318         return (numpy.exp(Lp)+1)*x_max
 319
 320     def model(self, params):
 321         """Extract the relavant arguments and use :func:`FJC_fn`.
 322
 323         Parameters
 324         ----------
 325         params : list of floats
 326             `[Lp, a]`, where the presence of `a` is optional.
 327         """
 328         # setup convenient aliases
 329         T = self.info['temperature (K)']
 330         x_data = self.info['x data (m)']
 331         L = self.L(params[0])
 332         if len(params) > 1:
 333             a = abs(params[1])
 334         else:
 335             a = self.info['Kuhn length (m)']
 336         # compute model data
 337         self._model_data[:] = FJC_fn(x_data, T, L, a)
 338         return self._model_data
 339
 340     def guess_initial_params(self, outqueue=None):
 341         """Guess initial fitting parameters.
 342
 343         Returns
 344         -------
 345         Lp : float
 346             A guess at the reparameterized contour length in meters.
 347         a : float (optional)
 348             A guess at the Kuhn length in meters.  If `.info` has a
 349             setting for `'Kuhn length (m)'`, `a` is not returned.
 350         """
 351         x_max = self.info['x data (m)'].max()
 352         a_given = 'Kuhn length (m)' in self.info
 353         if a_given:
 354             a = self.info['Kuhn length (m)']
 355         else:
 356             a = x_max / 10.0
 357         L = 1.5 * x_max
 358         Lp = self.Lp(L)
 359         if a_given:
 360             return [Lp]
 361         return [Lp, a]
 362
 363     def guess_scale(self, params, outqueue=None):
 364         """Guess parameter scales.
 365
 366         Returns
 367         -------
 368         Lp_scale : float
 369             A guess at the reparameterized contour length scale in meters.
 370         a_scale : float (optional)
 371             A guess at the Kuhn length in meters.  If the length of
 372             `params` is less than 2, `a_scale` is not returned.
 373         """
 374         Lp_scale = 1.0
 375         if len(params) == 1:
 376             return [Lp_scale]
 377         return [Lp_scale] + [x/10.0 for x in params[1:]]
 378
 379
 380 def inverse_FJC_PEG_fn(F_data, T=300, N=1, k=150., Lp=3.58e-10, Lh=2.8e-10, dG=3., a=7e-10):
 381     """Inverse poly(ethylene-glycol) adjusted extended FJC model.
 382
 383     Examples
 384     --------
 385     >>> kwargs = {'T':300.0, 'N':1, 'k':150.0, 'Lp':3.58e-10, 'Lh':2.8e-10,
 386     ...           'dG':3.0, 'a':7e-10}
 387     >>> inverse_FJC_PEG_fn(F_data=200e-12, **kwargs)  # doctest: +ELLIPSIS
 388     3.487...e-10
 389
 390     Notes
 391     -----
 392     The function is that proposed by F. Oesterhelt, et al. [#]_.
 393
 394     .. math::
 395       x(F) = N_\text{S} \cdot \left[
 396         \left(
 397             \frac{L_\text{planar}}{
 398                   \exp\left(-\Delta G/k_B T\right) + 1}
 399             + \frac{L_\text{helical}}{
 400                   \exp\left(+\Delta G/k_B T\right) + 1}
 401           \right) \cdot \left(
 402             \coth\left(\frac{Fa}{k_B T}\right)
 403             - \frac{k_B T}{Fa}
 404           \right)
 405         + \frac{F}{K_\text{S}}
 406
 407     where :math:`N_\text{S}` is the number of segments,
 408     :math:`K_\text{S}` is the segment elasticicty,
 409     :math:`L_\text{planar}` is the ttt contour length,
 410     :math:`L_\text{helical}` is the ttg contour length,
 411     :math:`\Delta G` is the Gibbs free energy difference
 412      :math:`G_\text{planar}-G_\text{helical}`,
 413     :math:`a` is the Kuhn length,
 414     and :math:`F` is the chain tension.
 415
 416     .. [#]: F. Oesterhelt, M. Rief, and H.E. Gaub.
 417       "Single molecule force spectroscopy by AFM indicates helical
 418       structure of poly(ethylene-glycol) in water."
 419       New Journal of Physics. 1999.
 420       doi: `10.1088/1367-2630/1/1/006 <http://dx.doi.org/10.1088/1367-2630/1/1/006>`
 421       (section 4.2)
 422     """
 423     kBT = kB * T
 424     g = (dG - F_data*(Lp-Lh)) / kBT
 425     z = F_data*a/kBT
 426     return N * ((Lp/(numpy.exp(-g)+1) + Lh/(numpy.exp(+g)+1)) * (coth(z)-1./z)
 427                 + F_data/k)
 428
 429 def FJC_PEG_fn(x_data, T, N, k, Lp, Lh, dG, a):
 430     """Poly(ethylene-glycol) adjusted extended FJC model.
 431
 432     Parameters
 433     ----------
 434     z : float or array_like
 435         object whose inverse Langevin will be returned
 436     extreme : float
 437         value (close to one) for which we assume the inverse is +/-
 438         infinity.  This avoids problems with Newton-Raphson
 439         convergence in regions with low slope.
 440
 441     Examples
 442     --------
 443     >>> kwargs = {'T':300.0, 'N':1, 'k':150.0, 'Lp':3.58e-10, 'Lh':2.8e-10,
 444     ...           'dG':3.0, 'a':7e-10}
 445     >>> FJC_PEG_fn(inverse_FJC_PEG_fn(1e-11, **kwargs), **kwargs)  # doctest: +ELLIPSIS
 446     9.999...e-12
 447     >>> FJC_PEG_fn(inverse_FJC_PEG_fn(3.4e-10, **kwargs), **kwargs)  # doctest: +ELLIPSIS
 448     3.400...e-10
 449     >>> FJC_PEG_fn(numpy.array([1e-10,2e-10,3e-10]), **kwargs)  # doctest: +ELLIPSIS
 450     array([  5.207...e-12,   1.257...e-11,   3.636...e-11])
 451     >>> kwargs['N'] = 123
 452     >>> FJC_PEG_fn(numpy.array([1e-8,2e-8,3e-8]), **kwargs)  # doctest: +ELLIPSIS
 453     array([  4.160...e-12,   9.302...e-12,   1.830...e-11])
 454
 455     Notes
 456     -----
 457     We approximate the PEG adjusted eJPG via bisection
 458     (:func:`scipy.optimize.bisect`) Because it is more stable than
 459     Newton's algorithm.  For the starting point, we use the standard
 460     JPG function with an averaged contour length.
 461     """
 462     kwargs = {'T':T, 'N':N, 'k':k, 'Lp':Lp, 'Lh':Lh, 'dG':dG, 'a':a}
 463     if is_iterable(x_data):
 464         ret = numpy.ndarray(shape=x_data.shape, dtype=x_data.dtype)
 465         for i,x in enumerate(x_data):
 466             ret[i] = FJC_PEG_fn(x, **kwargs)
 467         return ret
 468     if x_data == 0:
 469         return 0
 470     # Approximate with the simple FJC_fn()
 471     guess = numpy.inf
 472     L = N*max(Lp, Lh)
 473     while guess == numpy.inf:
 474         guess = FJC_fn(x_data, T=T, L=L, a=a)
 475         L *= 2.0
 476
 477     if False: # Bisection
 478         while inverse_FJC_PEG_fn(guess, **kwargs) < x_data:
 479             guess *= 2.0
 480         lower = guess
 481         while inverse_FJC_PEG_fn(lower, **kwargs) > x_data:
 482             lower /= 2.0
 483
 484         fn = lambda f: inverse_FJC_PEG_fn(guess*abs(f), **kwargs) - x_data
 485         return guess*abs(bisect(f=fn, a=lower/guess, b=1.0))
 486     else:  # Newton
 487         fn = lambda f: inverse_FJC_PEG_fn(guess*abs(f), **kwargs) - x_data
 488         ret = guess*abs(newton(func=fn, x0=1.0))
 489         return ret
 490
 491 class FJC_PEG (ModelFitter):
 492     """Fit the data to an poly(ethylene-glycol) adjusted extended FJC.
 493
 494     Examples
 495     --------
 496     Generate some example data
 497
 498     >>> kwargs = {'T':300.0, 'N':123, 'k':150.0, 'Lp':3.58e-10, 'Lh':2.8e-10,
 499     ...           'dG':3.0, 'a':7e-10}
 500     >>> x_data = numpy.linspace(10e-9, 30e-9, num=20)
 501     >>> d_data = FJC_PEG_fn(x_data, **kwargs)
 502
 503     Fit the example data with a two-parameter fit (`N` and `a`).
 504
 505     >>> info = {
 506     ...     'temperature (K)': kwargs['T'],
 507     ...     'x data (m)': x_data,
 508     ...     'section elasticity (N/m)': kwargs['k'],
 509     ...     'planar section length (m)': kwargs['Lp'],
 510     ...     'helical section length (m)': kwargs['Lh'],
 511     ...     'Gibbs free energy difference (Gp - Gh) (kBT)': kwargs['dG'],
 512     ...     }
 513     >>> model = FJC_PEG(d_data, info=info, rescale=True)
 514     >>> outqueue = Queue.Queue()
 515     >>> Nr,a = model.fit(outqueue=outqueue)
 516     >>> fit_info = outqueue.get(block=False)
 517     >>> model.L(Nr)  # doctest: +ELLIPSIS
 518     122.999...
 519     >>> a  # doctest: +ELLIPSIS
 520     6.999...e-10
 521
 522     Fit the example data with a one-parameter fit (`N`).  We introduce
 523     some error in our fixed Kuhn length for fun.
 524
 525     >>> info['Kuhn length (m)'] = 2*kwargs['a']
 526     >>> model = FJC_PEG(d_data, info=info, rescale=True)
 527     >>> Nr = model.fit(outqueue=outqueue)
 528     >>> fit_info = outqueue.get(block=False)
 529     >>> model.L(Nr)  # doctest: +ELLIPSIS
 530     96.931...
 531     """
 532     def Lr(self, L):
 533         """To avoid invalid values of `L`, we reparametrize with `Lr`.
 534
 535         Parameters
 536         ----------
 537         L : float
 538             contour length in meters
 539
 540         Returns
 541         -------
 542         Lr : float
 543             rescaled version of `L`.
 544
 545         Examples
 546         --------
 547         >>> x_data = numpy.linspace(1e-9, 20e-9, num=10)
 548         >>> model = FJC_PEG(data=x_data, info={'x data (m)':x_data})
 549         >>> model.Lr(20e-9)  # doctest: +ELLIPSIS
 550         0.0
 551         >>> model.Lr(19e-9)  # doctest: +ELLIPSIS
 552         -0.0512...
 553         >>> model.Lr(21e-9)  # doctest: +ELLIPSIS
 554         0.0487...
 555         >>> model.Lr(100e-9)  # doctest: +ELLIPSIS
 556         1.609...
 557
 558         Notes
 559         -----
 560         The rescaling is designed to limit `L` to values strictly
 561         greater than zero, so we use
 562
 563         .. math::
 564             L = \exp(L_p) * x_\text{max}
 565
 566         which has the inverse
 567
 568         .. math::
 569             L_p = \ln(L/x_\text{max})
 570
 571         This will obviously effect the interpretation of the fit's covariance matrix.
 572         """
 573         x_max = self.info['x data (m)'].max()
 574         return numpy.log(L/x_max)
 575
 576     def L(self, Lr):
 577         """Inverse of :meth:`Lr`.
 578
 579         Parameters
 580         ----------
 581         Lr : float
 582             rescaled version of `L`
 583
 584         Returns
 585         -------
 586         L : float
 587             contour length in meters
 588
 589         Examples
 590         --------
 591         >>> x_data = numpy.linspace(1e-9, 20e-9, num=10)
 592         >>> model = FJC_PEG(data=x_data, info={'x data (m)':x_data})
 593         >>> model.L(model.Lr(21e-9))  # doctest: +ELLIPSIS
 594         2.100...e-08
 595         >>> model.L(model.Lr(100e-9))  # doctest: +ELLIPSIS
 596         9.999...e-08
 597         """
 598         x_max = self.info['x data (m)'].max()
 599         return numpy.exp(Lr)*x_max
 600
 601     def _kwargs(self):
 602         return {
 603             'T': self.info['temperature (K)'],
 604             'k': self.info['section elasticity (N/m)'],
 605             'Lp': self.info['planar section length (m)'],
 606             'Lh': self.info['helical section length (m)'],
 607             'dG': self.info['Gibbs free energy difference (Gp - Gh) (kBT)'],
 608             }
 609
 610     def model(self, params):
 611         """Extract the relavant arguments and use :func:`FJC_PEG_fn`.
 612
 613         Parameters
 614         ----------
 615         params : list of floats
 616             `[N, a]`, where the presence of `a` is optional.
 617         """
 618         # setup convenient aliases
 619         T = self.info['temperature (K)']
 620         x_data = self.info['x data (m)']
 621         N = self.L(params[0])
 622         if len(params) > 1:
 623             a = abs(params[1])
 624         else:
 625             a = self.info['Kuhn length (m)']
 626         kwargs = self._kwargs()
 627         # compute model data
 628         self._model_data[:] = FJC_PEG_fn(x_data, N=N, a=a, **kwargs)
 629         return self._model_data
 630
 631     def guess_initial_params(self, outqueue=None):
 632         """Guess initial fitting parameters.
 633
 634         Returns
 635         -------
 636         N : float
 637             A guess at the section count.
 638         a : float (optional)
 639             A guess at the Kuhn length in meters.  If `.info` has a
 640             setting for `'Kuhn length (m)'`, `a` is not returned.
 641         """
 642         x_max = self.info['x data (m)'].max()
 643         a_given = 'Kuhn length (m)' in self.info
 644         if a_given:
 645             a = self.info['Kuhn length (m)']
 646         else:
 647             a = x_max / 10.0
 648         f_max = self._data.max()
 649         kwargs = self._kwargs()
 650         kwargs['a'] = a
 651         x_section = inverse_FJC_PEG_fn(F_data=f_max, **kwargs)
 652         N = x_max / x_section;
 653         Nr = self.Lr(N)
 654         if a_given:
 655             return [Nr]
 656         return [Nr, a]
 657
 658     def guess_scale(self, params, outqueue=None):
 659         """Guess parameter scales.
 660
 661         Returns
 662         -------
 663         N_scale : float
 664             A guess at the section count scale in meters.
 665         a_scale : float (optional)
 666             A guess at the Kuhn length in meters.  If the length of
 667             `params` is less than 2, `a_scale` is not returned.
 668         """
 669         return [x/10.0 for x in params]
 670
 671
 672 def WLC_fn(x_data, T, L, p):
 673     """The worm like chain interpolation model.
 674
 675     Parameters
 676     ----------
 677     x_data : array_like
 678         x values for which the WLC tension should be calculated.
 679     T : float
 680         temperature in Kelvin.
 681     L : float
 682         contour length in meters.
 683     p : float
 684         persistence length in meters.
 685
 686     Returns
 687     -------
 688     f_data : array
 689         `F(x_data)`.
 690
 691     Examples
 692     --------
 693     >>> WLC_fn(numpy.array([1e-9, 5e-9, 10e-9], dtype=numpy.float),
 694     ...        T=300, L=15e-9, p=2.5e-10) # doctest: +ELLIPSIS
 695     array([  1.717...e-12,   1.070...e-11,   4.418...e-11])
 696
 697     Notes
 698     -----
 699     The function is the simple polynomial worm like chain
 700     interpolation forumula proposed by C. Bustamante, et al. [#]_.
 701
 702     .. math::
 703       F(x) = \frac{k_B T}{p} \left[
 704         \frac{1}{4}\left( \frac{1}{(1-x/L)^2} - 1 \right)
 705         + \frac{x}{L}
 706                              \right]
 707
 708     .. [#] C. Bustamante, J.F. Marko, E.D. Siggia, and S.Smith.
 709     "Entropic elasticity of lambda-phage DNA."
 710     Science. 1994.
 711     doi: `10.1126/science.8079175 <http://dx.doi.org/10.1126/science.8079175>`
 712     """
 713     a = kB * T / p
 714     scaled_data = x_data / L
 715     return a * (0.25*((1-scaled_data)**-2 - 1) + scaled_data)
 716
 717 class WLC (ModelFitter):
 718     """Fit the data to a worm like chain.
 719
 720     Examples
 721     --------
 722     Generate some example data
 723
 724     >>> T = 300  # Kelvin
 725     >>> L = 35e-9  # meters
 726     >>> p = 2.5e-10  # meters
 727     >>> x_data = numpy.linspace(10e-9, 30e-9, num=20)
 728     >>> d_data = WLC_fn(x_data, T=T, L=L, p=p)
 729
 730     Fit the example data with a two-parameter fit (`L` and `p`).
 731
 732     >>> info = {
 733     ...     'temperature (K)': T,
 734     ...     'x data (m)': x_data,
 735     ...     }
 736     >>> model = WLC(d_data, info=info, rescale=True)
 737     >>> outqueue = Queue.Queue()
 738     >>> Lp,p = model.fit(outqueue=outqueue)
 739     >>> fit_info = outqueue.get(block=False)
 740     >>> model.L(Lp)  # doctest: +ELLIPSIS
 741     3.500...e-08
 742     >>> p  # doctest: +ELLIPSIS
 743     2.500...e-10
 744
 745     Fit the example data with a one-parameter fit (`L`).  We introduce
 746     some error in our fixed persistence length for fun.
 747
 748     >>> info['persistence length (m)'] = 2*p
 749     >>> model = WLC(d_data, info=info, rescale=True)
 750     >>> Lp = model.fit(outqueue=outqueue)
 751     >>> fit_info = outqueue.get(block=False)
 752     >>> model.L(Lp)  # doctest: +ELLIPSIS
 753     3.318...e-08
 754     """
 755     def Lp(self, L):
 756         """To avoid invalid values of `L`, we reparametrize with `Lp`.
 757
 758         Parameters
 759         ----------
 760         L : float
 761             contour length in meters
 762
 763         Returns
 764         -------
 765         Lp : float
 766             rescaled version of `L`.
 767
 768         Examples
 769         --------
 770         >>> x_data = numpy.linspace(1e-9, 20e-9, num=10)
 771         >>> model = WLC(data=x_data, info={'x data (m)':x_data})
 772         >>> model.Lp(20e-9)  # doctest: +ELLIPSIS
 773         -inf
 774         >>> model.Lp(19e-9)  # doctest: +ELLIPSIS
 775         nan
 776         >>> model.Lp(21e-9)  # doctest: +ELLIPSIS
 777         -2.99...
 778         >>> model.Lp(100e-9)  # doctest: +ELLIPSIS
 779         1.386...
 780
 781         Notes
 782         -----
 783         The rescaling is designed to limit `L` to values strictly
 784         greater than the maximum `x` data value, so we use
 785
 786         .. math::
 787             L = [\exp(L_p))+1] * x_\text{max}
 788
 789         which has the inverse
 790
 791         .. math::
 792             L_p = \ln(L/x_\text{max}-1)
 793
 794         This will obviously effect the interpretation of the fit's covariance matrix.
 795         """
 796         x_max = self.info['x data (m)'].max()
 797         return numpy.log(L/x_max - 1)
 798
 799     def L(self, Lp):
 800         """Inverse of :meth:`Lp`.
 801
 802         Parameters
 803         ----------
 804         Lp : float
 805             rescaled version of `L`
 806
 807         Returns
 808         -------
 809         L : float
 810             contour length in meters
 811
 812         Examples
 813         --------
 814         >>> x_data = numpy.linspace(1e-9, 20e-9, num=10)
 815         >>> model = WLC(data=x_data, info={'x data (m)':x_data})
 816         >>> model.L(model.Lp(21e-9))  # doctest: +ELLIPSIS
 817         2.100...e-08
 818         >>> model.L(model.Lp(100e-9))  # doctest: +ELLIPSIS
 819         9.999...e-08
 820         """
 821         x_max = self.info['x data (m)'].max()
 822         return (numpy.exp(Lp)+1)*x_max
 823
 824     def model(self, params):
 825         """Extract the relavant arguments and use :func:`WLC_fn`.
 826
 827         Parameters
 828         ----------
 829         params : list of floats
 830             `[Lp, p]`, where the presence of `p` is optional.
 831         """
 832         # setup convenient aliases
 833         T = self.info['temperature (K)']
 834         x_data = self.info['x data (m)']
 835         L = self.L(params[0])
 836         if len(params) > 1:
 837             p = abs(params[1])
 838         else:
 839             p = self.info['persistence length (m)']
 840         # compute model data
 841         self._model_data[:] = WLC_fn(x_data, T, L, p)
 842         return self._model_data
 843
 844     def guess_initial_params(self, outqueue=None):
 845         """Guess initial fitting parameters.
 846
 847         Returns
 848         -------
 849         Lp : float
 850             A guess at the reparameterized contour length in meters.
 851         p : float (optional)
 852             A guess at the persistence length in meters.  If `.info`
 853             has a setting for `'persistence length (m)'`, `p` is not
 854             returned.
 855         """
 856         x_max = self.info['x data (m)'].max()
 857         p_given = 'persistence length (m)' in self.info
 858         if p_given:
 859             p = self.info['persistence length (m)']
 860         else:
 861             p = x_max / 10.0
 862         L = 1.5 * x_max
 863         Lp = self.Lp(L)
 864         if p_given:
 865             return [Lp]
 866         return [Lp, p]
 867
 868     def guess_scale(self, params, outqueue=None):
 869         """Guess parameter scales.
 870
 871         Returns
 872         -------
 873         Lp_scale : float
 874             A guess at the reparameterized contour length scale in meters.
 875         p_scale : float (optional)
 876             A guess at the persistence length in meters.  If the
 877             length of `params` is less than 2, `p_scale` is not
 878             returned.
 879         """
 880         Lp_scale = 1.0
 881         if len(params) == 1:
 882             return [Lp_scale]
 883         return [Lp_scale] + [x/10.0 for x in params[1:]]
 884
 885
 886 class PolymerFitPlugin (Plugin):
 887     """Polymer model (WLC, FJC, etc.) fitting.
 888     """
 889     def __init__(self):
 890         super(PolymerFitPlugin, self).__init__(name='polymer_fit')
 891         self._commands = [PolymerFitCommand(self),]
 892
 893     def dependencies(self):
 894         return ['vclamp']
 895
 896     def default_settings(self):
 897         return [
 898             Setting(section=self.setting_section, help=self.__doc__),
 899             Setting(section=self.setting_section,
 900                     option='polymer model',
 901                     value='WLC',
 902                     help="Select the default polymer model for 'polymer fit'.  See the documentation for descriptions of available algorithms."),
 903             Setting(section=self.setting_section,
 904                     option='FJC Kuhn length',
 905                     value=4e-10, type='float',
 906                     help='Kuhn length in meters'),
 907             Setting(section=self.setting_section,
 908                     option='FJC-PEG Kuhn length',
 909                     value=4e-10, type='float',
 910                     help='Kuhn length in meters'),
 911             Setting(section=self.setting_section,
 912                     option='FJC-PEG elasticity',
 913                     value=150.0, type='float',
 914                     help='Elasticity of a PEG segment in Newtons per meter.'),
 915             Setting(section=self.setting_section,
 916                     option='FJC-PEG delta G',
 917                     value=3.0, type='float',
 918                     help='Gibbs free energy difference between trans-trans-trans (ttt) and trans-trans-gauche (ttg) PEG states in units of kBT.'),
 919             Setting(section=self.setting_section,
 920                     option='FJC-PEG L_helical',
 921                     value=2.8e-10, type='float',
 922                     help='Contour length of PEG in the ttg state.'),
 923             Setting(section=self.setting_section,
 924                     option='FJC-PEG L_planar',
 925                     value=3.58e-10, type='float',
 926                     help='Contour length of PEG in the ttt state.'),
 927             Setting(section=self.setting_section,
 928                     option='WLC persistence length',
 929                     value=4e-10, type='float',
 930                     help='Persistence length in meters'),
 931             ]
 932
 933
 934 class PolymerFitCommand (Command):
 935     """Polymer model (WLC, FJC, etc.) fitting.
 936
 937     Fits an entropic elasticity function to a given chunk of the
 938     curve.  Fit quality compares the residual with the thermal noise
 939     (a good fit should be 1 or less).
 940
 941     Because the models are complicated and varied, you should
 942     configure the command by setting configuration options instead of
 943     using command arguments.  TODO: argument_to_setting()
 944     """
 945     def __init__(self, plugin):
 946         super(PolymerFitCommand, self).__init__(
 947             name='polymer fit',
 948             arguments=[
 949                 CurveArgument,
 950                 Argument(name='block', aliases=['set'], type='int', default=0,
 951                          help="""
 952 Data block for which the fit should be calculated.  For an
 953 approach/retract force curve, `0` selects the approaching curve and
 954 `1` selects the retracting curve.
 955 """.strip()),
 956                 Argument(name='bounds', type='point', optional=False, count=2,
 957                          help="""
 958 Indicies of points bounding the selected data.
 959 """.strip()),
 960                 ],
 961             help=self.__doc__, plugin=plugin)
 962
 963     def _run(self, hooke, inqueue, outqueue, params):
 964         scale(hooke, params['curve'], params['block'])  # TODO: is autoscaling a good idea? (explicit is better than implicit)
 965         data = params['curve'].data[params['block']]
 966         # HACK? rely on params['curve'] being bound to the local hooke
 967         # playlist (i.e. not a copy, as you would get by passing a
 968         # curve through the queue).  Ugh.  Stupid queues.  As an
 969         # alternative, we could pass lookup information through the
 970         # queue...
 971         model = self.plugin.config['polymer model']
 972         new = Data((data.shape[0], data.shape[1]+1), dtype=data.dtype)
 973         new.info = copy.deepcopy(data.info)
 974         new[:,:-1] = data
 975         new.info['columns'].append('%s tension (N)' % model)  # TODO: WLC fit for each peak, etc.
 976         z_data = data[:,data.info['columns'].index('surface distance (m)')]
 977         d_data = data[:,data.info['columns'].index('deflection (N)')]
 978         start,stop = params['bounds']
 979         tension_data,ps = self.fit_polymer_model(
 980             params['curve'], z_data, d_data, start, stop, outqueue)
 981         new.info['%s polymer fit parameters' % model] = ps
 982         new[:,-1] = tension_data
 983         params['curve'].data[params['block']] = new
 984
 985     def fit_polymer_model(self, curve, z_data, d_data, start, stop,
 986                           outqueue=None):
 987         """Railyard for the `fit_*_model` family.
 988
 989         Uses the `polymer model` configuration setting to call the
 990         appropriate backend algorithm.
 991         """
 992         fn = getattr(self, 'fit_%s_model'
 993                      % self.plugin.config['polymer model'].replace('-','_'))
 994         return fn(curve, z_data, d_data, start, stop, outqueue)
 995
 996     def fit_FJC_model(self, curve, z_data, d_data, start, stop,
 997                       outqueue=None):
 998         """Fit the data with :class:`FJC`.
 999         """
1000         """Fit the data with :class:`WLC`.
1001         """
1002         info = {
1003             'temperature (K)': self.plugin.config['temperature'],
1004             'x data (m)': z_data,
1005             }
1006         if True:  # TODO: optionally free persistence length
1007             info['Kuhn length (m)'] = (
1008                 self.plugin.config['FJC Kuhn length'])
1009         model = FJC(d_data, info=info)
1010         queue = Queue.Queue()
1011         params = model.fit(outqueue=queue)
1012         if True:  # TODO: if Kuhn length fixed
1013             params = [params]
1014             a = info['Kuhn length (m)']
1015         else:
1016             a = params[1]
1017         Lp = params[0]
1018         L = model.L(Lp)
1019         T = info['temperature (K)']
1020         fit_info = queue.get(block=False)
1021         mask = numpy.zeros(z_data.shape, dtype=numpy.bool)
1022         mask[start:stop] = True
1023         return [FJC_fn(z_data, T=T, L=L, a=a) * mask,
1024                 fit_info]
1025
1026     def fit_FJC_PEG_model(self, curve, z_data, d_data, start, stop,
1027                           outqueue=None):
1028         """Fit the data with :class:`FJC_PEG`.
1029         """
1030         pass
1031
1032     def fit_WLC_model(self, curve, z_data, d_data, start, stop,
1033                       outqueue=None):
1034         """Fit the data with :class:`WLC`.
1035         """
1036         info = {
1037             'temperature (K)': self.plugin.config['temperature'],
1038             'x data (m)': z_data,
1039             }
1040         if True:  # TODO: optionally free persistence length
1041             info['persistence length (m)'] = (
1042                 self.plugin.config['WLC persistence length'])
1043         model = WLC(d_data, info=info)
1044         queue = Queue.Queue()
1045         params = model.fit(outqueue=queue)
1046         if True:  # TODO: if persistence length fixed
1047             params = [params]
1048             p = info['persistence length (m)']
1049         else:
1050             p = params[1]
1051         Lp = params[0]
1052         L = model.L(Lp)
1053         T = info['temperature (K)']
1054         fit_info = queue.get(block=False)
1055         mask = numpy.zeros(z_data.shape, dtype=numpy.bool)
1056         mask[start:stop] = True
1057         return [WLC_fn(z_data, T=T, L=L, p=p) * mask,
1058                 fit_info]
1059
1060
1061 # TODO:
1062 # def dist2fit(self):
1063 #     '''Calculates the average distance from data to fit, scaled by
1064 #     the standard deviation of the free cantilever area (thermal
1065 #     noise).
1066 #     '''