hooke/plugin/polymer_fit.py

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