From 4ed7e9602b43341002b3cf0cb97823f23accadef Mon Sep 17 00:00:00 2001 From: "W. Trevor King" Date: Tue, 10 Aug 2010 11:49:13 -0400 Subject: [PATCH] Added PolymerFitPeaksCommand to hooke.plugin.polymer_fit --- hooke/plugin/flatfilt.py | 7 +- hooke/plugin/polymer_fit.py | 181 +++++++++++++++++++++++------------- 2 files changed, 120 insertions(+), 68 deletions(-) diff --git a/hooke/plugin/flatfilt.py b/hooke/plugin/flatfilt.py index 757b6bb..a1857fb 100644 --- a/hooke/plugin/flatfilt.py +++ b/hooke/plugin/flatfilt.py @@ -33,7 +33,7 @@ See Also """ import copy -from multiprocessing import Queue +from Queue import Queue from numpy import diff from scipy.signal.signaltools import medfilt @@ -159,11 +159,10 @@ class FlatPeaksCommand (Command): deriv = diff(median) peaks = find_peaks(deriv, **_kwargs(params, find_peaks_arguments, argument_input_keys=True)) + d_name,d_unit = split_data_label(params['input deflection column']) for i,peak in enumerate(peaks): - peak.name = 'flat filter of %s' % ( - params['input deflection column']) + peak.name = 'flat filter peak %d of %s' % (i, d_name) peak.index += start_index - peak.count = i data.info[params['peak info name']] = peaks # Add a peak-masked deflection column. diff --git a/hooke/plugin/polymer_fit.py b/hooke/plugin/polymer_fit.py index 89962ba..66df3e6 100644 --- a/hooke/plugin/polymer_fit.py +++ b/hooke/plugin/polymer_fit.py @@ -27,7 +27,7 @@ velocity-clamp data to various polymer models (WLC, FJC, etc.). """ import copy -import Queue +from Queue import Queue import StringIO import sys @@ -37,7 +37,7 @@ from scipy.optimize import newton from ..command import Command, Argument, Failure from ..config import Setting from ..curve import Data -from ..plugin import Plugin +from ..plugin import Plugin, argument_to_setting from ..util.callback import is_iterable from ..util.fit import PoorFit, ModelFitter from ..util.si import join_data_label, split_data_label @@ -231,7 +231,7 @@ class FJC (ModelFitter): ... 'x data (m)': x_data, ... } >>> model = FJC(d_data, info=info, rescale=True) - >>> outqueue = Queue.Queue() + >>> outqueue = Queue() >>> Lp,a = model.fit(outqueue=outqueue) >>> fit_info = outqueue.get(block=False) >>> model.L(Lp) # doctest: +ELLIPSIS @@ -501,7 +501,7 @@ class FJC_PEG (ModelFitter): ... 'Gibbs free energy difference (Gp - Gh) (kBT)': kwargs['dG'], ... } >>> model = FJC_PEG(d_data, info=info, rescale=True) - >>> outqueue = Queue.Queue() + >>> outqueue = Queue() >>> Nr,a = model.fit(outqueue=outqueue) >>> fit_info = outqueue.get(block=False) >>> model.L(Nr) # doctest: +ELLIPSIS @@ -724,7 +724,7 @@ class WLC (ModelFitter): ... 'x data (m)': x_data, ... } >>> model = WLC(d_data, info=info, rescale=True) - >>> outqueue = Queue.Queue() + >>> outqueue = Queue() >>> Lp,p = model.fit(outqueue=outqueue) >>> fit_info = outqueue.get(block=False) >>> model.L(Lp) # doctest: +ELLIPSIS @@ -878,48 +878,56 @@ class PolymerFitPlugin (Plugin): """ def __init__(self): super(PolymerFitPlugin, self).__init__(name='polymer_fit') - self._commands = [PolymerFitCommand(self),] - - def dependencies(self): - return ['vclamp'] - - def default_settings(self): - return [ - Setting(section=self.setting_section, help=self.__doc__), - Setting(section=self.setting_section, - option='polymer model', - value='WLC', - help="Select the default polymer model for 'polymer fit'. See the documentation for descriptions of available algorithms."), - Setting(section=self.setting_section, - option='FJC Kuhn length', - value=4e-10, type='float', + self._arguments = [ # For Command initialization + Argument('polymer model', default='WLC', help=""" +Select the default polymer model for 'polymer fit'. See the +documentation for descriptions of available algorithms. +""".strip()), + Argument('FJC Kuhn length', type='float', default=4e-10, help='Kuhn length in meters'), - Setting(section=self.setting_section, - option='FJC-PEG Kuhn length', - value=4e-10, type='float', + Argument('FJC_PEG Kuhn length', type='float', default=4e-10, help='Kuhn length in meters'), - Setting(section=self.setting_section, - option='FJC-PEG elasticity', - value=150.0, type='float', + Argument('FJC_PEG elasticity', type='float', default=150.0, help='Elasticity of a PEG segment in Newtons per meter.'), - Setting(section=self.setting_section, - option='FJC-PEG delta G', - value=3.0, type='float', - help='Gibbs free energy difference between trans-trans-trans (ttt) and trans-trans-gauche (ttg) PEG states in units of kBT.'), - Setting(section=self.setting_section, - option='FJC-PEG L_helical', - value=2.8e-10, type='float', + Argument('FJC_PEG delta G', type='float', default=3.0, help=""" +Gibbs free energy difference between trans-trans-trans (ttt) and +trans-trans-gauche (ttg) PEG states in units of kBT. +""".strip()), + Argument('FJC_PEG L_helical', type='float', default=2.8e-10, help='Contour length of PEG in the ttg state.'), - Setting(section=self.setting_section, - option='FJC-PEG L_planar', - value=3.58e-10, type='float', + Argument('FJC_PEG L_planar', type='float', default=3.58e-10, help='Contour length of PEG in the ttt state.'), - Setting(section=self.setting_section, - option='WLC persistence length', - value=4e-10, type='float', - help='Persistence length in meters'), + Argument('WLC persistence length', type='float', default=4e-10, + help='Persistence length in meters'), + ] + self._settings = [ + Setting(section=self.setting_section, help=self.__doc__)] + for argument in self._arguments: + self._settings.append(argument_to_setting( + self.setting_section, argument)) + argument.default = None # if argument isn't given, use the config. + self._arguments.extend([ # Non-configurable arguments + Argument(name='input distance column', type='string', + default='cantilever adjusted extension (m)', + help=""" +Name of the column to use as the surface position input. +""".strip()), + Argument(name='input deflection column', type='string', + default='deflection (N)', + help=""" +Name of the column to use as the deflection input. +""".strip()), + ]) + self._commands = [ + PolymerFitCommand(self), PolymerFitPeaksCommand(self), ] + def dependencies(self): + return ['vclamp'] + + def default_settings(self): + return self._settings + class PolymerFitCommand (Command): """Polymer model (WLC, FJC, etc.) fitting. @@ -947,16 +955,7 @@ approach/retract force curve, `0` selects the approaching curve and help=""" Indicies of points bounding the selected data. """.strip()), - Argument(name='input distance column', type='string', - default='cantilever adjusted extension (m)', - help=""" -Name of the column to use as the surface position input. -""".strip()), - Argument(name='input deflection column', type='string', - default='deflection (N)', - help=""" -Name of the column to use as the deflection input. -""".strip()), + ] + plugin._arguments + [ Argument(name='output tension column', type='string', default='polymer tension', help=""" @@ -971,6 +970,10 @@ Name (without units) for storing the fit parameters in the `.info` dictionary. help=self.__doc__, plugin=plugin) def _run(self, hooke, inqueue, outqueue, params): + print params['curve'], id(params['curve']), 'PFC' + for key,value in params.items(): + if value == None: # Use configured default value. + params[key] = self.plugin.config[key] scale(hooke, params['curve'], params['block']) # TODO: is autoscaling a good idea? (explicit is better than implicit) data = params['curve'].data[params['block']] # HACK? rely on params['curve'] being bound to the local hooke @@ -978,25 +981,27 @@ Name (without units) for storing the fit parameters in the `.info` dictionary. # curve through the queue). Ugh. Stupid queues. As an # alternative, we could pass lookup information through the # queue... - model = self.plugin.config['polymer model'] + model = params['polymer model'] new = Data((data.shape[0], data.shape[1]+1), dtype=data.dtype) new.info = copy.deepcopy(data.info) new[:,:-1] = data new.info['columns'].append( join_data_label(params['output tension column'], 'N')) + print 'add column', id(params['curve']), new.info['columns'][-1], 'to', new.info['columns'][:-1] z_data = data[:,data.info['columns'].index( params['input distance column'])] d_data = data[:,data.info['columns'].index( params['input deflection column'])] start,stop = params['bounds'] tension_data,ps = self.fit_polymer_model( - params['curve'], z_data, d_data, start, stop, outqueue) + params, z_data, d_data, start, stop, outqueue) new.info[params['fit parameters info name']] = ps new.info[params['fit parameters info name']]['model'] = model + print 'add info', params['fit parameters info name'] new[:,-1] = tension_data params['curve'].data[params['block']] = new - def fit_polymer_model(self, curve, z_data, d_data, start, stop, + def fit_polymer_model(self, params, z_data, d_data, start, stop, outqueue=None): """Railyard for the `fit_*_model` family. @@ -1004,10 +1009,10 @@ Name (without units) for storing the fit parameters in the `.info` dictionary. appropriate backend algorithm. """ fn = getattr(self, 'fit_%s_model' - % self.plugin.config['polymer model'].replace('-','_')) - return fn(curve, z_data, d_data, start, stop, outqueue) + % params['polymer model'].replace('-','_')) + return fn(params, z_data, d_data, start, stop, outqueue) - def fit_FJC_model(self, curve, z_data, d_data, start, stop, + def fit_FJC_model(self, params, z_data, d_data, start, stop, outqueue=None): """Fit the data with :class:`FJC`. """ @@ -1017,9 +1022,9 @@ Name (without units) for storing the fit parameters in the `.info` dictionary. } if True: # TODO: optionally free persistence length info['Kuhn length (m)'] = ( - self.plugin.config['FJC Kuhn length']) + params['FJC Kuhn length']) model = FJC(d_data[start:stop], info=info, rescale=True) - queue = Queue.Queue() + queue = Queue() params = model.fit(outqueue=queue) if True: # TODO: if Kuhn length fixed params = [params] @@ -1035,7 +1040,7 @@ Name (without units) for storing the fit parameters in the `.info` dictionary. return [FJC_fn(z_data, T=T, L=L, a=a) * mask, fit_info] - def fit_FJC_PEG_model(self, curve, z_data, d_data, start, stop, + def fit_FJC_PEG_model(self, params, z_data, d_data, start, stop, outqueue=None): """Fit the data with :class:`FJC_PEG`. """ @@ -1046,9 +1051,9 @@ Name (without units) for storing the fit parameters in the `.info` dictionary. } if True: # TODO: optionally free persistence length info['Kuhn length (m)'] = ( - self.plugin.config['FJC Kuhn length']) + params['FJC Kuhn length']) model = FJC(d_data[start:stop], info=info, rescale=True) - queue = Queue.Queue() + queue = Queue() params = model.fit(outqueue=queue) if True: # TODO: if Kuhn length fixed params = [params] @@ -1064,7 +1069,7 @@ Name (without units) for storing the fit parameters in the `.info` dictionary. return [FJC_PEG_fn(z_data, **kwargs) * mask, fit_info] - def fit_WLC_model(self, curve, z_data, d_data, start, stop, + def fit_WLC_model(self, params, z_data, d_data, start, stop, outqueue=None): """Fit the data with :class:`WLC`. """ @@ -1074,9 +1079,9 @@ Name (without units) for storing the fit parameters in the `.info` dictionary. } if True: # TODO: optionally free persistence length info['persistence length (m)'] = ( - self.plugin.config['WLC persistence length']) + params['WLC persistence length']) model = WLC(d_data[start:stop], info=info, rescale=True) - queue = Queue.Queue() + queue = Queue() params = model.fit(outqueue=queue) if True: # TODO: if persistence length fixed params = [params] @@ -1093,6 +1098,54 @@ Name (without units) for storing the fit parameters in the `.info` dictionary. fit_info] + +class PolymerFitPeaksCommand (Command): + """Polymer model (WLC, FJC, etc.) fitting. + + Use :class:`PolymerFitCommand` to fit the each peak in a list of + previously determined peaks. + """ + def __init__(self, plugin): + super(PolymerFitPeaksCommand, self).__init__( + name='polymer fit peaks', + arguments=[ + CurveArgument, + Argument(name='block', aliases=['set'], type='int', default=0, + help=""" +Data block for which the fit should be calculated. For an +approach/retract force curve, `0` selects the approaching curve and +`1` selects the retracting curve. +""".strip()), + Argument(name='peak info name', type='string', + default='flat filter peaks', + help=""" +Name for storing the distance offset in the `.info` dictionary. +""".strip()), + Argument(name='peak index', type='int', count=-1, default=None, + help=""" +Index of the selected peak in the list of peaks. Use `None` to fit all peaks. +""".strip()), + ] + plugin._arguments, + help=self.__doc__, plugin=plugin) + + def _run(self, hooke, inqueue, outqueue, params): + print params['curve'], id(params['curve']), 'PFPC' + data = params['curve'].data[params['block']] + fit_command = [c for c in hooke.commands + if c.name=='polymer fit'][0] + inq = Queue() + outq = Queue() + p = copy.deepcopy(params) + p['curve'] = params['curve'] + del(p['peak info name']) + del(p['peak index']) + for i,peak in enumerate(data.info[params['peak info name']]): + if params['peak index'] == None or i in params['peak index']: + p['bounds'] = [peak.index, peak.post_index()] + p['output tension column'] = peak.name + p['fit parameters info name'] = peak.name + fit_command.run(hooke, inq, outq, **p) + # TODO: # def dist2fit(self): # '''Calculates the average distance from data to fit, scaled by -- 2.26.2