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