hooke/plugin/polymer_fit.py

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