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