1 # -*- coding: utf-8 -*-
3 # Copyright (C) 2008-2010 Alberto Gomez-Casado
5 # Massimo Sandal <devicerandom@gmail.com>
6 # W. Trevor King <wking@drexel.edu>
8 # This file is part of Hooke.
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.
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.
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/>.
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.).
30 from Queue import Queue
35 from scipy.optimize import newton
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
48 kB = 1.3806503e-23 # Joules/Kelvin
52 """Hyperbolic cotangent.
56 >>> coth(1.19967874) # doctest: +ELLIPSIS
64 \coth(z) \equiv \frac{\exp(z) + \exp(-z)}{\exp(z) - \exp(-z)}
68 http://mathworld.wolfram.com/HyperbolicCotangent.html
70 return 1.0/numpy.tanh(z)
73 """Inverse hyperbolic cotangent.
77 >>> arccoth(1.19967874) # doctest: +ELLIPSIS
79 >>> arccoth(coth(numpy.pi)) # doctest: +ELLIPSIS
84 Inverting from the :func:`definition <coth>`.
87 z \equiv \atanh\left( \frac{1}{\coth(z)} \right)
89 return numpy.arctanh(1.0/z)
96 >>> Langevin(numpy.pi) # doctest: +ELLIPSIS
101 From `Wikipedia`_ or Hatfield [#]_
104 L(z) \equiv \coth(z) - \frac{1}{z}
107 http://en.wikipedia.org/wiki/Brillouin_and_Langevin_functions#Langevin_Function
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>`
114 return coth(z) - 1.0/z
116 def inverse_Langevin(z, extreme=1.0 - 1e-8):
117 """Inverse Langevin function.
121 z : float or array_like
122 object whose inverse Langevin will be returned
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.
130 >>> inverse_Langevin(Langevin(numpy.pi)) # doctest: +ELLIPSIS
132 >>> inverse_Langevin(Langevin(numpy.array([1,2,3]))) # doctest: +ELLIPSIS
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`_).
142 L^{-1}(z) = 3z + \frac{9}{5}z^3 + \frac{297}{175}z^5 + \dots
145 http://en.wikipedia.org/wiki/Brillouin_and_Langevin_functions#Langevin_Function
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)
156 # Catch stdout since sometimes the newton function prints exit
157 # messages to stdout, e.g. "Tolerance of %s reached".
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,)
164 output = tmp_stdout.getvalue()
167 def FJC_fn(x_data, T, L, a):
168 """The freely jointed chain model.
173 x values for which the WLC tension should be calculated.
175 temperature in Kelvin.
177 contour length in meters.
179 Kuhn length in meters.
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])
194 The freely jointed chain model is
197 F(x) = \frac{k_B T}{a} L^{-1}\left( \frac{x}{L} \right)
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 [#]_.
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>`
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>`
212 return kB*T / a * inverse_Langevin(x_data/L)
214 class FJC (ModelFitter):
215 """Fit the data to a freely jointed chain.
219 Generate some example data
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)
227 Fit the example data with a two-parameter fit (`L` and `a`).
230 ... 'temperature (K)': T,
231 ... 'x data (m)': x_data,
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)
242 Fit the example data with a one-parameter fit (`L`). We introduce
243 some error in our fixed Kuhn length for fun.
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
253 """To avoid invalid values of `L`, we reparametrize with `Lp`.
258 contour length in meters
263 rescaled version of `L`.
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
271 >>> model.Lp(19e-9) # doctest: +ELLIPSIS
273 >>> model.Lp(21e-9) # doctest: +ELLIPSIS
275 >>> model.Lp(100e-9) # doctest: +ELLIPSIS
280 The rescaling is designed to limit `L` to values strictly
281 greater than the maximum `x` data value, so we use
284 L = [\exp(L_p))+1] * x_\text{max}
286 which has the inverse
289 L_p = \ln(L/x_\text{max}-1)
291 This will obviously effect the interpretation of the fit's covariance matrix.
293 x_max = self.info['x data (m)'].max()
294 return numpy.log(L/x_max - 1)
297 """Inverse of :meth:`Lp`.
302 rescaled version of `L`
307 contour length in meters
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
315 >>> model.L(model.Lp(100e-9)) # doctest: +ELLIPSIS
318 x_max = self.info['x data (m)'].max()
319 return (numpy.exp(Lp)+1)*x_max
321 def model(self, params):
322 """Extract the relavant arguments and use :func:`FJC_fn`.
326 params : list of floats
327 `[Lp, a]`, where the presence of `a` is optional.
329 # setup convenient aliases
330 T = self.info['temperature (K)']
331 x_data = self.info['x data (m)']
332 L = self.L(params[0])
336 a = self.info['Kuhn length (m)']
338 self._model_data[:] = FJC_fn(x_data, T, L, a)
339 return self._model_data
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
350 def guess_initial_params(self, outqueue=None):
351 """Guess initial fitting parameters.
356 A guess at the reparameterized contour length in meters.
358 A guess at the Kuhn length in meters. If `.info` has a
359 setting for `'Kuhn length (m)'`, `a` is not returned.
361 x_max = self.info['x data (m)'].max()
362 a_given = 'Kuhn length (m)' in self.info
364 a = self.info['Kuhn length (m)']
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.
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
386 The function is that proposed by F. Oesterhelt, et al. [#]_.
389 x(F) = N_\text{S} \cdot \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}
396 \coth\left(\frac{Fa}{k_B T}\right)
399 + \frac{F}{K_\text{S}}
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.
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>`
418 g = (dG - F_data*(Lp-Lh)) / kBT
420 return N * ((Lp/(numpy.exp(-g)+1) + Lh/(numpy.exp(+g)+1)) * (coth(z)-1./z)
423 def FJC_PEG_fn(x_data, T, N, k, Lp, Lh, dG, a):
424 """Poly(ethylene-glycol) adjusted extended FJC model.
428 z : float or array_like
429 object whose inverse Langevin will be returned
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.
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
441 >>> FJC_PEG_fn(inverse_FJC_PEG_fn(3.4e-10, **kwargs), **kwargs) # doctest: +ELLIPSIS
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])
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.
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)
463 # Approximate with the simple FJC_fn()
466 while guess == numpy.inf:
467 guess = FJC_fn(x_data, T=T, L=L, a=a)
470 fn = lambda f: inverse_FJC_PEG_fn(guess*abs(f), **kwargs) - x_data
471 ret = guess*abs(newton(func=fn, x0=1.0))
474 class FJC_PEG (ModelFitter):
475 """Fit the data to an poly(ethylene-glycol) adjusted extended FJC.
479 Generate some example data
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)
486 Fit the example data with a two-parameter fit (`N` and `a`).
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'],
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)
505 Fit the example data with a one-parameter fit (`N`). We introduce
506 some error in our fixed Kuhn length for fun.
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
516 """To avoid invalid values of `L`, we reparametrize with `Lr`.
521 contour length in meters
526 rescaled version of `L`.
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
534 >>> model.Lr(19e-9) # doctest: +ELLIPSIS
536 >>> model.Lr(21e-9) # doctest: +ELLIPSIS
538 >>> model.Lr(100e-9) # doctest: +ELLIPSIS
543 The rescaling is designed to limit `L` to values strictly
544 greater than zero, so we use
547 L = \exp(L_p) * x_\text{max}
549 which has the inverse
552 L_p = \ln(L/x_\text{max})
554 This will obviously effect the interpretation of the fit's covariance matrix.
556 x_max = self.info['x data (m)'].max()
557 return numpy.log(L/x_max)
560 """Inverse of :meth:`Lr`.
565 rescaled version of `L`
570 contour length in meters
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
578 >>> model.L(model.Lr(100e-9)) # doctest: +ELLIPSIS
581 x_max = self.info['x data (m)'].max()
582 return numpy.exp(Lr)*x_max
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)'],
593 def model(self, params):
594 """Extract the relavant arguments and use :func:`FJC_PEG_fn`.
598 params : list of floats
599 `[N, a]`, where the presence of `a` is optional.
601 # setup convenient aliases
602 T = self.info['temperature (K)']
603 x_data = self.info['x data (m)']
604 N = self.L(params[0])
608 a = self.info['Kuhn length (m)']
609 kwargs = self._kwargs()
611 self._model_data[:] = FJC_PEG_fn(x_data, N=N, a=a, **kwargs)
612 return self._model_data
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
623 def guess_initial_params(self, outqueue=None):
624 """Guess initial fitting parameters.
629 A guess at the section count.
631 A guess at the Kuhn length in meters. If `.info` has a
632 setting for `'Kuhn length (m)'`, `a` is not returned.
634 x_max = self.info['x data (m)'].max()
635 a_given = 'Kuhn length (m)' in self.info
637 a = self.info['Kuhn length (m)']
640 f_max = self._data.max()
641 kwargs = self._kwargs()
643 x_section = inverse_FJC_PEG_fn(F_data=f_max, **kwargs)
644 N = x_max / x_section;
651 def WLC_fn(x_data, T, L, p):
652 """The worm like chain interpolation model.
657 x values for which the WLC tension should be calculated.
659 temperature in Kelvin.
661 contour length in meters.
663 persistence length in meters.
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])
678 The function is the simple polynomial worm like chain
679 interpolation forumula proposed by C. Bustamante, et al. [#]_.
682 F(x) = \frac{k_B T}{p} \left[
683 \frac{1}{4}\left( \frac{1}{(1-x/L)^2} - 1 \right)
687 .. [#] C. Bustamante, J.F. Marko, E.D. Siggia, and S.Smith.
688 "Entropic elasticity of lambda-phage DNA."
690 doi: `10.1126/science.8079175 <http://dx.doi.org/10.1126/science.8079175>`
693 scaled_data = x_data / L
694 return a * (0.25*((1-scaled_data)**-2 - 1) + scaled_data)
696 class WLC (ModelFitter):
697 """Fit the data to a worm like chain.
701 Generate some example data
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)
709 Fit the example data with a two-parameter fit (`L` and `p`).
712 ... 'temperature (K)': T,
713 ... 'x data (m)': x_data,
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)
724 Fit the example data with a one-parameter fit (`L`). We introduce
725 some error in our fixed persistence length for fun.
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
735 """To avoid invalid values of `L`, we reparametrize with `Lp`.
740 contour length in meters
745 rescaled version of `L`.
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
753 >>> model.Lp(19e-9) # doctest: +ELLIPSIS
755 >>> model.Lp(21e-9) # doctest: +ELLIPSIS
757 >>> model.Lp(100e-9) # doctest: +ELLIPSIS
762 The rescaling is designed to limit `L` to values strictly
763 greater than the maximum `x` data value, so we use
766 L = [\exp(L_p))+1] * x_\text{max}
768 which has the inverse
771 L_p = \ln(L/x_\text{max}-1)
773 This will obviously effect the interpretation of the fit's covariance matrix.
775 x_max = self.info['x data (m)'].max()
776 return numpy.log(L/x_max - 1)
779 """Inverse of :meth:`Lp`.
784 rescaled version of `L`
789 contour length in meters
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
797 >>> model.L(model.Lp(100e-9)) # doctest: +ELLIPSIS
800 x_max = self.info['x data (m)'].max()
801 return (numpy.exp(Lp)+1)*x_max
803 def model(self, params):
804 """Extract the relavant arguments and use :func:`WLC_fn`.
808 params : list of floats
809 `[Lp, p]`, where the presence of `p` is optional.
811 # setup convenient aliases
812 T = self.info['temperature (K)']
813 x_data = self.info['x data (m)']
814 L = self.L(params[0])
818 p = self.info['persistence length (m)']
820 self._model_data[:] = WLC_fn(x_data, T, L, p)
821 return self._model_data
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
832 def guess_initial_params(self, outqueue=None):
833 """Guess initial fitting parameters.
838 A guess at the reparameterized contour length in meters.
840 A guess at the persistence length in meters. If `.info`
841 has a setting for `'persistence length (m)'`, `p` is not
844 x_max = self.info['x data (m)'].max()
845 p_given = 'persistence length (m)' in self.info
847 p = self.info['persistence length (m)']
857 class PolymerFitPlugin (Plugin):
858 """Polymer model (WLC, FJC, etc.) fitting.
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.
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.
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'),
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)',
893 Name of the column to use as the surface position input.
895 ('deflection column', 'deflection (N)',
897 Name of the column to use as the deflection input.
901 PolymerFitCommand(self), PolymerFitPeaksCommand(self),
902 TranslateFlatPeaksCommand(self),
905 def dependencies(self):
908 def default_settings(self):
909 return self._settings
912 class PolymerFitCommand (ColumnAddingCommand):
913 """Polymer model (WLC, FJC, etc.) fitting.
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).
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()
923 def __init__(self, plugin):
924 super(PolymerFitCommand, self).__init__(
926 columns=plugin._input_columns,
928 ('output tension column', 'polymer tension',
930 Name of the column (without units) to use as the polymer tension output.
934 Argument(name='fit parameters info name', type='string',
935 default='polymer fit',
937 Name (without units) for storing the fit parameters in the `.info` dictionary.
939 Argument(name='bounds', type='point', optional=False, count=2,
941 Indicies of points bounding the selected data.
943 ] + plugin._arguments,
944 help=self.__doc__, plugin=plugin)
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',
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)
972 def fit_polymer_model(self, params, dist_data, def_data, start, stop,
974 """Railyard for the `fit_*_model` family.
976 Uses the `polymer model` configuration setting to call the
977 appropriate backend algorithm.
979 fn = getattr(self, 'fit_%s_model'
980 % params['polymer model'].replace('-','_'))
981 return fn(params, dist_data, def_data, start, stop, outqueue)
983 def fit_FJC_model(self, params, dist_data, def_data, start, stop,
985 """Fit the data with :class:`FJC`.
988 'temperature (K)': self.plugin.config['temperature'],
989 'x data (m)': dist_data[start:stop],
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)
996 params = model.fit(outqueue=queue)
997 if True: # TODO: if Kuhn length fixed
999 a = info['Kuhn length (m)']
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]
1009 def fit_FJC_PEG_model(self, params, dist_data, def_data, start, stop,
1011 """Fit the data with :class:`FJC_PEG`.
1014 'temperature (K)': self.plugin.config['temperature'],
1015 'x data (m)': dist_data[start:stop],
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)
1023 params = model.fit(outqueue=queue)
1024 if True: # TODO: if Kuhn length fixed
1026 a = info['Kuhn length (m)']
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]
1036 def fit_WLC_model(self, params, dist_data, def_data, start, stop,
1038 """Fit the data with :class:`WLC`.
1041 'temperature (K)': self.plugin.config['temperature'],
1042 'x data (m)': dist_data[start:stop],
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)
1049 params = model.fit(outqueue=queue)
1050 if True: # TODO: if persistence length fixed
1052 p = info['persistence length (m)']
1055 info['persistence length (m)'] = p
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]
1067 class PolymerFitPeaksCommand (ColumnAccessCommand):
1068 """Polymer model (WLC, FJC, etc.) fitting.
1070 Use :class:`PolymerFitCommand` to fit the each peak in a list of
1071 previously determined peaks.
1073 def __init__(self, plugin):
1074 super(PolymerFitPeaksCommand, self).__init__(
1075 name='polymer fit peaks',
1076 columns=plugin._input_columns,
1078 Argument(name='peak info name', type='string',
1079 default='polymer peaks',
1081 Name for the list of peaks in the `.info` dictionary.
1083 Argument(name='peak index', type='int', count=-1, default=None,
1085 Index of the selected peak in the list of peaks. Use `None` to fit all peaks.
1087 ] + plugin._arguments,
1088 help=self.__doc__, plugin=plugin)
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]
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)
1109 if not isinstance(ret, Success):
1113 class TranslateFlatPeaksCommand (ColumnAccessCommand):
1114 """Translate flat filter peaks into polymer peaks for fitting.
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.
1124 def __init__(self, plugin):
1125 super(TranslateFlatPeaksCommand, self).__init__(
1126 name='flat peaks to polymer peaks',
1127 columns=plugin._input_columns,
1129 Argument(name='input peak info name', type='string',
1130 default='flat filter peaks',
1132 Name for the input peaks in the `.info` dictionary.
1134 Argument(name='output peak info name', type='string',
1135 default='polymer peaks',
1137 Name for the ouptput peaks in the `.info` dictionary.
1139 Argument(name='end offset', type='int', default=-1,
1141 Number of points between the end of a new peak and the start of the old.
1143 Argument(name='start fraction', type='float', default=0.2,
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
1150 ] + plugin._arguments,
1151 help=self.__doc__, plugin=plugin)
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()
1162 peaks = data.info[params['input peak info name']]
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']
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,
1177 values=def_data[start:stop])
1179 previous_old_stop = peak.post_index()
1180 data.info[params['output peak info name']] = new
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