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     >>> Lp,a = model.fit(outqueue=outqueue)
 237     >>> fit_info = outqueue.get(block=False)
 238     >>> model.L(Lp)  # doctest: +ELLIPSIS
 239     3.500...e-08
 240     >>> a  # doctest: +ELLIPSIS
 241     2.499...e-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     >>> Lp = model.fit(outqueue=outqueue)
 249     >>> fit_info = outqueue.get(block=False)
 250     >>> model.L(Lp)  # 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 guess_initial_params(self, outqueue=None):
 343         """Guess initial fitting parameters.
 344
 345         Returns
 346         -------
 347         Lp : float
 348             A guess at the reparameterized contour length in meters.
 349         a : float (optional)
 350             A guess at the Kuhn length in meters.  If `.info` has a
 351             setting for `'Kuhn length (m)'`, `a` is not returned.
 352         """
 353         x_max = self.info['x data (m)'].max()
 354         a_given = 'Kuhn length (m)' in self.info
 355         if a_given:
 356             a = self.info['Kuhn length (m)']
 357         else:
 358             a = x_max / 10.0
 359         L = 1.5 * x_max
 360         Lp = self.Lp(L)
 361         if a_given:
 362             return [Lp]
 363         return [Lp, a]
 364
 365     def guess_scale(self, params, outqueue=None):
 366         """Guess parameter scales.
 367
 368         Returns
 369         -------
 370         Lp_scale : float
 371             A guess at the reparameterized contour length scale in meters.
 372         a_scale : float (optional)
 373             A guess at the Kuhn length in meters.  If the length of
 374             `params` is less than 2, `a_scale` is not returned.
 375         """
 376         Lp_scale = 1.0
 377         if len(params) == 1:
 378             return [Lp_scale]
 379         return [Lp_scale] + [x/10.0 for x in params[1:]]
 380
 381
 382 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):
 383     """Inverse poly(ethylene-glycol) adjusted extended FJC model.
 384
 385     Examples
 386     --------
 387     >>> kwargs = {'T':300.0, 'N':1, 'k':150.0, 'Lp':3.58e-10, 'Lh':2.8e-10,
 388     ...           'dG':3.0, 'a':7e-10}
 389     >>> inverse_FJC_PEG_fn(F_data=200e-12, **kwargs)  # doctest: +ELLIPSIS
 390     3.487...e-10
 391
 392     Notes
 393     -----
 394     The function is that proposed by F. Oesterhelt, et al. [#]_.
 395
 396     .. math::
 397       x(F) = N_\text{S} \cdot \left[
 398         \left(
 399             \frac{L_\text{planar}}{
 400                   \exp\left(-\Delta G/k_B T\right) + 1}
 401             + \frac{L_\text{helical}}{
 402                   \exp\left(+\Delta G/k_B T\right) + 1}
 403           \right) \cdot \left(
 404             \coth\left(\frac{Fa}{k_B T}\right)
 405             - \frac{k_B T}{Fa}
 406           \right)
 407         + \frac{F}{K_\text{S}}
 408
 409     where :math:`N_\text{S}` is the number of segments,
 410     :math:`K_\text{S}` is the segment elasticicty,
 411     :math:`L_\text{planar}` is the ttt contour length,
 412     :math:`L_\text{helical}` is the ttg contour length,
 413     :math:`\Delta G` is the Gibbs free energy difference
 414      :math:`G_\text{planar}-G_\text{helical}`,
 415     :math:`a` is the Kuhn length,
 416     and :math:`F` is the chain tension.
 417
 418     .. [#]: F. Oesterhelt, M. Rief, and H.E. Gaub.
 419       "Single molecule force spectroscopy by AFM indicates helical
 420       structure of poly(ethylene-glycol) in water."
 421       New Journal of Physics. 1999.
 422       doi: `10.1088/1367-2630/1/1/006 <http://dx.doi.org/10.1088/1367-2630/1/1/006>`
 423       (section 4.2)
 424     """
 425     kBT = kB * T
 426     g = (dG - F_data*(Lp-Lh)) / kBT
 427     z = F_data*a/kBT
 428     return N * ((Lp/(numpy.exp(-g)+1) + Lh/(numpy.exp(+g)+1)) * (coth(z)-1./z)
 429                 + F_data/k)
 430
 431 def FJC_PEG_fn(x_data, T, N, k, Lp, Lh, dG, a):
 432     """Poly(ethylene-glycol) adjusted extended FJC model.
 433
 434     Parameters
 435     ----------
 436     z : float or array_like
 437         object whose inverse Langevin will be returned
 438     extreme : float
 439         value (close to one) for which we assume the inverse is +/-
 440         infinity.  This avoids problems with Newton-Raphson
 441         convergence in regions with low slope.
 442
 443     Examples
 444     --------
 445     >>> kwargs = {'T':300.0, 'N':1, 'k':150.0, 'Lp':3.58e-10, 'Lh':2.8e-10,
 446     ...           'dG':3.0, 'a':7e-10}
 447     >>> FJC_PEG_fn(inverse_FJC_PEG_fn(1e-11, **kwargs), **kwargs)  # doctest: +ELLIPSIS
 448     9.999...e-12
 449     >>> FJC_PEG_fn(inverse_FJC_PEG_fn(3.4e-10, **kwargs), **kwargs)  # doctest: +ELLIPSIS
 450     3.400...e-10
 451     >>> FJC_PEG_fn(numpy.array([1e-10,2e-10,3e-10]), **kwargs)  # doctest: +ELLIPSIS
 452     array([  5.207...e-12,   1.257...e-11,   3.636...e-11])
 453     >>> kwargs['N'] = 123
 454     >>> FJC_PEG_fn(numpy.array([1e-8,2e-8,3e-8]), **kwargs)  # doctest: +ELLIPSIS
 455     array([  4.160...e-12,   9.302...e-12,   1.830...e-11])
 456
 457     Notes
 458     -----
 459     We approximate the PEG adjusted eFJC via Newton's method
 460     (:func:`scipy.optimize.newton`).  For the starting point, we use
 461     the standard FJC function with an averaged contour length.
 462     """
 463     kwargs = {'T':T, 'N':N, 'k':k, 'Lp':Lp, 'Lh':Lh, 'dG':dG, 'a':a}
 464     if is_iterable(x_data):
 465         ret = numpy.ndarray(shape=x_data.shape, dtype=x_data.dtype)
 466         for i,x in enumerate(x_data):
 467             ret[i] = FJC_PEG_fn(x, **kwargs)
 468         return ret
 469     if x_data == 0:
 470         return 0
 471     # Approximate with the simple FJC_fn()
 472     guess = numpy.inf
 473     L = N*max(Lp, Lh)
 474     while guess == numpy.inf:
 475         guess = FJC_fn(x_data, T=T, L=L, a=a)
 476         L *= 2.0
 477
 478     fn = lambda f: inverse_FJC_PEG_fn(guess*abs(f), **kwargs) - x_data
 479     ret = guess*abs(newton(func=fn, x0=1.0))
 480     return ret
 481
 482 class FJC_PEG (ModelFitter):
 483     """Fit the data to an poly(ethylene-glycol) adjusted extended FJC.
 484
 485     Examples
 486     --------
 487     Generate some example data
 488
 489     >>> kwargs = {'T':300.0, 'N':123, 'k':150.0, 'Lp':3.58e-10, 'Lh':2.8e-10,
 490     ...           'dG':3.0, 'a':7e-10}
 491     >>> x_data = numpy.linspace(10e-9, 30e-9, num=20)
 492     >>> d_data = FJC_PEG_fn(x_data, **kwargs)
 493
 494     Fit the example data with a two-parameter fit (`N` and `a`).
 495
 496     >>> info = {
 497     ...     'temperature (K)': kwargs['T'],
 498     ...     'x data (m)': x_data,
 499     ...     'section elasticity (N/m)': kwargs['k'],
 500     ...     'planar section length (m)': kwargs['Lp'],
 501     ...     'helical section length (m)': kwargs['Lh'],
 502     ...     'Gibbs free energy difference (Gp - Gh) (kBT)': kwargs['dG'],
 503     ...     }
 504     >>> model = FJC_PEG(d_data, info=info, rescale=True)
 505     >>> outqueue = Queue()
 506     >>> Nr,a = model.fit(outqueue=outqueue)
 507     >>> fit_info = outqueue.get(block=False)
 508     >>> model.L(Nr)  # doctest: +ELLIPSIS
 509     122.999...
 510     >>> a  # doctest: +ELLIPSIS
 511     6.999...e-10
 512
 513     Fit the example data with a one-parameter fit (`N`).  We introduce
 514     some error in our fixed Kuhn length for fun.
 515
 516     >>> info['Kuhn length (m)'] = 2*kwargs['a']
 517     >>> model = FJC_PEG(d_data, info=info, rescale=True)
 518     >>> Nr = model.fit(outqueue=outqueue)
 519     >>> fit_info = outqueue.get(block=False)
 520     >>> model.L(Nr)  # doctest: +ELLIPSIS
 521     96.931...
 522     """
 523     def Lr(self, L):
 524         """To avoid invalid values of `L`, we reparametrize with `Lr`.
 525
 526         Parameters
 527         ----------
 528         L : float
 529             contour length in meters
 530
 531         Returns
 532         -------
 533         Lr : float
 534             rescaled version of `L`.
 535
 536         Examples
 537         --------
 538         >>> x_data = numpy.linspace(1e-9, 20e-9, num=10)
 539         >>> model = FJC_PEG(data=x_data, info={'x data (m)':x_data})
 540         >>> model.Lr(20e-9)  # doctest: +ELLIPSIS
 541         0.0
 542         >>> model.Lr(19e-9)  # doctest: +ELLIPSIS
 543         -0.0512...
 544         >>> model.Lr(21e-9)  # doctest: +ELLIPSIS
 545         0.0487...
 546         >>> model.Lr(100e-9)  # doctest: +ELLIPSIS
 547         1.609...
 548
 549         Notes
 550         -----
 551         The rescaling is designed to limit `L` to values strictly
 552         greater than zero, so we use
 553
 554         .. math::
 555             L = \exp(L_p) * x_\text{max}
 556
 557         which has the inverse
 558
 559         .. math::
 560             L_p = \ln(L/x_\text{max})
 561
 562         This will obviously effect the interpretation of the fit's covariance matrix.
 563         """
 564         x_max = self.info['x data (m)'].max()
 565         return numpy.log(L/x_max)
 566
 567     def L(self, Lr):
 568         """Inverse of :meth:`Lr`.
 569
 570         Parameters
 571         ----------
 572         Lr : float
 573             rescaled version of `L`
 574
 575         Returns
 576         -------
 577         L : float
 578             contour length in meters
 579
 580         Examples
 581         --------
 582         >>> x_data = numpy.linspace(1e-9, 20e-9, num=10)
 583         >>> model = FJC_PEG(data=x_data, info={'x data (m)':x_data})
 584         >>> model.L(model.Lr(21e-9))  # doctest: +ELLIPSIS
 585         2.100...e-08
 586         >>> model.L(model.Lr(100e-9))  # doctest: +ELLIPSIS
 587         9.999...e-08
 588         """
 589         x_max = self.info['x data (m)'].max()
 590         return numpy.exp(Lr)*x_max
 591
 592     def _kwargs(self):
 593         return {
 594             'T': self.info['temperature (K)'],
 595             'k': self.info['section elasticity (N/m)'],
 596             'Lp': self.info['planar section length (m)'],
 597             'Lh': self.info['helical section length (m)'],
 598             'dG': self.info['Gibbs free energy difference (Gp - Gh) (kBT)'],
 599             }
 600
 601     def model(self, params):
 602         """Extract the relavant arguments and use :func:`FJC_PEG_fn`.
 603
 604         Parameters
 605         ----------
 606         params : list of floats
 607             `[N, a]`, where the presence of `a` is optional.
 608         """
 609         # setup convenient aliases
 610         T = self.info['temperature (K)']
 611         x_data = self.info['x data (m)']
 612         N = self.L(params[0])
 613         if len(params) > 1:
 614             a = abs(params[1])
 615         else:
 616             a = self.info['Kuhn length (m)']
 617         kwargs = self._kwargs()
 618         # compute model data
 619         self._model_data[:] = FJC_PEG_fn(x_data, N=N, a=a, **kwargs)
 620         return self._model_data
 621
 622     def guess_initial_params(self, outqueue=None):
 623         """Guess initial fitting parameters.
 624
 625         Returns
 626         -------
 627         N : float
 628             A guess at the section count.
 629         a : float (optional)
 630             A guess at the Kuhn length in meters.  If `.info` has a
 631             setting for `'Kuhn length (m)'`, `a` is not returned.
 632         """
 633         x_max = self.info['x data (m)'].max()
 634         a_given = 'Kuhn length (m)' in self.info
 635         if a_given:
 636             a = self.info['Kuhn length (m)']
 637         else:
 638             a = x_max / 10.0
 639         f_max = self._data.max()
 640         kwargs = self._kwargs()
 641         kwargs['a'] = a
 642         x_section = inverse_FJC_PEG_fn(F_data=f_max, **kwargs)
 643         N = x_max / x_section;
 644         Nr = self.Lr(N)
 645         if a_given:
 646             return [Nr]
 647         return [Nr, a]
 648
 649     def guess_scale(self, params, outqueue=None):
 650         """Guess parameter scales.
 651
 652         Returns
 653         -------
 654         N_scale : float
 655             A guess at the section count scale in meters.
 656         a_scale : float (optional)
 657             A guess at the Kuhn length in meters.  If the length of
 658             `params` is less than 2, `a_scale` is not returned.
 659         """
 660         return [x/10.0 for x in params]
 661
 662
 663 def WLC_fn(x_data, T, L, p):
 664     """The worm like chain interpolation model.
 665
 666     Parameters
 667     ----------
 668     x_data : array_like
 669         x values for which the WLC tension should be calculated.
 670     T : float
 671         temperature in Kelvin.
 672     L : float
 673         contour length in meters.
 674     p : float
 675         persistence length in meters.
 676
 677     Returns
 678     -------
 679     f_data : array
 680         `F(x_data)`.
 681
 682     Examples
 683     --------
 684     >>> WLC_fn(numpy.array([1e-9, 5e-9, 10e-9], dtype=numpy.float),
 685     ...        T=300, L=15e-9, p=2.5e-10) # doctest: +ELLIPSIS
 686     array([  1.717...e-12,   1.070...e-11,   4.418...e-11])
 687
 688     Notes
 689     -----
 690     The function is the simple polynomial worm like chain
 691     interpolation forumula proposed by C. Bustamante, et al. [#]_.
 692
 693     .. math::
 694       F(x) = \frac{k_B T}{p} \left[
 695         \frac{1}{4}\left( \frac{1}{(1-x/L)^2} - 1 \right)
 696         + \frac{x}{L}
 697                              \right]
 698
 699     .. [#] C. Bustamante, J.F. Marko, E.D. Siggia, and S.Smith.
 700     "Entropic elasticity of lambda-phage DNA."
 701     Science. 1994.
 702     doi: `10.1126/science.8079175 <http://dx.doi.org/10.1126/science.8079175>`
 703     """
 704     a = kB * T / p
 705     scaled_data = x_data / L
 706     return a * (0.25*((1-scaled_data)**-2 - 1) + scaled_data)
 707
 708 class WLC (ModelFitter):
 709     """Fit the data to a worm like chain.
 710
 711     Examples
 712     --------
 713     Generate some example data
 714
 715     >>> T = 300  # Kelvin
 716     >>> L = 35e-9  # meters
 717     >>> p = 2.5e-10  # meters
 718     >>> x_data = numpy.linspace(10e-9, 30e-9, num=20)
 719     >>> d_data = WLC_fn(x_data, T=T, L=L, p=p)
 720
 721     Fit the example data with a two-parameter fit (`L` and `p`).
 722
 723     >>> info = {
 724     ...     'temperature (K)': T,
 725     ...     'x data (m)': x_data,
 726     ...     }
 727     >>> model = WLC(d_data, info=info, rescale=True)
 728     >>> outqueue = Queue()
 729     >>> Lp,p = model.fit(outqueue=outqueue)
 730     >>> fit_info = outqueue.get(block=False)
 731     >>> model.L(Lp)  # doctest: +ELLIPSIS
 732     3.500...e-08
 733     >>> p  # doctest: +ELLIPSIS
 734     2.500...e-10
 735
 736     Fit the example data with a one-parameter fit (`L`).  We introduce
 737     some error in our fixed persistence length for fun.
 738
 739     >>> info['persistence length (m)'] = 2*p
 740     >>> model = WLC(d_data, info=info, rescale=True)
 741     >>> Lp = model.fit(outqueue=outqueue)
 742     >>> fit_info = outqueue.get(block=False)
 743     >>> model.L(Lp)  # doctest: +ELLIPSIS
 744     3.318...e-08
 745     """
 746     def Lp(self, L):
 747         """To avoid invalid values of `L`, we reparametrize with `Lp`.
 748
 749         Parameters
 750         ----------
 751         L : float
 752             contour length in meters
 753
 754         Returns
 755         -------
 756         Lp : float
 757             rescaled version of `L`.
 758
 759         Examples
 760         --------
 761         >>> x_data = numpy.linspace(1e-9, 20e-9, num=10)
 762         >>> model = WLC(data=x_data, info={'x data (m)':x_data})
 763         >>> model.Lp(20e-9)  # doctest: +ELLIPSIS
 764         -inf
 765         >>> model.Lp(19e-9)  # doctest: +ELLIPSIS
 766         nan
 767         >>> model.Lp(21e-9)  # doctest: +ELLIPSIS
 768         -2.99...
 769         >>> model.Lp(100e-9)  # doctest: +ELLIPSIS
 770         1.386...
 771
 772         Notes
 773         -----
 774         The rescaling is designed to limit `L` to values strictly
 775         greater than the maximum `x` data value, so we use
 776
 777         .. math::
 778             L = [\exp(L_p))+1] * x_\text{max}
 779
 780         which has the inverse
 781
 782         .. math::
 783             L_p = \ln(L/x_\text{max}-1)
 784
 785         This will obviously effect the interpretation of the fit's covariance matrix.
 786         """
 787         x_max = self.info['x data (m)'].max()
 788         return numpy.log(L/x_max - 1)
 789
 790     def L(self, Lp):
 791         """Inverse of :meth:`Lp`.
 792
 793         Parameters
 794         ----------
 795         Lp : float
 796             rescaled version of `L`
 797
 798         Returns
 799         -------
 800         L : float
 801             contour length in meters
 802
 803         Examples
 804         --------
 805         >>> x_data = numpy.linspace(1e-9, 20e-9, num=10)
 806         >>> model = WLC(data=x_data, info={'x data (m)':x_data})
 807         >>> model.L(model.Lp(21e-9))  # doctest: +ELLIPSIS
 808         2.100...e-08
 809         >>> model.L(model.Lp(100e-9))  # doctest: +ELLIPSIS
 810         9.999...e-08
 811         """
 812         x_max = self.info['x data (m)'].max()
 813         return (numpy.exp(Lp)+1)*x_max
 814
 815     def model(self, params):
 816         """Extract the relavant arguments and use :func:`WLC_fn`.
 817
 818         Parameters
 819         ----------
 820         params : list of floats
 821             `[Lp, p]`, where the presence of `p` is optional.
 822         """
 823         # setup convenient aliases
 824         T = self.info['temperature (K)']
 825         x_data = self.info['x data (m)']
 826         L = self.L(params[0])
 827         if len(params) > 1:
 828             p = abs(params[1])
 829         else:
 830             p = self.info['persistence length (m)']
 831         # compute model data
 832         self._model_data[:] = WLC_fn(x_data, T, L, p)
 833         return self._model_data
 834
 835     def guess_initial_params(self, outqueue=None):
 836         """Guess initial fitting parameters.
 837
 838         Returns
 839         -------
 840         Lp : float
 841             A guess at the reparameterized contour length in meters.
 842         p : float (optional)
 843             A guess at the persistence length in meters.  If `.info`
 844             has a setting for `'persistence length (m)'`, `p` is not
 845             returned.
 846         """
 847         x_max = self.info['x data (m)'].max()
 848         p_given = 'persistence length (m)' in self.info
 849         if p_given:
 850             p = self.info['persistence length (m)']
 851         else:
 852             p = x_max / 10.0
 853         L = 1.5 * x_max
 854         Lp = self.Lp(L)
 855         if p_given:
 856             return [Lp]
 857         return [Lp, p]
 858
 859     def guess_scale(self, params, outqueue=None):
 860         """Guess parameter scales.
 861
 862         Returns
 863         -------
 864         Lp_scale : float
 865             A guess at the reparameterized contour length scale in meters.
 866         p_scale : float (optional)
 867             A guess at the persistence length in meters.  If the
 868             length of `params` is less than 2, `p_scale` is not
 869             returned.
 870         """
 871         Lp_scale = 1.0
 872         if len(params) == 1:
 873             return [Lp_scale]
 874         return [Lp_scale] + [x/10.0 for x in params[1:]]
 875
 876
 877 class PolymerFitPlugin (Plugin):
 878     """Polymer model (WLC, FJC, etc.) fitting.
 879     """
 880     def __init__(self):
 881         super(PolymerFitPlugin, self).__init__(name='polymer_fit')
 882         self._arguments = [  # For Command initialization
 883             Argument('polymer model', default='WLC', help="""
 884 Select the default polymer model for 'polymer fit'.  See the
 885 documentation for descriptions of available algorithms.
 886 """.strip()),
 887             Argument('FJC Kuhn length', type='float', default=4e-10,
 888                     help='Kuhn length in meters'),
 889             Argument('FJC_PEG Kuhn length', type='float', default=4e-10,
 890                     help='Kuhn length in meters'),
 891             Argument('FJC_PEG elasticity', type='float', default=150.0,
 892                     help='Elasticity of a PEG segment in Newtons per meter.'),
 893             Argument('FJC_PEG delta G', type='float', default=3.0, help="""
 894 Gibbs free energy difference between trans-trans-trans (ttt) and
 895 trans-trans-gauche (ttg) PEG states in units of kBT.
 896 """.strip()),
 897             Argument('FJC_PEG L_helical', type='float', default=2.8e-10,
 898                     help='Contour length of PEG in the ttg state.'),
 899             Argument('FJC_PEG L_planar', type='float', default=3.58e-10,
 900                     help='Contour length of PEG in the ttt state.'),
 901             Argument('WLC persistence length', type='float', default=4e-10,
 902                     help='Persistence length in meters'),
 903             ]
 904         self._settings = [
 905             Setting(section=self.setting_section, help=self.__doc__)]
 906         for argument in self._arguments:
 907             self._settings.append(argument_to_setting(
 908                     self.setting_section, argument))
 909             argument.default = None # if argument isn't given, use the config.
 910         self._arguments.extend([  # Non-configurable arguments
 911                 Argument(name='input distance column', type='string',
 912                          default='cantilever adjusted extension (m)',
 913                          help="""
 914 Name of the column to use as the surface position input.
 915 """.strip()),
 916                 Argument(name='input deflection column', type='string',
 917                          default='deflection (N)',
 918                          help="""
 919 Name of the column to use as the deflection input.
 920 """.strip()),
 921                 ])
 922         self._commands = [
 923             PolymerFitCommand(self), PolymerFitPeaksCommand(self),
 924             TranslateFlatPeaksCommand(self),
 925             ]
 926
 927     def dependencies(self):
 928         return ['vclamp']
 929
 930     def default_settings(self):
 931         return self._settings
 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                 ] + plugin._arguments + [
 961                 Argument(name='output tension column', type='string',
 962                          default='polymer tension',
 963                          help="""
 964 Name of the column (without units) to use as the polymer tension output.
 965 """.strip()),
 966                 Argument(name='fit parameters info name', type='string',
 967                          default='surface deflection offset',
 968                          help="""
 969 Name (without units) for storing the fit parameters in the `.info` dictionary.
 970 """.strip()),
 971                 ],
 972             help=self.__doc__, plugin=plugin)
 973
 974     def _run(self, hooke, inqueue, outqueue, params):
 975         for key,value in params.items():
 976             if value == None:  # Use configured default value.
 977                 params[key] = self.plugin.config[key]
 978         scale(hooke, params['curve'], params['block'])  # TODO: is autoscaling a good idea? (explicit is better than implicit)
 979         data = params['curve'].data[params['block']]
 980         # HACK? rely on params['curve'] being bound to the local hooke
 981         # playlist (i.e. not a copy, as you would get by passing a
 982         # curve through the queue).  Ugh.  Stupid queues.  As an
 983         # alternative, we could pass lookup information through the
 984         # queue...
 985         model = params['polymer model']
 986         new = Data((data.shape[0], data.shape[1]+1), dtype=data.dtype)
 987         new.info = copy.deepcopy(data.info)
 988         new[:,:-1] = data
 989         new.info['columns'].append(
 990             join_data_label(params['output tension column'], 'N'))
 991         z_data = data[:,data.info['columns'].index(
 992                 params['input distance column'])]
 993         d_data = data[:,data.info['columns'].index(
 994                 params['input deflection column'])]
 995         start,stop = params['bounds']
 996         tension_data,ps = self.fit_polymer_model(
 997             params, z_data, d_data, start, stop, outqueue)
 998         new.info[params['fit parameters info name']] = ps
 999         new.info[params['fit parameters info name']]['model'] = model
1000         new[:,-1] = tension_data
1001         params['curve'].data[params['block']] = new
1002
1003     def fit_polymer_model(self, params, z_data, d_data, start, stop,
1004                           outqueue=None):
1005         """Railyard for the `fit_*_model` family.
1006
1007         Uses the `polymer model` configuration setting to call the
1008         appropriate backend algorithm.
1009         """
1010         fn = getattr(self, 'fit_%s_model'
1011                      % params['polymer model'].replace('-','_'))
1012         return fn(params, z_data, d_data, start, stop, outqueue)
1013
1014     def fit_FJC_model(self, params, z_data, d_data, start, stop,
1015                       outqueue=None):
1016         """Fit the data with :class:`FJC`.
1017         """
1018         info = {
1019             'temperature (K)': self.plugin.config['temperature'],
1020             'x data (m)': z_data[start:stop],
1021             }
1022         if True:  # TODO: optionally free persistence length
1023             info['Kuhn length (m)'] = (
1024                 params['FJC Kuhn length'])
1025         model = FJC(d_data[start:stop], info=info, rescale=True)
1026         queue = Queue()
1027         params = model.fit(outqueue=queue)
1028         if True:  # TODO: if Kuhn length fixed
1029             params = [params]
1030             a = info['Kuhn length (m)']
1031         else:
1032             a = params[1]
1033         Lp = params[0]
1034         L = model.L(Lp)
1035         T = info['temperature (K)']
1036         fit_info = queue.get(block=False)
1037         f_data = numpy.ones(z_data.shape, dtype=z_data.dtype) * numpy.nan
1038         f_data[start:stop] = FJC_fn(z_data[start:stop], T=T, L=L, a=a)
1039         return [f_data, fit_info]
1040
1041     def fit_FJC_PEG_model(self, params, z_data, d_data, start, stop,
1042                           outqueue=None):
1043         """Fit the data with :class:`FJC_PEG`.
1044         """
1045         info = {
1046             'temperature (K)': self.plugin.config['temperature'],
1047             'x data (m)': z_data[start:stop],
1048             # TODO: more info
1049             }
1050         if True:  # TODO: optionally free persistence length
1051             info['Kuhn length (m)'] = (
1052                 params['FJC Kuhn length'])
1053         model = FJC_PEG(d_data[start:stop], info=info, rescale=True)
1054         queue = Queue()
1055         params = model.fit(outqueue=queue)
1056         if True:  # TODO: if Kuhn length fixed
1057             params = [params]
1058             a = info['Kuhn length (m)']
1059         else:
1060             a = params[1]
1061         Nr = params[0]
1062         N = model.L(Nr)
1063         T = info['temperature (K)']
1064         fit_info = queue.get(block=False)
1065         f_data = numpy.ones(z_data.shape, dtype=z_data.dtype) * numpy.nan
1066         f_data[start:stop] = FJC_PEG_fn(z_data[start:stop], **kwargs)
1067         return [f_data, fit_info]
1068
1069     def fit_WLC_model(self, params, z_data, d_data, start, stop,
1070                       outqueue=None):
1071         """Fit the data with :class:`WLC`.
1072         """
1073         info = {
1074             'temperature (K)': self.plugin.config['temperature'],
1075             'x data (m)': z_data[start:stop],
1076             }
1077         if True:  # TODO: optionally free persistence length
1078             info['persistence length (m)'] = (
1079                 params['WLC persistence length'])
1080         model = WLC(d_data[start:stop], info=info, rescale=True)
1081         queue = Queue()
1082         params = model.fit(outqueue=queue)
1083         if True:  # TODO: if persistence length fixed
1084             params = [params]
1085             p = info['persistence length (m)']
1086         else:
1087             p = params[1]
1088         Lp = params[0]
1089         L = model.L(Lp)
1090         T = info['temperature (K)']
1091         fit_info = queue.get(block=False)
1092         f_data = numpy.ones(z_data.shape, dtype=z_data.dtype) * numpy.nan
1093         f_data[start:stop] = WLC_fn(z_data[start:stop], T=T, L=L, p=p)
1094         return [f_data, fit_info]
1095
1096
1097 class PolymerFitPeaksCommand (Command):
1098     """Polymer model (WLC, FJC, etc.) fitting.
1099
1100     Use :class:`PolymerFitCommand` to fit the each peak in a list of
1101     previously determined peaks.
1102     """
1103     def __init__(self, plugin):
1104         super(PolymerFitPeaksCommand, self).__init__(
1105             name='polymer fit peaks',
1106             arguments=[
1107                 CurveArgument,
1108                 Argument(name='block', aliases=['set'], type='int', default=0,
1109                          help="""
1110 Data block for which the fit should be calculated.  For an
1111 approach/retract force curve, `0` selects the approaching curve and
1112 `1` selects the retracting curve.
1113 """.strip()),
1114                 Argument(name='peak info name', type='string',
1115                          default='polymer peaks',
1116                          help="""
1117 Name for the list of peaks in the `.info` dictionary.
1118 """.strip()),
1119                 Argument(name='peak index', type='int', count=-1, default=None,
1120                          help="""
1121 Index of the selected peak in the list of peaks.  Use `None` to fit all peaks.
1122 """.strip()),
1123                 ] + plugin._arguments,
1124             help=self.__doc__, plugin=plugin)
1125
1126     def _run(self, hooke, inqueue, outqueue, params):
1127         data = params['curve'].data[params['block']]
1128         fit_command = [c for c in hooke.commands
1129                        if c.name=='polymer fit'][0]
1130         inq = Queue()
1131         outq = Queue()
1132         p = copy.deepcopy(params)
1133         p['curve'] = params['curve']
1134         del(p['peak info name'])
1135         del(p['peak index'])
1136         for i,peak in enumerate(data.info[params['peak info name']]):
1137             if params['peak index'] == None or i in params['peak index']:
1138                 p['bounds'] = [peak.index, peak.post_index()]
1139                 p['output tension column'] = peak.name
1140                 p['fit parameters info name'] = peak.name
1141                 fit_command.run(hooke, inq, outq, **p)
1142             ret = outq.get()
1143             if not isinstance(ret, Success):
1144                 raise ret
1145
1146
1147 class TranslateFlatPeaksCommand (Command):
1148     """Translate flat filter peaks into polymer peaks for fitting.
1149
1150     Use :class:`~hooke.plugin.flatfilt.FlatPeaksCommand` creates a
1151     list of peaks for regions with large derivatives.  For velocity
1152     clamp measurements, these regions are usually the rebound phase
1153     after a protein domain unfolds, the cantilever detaches, etc.
1154     Because these features occur after the polymer loading phase, we
1155     need to shift the selected regions back to align them with the
1156     polymer loading regions.
1157     """
1158     def __init__(self, plugin):
1159         plugin_arguments = [a for a in plugin._arguments
1160                             if a.name in ['input distance column',
1161                                           'input deflection column']]
1162         arguments = [
1163             CurveArgument,
1164             Argument(name='block', aliases=['set'], type='int', default=0,
1165                      help="""
1166 Data block for which the fit should be calculated.  For an
1167 approach/retract force curve, `0` selects the approaching curve and
1168 `1` selects the retracting curve.
1169 """.strip()),
1170             ] + plugin_arguments + [
1171             Argument(name='input peak info name', type='string',
1172                      default='flat filter peaks',
1173                      help="""
1174 Name for the input peaks in the `.info` dictionary.
1175 """.strip()),
1176             Argument(name='output peak info name', type='string',
1177                      default='polymer peaks',
1178                      help="""
1179 Name for the ouptput peaks in the `.info` dictionary.
1180 """.strip()),
1181             Argument(name='end offset', type='int', default=-1,
1182                      help="""
1183 Number of points between the end of a new peak and the start of the old.
1184 """.strip()),
1185             Argument(name='start fraction', type='float', default=0.2,
1186                      help="""
1187 Place the start of the new peak at `start_fraction` from the end of
1188 the previous old peak to the end of the new peak.  Because the first
1189 new peak will have no previous old peak, it uses a distance of zero
1190 instead.
1191 """.strip()),
1192             ]
1193         super(TranslateFlatPeaksCommand, self).__init__(
1194             name='flat peaks to polymer peaks',
1195             arguments=arguments,
1196             help=self.__doc__, plugin=plugin)
1197
1198     def _run(self, hooke, inqueue, outqueue, params):
1199         data = params['curve'].data[params['block']]
1200         z_data = data[:,data.info['columns'].index(
1201                 params['input distance column'])]
1202         d_data = data[:,data.info['columns'].index(
1203                 params['input deflection column'])]
1204         previous_old_stop = numpy.absolute(z_data).argmin()
1205         new = []
1206         for i,peak in enumerate(data.info[params['input peak info name']]):
1207             next_old_start = peak.index
1208             stop = next_old_start + params['end offset']
1209             z_start = z_data[previous_old_stop] + params['start fraction']*(
1210                 z_data[stop] - z_data[previous_old_stop])
1211             start = numpy.absolute(z_data - z_start).argmin()
1212             p = Peak('polymer peak %d' % i,
1213                      index=start,
1214                      values=d_data[start:stop])
1215             new.append(p)
1216             previous_old_stop = peak.post_index()
1217         data.info[params['output peak info name']] = new
1218
1219
1220 # TODO:
1221 # def dist2fit(self):
1222 #     '''Calculates the average distance from data to fit, scaled by
1223 #     the standard deviation of the free cantilever area (thermal
1224 #     noise).
1225 #     '''