# -*- coding: utf-8 -*-
#
-# Copyright (C) 2008-2010 Alberto Gomez-Casado
-# Fabrizio Benedetti
-# Massimo Sandal <devicerandom@gmail.com>
-# W. Trevor King <wking@drexel.edu>
+# Copyright (C) 2010-2012 W. Trevor King <wking@tremily.us>
#
# This file is part of Hooke.
#
-# Hooke is free software: you can redistribute it and/or modify it
-# under the terms of the GNU Lesser General Public License as
-# published by the Free Software Foundation, either version 3 of the
-# License, or (at your option) any later version.
+# Hooke is free software: you can redistribute it and/or modify it under the
+# terms of the GNU Lesser General Public License as published by the Free
+# Software Foundation, either version 3 of the License, or (at your option) any
+# later version.
#
-# Hooke is distributed in the hope that it will be useful, but WITHOUT
-# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
-# or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General
-# Public License for more details.
+# Hooke is distributed in the hope that it will be useful, but WITHOUT ANY
+# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
+# A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
+# details.
#
-# You should have received a copy of the GNU Lesser General Public
-# License along with Hooke. If not, see
-# <http://www.gnu.org/licenses/>.
+# You should have received a copy of the GNU Lesser General Public License
+# along with Hooke. If not, see <http://www.gnu.org/licenses/>.
"""The ``polymer_fit`` module proviews :class:`PolymerFitPlugin` and
serveral associated :class:`hooke.command.Command`\s for Fitting
import copy
from Queue import Queue
+import logging
import StringIO
import sys
from ..util.callback import is_iterable
from ..util.fit import PoorFit, ModelFitter
from ..util.peak import Peak
+from ..util.quickhull import qhull, points_inside_hull
from ..util.si import join_data_label, split_data_label
from . import Plugin, argument_to_setting
from .curve import ColumnAccessCommand, ColumnAddingCommand
>>> outqueue = Queue()
>>> L,a = model.fit(outqueue=outqueue)
>>> fit_info = outqueue.get(block=False)
- >>> print L
+ >>> print(L)
3.5e-08
- >>> print a
+ >>> print(a)
2.5e-10
Fit the example data with a one-parameter fit (`L`). We introduce
>>> info['Kuhn length (m)'] = 2*a
>>> model = FJC(d_data, info=info, rescale=True)
- >>> L = model.fit(outqueue=outqueue)
+ >>> L, = model.fit(outqueue=outqueue)
>>> fit_info = outqueue.get(block=False)
- >>> print L # doctest: +ELLIPSIS
+ >>> print(L) # doctest: +ELLIPSIS
3.199...e-08
"""
def Lp(self, L):
def fit(self, *args, **kwargs):
params = super(FJC, self).fit(*args, **kwargs)
- if is_iterable(params):
- params[0] = self.L(params[0]) # convert Lp -> L
+ params[0] = self.L(params[0]) # convert Lp -> L
+ if len(params) > 1:
params[1] = abs(params[1]) # take the absolute value of `a`
- else: # params is a float
- params = self.L(params) # convert Lp -> L
return params
def guess_initial_params(self, outqueue=None):
>>> outqueue = Queue()
>>> N,a = model.fit(outqueue=outqueue)
>>> fit_info = outqueue.get(block=False)
- >>> print N
+ >>> print(N)
123.0
- >>> print a
+ >>> print(a)
7e-10
Fit the example data with a one-parameter fit (`N`). We introduce
>>> info['Kuhn length (m)'] = 2*kwargs['a']
>>> model = FJC_PEG(d_data, info=info, rescale=True)
- >>> N = model.fit(outqueue=outqueue)
+ >>> N, = model.fit(outqueue=outqueue)
>>> fit_info = outqueue.get(block=False)
- >>> print N # doctest: +ELLIPSIS
+ >>> print(N) # doctest: +ELLIPSIS
96.931...
"""
def Lr(self, L):
def fit(self, *args, **kwargs):
params = super(FJC_PEG, self).fit(*args, **kwargs)
- if is_iterable(params):
- params[0] = self.L(params[0]) # convert Nr -> N
+ params[0] = self.L(params[0]) # convert Nr -> N
+ if len(params) > 1:
params[1] = abs(params[1]) # take the absolute value of `a`
- else: # params is a float
- params = self.L(params) # convert Nr -> N
return params
def guess_initial_params(self, outqueue=None):
>>> outqueue = Queue()
>>> L,p = model.fit(outqueue=outqueue)
>>> fit_info = outqueue.get(block=False)
- >>> print L
+ >>> print(L)
3.5e-08
- >>> print p
+ >>> print(p)
2.5e-10
Fit the example data with a one-parameter fit (`L`). We introduce
>>> info['persistence length (m)'] = 2*p
>>> model = WLC(d_data, info=info, rescale=True)
- >>> L = model.fit(outqueue=outqueue)
+ >>> L, = model.fit(outqueue=outqueue)
>>> fit_info = outqueue.get(block=False)
- >>> print L # doctest: +ELLIPSIS
+ >>> print(L) # doctest: +ELLIPSIS
3.318...e-08
"""
def Lp(self, L):
def fit(self, *args, **kwargs):
params = super(WLC, self).fit(*args, **kwargs)
- if is_iterable(params):
- params[0] = self.L(params[0]) # convert Lp -> L
+ params[0] = self.L(params[0]) # convert Lp -> L
+ if len(params) > 1:
params[1] = abs(params[1]) # take the absolute value of `p`
- else: # params is a float
- params = self.L(params) # convert Lp -> L
return params
def guess_initial_params(self, outqueue=None):
help="""
Indicies of points bounding the selected data.
""".strip()),
+ Argument(name='peak extraction method', default='convex hull',
+ help=(
+ 'Select the method used to extract the peak '
+ 'deflection from the fitted model. `convex hull`, '
+ '`peak data`, or `peak model`.')),
] + plugin._arguments,
help=self.__doc__, plugin=plugin)
def _run(self, hooke, inqueue, outqueue, params):
- params = self.__setup_params(hooke, params)
+ self._add_to_command_stack(params)
+ log = logging.getLogger('hooke')
+ params = self._setup_params(hooke, params)
data = self._block(hooke, params)
model = params['polymer model']
dist_data = self._get_column(
start,stop = params['bounds']
tension_data,ps = self.fit_polymer_model(
params, dist_data, def_data, start, stop, outqueue)
+ if params['peak extraction method'] == 'convex hull':
+ peak_def,hull = self._convex_hull_deflection(
+ distance_data=dist_data, deflection_data=def_data,
+ deflection_model=tension_data, start=start, stop=stop,
+ outqueue=outqueue)
+ elif params['peak extraction method'] == 'peak data':
+ peak_def = numpy.nanmax(def_data[start:stop])
+ elif params['peak extraction method'] == 'peak model':
+ peak_def = numpy.nanmax(tension_data[start:stop])
+ else:
+ raise ValueError(params['peak extraction method'])
+ log.debug('{}: {}'.format(params['peak extraction method'], peak_def))
+ ps['peak deflection'] = peak_def
+ ps['peak extraction method'] = params['peak extraction method']
+ ps['model'] = model
data.info[params['fit parameters info name']] = ps
- data.info[params['fit parameters info name']]['model'] = model
self._set_column(hooke=hooke, params=params,
column_name='output tension column',
values=tension_data)
- def __setup_params(self, hooke, params):
+ def _setup_params(self, hooke, params):
for key,value in params.items():
if value == None: # Use configured default value.
params[key] = self.plugin.config[key]
queue = Queue()
params = model.fit(outqueue=queue)
if True: # TODO: if Kuhn length fixed
- params = [params]
a = info['Kuhn length (m)']
else:
a = params[1]
queue = Queue()
params = model.fit(outqueue=queue)
if True: # TODO: if Kuhn length fixed
- params = [params]
a = info['Kuhn length (m)']
else:
a = params[1]
queue = Queue()
params = model.fit(outqueue=queue)
if True: # TODO: if persistence length fixed
- params = [params]
p = info['persistence length (m)']
else:
p = params[1]
f_data[start:stop] = WLC_fn(dist_data[start:stop], T=T, L=L, p=p)
return [f_data, info]
+ def _convex_hull_deflection(self, distance_data, deflection_data,
+ deflection_model, start, stop, outqueue=None):
+ """Return the last model deflection point inside the data hull.
+
+ If the `stop` point is a bit past it's ideal location, the
+ rapid rise of some polymer models can lead to artificially
+ high peak deflections if the largest distance value is plugged
+ into the polymer model directly. As a more robust estimation
+ of the peak deflection, we calculate the convex hull
+ surrounding the distance-deflection data and return the last
+ model deflection point that is inside that hull.
+ """
+ distance_data = distance_data[start:stop]
+ deflection_data = deflection_data[start:stop]
+ deflection_model = deflection_model[start:stop]
+ unfolds = []
+ hulls = []
+ data = numpy.array(
+ [[x,y] for x,y in zip(distance_data, deflection_data)])
+ model = numpy.array(
+ [[x,y] for x,y in zip(distance_data, deflection_model)])
+ hull = qhull(data)
+ inside = points_inside_hull(hull, model)
+ # distance data is not necessarily monatonic
+ index_inside = [j for j,i in enumerate(inside) if i is True]
+ distance_inside = [(distance_data[i],i) for i in index_inside]
+ if distance_inside:
+ peak_dist,peak_index = max(distance_inside)
+ unfold = deflection_model[peak_index]
+ else: # all model points are outside the hull?!
+ unfold = None
+ return (unfold, hull)
+
class PolymerFitPeaksCommand (ColumnAccessCommand):
"""Polymer model (WLC, FJC, etc.) fitting.
help=self.__doc__, plugin=plugin)
def _run(self, hooke, inqueue, outqueue, params):
+ self._add_to_command_stack(params)
data = self._block(hooke, params)
- fit_command = [c for c in hooke.commands
- if c.name=='polymer fit'][0]
+ fit_command = hooke.command_by_name['polymer fit']
inq = Queue()
outq = Queue()
curve = params['curve']
p['bounds'] = [peak.index, peak.post_index()]
p['output tension column'] = peak.name
p['fit parameters info name'] = peak.name
+ p['stack'] = False
fit_command.run(hooke, inq, outq, **p)
ret = outq.get()
if not isinstance(ret, Success):
help=self.__doc__, plugin=plugin)
def _run(self, hooke, inqueue, outqueue, params):
+ self._add_to_command_stack(params)
data = self._block(hooke, params)
dist_data = self._get_column(
hooke=hooke, params=params, column_name='distance column')
projects / hooke.git / blobdiff