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.fit import PoorFit, ModelFitter
42 from ..util.callback import is_iterable
43 from .curve import CurveArgument
44 from .vclamp import scale
47 kB = 1.3806503e-23 # Joules/Kelvin
51 """Hyperbolic cotangent.
55 >>> coth(1.19967874) # doctest: +ELLIPSIS
63 \coth(z) \equiv \frac{\exp(z) + \exp(-z)}{\exp(z) - \exp(-z)}
67 http://mathworld.wolfram.com/HyperbolicCotangent.html
69 return 1.0/numpy.tanh(z)
72 """Inverse hyperbolic cotangent.
76 >>> arccoth(1.19967874) # doctest: +ELLIPSIS
78 >>> arccoth(coth(numpy.pi)) # doctest: +ELLIPSIS
83 Inverting from the :func:`definition <coth>`.
86 z \equiv \atanh\left( \frac{1}{\coth(z)} \right)
88 return numpy.arctanh(1.0/z)
95 >>> Langevin(numpy.pi) # doctest: +ELLIPSIS
100 From `Wikipedia`_ or Hatfield [#]_
103 L(z) \equiv \coth(z) - \frac{1}{z}
106 http://en.wikipedia.org/wiki/Brillouin_and_Langevin_functions#Langevin_Function
108 .. [#]: J.W. Hatfield and S.R. Quake
109 "Dynamic Properties of an Extended Polymer in Solution."
110 Phys. Rev. Lett. 1999.
111 doi: `10.1103/PhysRevLett.82.3548 <http://dx.doi.org/10.1103/PhysRevLett.82.3548>`
113 return coth(z) - 1.0/z
115 def inverse_Langevin(z, extreme=1.0 - 1e-8):
116 """Inverse Langevin function.
120 z : float or array_like
121 object whose inverse Langevin will be returned
123 value (close to one) for which we assume the inverse is +/-
124 infinity. This avoids problems with Newton-Raphson
125 convergence in regions with low slope.
129 >>> inverse_Langevin(Langevin(numpy.pi)) # doctest: +ELLIPSIS
131 >>> inverse_Langevin(Langevin(numpy.array([1,2,3]))) # doctest: +ELLIPSIS
136 We approximate the inverse Langevin via Newton's method
137 (:func:`scipy.optimize.newton`). For the starting point, we take
138 the first three terms from the Taylor series (from `Wikipedia`_).
141 L^{-1}(z) = 3z + \frac{9}{5}z^3 + \frac{297}{175}z^5 + \dots
144 http://en.wikipedia.org/wiki/Brillouin_and_Langevin_functions#Langevin_Function
147 ret = numpy.ndarray(shape=z.shape, dtype=z.dtype)
148 for i,z_i in enumerate(z):
149 ret[i] = inverse_Langevin(z_i)
155 # Catch stdout since sometimes the newton function prints exit
156 # messages to stdout, e.g. "Tolerance of %s reached".
158 tmp_stdout = StringIO.StringIO()
159 sys.stdout = tmp_stdout
160 ret = newton(func=lambda x: Langevin(x)-z,
161 x0=3*z + 9./5.*z**3 + 297./175.*z**5,)
163 output = tmp_stdout.getvalue()
166 def FJC_fn(x_data, T, L, a):
167 """The freely jointed chain model.
172 x values for which the WLC tension should be calculated.
174 temperature in Kelvin.
176 contour length in meters.
178 Kuhn length in meters.
187 >>> FJC_fn(numpy.array([1e-9, 5e-9, 10e-9], dtype=numpy.float),
188 ... T=300, L=15e-9, a=2.5e-10) # doctest: +ELLIPSIS
189 array([ 3.322...-12, 1.78...e-11, 4.889...e-11])
193 The freely jointed chain model is
196 F(x) = \frac{k_B T}{a} L^{-1}\left( \frac{x}{L} \right)
198 where :math:`L^{-1}` is the :func:`inverse_Langevin`, :math:`a` is
199 the Kuhn length, and :math:`L` is the contour length [#]_. For
200 the inverse formulation (:math:`x(F)`), see Ray and Akhremitchev [#]_.
203 .. [#]: J.W. Hatfield and S.R. Quake
204 "Dynamic Properties of an Extended Polymer in Solution."
205 Phys. Rev. Lett. 1999.
206 doi: `10.1103/PhysRevLett.82.3548 <http://dx.doi.org/10.1103/PhysRevLett.82.3548>`
208 .. [#] C. Ray and B.B. Akhremitchev.
209 `"Fitting Force Curves by the Extended Freely Jointed Chain Model" <http://www.chem.duke.edu/~boris/research/force_spectroscopy/fit_efjc.pdf>`
211 return kB*T / a * inverse_Langevin(x_data/L)
213 class FJC (ModelFitter):
214 """Fit the data to a freely jointed chain.
218 Generate some example data
221 >>> L = 35e-9 # meters
222 >>> a = 2.5e-10 # meters
223 >>> x_data = numpy.linspace(10e-9, 30e-9, num=20)
224 >>> d_data = FJC_fn(x_data, T=T, L=L, a=a)
226 Fit the example data with a two-parameter fit (`L` and `a`).
229 ... 'temperature (K)': T,
230 ... 'x data (m)': x_data,
232 >>> model = FJC(d_data, info=info, rescale=True)
233 >>> outqueue = Queue.Queue()
234 >>> Lp,a = model.fit(outqueue=outqueue)
235 >>> fit_info = outqueue.get(block=False)
236 >>> model.L(Lp) # doctest: +ELLIPSIS
238 >>> a # doctest: +ELLIPSIS
241 Fit the example data with a one-parameter fit (`L`). We introduce
242 some error in our fixed Kuhn length for fun.
244 >>> info['Kuhn length (m)'] = 2*a
245 >>> model = FJC(d_data, info=info, rescale=True)
246 >>> Lp = model.fit(outqueue=outqueue)
247 >>> fit_info = outqueue.get(block=False)
248 >>> model.L(Lp) # doctest: +ELLIPSIS
252 """To avoid invalid values of `L`, we reparametrize with `Lp`.
257 contour length in meters
262 rescaled version of `L`.
266 >>> x_data = numpy.linspace(1e-9, 20e-9, num=10)
267 >>> model = FJC(data=x_data, info={'x data (m)':x_data})
268 >>> model.Lp(20e-9) # doctest: +ELLIPSIS
270 >>> model.Lp(19e-9) # doctest: +ELLIPSIS
272 >>> model.Lp(21e-9) # doctest: +ELLIPSIS
274 >>> model.Lp(100e-9) # doctest: +ELLIPSIS
279 The rescaling is designed to limit `L` to values strictly
280 greater than the maximum `x` data value, so we use
283 L = [\exp(L_p))+1] * x_\text{max}
285 which has the inverse
288 L_p = \ln(L/x_\text{max}-1)
290 This will obviously effect the interpretation of the fit's covariance matrix.
292 x_max = self.info['x data (m)'].max()
293 return numpy.log(L/x_max - 1)
296 """Inverse of :meth:`Lp`.
301 rescaled version of `L`
306 contour length in meters
310 >>> x_data = numpy.linspace(1e-9, 20e-9, num=10)
311 >>> model = FJC(data=x_data, info={'x data (m)':x_data})
312 >>> model.L(model.Lp(21e-9)) # doctest: +ELLIPSIS
314 >>> model.L(model.Lp(100e-9)) # doctest: +ELLIPSIS
317 x_max = self.info['x data (m)'].max()
318 return (numpy.exp(Lp)+1)*x_max
320 def model(self, params):
321 """Extract the relavant arguments and use :func:`FJC_fn`.
325 params : list of floats
326 `[Lp, a]`, where the presence of `a` is optional.
328 # setup convenient aliases
329 T = self.info['temperature (K)']
330 x_data = self.info['x data (m)']
331 L = self.L(params[0])
335 a = self.info['Kuhn length (m)']
337 self._model_data[:] = FJC_fn(x_data, T, L, a)
338 return self._model_data
340 def guess_initial_params(self, outqueue=None):
341 """Guess initial fitting parameters.
346 A guess at the reparameterized contour length in meters.
348 A guess at the Kuhn length in meters. If `.info` has a
349 setting for `'Kuhn length (m)'`, `a` is not returned.
351 x_max = self.info['x data (m)'].max()
352 a_given = 'Kuhn length (m)' in self.info
354 a = self.info['Kuhn length (m)']
363 def guess_scale(self, params, outqueue=None):
364 """Guess parameter scales.
369 A guess at the reparameterized contour length scale in meters.
370 a_scale : float (optional)
371 A guess at the Kuhn length in meters. If the length of
372 `params` is less than 2, `a_scale` is not returned.
377 return [Lp_scale] + [x/10.0 for x in params[1:]]
380 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):
381 """Inverse poly(ethylene-glycol) adjusted extended FJC model.
385 >>> kwargs = {'T':300.0, 'N':1, 'k':150.0, 'Lp':3.58e-10, 'Lh':2.8e-10,
386 ... 'dG':3.0, 'a':7e-10}
387 >>> inverse_FJC_PEG_fn(F_data=200e-12, **kwargs) # doctest: +ELLIPSIS
392 The function is that proposed by F. Oesterhelt, et al. [#]_.
395 x(F) = N_\text{S} \cdot \left[
397 \frac{L_\text{planar}}{
398 \exp\left(-\Delta G/k_B T\right) + 1}
399 + \frac{L_\text{helical}}{
400 \exp\left(+\Delta G/k_B T\right) + 1}
402 \coth\left(\frac{Fa}{k_B T}\right)
405 + \frac{F}{K_\text{S}}
407 where :math:`N_\text{S}` is the number of segments,
408 :math:`K_\text{S}` is the segment elasticicty,
409 :math:`L_\text{planar}` is the ttt contour length,
410 :math:`L_\text{helical}` is the ttg contour length,
411 :math:`\Delta G` is the Gibbs free energy difference
412 :math:`G_\text{planar}-G_\text{helical}`,
413 :math:`a` is the Kuhn length,
414 and :math:`F` is the chain tension.
416 .. [#]: F. Oesterhelt, M. Rief, and H.E. Gaub.
417 "Single molecule force spectroscopy by AFM indicates helical
418 structure of poly(ethylene-glycol) in water."
419 New Journal of Physics. 1999.
420 doi: `10.1088/1367-2630/1/1/006 <http://dx.doi.org/10.1088/1367-2630/1/1/006>`
424 g = (dG - F_data*(Lp-Lh)) / kBT
426 return N * ((Lp/(numpy.exp(-g)+1) + Lh/(numpy.exp(+g)+1)) * (coth(z)-1./z)
429 def FJC_PEG_fn(x_data, T, N, k, Lp, Lh, dG, a):
430 """Poly(ethylene-glycol) adjusted extended FJC model.
434 z : float or array_like
435 object whose inverse Langevin will be returned
437 value (close to one) for which we assume the inverse is +/-
438 infinity. This avoids problems with Newton-Raphson
439 convergence in regions with low slope.
443 >>> kwargs = {'T':300.0, 'N':1, 'k':150.0, 'Lp':3.58e-10, 'Lh':2.8e-10,
444 ... 'dG':3.0, 'a':7e-10}
445 >>> FJC_PEG_fn(inverse_FJC_PEG_fn(1e-11, **kwargs), **kwargs) # doctest: +ELLIPSIS
447 >>> FJC_PEG_fn(inverse_FJC_PEG_fn(3.4e-10, **kwargs), **kwargs) # doctest: +ELLIPSIS
449 >>> FJC_PEG_fn(numpy.array([1e-10,2e-10,3e-10]), **kwargs) # doctest: +ELLIPSIS
450 array([ 5.207...e-12, 1.257...e-11, 3.636...e-11])
451 >>> kwargs['N'] = 123
452 >>> FJC_PEG_fn(numpy.array([1e-8,2e-8,3e-8]), **kwargs) # doctest: +ELLIPSIS
453 array([ 4.160...e-12, 9.302...e-12, 1.830...e-11])
457 We approximate the PEG adjusted eFJC via Newton's method
458 (:func:`scipy.optimize.newton`). For the starting point, we use
459 the standard FJC function with an averaged contour length.
461 kwargs = {'T':T, 'N':N, 'k':k, 'Lp':Lp, 'Lh':Lh, 'dG':dG, 'a':a}
462 if is_iterable(x_data):
463 ret = numpy.ndarray(shape=x_data.shape, dtype=x_data.dtype)
464 for i,x in enumerate(x_data):
465 ret[i] = FJC_PEG_fn(x, **kwargs)
469 # Approximate with the simple FJC_fn()
472 while guess == numpy.inf:
473 guess = FJC_fn(x_data, T=T, L=L, a=a)
476 fn = lambda f: inverse_FJC_PEG_fn(guess*abs(f), **kwargs) - x_data
477 ret = guess*abs(newton(func=fn, x0=1.0))
480 class FJC_PEG (ModelFitter):
481 """Fit the data to an poly(ethylene-glycol) adjusted extended FJC.
485 Generate some example data
487 >>> kwargs = {'T':300.0, 'N':123, 'k':150.0, 'Lp':3.58e-10, 'Lh':2.8e-10,
488 ... 'dG':3.0, 'a':7e-10}
489 >>> x_data = numpy.linspace(10e-9, 30e-9, num=20)
490 >>> d_data = FJC_PEG_fn(x_data, **kwargs)
492 Fit the example data with a two-parameter fit (`N` and `a`).
495 ... 'temperature (K)': kwargs['T'],
496 ... 'x data (m)': x_data,
497 ... 'section elasticity (N/m)': kwargs['k'],
498 ... 'planar section length (m)': kwargs['Lp'],
499 ... 'helical section length (m)': kwargs['Lh'],
500 ... 'Gibbs free energy difference (Gp - Gh) (kBT)': kwargs['dG'],
502 >>> model = FJC_PEG(d_data, info=info, rescale=True)
503 >>> outqueue = Queue.Queue()
504 >>> Nr,a = model.fit(outqueue=outqueue)
505 >>> fit_info = outqueue.get(block=False)
506 >>> model.L(Nr) # doctest: +ELLIPSIS
508 >>> a # doctest: +ELLIPSIS
511 Fit the example data with a one-parameter fit (`N`). We introduce
512 some error in our fixed Kuhn length for fun.
514 >>> info['Kuhn length (m)'] = 2*kwargs['a']
515 >>> model = FJC_PEG(d_data, info=info, rescale=True)
516 >>> Nr = model.fit(outqueue=outqueue)
517 >>> fit_info = outqueue.get(block=False)
518 >>> model.L(Nr) # doctest: +ELLIPSIS
522 """To avoid invalid values of `L`, we reparametrize with `Lr`.
527 contour length in meters
532 rescaled version of `L`.
536 >>> x_data = numpy.linspace(1e-9, 20e-9, num=10)
537 >>> model = FJC_PEG(data=x_data, info={'x data (m)':x_data})
538 >>> model.Lr(20e-9) # doctest: +ELLIPSIS
540 >>> model.Lr(19e-9) # doctest: +ELLIPSIS
542 >>> model.Lr(21e-9) # doctest: +ELLIPSIS
544 >>> model.Lr(100e-9) # doctest: +ELLIPSIS
549 The rescaling is designed to limit `L` to values strictly
550 greater than zero, so we use
553 L = \exp(L_p) * x_\text{max}
555 which has the inverse
558 L_p = \ln(L/x_\text{max})
560 This will obviously effect the interpretation of the fit's covariance matrix.
562 x_max = self.info['x data (m)'].max()
563 return numpy.log(L/x_max)
566 """Inverse of :meth:`Lr`.
571 rescaled version of `L`
576 contour length in meters
580 >>> x_data = numpy.linspace(1e-9, 20e-9, num=10)
581 >>> model = FJC_PEG(data=x_data, info={'x data (m)':x_data})
582 >>> model.L(model.Lr(21e-9)) # doctest: +ELLIPSIS
584 >>> model.L(model.Lr(100e-9)) # doctest: +ELLIPSIS
587 x_max = self.info['x data (m)'].max()
588 return numpy.exp(Lr)*x_max
592 'T': self.info['temperature (K)'],
593 'k': self.info['section elasticity (N/m)'],
594 'Lp': self.info['planar section length (m)'],
595 'Lh': self.info['helical section length (m)'],
596 'dG': self.info['Gibbs free energy difference (Gp - Gh) (kBT)'],
599 def model(self, params):
600 """Extract the relavant arguments and use :func:`FJC_PEG_fn`.
604 params : list of floats
605 `[N, a]`, where the presence of `a` is optional.
607 # setup convenient aliases
608 T = self.info['temperature (K)']
609 x_data = self.info['x data (m)']
610 N = self.L(params[0])
614 a = self.info['Kuhn length (m)']
615 kwargs = self._kwargs()
617 self._model_data[:] = FJC_PEG_fn(x_data, N=N, a=a, **kwargs)
618 return self._model_data
620 def guess_initial_params(self, outqueue=None):
621 """Guess initial fitting parameters.
626 A guess at the section count.
628 A guess at the Kuhn length in meters. If `.info` has a
629 setting for `'Kuhn length (m)'`, `a` is not returned.
631 x_max = self.info['x data (m)'].max()
632 a_given = 'Kuhn length (m)' in self.info
634 a = self.info['Kuhn length (m)']
637 f_max = self._data.max()
638 kwargs = self._kwargs()
640 x_section = inverse_FJC_PEG_fn(F_data=f_max, **kwargs)
641 N = x_max / x_section;
647 def guess_scale(self, params, outqueue=None):
648 """Guess parameter scales.
653 A guess at the section count scale in meters.
654 a_scale : float (optional)
655 A guess at the Kuhn length in meters. If the length of
656 `params` is less than 2, `a_scale` is not returned.
658 return [x/10.0 for x in params]
661 def WLC_fn(x_data, T, L, p):
662 """The worm like chain interpolation model.
667 x values for which the WLC tension should be calculated.
669 temperature in Kelvin.
671 contour length in meters.
673 persistence length in meters.
682 >>> WLC_fn(numpy.array([1e-9, 5e-9, 10e-9], dtype=numpy.float),
683 ... T=300, L=15e-9, p=2.5e-10) # doctest: +ELLIPSIS
684 array([ 1.717...e-12, 1.070...e-11, 4.418...e-11])
688 The function is the simple polynomial worm like chain
689 interpolation forumula proposed by C. Bustamante, et al. [#]_.
692 F(x) = \frac{k_B T}{p} \left[
693 \frac{1}{4}\left( \frac{1}{(1-x/L)^2} - 1 \right)
697 .. [#] C. Bustamante, J.F. Marko, E.D. Siggia, and S.Smith.
698 "Entropic elasticity of lambda-phage DNA."
700 doi: `10.1126/science.8079175 <http://dx.doi.org/10.1126/science.8079175>`
703 scaled_data = x_data / L
704 return a * (0.25*((1-scaled_data)**-2 - 1) + scaled_data)
706 class WLC (ModelFitter):
707 """Fit the data to a worm like chain.
711 Generate some example data
714 >>> L = 35e-9 # meters
715 >>> p = 2.5e-10 # meters
716 >>> x_data = numpy.linspace(10e-9, 30e-9, num=20)
717 >>> d_data = WLC_fn(x_data, T=T, L=L, p=p)
719 Fit the example data with a two-parameter fit (`L` and `p`).
722 ... 'temperature (K)': T,
723 ... 'x data (m)': x_data,
725 >>> model = WLC(d_data, info=info, rescale=True)
726 >>> outqueue = Queue.Queue()
727 >>> Lp,p = model.fit(outqueue=outqueue)
728 >>> fit_info = outqueue.get(block=False)
729 >>> model.L(Lp) # doctest: +ELLIPSIS
731 >>> p # doctest: +ELLIPSIS
734 Fit the example data with a one-parameter fit (`L`). We introduce
735 some error in our fixed persistence length for fun.
737 >>> info['persistence length (m)'] = 2*p
738 >>> model = WLC(d_data, info=info, rescale=True)
739 >>> Lp = model.fit(outqueue=outqueue)
740 >>> fit_info = outqueue.get(block=False)
741 >>> model.L(Lp) # doctest: +ELLIPSIS
745 """To avoid invalid values of `L`, we reparametrize with `Lp`.
750 contour length in meters
755 rescaled version of `L`.
759 >>> x_data = numpy.linspace(1e-9, 20e-9, num=10)
760 >>> model = WLC(data=x_data, info={'x data (m)':x_data})
761 >>> model.Lp(20e-9) # doctest: +ELLIPSIS
763 >>> model.Lp(19e-9) # doctest: +ELLIPSIS
765 >>> model.Lp(21e-9) # doctest: +ELLIPSIS
767 >>> model.Lp(100e-9) # doctest: +ELLIPSIS
772 The rescaling is designed to limit `L` to values strictly
773 greater than the maximum `x` data value, so we use
776 L = [\exp(L_p))+1] * x_\text{max}
778 which has the inverse
781 L_p = \ln(L/x_\text{max}-1)
783 This will obviously effect the interpretation of the fit's covariance matrix.
785 x_max = self.info['x data (m)'].max()
786 return numpy.log(L/x_max - 1)
789 """Inverse of :meth:`Lp`.
794 rescaled version of `L`
799 contour length in meters
803 >>> x_data = numpy.linspace(1e-9, 20e-9, num=10)
804 >>> model = WLC(data=x_data, info={'x data (m)':x_data})
805 >>> model.L(model.Lp(21e-9)) # doctest: +ELLIPSIS
807 >>> model.L(model.Lp(100e-9)) # doctest: +ELLIPSIS
810 x_max = self.info['x data (m)'].max()
811 return (numpy.exp(Lp)+1)*x_max
813 def model(self, params):
814 """Extract the relavant arguments and use :func:`WLC_fn`.
818 params : list of floats
819 `[Lp, p]`, where the presence of `p` is optional.
821 # setup convenient aliases
822 T = self.info['temperature (K)']
823 x_data = self.info['x data (m)']
824 L = self.L(params[0])
828 p = self.info['persistence length (m)']
830 self._model_data[:] = WLC_fn(x_data, T, L, p)
831 return self._model_data
833 def guess_initial_params(self, outqueue=None):
834 """Guess initial fitting parameters.
839 A guess at the reparameterized contour length in meters.
841 A guess at the persistence length in meters. If `.info`
842 has a setting for `'persistence length (m)'`, `p` is not
845 x_max = self.info['x data (m)'].max()
846 p_given = 'persistence length (m)' in self.info
848 p = self.info['persistence length (m)']
857 def guess_scale(self, params, outqueue=None):
858 """Guess parameter scales.
863 A guess at the reparameterized contour length scale in meters.
864 p_scale : float (optional)
865 A guess at the persistence length in meters. If the
866 length of `params` is less than 2, `p_scale` is not
872 return [Lp_scale] + [x/10.0 for x in params[1:]]
875 class PolymerFitPlugin (Plugin):
876 """Polymer model (WLC, FJC, etc.) fitting.
879 super(PolymerFitPlugin, self).__init__(name='polymer_fit')
880 self._commands = [PolymerFitCommand(self),]
882 def dependencies(self):
885 def default_settings(self):
887 Setting(section=self.setting_section, help=self.__doc__),
888 Setting(section=self.setting_section,
889 option='polymer model',
891 help="Select the default polymer model for 'polymer fit'. See the documentation for descriptions of available algorithms."),
892 Setting(section=self.setting_section,
893 option='FJC Kuhn length',
894 value=4e-10, type='float',
895 help='Kuhn length in meters'),
896 Setting(section=self.setting_section,
897 option='FJC-PEG Kuhn length',
898 value=4e-10, type='float',
899 help='Kuhn length in meters'),
900 Setting(section=self.setting_section,
901 option='FJC-PEG elasticity',
902 value=150.0, type='float',
903 help='Elasticity of a PEG segment in Newtons per meter.'),
904 Setting(section=self.setting_section,
905 option='FJC-PEG delta G',
906 value=3.0, type='float',
907 help='Gibbs free energy difference between trans-trans-trans (ttt) and trans-trans-gauche (ttg) PEG states in units of kBT.'),
908 Setting(section=self.setting_section,
909 option='FJC-PEG L_helical',
910 value=2.8e-10, type='float',
911 help='Contour length of PEG in the ttg state.'),
912 Setting(section=self.setting_section,
913 option='FJC-PEG L_planar',
914 value=3.58e-10, type='float',
915 help='Contour length of PEG in the ttt state.'),
916 Setting(section=self.setting_section,
917 option='WLC persistence length',
918 value=4e-10, type='float',
919 help='Persistence length in meters'),
923 class PolymerFitCommand (Command):
924 """Polymer model (WLC, FJC, etc.) fitting.
926 Fits an entropic elasticity function to a given chunk of the
927 curve. Fit quality compares the residual with the thermal noise
928 (a good fit should be 1 or less).
930 Because the models are complicated and varied, you should
931 configure the command by setting configuration options instead of
932 using command arguments. TODO: argument_to_setting()
934 def __init__(self, plugin):
935 super(PolymerFitCommand, self).__init__(
939 Argument(name='block', aliases=['set'], type='int', default=0,
941 Data block for which the fit should be calculated. For an
942 approach/retract force curve, `0` selects the approaching curve and
943 `1` selects the retracting curve.
945 Argument(name='bounds', type='point', optional=False, count=2,
947 Indicies of points bounding the selected data.
950 help=self.__doc__, plugin=plugin)
952 def _run(self, hooke, inqueue, outqueue, params):
953 scale(hooke, params['curve'], params['block']) # TODO: is autoscaling a good idea? (explicit is better than implicit)
954 data = params['curve'].data[params['block']]
955 # HACK? rely on params['curve'] being bound to the local hooke
956 # playlist (i.e. not a copy, as you would get by passing a
957 # curve through the queue). Ugh. Stupid queues. As an
958 # alternative, we could pass lookup information through the
960 model = self.plugin.config['polymer model']
961 new = Data((data.shape[0], data.shape[1]+1), dtype=data.dtype)
962 new.info = copy.deepcopy(data.info)
964 new.info['columns'].append('%s tension (N)' % model) # TODO: WLC fit for each peak, etc.
965 z_data = data[:,data.info['columns'].index(
966 'cantilever adjusted extension (m)')]
967 d_data = data[:,data.info['columns'].index('deflection (N)')]
968 start,stop = params['bounds']
969 tension_data,ps = self.fit_polymer_model(
970 params['curve'], z_data, d_data, start, stop, outqueue)
971 new.info['%s polymer fit parameters' % model] = ps
972 new[:,-1] = tension_data
973 params['curve'].data[params['block']] = new
975 def fit_polymer_model(self, curve, z_data, d_data, start, stop,
977 """Railyard for the `fit_*_model` family.
979 Uses the `polymer model` configuration setting to call the
980 appropriate backend algorithm.
982 fn = getattr(self, 'fit_%s_model'
983 % self.plugin.config['polymer model'].replace('-','_'))
984 return fn(curve, z_data, d_data, start, stop, outqueue)
986 def fit_FJC_model(self, curve, z_data, d_data, start, stop,
988 """Fit the data with :class:`FJC`.
991 'temperature (K)': self.plugin.config['temperature'],
992 'x data (m)': z_data[start:stop],
994 if True: # TODO: optionally free persistence length
995 info['Kuhn length (m)'] = (
996 self.plugin.config['FJC Kuhn length'])
997 model = FJC(d_data[start:stop], info=info, rescale=True)
998 queue = Queue.Queue()
999 params = model.fit(outqueue=queue)
1000 if True: # TODO: if Kuhn length fixed
1002 a = info['Kuhn length (m)']
1007 T = info['temperature (K)']
1008 fit_info = queue.get(block=False)
1009 mask = numpy.zeros(z_data.shape, dtype=numpy.bool)
1010 mask[start:stop] = True
1011 return [FJC_fn(z_data, T=T, L=L, a=a) * mask,
1014 def fit_FJC_PEG_model(self, curve, z_data, d_data, start, stop,
1016 """Fit the data with :class:`FJC_PEG`.
1019 'temperature (K)': self.plugin.config['temperature'],
1020 'x data (m)': z_data[start:stop],
1023 if True: # TODO: optionally free persistence length
1024 info['Kuhn length (m)'] = (
1025 self.plugin.config['FJC Kuhn length'])
1026 model = FJC(d_data[start:stop], info=info, rescale=True)
1027 queue = Queue.Queue()
1028 params = model.fit(outqueue=queue)
1029 if True: # TODO: if Kuhn length fixed
1031 a = info['Kuhn length (m)']
1036 T = info['temperature (K)']
1037 fit_info = queue.get(block=False)
1038 mask = numpy.zeros(z_data.shape, dtype=numpy.bool)
1039 mask[start:stop] = True
1040 return [FJC_PEG_fn(z_data, **kwargs) * mask,
1043 def fit_WLC_model(self, curve, z_data, d_data, start, stop,
1045 """Fit the data with :class:`WLC`.
1048 'temperature (K)': self.plugin.config['temperature'],
1049 'x data (m)': z_data[start:stop],
1051 if True: # TODO: optionally free persistence length
1052 info['persistence length (m)'] = (
1053 self.plugin.config['WLC persistence length'])
1054 model = WLC(d_data[start:stop], info=info, rescale=True)
1055 queue = Queue.Queue()
1056 params = model.fit(outqueue=queue)
1057 if True: # TODO: if persistence length fixed
1059 p = info['persistence length (m)']
1064 T = info['temperature (K)']
1065 fit_info = queue.get(block=False)
1066 mask = numpy.zeros(z_data.shape, dtype=numpy.bool)
1067 mask[start:stop] = True
1068 return [WLC_fn(z_data, T=T, L=L, p=p) * mask,
1073 # def dist2fit(self):
1074 # '''Calculates the average distance from data to fit, scaled by
1075 # the standard deviation of the free cantilever area (thermal