hooke/plugin/polymer_fit.py

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