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.).
35 from scipy.optimize import newton
37 from ..command import Command, Argument, Failure
38 from ..config import Setting
39 from ..curve import Data
40 from ..plugin import Plugin
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
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.Queue()
235 >>> Lp,a = model.fit(outqueue=outqueue)
236 >>> fit_info = outqueue.get(block=False)
237 >>> model.L(Lp) # doctest: +ELLIPSIS
239 >>> a # doctest: +ELLIPSIS
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 >>> Lp = model.fit(outqueue=outqueue)
248 >>> fit_info = outqueue.get(block=False)
249 >>> model.L(Lp) # 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 guess_initial_params(self, outqueue=None):
342 """Guess initial fitting parameters.
347 A guess at the reparameterized contour length in meters.
349 A guess at the Kuhn length in meters. If `.info` has a
350 setting for `'Kuhn length (m)'`, `a` is not returned.
352 x_max = self.info['x data (m)'].max()
353 a_given = 'Kuhn length (m)' in self.info
355 a = self.info['Kuhn length (m)']
364 def guess_scale(self, params, outqueue=None):
365 """Guess parameter scales.
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.
378 return [Lp_scale] + [x/10.0 for x in params[1:]]
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.
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
393 The function is that proposed by F. Oesterhelt, et al. [#]_.
396 x(F) = N_\text{S} \cdot \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}
403 \coth\left(\frac{Fa}{k_B T}\right)
406 + \frac{F}{K_\text{S}}
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.
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>`
425 g = (dG - F_data*(Lp-Lh)) / kBT
427 return N * ((Lp/(numpy.exp(-g)+1) + Lh/(numpy.exp(+g)+1)) * (coth(z)-1./z)
430 def FJC_PEG_fn(x_data, T, N, k, Lp, Lh, dG, a):
431 """Poly(ethylene-glycol) adjusted extended FJC model.
435 z : float or array_like
436 object whose inverse Langevin will be returned
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.
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
448 >>> FJC_PEG_fn(inverse_FJC_PEG_fn(3.4e-10, **kwargs), **kwargs) # doctest: +ELLIPSIS
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])
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.
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)
470 # Approximate with the simple FJC_fn()
473 while guess == numpy.inf:
474 guess = FJC_fn(x_data, T=T, L=L, a=a)
477 fn = lambda f: inverse_FJC_PEG_fn(guess*abs(f), **kwargs) - x_data
478 ret = guess*abs(newton(func=fn, x0=1.0))
481 class FJC_PEG (ModelFitter):
482 """Fit the data to an poly(ethylene-glycol) adjusted extended FJC.
486 Generate some example data
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)
493 Fit the example data with a two-parameter fit (`N` and `a`).
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'],
503 >>> model = FJC_PEG(d_data, info=info, rescale=True)
504 >>> outqueue = Queue.Queue()
505 >>> Nr,a = model.fit(outqueue=outqueue)
506 >>> fit_info = outqueue.get(block=False)
507 >>> model.L(Nr) # doctest: +ELLIPSIS
509 >>> a # doctest: +ELLIPSIS
512 Fit the example data with a one-parameter fit (`N`). We introduce
513 some error in our fixed Kuhn length for fun.
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
523 """To avoid invalid values of `L`, we reparametrize with `Lr`.
528 contour length in meters
533 rescaled version of `L`.
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
541 >>> model.Lr(19e-9) # doctest: +ELLIPSIS
543 >>> model.Lr(21e-9) # doctest: +ELLIPSIS
545 >>> model.Lr(100e-9) # doctest: +ELLIPSIS
550 The rescaling is designed to limit `L` to values strictly
551 greater than zero, so we use
554 L = \exp(L_p) * x_\text{max}
556 which has the inverse
559 L_p = \ln(L/x_\text{max})
561 This will obviously effect the interpretation of the fit's covariance matrix.
563 x_max = self.info['x data (m)'].max()
564 return numpy.log(L/x_max)
567 """Inverse of :meth:`Lr`.
572 rescaled version of `L`
577 contour length in meters
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
585 >>> model.L(model.Lr(100e-9)) # doctest: +ELLIPSIS
588 x_max = self.info['x data (m)'].max()
589 return numpy.exp(Lr)*x_max
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)'],
600 def model(self, params):
601 """Extract the relavant arguments and use :func:`FJC_PEG_fn`.
605 params : list of floats
606 `[N, a]`, where the presence of `a` is optional.
608 # setup convenient aliases
609 T = self.info['temperature (K)']
610 x_data = self.info['x data (m)']
611 N = self.L(params[0])
615 a = self.info['Kuhn length (m)']
616 kwargs = self._kwargs()
618 self._model_data[:] = FJC_PEG_fn(x_data, N=N, a=a, **kwargs)
619 return self._model_data
621 def guess_initial_params(self, outqueue=None):
622 """Guess initial fitting parameters.
627 A guess at the section count.
629 A guess at the Kuhn length in meters. If `.info` has a
630 setting for `'Kuhn length (m)'`, `a` is not returned.
632 x_max = self.info['x data (m)'].max()
633 a_given = 'Kuhn length (m)' in self.info
635 a = self.info['Kuhn length (m)']
638 f_max = self._data.max()
639 kwargs = self._kwargs()
641 x_section = inverse_FJC_PEG_fn(F_data=f_max, **kwargs)
642 N = x_max / x_section;
648 def guess_scale(self, params, outqueue=None):
649 """Guess parameter scales.
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.
659 return [x/10.0 for x in params]
662 def WLC_fn(x_data, T, L, p):
663 """The worm like chain interpolation model.
668 x values for which the WLC tension should be calculated.
670 temperature in Kelvin.
672 contour length in meters.
674 persistence length in meters.
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])
689 The function is the simple polynomial worm like chain
690 interpolation forumula proposed by C. Bustamante, et al. [#]_.
693 F(x) = \frac{k_B T}{p} \left[
694 \frac{1}{4}\left( \frac{1}{(1-x/L)^2} - 1 \right)
698 .. [#] C. Bustamante, J.F. Marko, E.D. Siggia, and S.Smith.
699 "Entropic elasticity of lambda-phage DNA."
701 doi: `10.1126/science.8079175 <http://dx.doi.org/10.1126/science.8079175>`
704 scaled_data = x_data / L
705 return a * (0.25*((1-scaled_data)**-2 - 1) + scaled_data)
707 class WLC (ModelFitter):
708 """Fit the data to a worm like chain.
712 Generate some example data
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)
720 Fit the example data with a two-parameter fit (`L` and `p`).
723 ... 'temperature (K)': T,
724 ... 'x data (m)': x_data,
726 >>> model = WLC(d_data, info=info, rescale=True)
727 >>> outqueue = Queue.Queue()
728 >>> Lp,p = model.fit(outqueue=outqueue)
729 >>> fit_info = outqueue.get(block=False)
730 >>> model.L(Lp) # doctest: +ELLIPSIS
732 >>> p # doctest: +ELLIPSIS
735 Fit the example data with a one-parameter fit (`L`). We introduce
736 some error in our fixed persistence length for fun.
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
746 """To avoid invalid values of `L`, we reparametrize with `Lp`.
751 contour length in meters
756 rescaled version of `L`.
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
764 >>> model.Lp(19e-9) # doctest: +ELLIPSIS
766 >>> model.Lp(21e-9) # doctest: +ELLIPSIS
768 >>> model.Lp(100e-9) # doctest: +ELLIPSIS
773 The rescaling is designed to limit `L` to values strictly
774 greater than the maximum `x` data value, so we use
777 L = [\exp(L_p))+1] * x_\text{max}
779 which has the inverse
782 L_p = \ln(L/x_\text{max}-1)
784 This will obviously effect the interpretation of the fit's covariance matrix.
786 x_max = self.info['x data (m)'].max()
787 return numpy.log(L/x_max - 1)
790 """Inverse of :meth:`Lp`.
795 rescaled version of `L`
800 contour length in meters
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
808 >>> model.L(model.Lp(100e-9)) # doctest: +ELLIPSIS
811 x_max = self.info['x data (m)'].max()
812 return (numpy.exp(Lp)+1)*x_max
814 def model(self, params):
815 """Extract the relavant arguments and use :func:`WLC_fn`.
819 params : list of floats
820 `[Lp, p]`, where the presence of `p` is optional.
822 # setup convenient aliases
823 T = self.info['temperature (K)']
824 x_data = self.info['x data (m)']
825 L = self.L(params[0])
829 p = self.info['persistence length (m)']
831 self._model_data[:] = WLC_fn(x_data, T, L, p)
832 return self._model_data
834 def guess_initial_params(self, outqueue=None):
835 """Guess initial fitting parameters.
840 A guess at the reparameterized contour length in meters.
842 A guess at the persistence length in meters. If `.info`
843 has a setting for `'persistence length (m)'`, `p` is not
846 x_max = self.info['x data (m)'].max()
847 p_given = 'persistence length (m)' in self.info
849 p = self.info['persistence length (m)']
858 def guess_scale(self, params, outqueue=None):
859 """Guess parameter scales.
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
873 return [Lp_scale] + [x/10.0 for x in params[1:]]
876 class PolymerFitPlugin (Plugin):
877 """Polymer model (WLC, FJC, etc.) fitting.
880 super(PolymerFitPlugin, self).__init__(name='polymer_fit')
881 self._commands = [PolymerFitCommand(self),]
883 def dependencies(self):
886 def default_settings(self):
888 Setting(section=self.setting_section, help=self.__doc__),
889 Setting(section=self.setting_section,
890 option='polymer model',
892 help="Select the default polymer model for 'polymer fit'. See the documentation for descriptions of available algorithms."),
893 Setting(section=self.setting_section,
894 option='FJC Kuhn length',
895 value=4e-10, type='float',
896 help='Kuhn length in meters'),
897 Setting(section=self.setting_section,
898 option='FJC-PEG Kuhn length',
899 value=4e-10, type='float',
900 help='Kuhn length in meters'),
901 Setting(section=self.setting_section,
902 option='FJC-PEG elasticity',
903 value=150.0, type='float',
904 help='Elasticity of a PEG segment in Newtons per meter.'),
905 Setting(section=self.setting_section,
906 option='FJC-PEG delta G',
907 value=3.0, type='float',
908 help='Gibbs free energy difference between trans-trans-trans (ttt) and trans-trans-gauche (ttg) PEG states in units of kBT.'),
909 Setting(section=self.setting_section,
910 option='FJC-PEG L_helical',
911 value=2.8e-10, type='float',
912 help='Contour length of PEG in the ttg state.'),
913 Setting(section=self.setting_section,
914 option='FJC-PEG L_planar',
915 value=3.58e-10, type='float',
916 help='Contour length of PEG in the ttt state.'),
917 Setting(section=self.setting_section,
918 option='WLC persistence length',
919 value=4e-10, type='float',
920 help='Persistence length in meters'),
924 class PolymerFitCommand (Command):
925 """Polymer model (WLC, FJC, etc.) fitting.
927 Fits an entropic elasticity function to a given chunk of the
928 curve. Fit quality compares the residual with the thermal noise
929 (a good fit should be 1 or less).
931 Because the models are complicated and varied, you should
932 configure the command by setting configuration options instead of
933 using command arguments. TODO: argument_to_setting()
935 def __init__(self, plugin):
936 super(PolymerFitCommand, self).__init__(
940 Argument(name='block', aliases=['set'], type='int', default=0,
942 Data block for which the fit should be calculated. For an
943 approach/retract force curve, `0` selects the approaching curve and
944 `1` selects the retracting curve.
946 Argument(name='bounds', type='point', optional=False, count=2,
948 Indicies of points bounding the selected data.
950 Argument(name='input distance column', type='string',
951 default='cantilever adjusted extension (m)',
953 Name of the column to use as the surface positioning input.
955 Argument(name='input deflection column', type='string',
956 default='deflection (N)',
958 Name of the column to use as the deflection input.
960 Argument(name='output tension column', type='string',
961 default='polymer tension',
963 Name of the column (without units) to use as the polymer tension output.
965 Argument(name='fit parameters info name', type='string',
966 default='surface deflection offset',
968 Name (without units) for storing the fit parameters in the `.info` dictionary.
971 help=self.__doc__, plugin=plugin)
973 def _run(self, hooke, inqueue, outqueue, params):
974 scale(hooke, params['curve'], params['block']) # TODO: is autoscaling a good idea? (explicit is better than implicit)
975 data = params['curve'].data[params['block']]
976 # HACK? rely on params['curve'] being bound to the local hooke
977 # playlist (i.e. not a copy, as you would get by passing a
978 # curve through the queue). Ugh. Stupid queues. As an
979 # alternative, we could pass lookup information through the
981 model = self.plugin.config['polymer model']
982 new = Data((data.shape[0], data.shape[1]+1), dtype=data.dtype)
983 new.info = copy.deepcopy(data.info)
985 new.info['columns'].append(
986 join_data_label(params['output tension column'], 'N'))
987 z_data = data[:,data.info['columns'].index(
988 params['input distance column'])]
989 d_data = data[:,data.info['columns'].index(
990 params['input deflection column'])]
991 start,stop = params['bounds']
992 tension_data,ps = self.fit_polymer_model(
993 params['curve'], z_data, d_data, start, stop, outqueue)
994 new.info[params['fit parameters info name']] = ps
995 new.info[params['fit parameters info name']]['model'] = model
996 new[:,-1] = tension_data
997 params['curve'].data[params['block']] = new
999 def fit_polymer_model(self, curve, z_data, d_data, start, stop,
1001 """Railyard for the `fit_*_model` family.
1003 Uses the `polymer model` configuration setting to call the
1004 appropriate backend algorithm.
1006 fn = getattr(self, 'fit_%s_model'
1007 % self.plugin.config['polymer model'].replace('-','_'))
1008 return fn(curve, z_data, d_data, start, stop, outqueue)
1010 def fit_FJC_model(self, curve, z_data, d_data, start, stop,
1012 """Fit the data with :class:`FJC`.
1015 'temperature (K)': self.plugin.config['temperature'],
1016 'x data (m)': z_data[start:stop],
1018 if True: # TODO: optionally free persistence length
1019 info['Kuhn length (m)'] = (
1020 self.plugin.config['FJC Kuhn length'])
1021 model = FJC(d_data[start:stop], info=info, rescale=True)
1022 queue = Queue.Queue()
1023 params = model.fit(outqueue=queue)
1024 if True: # TODO: if Kuhn length fixed
1026 a = info['Kuhn length (m)']
1031 T = info['temperature (K)']
1032 fit_info = queue.get(block=False)
1033 mask = numpy.zeros(z_data.shape, dtype=numpy.bool)
1034 mask[start:stop] = True
1035 return [FJC_fn(z_data, T=T, L=L, a=a) * mask,
1038 def fit_FJC_PEG_model(self, curve, z_data, d_data, start, stop,
1040 """Fit the data with :class:`FJC_PEG`.
1043 'temperature (K)': self.plugin.config['temperature'],
1044 'x data (m)': z_data[start:stop],
1047 if True: # TODO: optionally free persistence length
1048 info['Kuhn length (m)'] = (
1049 self.plugin.config['FJC Kuhn length'])
1050 model = FJC(d_data[start:stop], info=info, rescale=True)
1051 queue = Queue.Queue()
1052 params = model.fit(outqueue=queue)
1053 if True: # TODO: if Kuhn length fixed
1055 a = info['Kuhn length (m)']
1060 T = info['temperature (K)']
1061 fit_info = queue.get(block=False)
1062 mask = numpy.zeros(z_data.shape, dtype=numpy.bool)
1063 mask[start:stop] = True
1064 return [FJC_PEG_fn(z_data, **kwargs) * mask,
1067 def fit_WLC_model(self, curve, z_data, d_data, start, stop,
1069 """Fit the data with :class:`WLC`.
1072 'temperature (K)': self.plugin.config['temperature'],
1073 'x data (m)': z_data[start:stop],
1075 if True: # TODO: optionally free persistence length
1076 info['persistence length (m)'] = (
1077 self.plugin.config['WLC persistence length'])
1078 model = WLC(d_data[start:stop], info=info, rescale=True)
1079 queue = Queue.Queue()
1080 params = model.fit(outqueue=queue)
1081 if True: # TODO: if persistence length fixed
1083 p = info['persistence length (m)']
1088 T = info['temperature (K)']
1089 fit_info = queue.get(block=False)
1090 mask = numpy.zeros(z_data.shape, dtype=numpy.bool)
1091 mask[start:stop] = True
1092 return [WLC_fn(z_data, T=T, L=L, p=p) * mask,
1097 # def dist2fit(self):
1098 # '''Calculates the average distance from data to fit, scaled by
1099 # the standard deviation of the free cantilever area (thermal