From ae09e7a6081e691c8cecf57fdb7a208c07fc663d Mon Sep 17 00:00:00 2001 From: "W. Trevor King" Date: Thu, 8 Mar 2012 10:35:22 -0500 Subject: [PATCH] Add `peak extraction method` option to `polymer fit` command. --- .update-copyright.conf | 2 +- hooke/plugin/polymer_fit.py | 57 +++++++++++++++- hooke/util/quickhull.py | 128 ++++++++++++++++++++++++++++++++++++ 3 files changed, 185 insertions(+), 2 deletions(-) create mode 100644 hooke/util/quickhull.py diff --git a/.update-copyright.conf b/.update-copyright.conf index 8238a4f..f064d46 100644 --- a/.update-copyright.conf +++ b/.update-copyright.conf @@ -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] diff --git a/hooke/plugin/polymer_fit.py b/hooke/plugin/polymer_fit.py index 45ead43..60a1284 100644 --- 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. diff --git a/hooke/util/quickhull.py b/hooke/util/quickhull.py new file mode 100644 index 0000000..1b54165 --- /dev/null +++ b/hooke/util/quickhull.py @@ -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) -- 2.26.2