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
  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 eFJC via Newton's method
 458     (:func:`scipy.optimize.newton`).  For the starting point, we use
 459     the standard FJC function with an averaged contour length.
 460     """
 461     kwargs = {'T':T, 'N':N, 'k':k, 'Lp':Lp, 'Lh':Lh, 'dG':dG, 'a':a}
 462     if is_iterable(x_data):
 463         ret = numpy.ndarray(shape=x_data.shape, dtype=x_data.dtype)
 464         for i,x in enumerate(x_data):
 465             ret[i] = FJC_PEG_fn(x, **kwargs)
 466         return ret
 467     if x_data == 0:
 468         return 0
 469     # Approximate with the simple FJC_fn()
 470     guess = numpy.inf
 471     L = N*max(Lp, Lh)
 472     while guess == numpy.inf:
 473         guess = FJC_fn(x_data, T=T, L=L, a=a)
 474         L *= 2.0
 475
 476     fn = lambda f: inverse_FJC_PEG_fn(guess*abs(f), **kwargs) - x_data
 477     ret = guess*abs(newton(func=fn, x0=1.0))
 478     return ret
 479
 480 class FJC_PEG (ModelFitter):
 481     """Fit the data to an poly(ethylene-glycol) adjusted extended FJC.
 482
 483     Examples
 484     --------
 485     Generate some example data
 486
 487     >>> kwargs = {'T':300.0, 'N':123, 'k':150.0, 'Lp':3.58e-10, 'Lh':2.8e-10,
 488     ...           'dG':3.0, 'a':7e-10}
 489     >>> x_data = numpy.linspace(10e-9, 30e-9, num=20)
 490     >>> d_data = FJC_PEG_fn(x_data, **kwargs)
 491
 492     Fit the example data with a two-parameter fit (`N` and `a`).
 493
 494     >>> info = {
 495     ...     'temperature (K)': kwargs['T'],
 496     ...     'x data (m)': x_data,
 497     ...     'section elasticity (N/m)': kwargs['k'],
 498     ...     'planar section length (m)': kwargs['Lp'],
 499     ...     'helical section length (m)': kwargs['Lh'],
 500     ...     'Gibbs free energy difference (Gp - Gh) (kBT)': kwargs['dG'],
 501     ...     }
 502     >>> model = FJC_PEG(d_data, info=info, rescale=True)
 503     >>> outqueue = Queue.Queue()
 504     >>> Nr,a = model.fit(outqueue=outqueue)
 505     >>> fit_info = outqueue.get(block=False)
 506     >>> model.L(Nr)  # doctest: +ELLIPSIS
 507     122.999...
 508     >>> a  # doctest: +ELLIPSIS
 509     6.999...e-10
 510
 511     Fit the example data with a one-parameter fit (`N`).  We introduce
 512     some error in our fixed Kuhn length for fun.
 513
 514     >>> info['Kuhn length (m)'] = 2*kwargs['a']
 515     >>> model = FJC_PEG(d_data, info=info, rescale=True)
 516     >>> Nr = model.fit(outqueue=outqueue)
 517     >>> fit_info = outqueue.get(block=False)
 518     >>> model.L(Nr)  # doctest: +ELLIPSIS
 519     96.931...
 520     """
 521     def Lr(self, L):
 522         """To avoid invalid values of `L`, we reparametrize with `Lr`.
 523
 524         Parameters
 525         ----------
 526         L : float
 527             contour length in meters
 528
 529         Returns
 530         -------
 531         Lr : float
 532             rescaled version of `L`.
 533
 534         Examples
 535         --------
 536         >>> x_data = numpy.linspace(1e-9, 20e-9, num=10)
 537         >>> model = FJC_PEG(data=x_data, info={'x data (m)':x_data})
 538         >>> model.Lr(20e-9)  # doctest: +ELLIPSIS
 539         0.0
 540         >>> model.Lr(19e-9)  # doctest: +ELLIPSIS
 541         -0.0512...
 542         >>> model.Lr(21e-9)  # doctest: +ELLIPSIS
 543         0.0487...
 544         >>> model.Lr(100e-9)  # doctest: +ELLIPSIS
 545         1.609...
 546
 547         Notes
 548         -----
 549         The rescaling is designed to limit `L` to values strictly
 550         greater than zero, so we use
 551
 552         .. math::
 553             L = \exp(L_p) * x_\text{max}
 554
 555         which has the inverse
 556
 557         .. math::
 558             L_p = \ln(L/x_\text{max})
 559
 560         This will obviously effect the interpretation of the fit's covariance matrix.
 561         """
 562         x_max = self.info['x data (m)'].max()
 563         return numpy.log(L/x_max)
 564
 565     def L(self, Lr):
 566         """Inverse of :meth:`Lr`.
 567
 568         Parameters
 569         ----------
 570         Lr : float
 571             rescaled version of `L`
 572
 573         Returns
 574         -------
 575         L : float
 576             contour length in meters
 577
 578         Examples
 579         --------
 580         >>> x_data = numpy.linspace(1e-9, 20e-9, num=10)
 581         >>> model = FJC_PEG(data=x_data, info={'x data (m)':x_data})
 582         >>> model.L(model.Lr(21e-9))  # doctest: +ELLIPSIS
 583         2.100...e-08
 584         >>> model.L(model.Lr(100e-9))  # doctest: +ELLIPSIS
 585         9.999...e-08
 586         """
 587         x_max = self.info['x data (m)'].max()
 588         return numpy.exp(Lr)*x_max
 589
 590     def _kwargs(self):
 591         return {
 592             'T': self.info['temperature (K)'],
 593             'k': self.info['section elasticity (N/m)'],
 594             'Lp': self.info['planar section length (m)'],
 595             'Lh': self.info['helical section length (m)'],
 596             'dG': self.info['Gibbs free energy difference (Gp - Gh) (kBT)'],
 597             }
 598
 599     def model(self, params):
 600         """Extract the relavant arguments and use :func:`FJC_PEG_fn`.
 601
 602         Parameters
 603         ----------
 604         params : list of floats
 605             `[N, a]`, where the presence of `a` is optional.
 606         """
 607         # setup convenient aliases
 608         T = self.info['temperature (K)']
 609         x_data = self.info['x data (m)']
 610         N = self.L(params[0])
 611         if len(params) > 1:
 612             a = abs(params[1])
 613         else:
 614             a = self.info['Kuhn length (m)']
 615         kwargs = self._kwargs()
 616         # compute model data
 617         self._model_data[:] = FJC_PEG_fn(x_data, N=N, a=a, **kwargs)
 618         return self._model_data
 619
 620     def guess_initial_params(self, outqueue=None):
 621         """Guess initial fitting parameters.
 622
 623         Returns
 624         -------
 625         N : float
 626             A guess at the section count.
 627         a : float (optional)
 628             A guess at the Kuhn length in meters.  If `.info` has a
 629             setting for `'Kuhn length (m)'`, `a` is not returned.
 630         """
 631         x_max = self.info['x data (m)'].max()
 632         a_given = 'Kuhn length (m)' in self.info
 633         if a_given:
 634             a = self.info['Kuhn length (m)']
 635         else:
 636             a = x_max / 10.0
 637         f_max = self._data.max()
 638         kwargs = self._kwargs()
 639         kwargs['a'] = a
 640         x_section = inverse_FJC_PEG_fn(F_data=f_max, **kwargs)
 641         N = x_max / x_section;
 642         Nr = self.Lr(N)
 643         if a_given:
 644             return [Nr]
 645         return [Nr, a]
 646
 647     def guess_scale(self, params, outqueue=None):
 648         """Guess parameter scales.
 649
 650         Returns
 651         -------
 652         N_scale : float
 653             A guess at the section count scale in meters.
 654         a_scale : float (optional)
 655             A guess at the Kuhn length in meters.  If the length of
 656             `params` is less than 2, `a_scale` is not returned.
 657         """
 658         return [x/10.0 for x in params]
 659
 660
 661 def WLC_fn(x_data, T, L, p):
 662     """The worm like chain interpolation model.
 663
 664     Parameters
 665     ----------
 666     x_data : array_like
 667         x values for which the WLC tension should be calculated.
 668     T : float
 669         temperature in Kelvin.
 670     L : float
 671         contour length in meters.
 672     p : float
 673         persistence length in meters.
 674
 675     Returns
 676     -------
 677     f_data : array
 678         `F(x_data)`.
 679
 680     Examples
 681     --------
 682     >>> WLC_fn(numpy.array([1e-9, 5e-9, 10e-9], dtype=numpy.float),
 683     ...        T=300, L=15e-9, p=2.5e-10) # doctest: +ELLIPSIS
 684     array([  1.717...e-12,   1.070...e-11,   4.418...e-11])
 685
 686     Notes
 687     -----
 688     The function is the simple polynomial worm like chain
 689     interpolation forumula proposed by C. Bustamante, et al. [#]_.
 690
 691     .. math::
 692       F(x) = \frac{k_B T}{p} \left[
 693         \frac{1}{4}\left( \frac{1}{(1-x/L)^2} - 1 \right)
 694         + \frac{x}{L}
 695                              \right]
 696
 697     .. [#] C. Bustamante, J.F. Marko, E.D. Siggia, and S.Smith.
 698     "Entropic elasticity of lambda-phage DNA."
 699     Science. 1994.
 700     doi: `10.1126/science.8079175 <http://dx.doi.org/10.1126/science.8079175>`
 701     """
 702     a = kB * T / p
 703     scaled_data = x_data / L
 704     return a * (0.25*((1-scaled_data)**-2 - 1) + scaled_data)
 705
 706 class WLC (ModelFitter):
 707     """Fit the data to a worm like chain.
 708
 709     Examples
 710     --------
 711     Generate some example data
 712
 713     >>> T = 300  # Kelvin
 714     >>> L = 35e-9  # meters
 715     >>> p = 2.5e-10  # meters
 716     >>> x_data = numpy.linspace(10e-9, 30e-9, num=20)
 717     >>> d_data = WLC_fn(x_data, T=T, L=L, p=p)
 718
 719     Fit the example data with a two-parameter fit (`L` and `p`).
 720
 721     >>> info = {
 722     ...     'temperature (K)': T,
 723     ...     'x data (m)': x_data,
 724     ...     }
 725     >>> model = WLC(d_data, info=info, rescale=True)
 726     >>> outqueue = Queue.Queue()
 727     >>> Lp,p = model.fit(outqueue=outqueue)
 728     >>> fit_info = outqueue.get(block=False)
 729     >>> model.L(Lp)  # doctest: +ELLIPSIS
 730     3.500...e-08
 731     >>> p  # doctest: +ELLIPSIS
 732     2.500...e-10
 733
 734     Fit the example data with a one-parameter fit (`L`).  We introduce
 735     some error in our fixed persistence length for fun.
 736
 737     >>> info['persistence length (m)'] = 2*p
 738     >>> model = WLC(d_data, info=info, rescale=True)
 739     >>> Lp = model.fit(outqueue=outqueue)
 740     >>> fit_info = outqueue.get(block=False)
 741     >>> model.L(Lp)  # doctest: +ELLIPSIS
 742     3.318...e-08
 743     """
 744     def Lp(self, L):
 745         """To avoid invalid values of `L`, we reparametrize with `Lp`.
 746
 747         Parameters
 748         ----------
 749         L : float
 750             contour length in meters
 751
 752         Returns
 753         -------
 754         Lp : float
 755             rescaled version of `L`.
 756
 757         Examples
 758         --------
 759         >>> x_data = numpy.linspace(1e-9, 20e-9, num=10)
 760         >>> model = WLC(data=x_data, info={'x data (m)':x_data})
 761         >>> model.Lp(20e-9)  # doctest: +ELLIPSIS
 762         -inf
 763         >>> model.Lp(19e-9)  # doctest: +ELLIPSIS
 764         nan
 765         >>> model.Lp(21e-9)  # doctest: +ELLIPSIS
 766         -2.99...
 767         >>> model.Lp(100e-9)  # doctest: +ELLIPSIS
 768         1.386...
 769
 770         Notes
 771         -----
 772         The rescaling is designed to limit `L` to values strictly
 773         greater than the maximum `x` data value, so we use
 774
 775         .. math::
 776             L = [\exp(L_p))+1] * x_\text{max}
 777
 778         which has the inverse
 779
 780         .. math::
 781             L_p = \ln(L/x_\text{max}-1)
 782
 783         This will obviously effect the interpretation of the fit's covariance matrix.
 784         """
 785         x_max = self.info['x data (m)'].max()
 786         return numpy.log(L/x_max - 1)
 787
 788     def L(self, Lp):
 789         """Inverse of :meth:`Lp`.
 790
 791         Parameters
 792         ----------
 793         Lp : float
 794             rescaled version of `L`
 795
 796         Returns
 797         -------
 798         L : float
 799             contour length in meters
 800
 801         Examples
 802         --------
 803         >>> x_data = numpy.linspace(1e-9, 20e-9, num=10)
 804         >>> model = WLC(data=x_data, info={'x data (m)':x_data})
 805         >>> model.L(model.Lp(21e-9))  # doctest: +ELLIPSIS
 806         2.100...e-08
 807         >>> model.L(model.Lp(100e-9))  # doctest: +ELLIPSIS
 808         9.999...e-08
 809         """
 810         x_max = self.info['x data (m)'].max()
 811         return (numpy.exp(Lp)+1)*x_max
 812
 813     def model(self, params):
 814         """Extract the relavant arguments and use :func:`WLC_fn`.
 815
 816         Parameters
 817         ----------
 818         params : list of floats
 819             `[Lp, p]`, where the presence of `p` is optional.
 820         """
 821         # setup convenient aliases
 822         T = self.info['temperature (K)']
 823         x_data = self.info['x data (m)']
 824         L = self.L(params[0])
 825         if len(params) > 1:
 826             p = abs(params[1])
 827         else:
 828             p = self.info['persistence length (m)']
 829         # compute model data
 830         self._model_data[:] = WLC_fn(x_data, T, L, p)
 831         return self._model_data
 832
 833     def guess_initial_params(self, outqueue=None):
 834         """Guess initial fitting parameters.
 835
 836         Returns
 837         -------
 838         Lp : float
 839             A guess at the reparameterized contour length in meters.
 840         p : float (optional)
 841             A guess at the persistence length in meters.  If `.info`
 842             has a setting for `'persistence length (m)'`, `p` is not
 843             returned.
 844         """
 845         x_max = self.info['x data (m)'].max()
 846         p_given = 'persistence length (m)' in self.info
 847         if p_given:
 848             p = self.info['persistence length (m)']
 849         else:
 850             p = x_max / 10.0
 851         L = 1.5 * x_max
 852         Lp = self.Lp(L)
 853         if p_given:
 854             return [Lp]
 855         return [Lp, p]
 856
 857     def guess_scale(self, params, outqueue=None):
 858         """Guess parameter scales.
 859
 860         Returns
 861         -------
 862         Lp_scale : float
 863             A guess at the reparameterized contour length scale in meters.
 864         p_scale : float (optional)
 865             A guess at the persistence length in meters.  If the
 866             length of `params` is less than 2, `p_scale` is not
 867             returned.
 868         """
 869         Lp_scale = 1.0
 870         if len(params) == 1:
 871             return [Lp_scale]
 872         return [Lp_scale] + [x/10.0 for x in params[1:]]
 873
 874
 875 class PolymerFitPlugin (Plugin):
 876     """Polymer model (WLC, FJC, etc.) fitting.
 877     """
 878     def __init__(self):
 879         super(PolymerFitPlugin, self).__init__(name='polymer_fit')
 880         self._commands = [PolymerFitCommand(self),]
 881
 882     def dependencies(self):
 883         return ['vclamp']
 884
 885     def default_settings(self):
 886         return [
 887             Setting(section=self.setting_section, help=self.__doc__),
 888             Setting(section=self.setting_section,
 889                     option='polymer model',
 890                     value='WLC',
 891                     help="Select the default polymer model for 'polymer fit'.  See the documentation for descriptions of available algorithms."),
 892             Setting(section=self.setting_section,
 893                     option='FJC Kuhn length',
 894                     value=4e-10, type='float',
 895                     help='Kuhn length in meters'),
 896             Setting(section=self.setting_section,
 897                     option='FJC-PEG Kuhn length',
 898                     value=4e-10, type='float',
 899                     help='Kuhn length in meters'),
 900             Setting(section=self.setting_section,
 901                     option='FJC-PEG elasticity',
 902                     value=150.0, type='float',
 903                     help='Elasticity of a PEG segment in Newtons per meter.'),
 904             Setting(section=self.setting_section,
 905                     option='FJC-PEG delta G',
 906                     value=3.0, type='float',
 907                     help='Gibbs free energy difference between trans-trans-trans (ttt) and trans-trans-gauche (ttg) PEG states in units of kBT.'),
 908             Setting(section=self.setting_section,
 909                     option='FJC-PEG L_helical',
 910                     value=2.8e-10, type='float',
 911                     help='Contour length of PEG in the ttg state.'),
 912             Setting(section=self.setting_section,
 913                     option='FJC-PEG L_planar',
 914                     value=3.58e-10, type='float',
 915                     help='Contour length of PEG in the ttt state.'),
 916             Setting(section=self.setting_section,
 917                     option='WLC persistence length',
 918                     value=4e-10, type='float',
 919                     help='Persistence length in meters'),
 920             ]
 921
 922
 923 class PolymerFitCommand (Command):
 924     """Polymer model (WLC, FJC, etc.) fitting.
 925
 926     Fits an entropic elasticity function to a given chunk of the
 927     curve.  Fit quality compares the residual with the thermal noise
 928     (a good fit should be 1 or less).
 929
 930     Because the models are complicated and varied, you should
 931     configure the command by setting configuration options instead of
 932     using command arguments.  TODO: argument_to_setting()
 933     """
 934     def __init__(self, plugin):
 935         super(PolymerFitCommand, self).__init__(
 936             name='polymer fit',
 937             arguments=[
 938                 CurveArgument,
 939                 Argument(name='block', aliases=['set'], type='int', default=0,
 940                          help="""
 941 Data block for which the fit should be calculated.  For an
 942 approach/retract force curve, `0` selects the approaching curve and
 943 `1` selects the retracting curve.
 944 """.strip()),
 945                 Argument(name='bounds', type='point', optional=False, count=2,
 946                          help="""
 947 Indicies of points bounding the selected data.
 948 """.strip()),
 949                 ],
 950             help=self.__doc__, plugin=plugin)
 951
 952     def _run(self, hooke, inqueue, outqueue, params):
 953         scale(hooke, params['curve'], params['block'])  # TODO: is autoscaling a good idea? (explicit is better than implicit)
 954         data = params['curve'].data[params['block']]
 955         # HACK? rely on params['curve'] being bound to the local hooke
 956         # playlist (i.e. not a copy, as you would get by passing a
 957         # curve through the queue).  Ugh.  Stupid queues.  As an
 958         # alternative, we could pass lookup information through the
 959         # queue...
 960         model = self.plugin.config['polymer model']
 961         new = Data((data.shape[0], data.shape[1]+1), dtype=data.dtype)
 962         new.info = copy.deepcopy(data.info)
 963         new[:,:-1] = data
 964         new.info['columns'].append('%s tension (N)' % model)  # TODO: WLC fit for each peak, etc.
 965         z_data = data[:,data.info['columns'].index(
 966                 'cantilever adjusted extension (m)')]
 967         d_data = data[:,data.info['columns'].index('deflection (N)')]
 968         start,stop = params['bounds']
 969         tension_data,ps = self.fit_polymer_model(
 970             params['curve'], z_data, d_data, start, stop, outqueue)
 971         new.info['%s polymer fit parameters' % model] = ps
 972         new[:,-1] = tension_data
 973         params['curve'].data[params['block']] = new
 974
 975     def fit_polymer_model(self, curve, z_data, d_data, start, stop,
 976                           outqueue=None):
 977         """Railyard for the `fit_*_model` family.
 978
 979         Uses the `polymer model` configuration setting to call the
 980         appropriate backend algorithm.
 981         """
 982         fn = getattr(self, 'fit_%s_model'
 983                      % self.plugin.config['polymer model'].replace('-','_'))
 984         return fn(curve, z_data, d_data, start, stop, outqueue)
 985
 986     def fit_FJC_model(self, curve, z_data, d_data, start, stop,
 987                       outqueue=None):
 988         """Fit the data with :class:`FJC`.
 989         """
 990         info = {
 991             'temperature (K)': self.plugin.config['temperature'],
 992             'x data (m)': z_data[start:stop],
 993             }
 994         if True:  # TODO: optionally free persistence length
 995             info['Kuhn length (m)'] = (
 996                 self.plugin.config['FJC Kuhn length'])
 997         model = FJC(d_data[start:stop], info=info, rescale=True)
 998         queue = Queue.Queue()
 999         params = model.fit(outqueue=queue)
1000         if True:  # TODO: if Kuhn length fixed
1001             params = [params]
1002             a = info['Kuhn length (m)']
1003         else:
1004             a = params[1]
1005         Lp = params[0]
1006         L = model.L(Lp)
1007         T = info['temperature (K)']
1008         fit_info = queue.get(block=False)
1009         mask = numpy.zeros(z_data.shape, dtype=numpy.bool)
1010         mask[start:stop] = True
1011         return [FJC_fn(z_data, T=T, L=L, a=a) * mask,
1012                 fit_info]
1013
1014     def fit_FJC_PEG_model(self, curve, z_data, d_data, start, stop,
1015                           outqueue=None):
1016         """Fit the data with :class:`FJC_PEG`.
1017         """
1018         info = {
1019             'temperature (K)': self.plugin.config['temperature'],
1020             'x data (m)': z_data[start:stop],
1021             # TODO: more info
1022             }
1023         if True:  # TODO: optionally free persistence length
1024             info['Kuhn length (m)'] = (
1025                 self.plugin.config['FJC Kuhn length'])
1026         model = FJC(d_data[start:stop], info=info, rescale=True)
1027         queue = Queue.Queue()
1028         params = model.fit(outqueue=queue)
1029         if True:  # TODO: if Kuhn length fixed
1030             params = [params]
1031             a = info['Kuhn length (m)']
1032         else:
1033             a = params[1]
1034         Nr = params[0]
1035         N = model.L(Nr)
1036         T = info['temperature (K)']
1037         fit_info = queue.get(block=False)
1038         mask = numpy.zeros(z_data.shape, dtype=numpy.bool)
1039         mask[start:stop] = True
1040         return [FJC_PEG_fn(z_data, **kwargs) * mask,
1041                 fit_info]
1042
1043     def fit_WLC_model(self, curve, z_data, d_data, start, stop,
1044                       outqueue=None):
1045         """Fit the data with :class:`WLC`.
1046         """
1047         info = {
1048             'temperature (K)': self.plugin.config['temperature'],
1049             'x data (m)': z_data[start:stop],
1050             }
1051         if True:  # TODO: optionally free persistence length
1052             info['persistence length (m)'] = (
1053                 self.plugin.config['WLC persistence length'])
1054         model = WLC(d_data[start:stop], info=info, rescale=True)
1055         queue = Queue.Queue()
1056         params = model.fit(outqueue=queue)
1057         if True:  # TODO: if persistence length fixed
1058             params = [params]
1059             p = info['persistence length (m)']
1060         else:
1061             p = params[1]
1062         Lp = params[0]
1063         L = model.L(Lp)
1064         T = info['temperature (K)']
1065         fit_info = queue.get(block=False)
1066         mask = numpy.zeros(z_data.shape, dtype=numpy.bool)
1067         mask[start:stop] = True
1068         return [WLC_fn(z_data, T=T, L=L, p=p) * mask,
1069                 fit_info]
1070
1071
1072 # TODO:
1073 # def dist2fit(self):
1074 #     '''Calculates the average distance from data to fit, scaled by
1075 #     the standard deviation of the free cantilever area (thermal
1076 #     noise).
1077 #     '''