Add `peak extraction method` option to `polymer fit` command.

[hooke.git] / hooke / plugin / polymer_fit.py

diff --git a/hooke/plugin/polymer_fit.py b/hooke/plugin/polymer_fit.py
index 45ead438c884bdff797e937f4a391ac335e6bebc..60a1284e28d8f1101dcff12d0825eb68393dfbfd 100644 (file)
--- a/hooke/plugin/polymer_fit.py
+++ b/hooke/plugin/polymer_fit.py
@@ -24,6 +24,7 @@ velocity-clamp data to various polymer models (WLC, FJC, etc.).
 
 import copy
 from Queue import Queue
+import logging
 import StringIO
 import sys
 
@@ -36,6 +37,7 @@ from ..curve import Data
 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
@@ -930,10 +932,16 @@ Name (without units) for storing the fit parameters in the `.info` dictionary.
                          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):
+        log = logging.getLogger('hooke')
         params = self._setup_params(hooke, params)
         data = self._block(hooke, params)
         model = params['polymer model']
@@ -944,8 +952,22 @@ Indicies of points bounding the selected data.
         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)
@@ -1050,6 +1072,39 @@ Indicies of points bounding the selected data.
         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.