Added PolymerFitPeaksCommand to hooke.plugin.polymer_fit
authorW. Trevor King <wking@drexel.edu>
Tue, 10 Aug 2010 15:49:13 +0000 (11:49 -0400)
committerW. Trevor King <wking@drexel.edu>
Tue, 10 Aug 2010 15:49:13 +0000 (11:49 -0400)
hooke/plugin/flatfilt.py
hooke/plugin/polymer_fit.py

index 757b6bba940838aa7ca63f8ee3866cd8c9919768..a1857fbfd66f648a6ac370c0102fd98149ac11b8 100644 (file)
@@ -33,7 +33,7 @@ See Also
 """
 
 import copy
 """
 
 import copy
-from multiprocessing import Queue
+from Queue import Queue
 
 from numpy import diff
 from scipy.signal.signaltools import medfilt
 
 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))
         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):
         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.index += start_index
-            peak.count = i
         data.info[params['peak info name']] = peaks
 
         # Add a peak-masked deflection column.
         data.info[params['peak info name']] = peaks
 
         # Add a peak-masked deflection column.
index 89962ba0a6091b1e853d00f8f7f7ee437976a3b3..66df3e64a125fd9f67753589ffc93bd7f56ccb93 100644 (file)
@@ -27,7 +27,7 @@ velocity-clamp data to various polymer models (WLC, FJC, etc.).
 """
 
 import copy
 """
 
 import copy
-import Queue
+from Queue import Queue
 import StringIO
 import sys
 
 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 ..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
 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)
     ...     '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
     >>> 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)
     ...     '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
     >>> 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)
     ...     '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
     >>> 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')
     """
     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'),
                     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'),
                     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.'),
                     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.'),
                     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.'),
                     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.
 
 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()),
                          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="""
                 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):
             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
         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...
         # 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'))
         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(
         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
         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
 
         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.
 
                           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'
         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`.
         """
                       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)'] = (
             }
         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)
         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]
         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]
 
         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`.
         """
                           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)'] = (
             }
         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)
         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]
         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]
 
         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`.
         """
                       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)'] = (
             }
         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)
         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]
         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]
 
 
                 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
 # TODO:
 # def dist2fit(self):
 #     '''Calculates the average distance from data to fit, scaled by