Add `peak extraction method` option to `polymer fit` command.
authorW. Trevor King <wking@drexel.edu>
Thu, 8 Mar 2012 15:35:22 +0000 (10:35 -0500)
committerW. Trevor King <wking@drexel.edu>
Thu, 8 Mar 2012 15:35:22 +0000 (10:35 -0500)
.update-copyright.conf
hooke/plugin/polymer_fit.py
hooke/util/quickhull.py [new file with mode: 0644]

index 8238a4f..f064d46 100644 (file)
@@ -5,7 +5,7 @@ vcs: Git
 [files]
 authors: yes
 files: yes
-ignored: COPYING, README, .update-copyright.conf, .git*, conf/*, test/data/*, doc/*
+ignored: COPYING, README, .update-copyright.conf, .git*, */quickhull.py, conf/*, test/data/*, doc/*
 pyfile: hooke/license.py
 
 [author-hacks]
index 45ead43..60a1284 100644 (file)
@@ -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.
diff --git a/hooke/util/quickhull.py b/hooke/util/quickhull.py
new file mode 100644 (file)
index 0000000..1b54165
--- /dev/null
@@ -0,0 +1,128 @@
+# Copyright (c) 2009 W. Trevor King and the authors listed at the
+# following URL, and/or the authors of referenced articles or
+# incorporated external code:
+# http://en.literateprograms.org/Quickhull_(Python,_arrays)?action=history&offset=20090125100840
+# 
+# Permission is hereby granted, free of charge, to any person obtaining
+# a copy of this software and associated documentation files (the
+# "Software"), to deal in the Software without restriction, including
+# without limitation the rights to use, copy, modify, merge, publish,
+# distribute, sublicense, and/or sell copies of the Software, and to
+# permit persons to whom the Software is furnished to do so, subject to
+# the following conditions:
+# 
+# The above copyright notice and this permission notice shall be
+# included in all copies or substantial portions of the Software.
+# 
+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+# EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+# MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
+# IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
+# CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
+# TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
+# SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+# 
+# Retrieved from: http://en.literateprograms.org/Quickhull_(Python,_arrays)?oldid=16036
+# Modified by W. Trevor King to add Simplex and point_inside_hull
+
+from numpy import *
+
+link = lambda a,b: concatenate((a,b[1:]))
+edge = lambda a,b: concatenate(([a],[b]))
+
+def qhull(sample):
+    def dome(sample,base): 
+        h, t = base
+        dists = dot(sample-h, dot(((0,-1),(1,0)),(t-h)))
+        outer = repeat(sample, dists>0, 0)
+
+        if len(outer):
+            pivot = sample[argmax(dists)]
+            return link(dome(outer, edge(h, pivot)),
+                        dome(outer, edge(pivot, t)))
+        else:
+            return base
+
+    if len(sample) > 2:
+       axis = sample[:,0]
+       base = take(sample, [argmin(axis), argmax(axis)], 0)
+       return link(dome(sample, base),
+                    dome(sample, base[::-1]))
+    else:
+       return sample
+
+class Simplex(object): # N+1 points in N dimensions
+    def __init__(self, vertices):
+        self.vertices = vertices
+        self.origin = self.vertices[0]
+        self.basis = array([v-self.origin for v in self.vertices[1:]],
+                           dtype=double).transpose()
+        self.inverse_basis = linalg.inv(self.basis)
+    def simplex_coordinates(self, points):
+        rel_coords = (points-self.origin).transpose()
+        return dot(self.inverse_basis, rel_coords).transpose()
+    def inside_single_point_simplex_coords(self, point, tol=1e-15):
+        # check that we're in the positive coordinate
+        for coordinate in point:
+            if coordinate < -tol:
+                return False
+        # check that we're under the (1,1,1,...) plane
+        if point.sum() > 1+tol:
+            return False
+        return True
+    def inside(self, points):
+        ps = self.simplex_coordinates(points)
+        if not hasattr(ps[0], "__len__"): # single point
+            return self.inside_single_point_simplex_coords(ps)
+        else:
+            insides = []
+            for p in ps:
+                insides.append(self.inside_single_point_simplex_coords(p))
+            return insides
+
+def points_inside_hull(hull, points):
+    if not hasattr(points[0], "__len__"): # single point
+        points = numpy.array([points])
+    inside = [False]*len(points)
+    pivot = hull[0]
+    for a,b in zip(hull[1:], hull[2:]):
+        rb = b-pivot
+        if vdot(rb, rb) == 0: # hull is a closed loop
+            continue
+        simplex = Simplex([pivot,a,b])
+        simplex_inside = simplex.inside(points)
+        if any(simplex_inside):
+            for i,value in enumerate(simplex_inside):
+                if value == True:
+                    inside[i] = True
+    return inside
+
+def print_postscript(sample, hull):
+    print "%!"
+    print "100 500 translate 2 2 scale 0 0 moveto"
+    print "/tick {moveto 0 2 rlineto 0 -4 rlineto 0 2 rlineto"
+    print "              2 0 rlineto -4 0 rlineto 2 0 rlineto} def"
+    for (x,y) in sample:
+       print x, y, "tick"
+    print "stroke"
+    print hull[0,0], hull[0,1], "moveto"
+    for (x,y) in hull[1:]:
+       print x, y, "lineto"
+    print "closepath stroke showpage"
+
+if __name__ == "__main__":
+    #sample = 10*array([(x,y) for x in arange(10) for y in arange(10)])
+    sample = 100*random.random((32,2))
+    hull = qhull(sample)
+    if all(point_inside_hull(hull, sample)) != True:
+        for point in sample:
+            assert point_inside_hull(hull, point) == True, \
+                "point %s should be in hull\n%s" % (sample, hull)
+    apoints = hull[0:-1]
+    bpoints = hull[1:]
+    if any(point_inside_hull(hull, apoints - 0.1*(bpoints-apoints))) == True:
+        for a,b in zip(apoints, bpoints):
+            assert point_inside_hull(hull, a - 0.1*(b-a)) == False, \
+                "point %s (a=%s, b=%s) should be outside hull\n%s" \
+                % (point, a, b, hull)
+    print_postscript(sample, hull)