1 # -*- coding: utf-8 -*-
3 # Copyright (C) 2010-2012 W. Trevor King <wking@drexel.edu>
5 # This file is part of Hooke.
7 # Hooke is free software: you can redistribute it and/or modify it
8 # under the terms of the GNU Lesser General Public License as
9 # published by the Free Software Foundation, either version 3 of the
10 # License, or (at your option) any later version.
12 # Hooke is distributed in the hope that it will be useful, but WITHOUT
13 # ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
14 # or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General
15 # Public License for more details.
17 # You should have received a copy of the GNU Lesser General Public
18 # License along with Hooke. If not, see
19 # <http://www.gnu.org/licenses/>.
21 """The ``polymer_fit`` module proviews :class:`PolymerFitPlugin` and
22 serveral associated :class:`hooke.command.Command`\s for Fitting
23 velocity-clamp data to various polymer models (WLC, FJC, etc.).
27 from Queue import Queue
32 from scipy.optimize import newton
34 from ..command import Command, Argument, Success, Failure
35 from ..config import Setting
36 from ..curve import Data
37 from ..util.callback import is_iterable
38 from ..util.fit import PoorFit, ModelFitter
39 from ..util.peak import Peak
40 from ..util.si import join_data_label, split_data_label
41 from . import Plugin, argument_to_setting
42 from .curve import ColumnAccessCommand, ColumnAddingCommand
45 kB = 1.3806503e-23 # Joules/Kelvin
49 """Hyperbolic cotangent.
53 >>> coth(1.19967874) # doctest: +ELLIPSIS
61 \coth(z) \equiv \frac{\exp(z) + \exp(-z)}{\exp(z) - \exp(-z)}
65 http://mathworld.wolfram.com/HyperbolicCotangent.html
67 return 1.0/numpy.tanh(z)
70 """Inverse hyperbolic cotangent.
74 >>> arccoth(1.19967874) # doctest: +ELLIPSIS
76 >>> arccoth(coth(numpy.pi)) # doctest: +ELLIPSIS
81 Inverting from the :func:`definition <coth>`.
84 z \equiv \atanh\left( \frac{1}{\coth(z)} \right)
86 return numpy.arctanh(1.0/z)
93 >>> Langevin(numpy.pi) # doctest: +ELLIPSIS
98 From `Wikipedia`_ or Hatfield [#]_
101 L(z) \equiv \coth(z) - \frac{1}{z}
104 http://en.wikipedia.org/wiki/Brillouin_and_Langevin_functions#Langevin_Function
106 .. [#]: J.W. Hatfield and S.R. Quake
107 "Dynamic Properties of an Extended Polymer in Solution."
108 Phys. Rev. Lett. 1999.
109 doi: `10.1103/PhysRevLett.82.3548 <http://dx.doi.org/10.1103/PhysRevLett.82.3548>`
111 return coth(z) - 1.0/z
113 def inverse_Langevin(z, extreme=1.0 - 1e-8):
114 """Inverse Langevin function.
118 z : float or array_like
119 object whose inverse Langevin will be returned
121 value (close to one) for which we assume the inverse is +/-
122 infinity. This avoids problems with Newton-Raphson
123 convergence in regions with low slope.
127 >>> inverse_Langevin(Langevin(numpy.pi)) # doctest: +ELLIPSIS
129 >>> inverse_Langevin(Langevin(numpy.array([1,2,3]))) # doctest: +ELLIPSIS
134 We approximate the inverse Langevin via Newton's method
135 (:func:`scipy.optimize.newton`). For the starting point, we take
136 the first three terms from the Taylor series (from `Wikipedia`_).
139 L^{-1}(z) = 3z + \frac{9}{5}z^3 + \frac{297}{175}z^5 + \dots
142 http://en.wikipedia.org/wiki/Brillouin_and_Langevin_functions#Langevin_Function
145 ret = numpy.ndarray(shape=z.shape, dtype=z.dtype)
146 for i,z_i in enumerate(z):
147 ret[i] = inverse_Langevin(z_i)
153 # Catch stdout since sometimes the newton function prints exit
154 # messages to stdout, e.g. "Tolerance of %s reached".
156 tmp_stdout = StringIO.StringIO()
157 sys.stdout = tmp_stdout
158 ret = newton(func=lambda x: Langevin(x)-z,
159 x0=3*z + 9./5.*z**3 + 297./175.*z**5,)
161 output = tmp_stdout.getvalue()
164 def FJC_fn(x_data, T, L, a):
165 """The freely jointed chain model.
170 x values for which the WLC tension should be calculated.
172 temperature in Kelvin.
174 contour length in meters.
176 Kuhn length in meters.
185 >>> FJC_fn(numpy.array([1e-9, 5e-9, 10e-9], dtype=numpy.float),
186 ... T=300, L=15e-9, a=2.5e-10) # doctest: +ELLIPSIS
187 array([ 3.322...-12, 1.78...e-11, 4.889...e-11])
191 The freely jointed chain model is
194 F(x) = \frac{k_B T}{a} L^{-1}\left( \frac{x}{L} \right)
196 where :math:`L^{-1}` is the :func:`inverse_Langevin`, :math:`a` is
197 the Kuhn length, and :math:`L` is the contour length [#]_. For
198 the inverse formulation (:math:`x(F)`), see Ray and Akhremitchev [#]_.
201 .. [#]: J.W. Hatfield and S.R. Quake
202 "Dynamic Properties of an Extended Polymer in Solution."
203 Phys. Rev. Lett. 1999.
204 doi: `10.1103/PhysRevLett.82.3548 <http://dx.doi.org/10.1103/PhysRevLett.82.3548>`
206 .. [#] C. Ray and B.B. Akhremitchev.
207 `"Fitting Force Curves by the Extended Freely Jointed Chain Model" <http://www.chem.duke.edu/~boris/research/force_spectroscopy/fit_efjc.pdf>`
209 return kB*T / a * inverse_Langevin(x_data/L)
211 class FJC (ModelFitter):
212 """Fit the data to a freely jointed chain.
216 Generate some example data
219 >>> L = 35e-9 # meters
220 >>> a = 2.5e-10 # meters
221 >>> x_data = numpy.linspace(10e-9, 30e-9, num=20)
222 >>> d_data = FJC_fn(x_data, T=T, L=L, a=a)
224 Fit the example data with a two-parameter fit (`L` and `a`).
227 ... 'temperature (K)': T,
228 ... 'x data (m)': x_data,
230 >>> model = FJC(d_data, info=info, rescale=True)
231 >>> outqueue = Queue()
232 >>> L,a = model.fit(outqueue=outqueue)
233 >>> fit_info = outqueue.get(block=False)
239 Fit the example data with a one-parameter fit (`L`). We introduce
240 some error in our fixed Kuhn length for fun.
242 >>> info['Kuhn length (m)'] = 2*a
243 >>> model = FJC(d_data, info=info, rescale=True)
244 >>> L, = model.fit(outqueue=outqueue)
245 >>> fit_info = outqueue.get(block=False)
246 >>> print L # doctest: +ELLIPSIS
250 """To avoid invalid values of `L`, we reparametrize with `Lp`.
255 contour length in meters
260 rescaled version of `L`.
264 >>> x_data = numpy.linspace(1e-9, 20e-9, num=10)
265 >>> model = FJC(data=x_data, info={'x data (m)':x_data})
266 >>> model.Lp(20e-9) # doctest: +ELLIPSIS
268 >>> model.Lp(19e-9) # doctest: +ELLIPSIS
270 >>> model.Lp(21e-9) # doctest: +ELLIPSIS
272 >>> model.Lp(100e-9) # doctest: +ELLIPSIS
277 The rescaling is designed to limit `L` to values strictly
278 greater than the maximum `x` data value, so we use
281 L = [\exp(L_p))+1] * x_\text{max}
283 which has the inverse
286 L_p = \ln(L/x_\text{max}-1)
288 This will obviously effect the interpretation of the fit's covariance matrix.
290 x_max = self.info['x data (m)'].max()
291 return numpy.log(L/x_max - 1)
294 """Inverse of :meth:`Lp`.
299 rescaled version of `L`
304 contour length in meters
308 >>> x_data = numpy.linspace(1e-9, 20e-9, num=10)
309 >>> model = FJC(data=x_data, info={'x data (m)':x_data})
310 >>> model.L(model.Lp(21e-9)) # doctest: +ELLIPSIS
312 >>> model.L(model.Lp(100e-9)) # doctest: +ELLIPSIS
315 x_max = self.info['x data (m)'].max()
316 return (numpy.exp(Lp)+1)*x_max
318 def model(self, params):
319 """Extract the relavant arguments and use :func:`FJC_fn`.
323 params : list of floats
324 `[Lp, a]`, where the presence of `a` is optional.
326 # setup convenient aliases
327 T = self.info['temperature (K)']
328 x_data = self.info['x data (m)']
329 L = self.L(params[0])
333 a = self.info['Kuhn length (m)']
335 self._model_data[:] = FJC_fn(x_data, T, L, a)
336 return self._model_data
338 def fit(self, *args, **kwargs):
339 params = super(FJC, self).fit(*args, **kwargs)
340 params[0] = self.L(params[0]) # convert Lp -> L
342 params[1] = abs(params[1]) # take the absolute value of `a`
345 def guess_initial_params(self, outqueue=None):
346 """Guess initial fitting parameters.
351 A guess at the reparameterized contour length in meters.
353 A guess at the Kuhn length in meters. If `.info` has a
354 setting for `'Kuhn length (m)'`, `a` is not returned.
356 x_max = self.info['x data (m)'].max()
357 a_given = 'Kuhn length (m)' in self.info
359 a = self.info['Kuhn length (m)']
369 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):
370 """Inverse poly(ethylene-glycol) adjusted extended FJC model.
374 >>> kwargs = {'T':300.0, 'N':1, 'k':150.0, 'Lp':3.58e-10, 'Lh':2.8e-10,
375 ... 'dG':3.0, 'a':7e-10}
376 >>> inverse_FJC_PEG_fn(F_data=200e-12, **kwargs) # doctest: +ELLIPSIS
381 The function is that proposed by F. Oesterhelt, et al. [#]_.
384 x(F) = N_\text{S} \cdot \left[
386 \frac{L_\text{planar}}{
387 \exp\left(-\Delta G/k_B T\right) + 1}
388 + \frac{L_\text{helical}}{
389 \exp\left(+\Delta G/k_B T\right) + 1}
391 \coth\left(\frac{Fa}{k_B T}\right)
394 + \frac{F}{K_\text{S}}
396 where :math:`N_\text{S}` is the number of segments,
397 :math:`K_\text{S}` is the segment elasticicty,
398 :math:`L_\text{planar}` is the ttt contour length,
399 :math:`L_\text{helical}` is the ttg contour length,
400 :math:`\Delta G` is the Gibbs free energy difference
401 :math:`G_\text{planar}-G_\text{helical}`,
402 :math:`a` is the Kuhn length,
403 and :math:`F` is the chain tension.
405 .. [#]: F. Oesterhelt, M. Rief, and H.E. Gaub.
406 "Single molecule force spectroscopy by AFM indicates helical
407 structure of poly(ethylene-glycol) in water."
408 New Journal of Physics. 1999.
409 doi: `10.1088/1367-2630/1/1/006 <http://dx.doi.org/10.1088/1367-2630/1/1/006>`
413 g = (dG - F_data*(Lp-Lh)) / kBT
415 return N * ((Lp/(numpy.exp(-g)+1) + Lh/(numpy.exp(+g)+1)) * (coth(z)-1./z)
418 def FJC_PEG_fn(x_data, T, N, k, Lp, Lh, dG, a):
419 """Poly(ethylene-glycol) adjusted extended FJC model.
423 z : float or array_like
424 object whose inverse Langevin will be returned
426 value (close to one) for which we assume the inverse is +/-
427 infinity. This avoids problems with Newton-Raphson
428 convergence in regions with low slope.
432 >>> kwargs = {'T':300.0, 'N':1, 'k':150.0, 'Lp':3.58e-10, 'Lh':2.8e-10,
433 ... 'dG':3.0, 'a':7e-10}
434 >>> FJC_PEG_fn(inverse_FJC_PEG_fn(1e-11, **kwargs), **kwargs) # doctest: +ELLIPSIS
436 >>> FJC_PEG_fn(inverse_FJC_PEG_fn(3.4e-10, **kwargs), **kwargs) # doctest: +ELLIPSIS
438 >>> FJC_PEG_fn(numpy.array([1e-10,2e-10,3e-10]), **kwargs) # doctest: +ELLIPSIS
439 array([ 5.207...e-12, 1.257...e-11, 3.636...e-11])
440 >>> kwargs['N'] = 123
441 >>> FJC_PEG_fn(numpy.array([1e-8,2e-8,3e-8]), **kwargs) # doctest: +ELLIPSIS
442 array([ 4.160...e-12, 9.302...e-12, 1.830...e-11])
446 We approximate the PEG adjusted eFJC via Newton's method
447 (:func:`scipy.optimize.newton`). For the starting point, we use
448 the standard FJC function with an averaged contour length.
450 kwargs = {'T':T, 'N':N, 'k':k, 'Lp':Lp, 'Lh':Lh, 'dG':dG, 'a':a}
451 if is_iterable(x_data):
452 ret = numpy.ndarray(shape=x_data.shape, dtype=x_data.dtype)
453 for i,x in enumerate(x_data):
454 ret[i] = FJC_PEG_fn(x, **kwargs)
458 # Approximate with the simple FJC_fn()
461 while guess == numpy.inf:
462 guess = FJC_fn(x_data, T=T, L=L, a=a)
465 fn = lambda f: inverse_FJC_PEG_fn(guess*abs(f), **kwargs) - x_data
466 ret = guess*abs(newton(func=fn, x0=1.0))
469 class FJC_PEG (ModelFitter):
470 """Fit the data to an poly(ethylene-glycol) adjusted extended FJC.
474 Generate some example data
476 >>> kwargs = {'T':300.0, 'N':123, 'k':150.0, 'Lp':3.58e-10, 'Lh':2.8e-10,
477 ... 'dG':3.0, 'a':7e-10}
478 >>> x_data = numpy.linspace(10e-9, 30e-9, num=20)
479 >>> d_data = FJC_PEG_fn(x_data, **kwargs)
481 Fit the example data with a two-parameter fit (`N` and `a`).
484 ... 'temperature (K)': kwargs['T'],
485 ... 'x data (m)': x_data,
486 ... 'section elasticity (N/m)': kwargs['k'],
487 ... 'planar section length (m)': kwargs['Lp'],
488 ... 'helical section length (m)': kwargs['Lh'],
489 ... 'Gibbs free energy difference (Gp - Gh) (kBT)': kwargs['dG'],
491 >>> model = FJC_PEG(d_data, info=info, rescale=True)
492 >>> outqueue = Queue()
493 >>> N,a = model.fit(outqueue=outqueue)
494 >>> fit_info = outqueue.get(block=False)
500 Fit the example data with a one-parameter fit (`N`). We introduce
501 some error in our fixed Kuhn length for fun.
503 >>> info['Kuhn length (m)'] = 2*kwargs['a']
504 >>> model = FJC_PEG(d_data, info=info, rescale=True)
505 >>> N, = model.fit(outqueue=outqueue)
506 >>> fit_info = outqueue.get(block=False)
507 >>> print N # doctest: +ELLIPSIS
511 """To avoid invalid values of `L`, we reparametrize with `Lr`.
516 contour length in meters
521 rescaled version of `L`.
525 >>> x_data = numpy.linspace(1e-9, 20e-9, num=10)
526 >>> model = FJC_PEG(data=x_data, info={'x data (m)':x_data})
527 >>> model.Lr(20e-9) # doctest: +ELLIPSIS
529 >>> model.Lr(19e-9) # doctest: +ELLIPSIS
531 >>> model.Lr(21e-9) # doctest: +ELLIPSIS
533 >>> model.Lr(100e-9) # doctest: +ELLIPSIS
538 The rescaling is designed to limit `L` to values strictly
539 greater than zero, so we use
542 L = \exp(L_p) * x_\text{max}
544 which has the inverse
547 L_p = \ln(L/x_\text{max})
549 This will obviously effect the interpretation of the fit's covariance matrix.
551 x_max = self.info['x data (m)'].max()
552 return numpy.log(L/x_max)
555 """Inverse of :meth:`Lr`.
560 rescaled version of `L`
565 contour length in meters
569 >>> x_data = numpy.linspace(1e-9, 20e-9, num=10)
570 >>> model = FJC_PEG(data=x_data, info={'x data (m)':x_data})
571 >>> model.L(model.Lr(21e-9)) # doctest: +ELLIPSIS
573 >>> model.L(model.Lr(100e-9)) # doctest: +ELLIPSIS
576 x_max = self.info['x data (m)'].max()
577 return numpy.exp(Lr)*x_max
581 'T': self.info['temperature (K)'],
582 'k': self.info['section elasticity (N/m)'],
583 'Lp': self.info['planar section length (m)'],
584 'Lh': self.info['helical section length (m)'],
585 'dG': self.info['Gibbs free energy difference (Gp - Gh) (kBT)'],
588 def model(self, params):
589 """Extract the relavant arguments and use :func:`FJC_PEG_fn`.
593 params : list of floats
594 `[N, a]`, where the presence of `a` is optional.
596 # setup convenient aliases
597 T = self.info['temperature (K)']
598 x_data = self.info['x data (m)']
599 N = self.L(params[0])
603 a = self.info['Kuhn length (m)']
604 kwargs = self._kwargs()
606 self._model_data[:] = FJC_PEG_fn(x_data, N=N, a=a, **kwargs)
607 return self._model_data
609 def fit(self, *args, **kwargs):
610 params = super(FJC_PEG, self).fit(*args, **kwargs)
611 params[0] = self.L(params[0]) # convert Nr -> N
613 params[1] = abs(params[1]) # take the absolute value of `a`
616 def guess_initial_params(self, outqueue=None):
617 """Guess initial fitting parameters.
622 A guess at the section count.
624 A guess at the Kuhn length in meters. If `.info` has a
625 setting for `'Kuhn length (m)'`, `a` is not returned.
627 x_max = self.info['x data (m)'].max()
628 a_given = 'Kuhn length (m)' in self.info
630 a = self.info['Kuhn length (m)']
633 f_max = self._data.max()
634 kwargs = self._kwargs()
636 x_section = inverse_FJC_PEG_fn(F_data=f_max, **kwargs)
637 N = x_max / x_section;
644 def WLC_fn(x_data, T, L, p):
645 """The worm like chain interpolation model.
650 x values for which the WLC tension should be calculated.
652 temperature in Kelvin.
654 contour length in meters.
656 persistence length in meters.
665 >>> WLC_fn(numpy.array([1e-9, 5e-9, 10e-9], dtype=numpy.float),
666 ... T=300, L=15e-9, p=2.5e-10) # doctest: +ELLIPSIS
667 array([ 1.717...e-12, 1.070...e-11, 4.418...e-11])
671 The function is the simple polynomial worm like chain
672 interpolation forumula proposed by C. Bustamante, et al. [#]_.
675 F(x) = \frac{k_B T}{p} \left[
676 \frac{1}{4}\left( \frac{1}{(1-x/L)^2} - 1 \right)
680 .. [#] C. Bustamante, J.F. Marko, E.D. Siggia, and S.Smith.
681 "Entropic elasticity of lambda-phage DNA."
683 doi: `10.1126/science.8079175 <http://dx.doi.org/10.1126/science.8079175>`
686 scaled_data = x_data / L
687 return a * (0.25*((1-scaled_data)**-2 - 1) + scaled_data)
689 class WLC (ModelFitter):
690 """Fit the data to a worm like chain.
694 Generate some example data
697 >>> L = 35e-9 # meters
698 >>> p = 2.5e-10 # meters
699 >>> x_data = numpy.linspace(10e-9, 30e-9, num=20)
700 >>> d_data = WLC_fn(x_data, T=T, L=L, p=p)
702 Fit the example data with a two-parameter fit (`L` and `p`).
705 ... 'temperature (K)': T,
706 ... 'x data (m)': x_data,
708 >>> model = WLC(d_data, info=info, rescale=True)
709 >>> outqueue = Queue()
710 >>> L,p = model.fit(outqueue=outqueue)
711 >>> fit_info = outqueue.get(block=False)
717 Fit the example data with a one-parameter fit (`L`). We introduce
718 some error in our fixed persistence length for fun.
720 >>> info['persistence length (m)'] = 2*p
721 >>> model = WLC(d_data, info=info, rescale=True)
722 >>> L, = model.fit(outqueue=outqueue)
723 >>> fit_info = outqueue.get(block=False)
724 >>> print L # doctest: +ELLIPSIS
728 """To avoid invalid values of `L`, we reparametrize with `Lp`.
733 contour length in meters
738 rescaled version of `L`.
742 >>> x_data = numpy.linspace(1e-9, 20e-9, num=10)
743 >>> model = WLC(data=x_data, info={'x data (m)':x_data})
744 >>> model.Lp(20e-9) # doctest: +ELLIPSIS
746 >>> model.Lp(19e-9) # doctest: +ELLIPSIS
748 >>> model.Lp(21e-9) # doctest: +ELLIPSIS
750 >>> model.Lp(100e-9) # doctest: +ELLIPSIS
755 The rescaling is designed to limit `L` to values strictly
756 greater than the maximum `x` data value, so we use
759 L = [\exp(L_p))+1] * x_\text{max}
761 which has the inverse
764 L_p = \ln(L/x_\text{max}-1)
766 This will obviously effect the interpretation of the fit's covariance matrix.
768 x_max = self.info['x data (m)'].max()
769 return numpy.log(L/x_max - 1)
772 """Inverse of :meth:`Lp`.
777 rescaled version of `L`
782 contour length in meters
786 >>> x_data = numpy.linspace(1e-9, 20e-9, num=10)
787 >>> model = WLC(data=x_data, info={'x data (m)':x_data})
788 >>> model.L(model.Lp(21e-9)) # doctest: +ELLIPSIS
790 >>> model.L(model.Lp(100e-9)) # doctest: +ELLIPSIS
793 x_max = self.info['x data (m)'].max()
794 return (numpy.exp(Lp)+1)*x_max
796 def model(self, params):
797 """Extract the relavant arguments and use :func:`WLC_fn`.
801 params : list of floats
802 `[Lp, p]`, where the presence of `p` is optional.
804 # setup convenient aliases
805 T = self.info['temperature (K)']
806 x_data = self.info['x data (m)']
807 L = self.L(params[0])
811 p = self.info['persistence length (m)']
813 self._model_data[:] = WLC_fn(x_data, T, L, p)
814 return self._model_data
816 def fit(self, *args, **kwargs):
817 params = super(WLC, self).fit(*args, **kwargs)
818 params[0] = self.L(params[0]) # convert Lp -> L
820 params[1] = abs(params[1]) # take the absolute value of `p`
823 def guess_initial_params(self, outqueue=None):
824 """Guess initial fitting parameters.
829 A guess at the reparameterized contour length in meters.
831 A guess at the persistence length in meters. If `.info`
832 has a setting for `'persistence length (m)'`, `p` is not
835 x_max = self.info['x data (m)'].max()
836 p_given = 'persistence length (m)' in self.info
838 p = self.info['persistence length (m)']
848 class PolymerFitPlugin (Plugin):
849 """Polymer model (WLC, FJC, etc.) fitting.
852 super(PolymerFitPlugin, self).__init__(name='polymer_fit')
853 self._arguments = [ # For Command initialization
854 Argument('polymer model', default='WLC', help="""
855 Select the default polymer model for 'polymer fit'. See the
856 documentation for descriptions of available algorithms.
858 Argument('FJC Kuhn length', type='float', default=4e-10,
859 help='Kuhn length in meters'),
860 Argument('FJC_PEG Kuhn length', type='float', default=4e-10,
861 help='Kuhn length in meters'),
862 Argument('FJC_PEG elasticity', type='float', default=150.0,
863 help='Elasticity of a PEG segment in Newtons per meter.'),
864 Argument('FJC_PEG delta G', type='float', default=3.0, help="""
865 Gibbs free energy difference between trans-trans-trans (ttt) and
866 trans-trans-gauche (ttg) PEG states in units of kBT.
868 Argument('FJC_PEG L_helical', type='float', default=2.8e-10,
869 help='Contour length of PEG in the ttg state.'),
870 Argument('FJC_PEG L_planar', type='float', default=3.58e-10,
871 help='Contour length of PEG in the ttt state.'),
872 Argument('WLC persistence length', type='float', default=4e-10,
873 help='Persistence length in meters'),
876 Setting(section=self.setting_section, help=self.__doc__)]
877 for argument in self._arguments:
878 self._settings.append(argument_to_setting(
879 self.setting_section, argument))
880 argument.default = None # if argument isn't given, use the config.
881 self._input_columns = [
882 ('distance column', 'cantilever adjusted extension (m)',
884 Name of the column to use as the surface position input.
886 ('deflection column', 'deflection (N)',
888 Name of the column to use as the deflection input.
892 PolymerFitCommand(self), PolymerFitPeaksCommand(self),
893 TranslateFlatPeaksCommand(self),
896 def dependencies(self):
899 def default_settings(self):
900 return self._settings
903 class PolymerFitCommand (ColumnAddingCommand):
904 """Polymer model (WLC, FJC, etc.) fitting.
906 Fits an entropic elasticity function to a given chunk of the
907 curve. Fit quality compares the residual with the thermal noise
908 (a good fit should be 1 or less).
910 Because the models are complicated and varied, you should
911 configure the command by setting configuration options instead of
912 using command arguments. TODO: argument_to_setting()
914 def __init__(self, plugin):
915 super(PolymerFitCommand, self).__init__(
917 columns=plugin._input_columns,
919 ('output tension column', 'polymer tension',
921 Name of the column (without units) to use as the polymer tension output.
925 Argument(name='fit parameters info name', type='string',
926 default='polymer fit',
928 Name (without units) for storing the fit parameters in the `.info` dictionary.
930 Argument(name='bounds', type='point', optional=False, count=2,
932 Indicies of points bounding the selected data.
934 ] + plugin._arguments,
935 help=self.__doc__, plugin=plugin)
937 def _run(self, hooke, inqueue, outqueue, params):
938 params = self._setup_params(hooke, params)
939 data = self._block(hooke, params)
940 model = params['polymer model']
941 dist_data = self._get_column(
942 hooke=hooke, params=params, column_name='distance column')
943 def_data = self._get_column(
944 hooke=hooke, params=params, column_name='deflection column')
945 start,stop = params['bounds']
946 tension_data,ps = self.fit_polymer_model(
947 params, dist_data, def_data, start, stop, outqueue)
948 data.info[params['fit parameters info name']] = ps
949 data.info[params['fit parameters info name']]['model'] = model
950 self._set_column(hooke=hooke, params=params,
951 column_name='output tension column',
954 def _setup_params(self, hooke, params):
955 for key,value in params.items():
956 if value == None: # Use configured default value.
957 params[key] = self.plugin.config[key]
958 name,def_unit = split_data_label(params['deflection column'])
959 params['output tension column'] = join_data_label(
960 params['output tension column'], def_unit)
963 def fit_polymer_model(self, params, dist_data, def_data, start, stop,
965 """Railyard for the `fit_*_model` family.
967 Uses the `polymer model` configuration setting to call the
968 appropriate backend algorithm.
970 fn = getattr(self, 'fit_%s_model'
971 % params['polymer model'].replace('-','_'))
972 return fn(params, dist_data, def_data, start, stop, outqueue)
974 def fit_FJC_model(self, params, dist_data, def_data, start, stop,
976 """Fit the data with :class:`FJC`.
979 'temperature (K)': self.plugin.config['temperature'],
980 'x data (m)': dist_data[start:stop],
982 if True: # TODO: optionally free persistence length
983 info['Kuhn length (m)'] = (
984 params['FJC Kuhn length'])
985 model = FJC(def_data[start:stop], info=info, rescale=True)
987 params = model.fit(outqueue=queue)
988 if True: # TODO: if Kuhn length fixed
989 a = info['Kuhn length (m)']
993 T = info['temperature (K)']
994 fit_info = queue.get(block=False)
995 f_data = numpy.ones(dist_data.shape, dtype=dist_data.dtype) * numpy.nan
996 f_data[start:stop] = FJC_fn(dist_data[start:stop], T=T, L=L, a=a)
997 return [f_data, fit_info]
999 def fit_FJC_PEG_model(self, params, dist_data, def_data, start, stop,
1001 """Fit the data with :class:`FJC_PEG`.
1004 'temperature (K)': self.plugin.config['temperature'],
1005 'x data (m)': dist_data[start:stop],
1008 if True: # TODO: optionally free persistence length
1009 info['Kuhn length (m)'] = (
1010 params['FJC Kuhn length'])
1011 model = FJC_PEG(def_data[start:stop], info=info, rescale=True)
1013 params = model.fit(outqueue=queue)
1014 if True: # TODO: if Kuhn length fixed
1015 a = info['Kuhn length (m)']
1019 T = info['temperature (K)']
1020 fit_info = queue.get(block=False)
1021 f_data = numpy.ones(dist_data.shape, dtype=dist_data.dtype) * numpy.nan
1022 f_data[start:stop] = FJC_PEG_fn(dist_data[start:stop], **kwargs)
1023 return [f_data, fit_info]
1025 def fit_WLC_model(self, params, dist_data, def_data, start, stop,
1027 """Fit the data with :class:`WLC`.
1030 'temperature (K)': self.plugin.config['temperature'],
1031 'x data (m)': dist_data[start:stop],
1033 if True: # TODO: optionally free persistence length
1034 info['persistence length (m)'] = (
1035 params['WLC persistence length'])
1036 model = WLC(def_data[start:stop], info=info, rescale=True)
1038 params = model.fit(outqueue=queue)
1039 if True: # TODO: if persistence length fixed
1040 p = info['persistence length (m)']
1043 info['persistence length (m)'] = p
1045 info['contour length (m)'] = L
1046 T = info['temperature (K)']
1047 fit_info = queue.get(block=False)
1048 info['fit'] = fit_info
1049 info.pop('x data (m)')
1050 f_data = numpy.ones(dist_data.shape, dtype=dist_data.dtype) * numpy.nan
1051 f_data[start:stop] = WLC_fn(dist_data[start:stop], T=T, L=L, p=p)
1052 return [f_data, info]
1055 class PolymerFitPeaksCommand (ColumnAccessCommand):
1056 """Polymer model (WLC, FJC, etc.) fitting.
1058 Use :class:`PolymerFitCommand` to fit the each peak in a list of
1059 previously determined peaks.
1061 def __init__(self, plugin):
1062 super(PolymerFitPeaksCommand, self).__init__(
1063 name='polymer fit peaks',
1064 columns=plugin._input_columns,
1066 Argument(name='peak info name', type='string',
1067 default='polymer peaks',
1069 Name for the list of peaks in the `.info` dictionary.
1071 Argument(name='peak index', type='int', count=-1, default=None,
1073 Index of the selected peak in the list of peaks. Use `None` to fit all peaks.
1075 ] + plugin._arguments,
1076 help=self.__doc__, plugin=plugin)
1078 def _run(self, hooke, inqueue, outqueue, params):
1079 data = self._block(hooke, params)
1080 fit_command = hooke.command_by_name['polymer fit']
1083 curve = params['curve']
1084 params['curve'] = None
1085 p = copy.deepcopy(params)
1086 p['curve'] = params['curve']
1087 del(p['peak info name'])
1088 del(p['peak index'])
1089 for i,peak in enumerate(data.info[params['peak info name']]):
1090 if params['peak index'] == None or i in params['peak index']:
1091 p['bounds'] = [peak.index, peak.post_index()]
1092 p['output tension column'] = peak.name
1093 p['fit parameters info name'] = peak.name
1094 fit_command.run(hooke, inq, outq, **p)
1096 if not isinstance(ret, Success):
1100 class TranslateFlatPeaksCommand (ColumnAccessCommand):
1101 """Translate flat filter peaks into polymer peaks for fitting.
1103 Use :class:`~hooke.plugin.flatfilt.FlatPeaksCommand` creates a
1104 list of peaks for regions with large derivatives. For velocity
1105 clamp measurements, these regions are usually the rebound phase
1106 after a protein domain unfolds, the cantilever detaches, etc.
1107 Because these features occur after the polymer loading phase, we
1108 need to shift the selected regions back to align them with the
1109 polymer loading regions.
1111 def __init__(self, plugin):
1112 super(TranslateFlatPeaksCommand, self).__init__(
1113 name='flat peaks to polymer peaks',
1114 columns=plugin._input_columns,
1116 Argument(name='input peak info name', type='string',
1117 default='flat filter peaks',
1119 Name for the input peaks in the `.info` dictionary.
1121 Argument(name='output peak info name', type='string',
1122 default='polymer peaks',
1124 Name for the ouptput peaks in the `.info` dictionary.
1126 Argument(name='end offset', type='int', default=-1,
1128 Number of points between the end of a new peak and the start of the old.
1130 Argument(name='start fraction', type='float', default=0.2,
1132 Place the start of the new peak at `start_fraction` from the end of
1133 the previous old peak to the end of the new peak. Because the first
1134 new peak will have no previous old peak, it uses a distance of zero
1137 ] + plugin._arguments,
1138 help=self.__doc__, plugin=plugin)
1140 def _run(self, hooke, inqueue, outqueue, params):
1141 data = self._block(hooke, params)
1142 dist_data = self._get_column(
1143 hooke=hooke, params=params, column_name='distance column')
1144 def_data = self._get_column(
1145 hooke=hooke, params=params, column_name='deflection column')
1146 previous_old_stop = numpy.absolute(dist_data).argmin()
1149 peaks = data.info[params['input peak info name']]
1151 raise Failure('No value for %s in %s.info (%s): %s'
1152 % (params['input peak info name'], data.info['name'],
1153 sorted(data.info.keys()), e))
1154 for i,peak in enumerate(peaks):
1155 next_old_start = peak.index
1156 stop = next_old_start + params['end offset']
1158 dist_data[previous_old_stop]
1159 + params['start fraction']*(
1160 dist_data[stop] - dist_data[previous_old_stop]))
1161 start = numpy.absolute(dist_data - dist_start).argmin()
1162 p = Peak('polymer peak %d' % i,
1164 values=def_data[start:stop])
1166 previous_old_stop = peak.post_index()
1167 data.info[params['output peak info name']] = new
1171 # def dist2fit(self):
1172 # '''Calculates the average distance from data to fit, scaled by
1173 # the standard deviation of the free cantilever area (thermal