Remove the old multiplier plotmanip from the vclamp plugin.
[hooke.git] / hooke / plugin / vclamp.py
index 2ab03727a7376387081c9c5d37a4d771ff90dd6f..fbd6335aec334c78950ad04a19361c1f32aad65b 100644 (file)
@@ -6,15 +6,15 @@
 #
 # This file is part of Hooke.
 #
-# Hooke is free software: you can redistribute it and/or
-# modify it under the terms of the GNU Lesser General Public
-# License as published by the Free Software Foundation, either
-# version 3 of the License, or (at your option) any later version.
+# Hooke is free software: you can redistribute it and/or modify it
+# under the terms of the GNU Lesser General Public License as
+# published by the Free Software Foundation, either version 3 of the
+# License, or (at your option) any later version.
 #
-# Hooke is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-# GNU Lesser General Public License for more details.
+# Hooke is distributed in the hope that it will be useful, but WITHOUT
+# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
+# or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General
+# Public License for more details.
 #
 # You should have received a copy of the GNU Lesser General Public
 # License along with Hooke.  If not, see
@@ -28,30 +28,55 @@ common velocity clamp analysis tasks.
 import copy
 
 import numpy
+import scipy
 
 from ..command import Command, Argument, Failure, NullQueue
 from ..config import Setting
 from ..curve import Data
-from ..plugin import Builtin
+from ..plugin import Plugin
 from ..util.fit import PoorFit, ModelFitter
+from ..util.si import join_data_label, split_data_label
 from .curve import CurveArgument
 
 
-def scale(hooke, curve):
+def scale(hooke, curve, block=None):
+    """Run 'add block force array' on `block`.
+
+    Runs 'zero block surface contact point' first, if necessary.  Does
+    not run either command if the columns they add to the block are
+    already present.
+
+    If `block` is `None`, scale all blocks in `curve`.
+    """
     commands = hooke.commands
     contact = [c for c in hooke.commands
                if c.name == 'zero block surface contact point'][0]
     force = [c for c in hooke.commands if c.name == 'add block force array'][0]
+    cant_adjust = [c for c in hooke.commands
+               if c.name == 'add block cantilever adjusted extension array'][0]
     inqueue = None
     outqueue = NullQueue()
-    for i,block in enumerate(curve.data):
-        params = {'curve':curve, 'block':i}
-        contact._run(hooke, inqueue, outqueue, params)
-        force._run(hooke, inqueue, outqueue, params)
+    if block == None:
+        for i in range(len(curve.data)):
+            scale(hooke, curve, block=i)
+    else:
+        params = {'curve':curve, 'block':block}
+        b = curve.data[block]
+        if ('surface distance (m)' not in b.info['columns']
+            or 'surface deflection (m)' not in b.info['columns']):
+            try:
+                contact.run(hooke, inqueue, outqueue, **params)
+            except PoorFit, e:
+                raise PoorFit('Could not fit %s %s: %s'
+                              % (curve.path, block, str(e)))
+        if ('deflection (N)' not in b.info['columns']):
+            force.run(hooke, inqueue, outqueue, **params)
+        if ('cantilever adjusted extension (m)' not in b.info['columns']):
+            cant_adjust.run(hooke, inqueue, outqueue, **params)
     return curve
 
 class SurfacePositionModel (ModelFitter):
-    """
+    """Bilinear surface position model.
 
     The bilinear model is symmetric, but the parameter guessing and
     sanity checks assume the contact region occurs for lower indicies
@@ -72,9 +97,6 @@ class SurfacePositionModel (ModelFitter):
     should be enough data in the off-surface region that interactions
     due to proteins, etc. will not seriously skew the fit in the
     off-surface region.
-
-
-    We guess
     """
     def model(self, params):
         """A continuous, bilinear model.
@@ -96,11 +118,12 @@ class SurfacePositionModel (ModelFitter):
         if self.info.get('force zero non-contact slope', None) == True:
             p = list(p)
             p.append(0.)  # restore the non-contact slope parameter
-        r2 = numpy.round(p[2])
+        r2 = numpy.round(abs(p[2]))
         if r2 >= 1:
             self._model_data[:r2] = p[0] + p[1] * numpy.arange(r2)
         if r2 < len(self._data)-1:
-            self._model_data[r2:] = p[0] + p[1]*p[2] + p[3] * numpy.arange(len(self._data)-r2)
+            self._model_data[r2:] = \
+                p[0] + p[1]*p[2] + p[3] * numpy.arange(len(self._data)-r2)
         return self._model_data
 
     def set_data(self, data, info=None):
@@ -126,9 +149,9 @@ class SurfacePositionModel (ModelFitter):
         produces the maximum deflection at the left-most point, and the
         final (non-contact) slope (:math:`p_3`) is zero.
         """
-        left_offset = self.info['max deflection']
-        left_slope = 2*(-self.info['deflection range']
-                         /self.info['position range'])
+        left_offset = self.info['min deflection']
+        left_slope = 2*(self.info['deflection range']
+                        /self.info['position range'])
         kink_position = (self.info['max position']
                          +self.info['min position'])/2.0
         right_slope = 0
@@ -143,7 +166,6 @@ class SurfacePositionModel (ModelFitter):
 
         Notes
         -----
-
         We guess offset scale (:math:`p_0`) as one tenth of the total
         deflection range, the kink scale (:math:`p_2`) as one tenth of
         the total index range, the initial (contact) slope scale
@@ -163,16 +185,17 @@ class SurfacePositionModel (ModelFitter):
     def fit(self, *args, **kwargs):
         self.info['guessed contact slope'] = None
         params = super(SurfacePositionModel, self).fit(*args, **kwargs)
+        params[2] = abs(params[2])
         if self.info.get('force zero non-contact slope', None) == True:
             params = list(params)
             params.append(0.)  # restore the non-contact slope parameter
 
         # check that the fit is reasonable, see the :meth:`model` docstring
-        # for parameter descriptions.
+        # for parameter descriptions.  HACK: hardcoded cutoffs.
         if abs(params[3]*10) > abs(params[1]) :
             raise PoorFit('Slope in non-contact region, or no slope in contact')
         if params[2] < self.info['min position']+0.02*self.info['position range']:
-            raise poorFit(
+            raise PoorFit(
                 'No kink (kink %g less than %g, need more space to left)'
                 % (params[2],
                    self.info['min position']+0.02*self.info['position range']))
@@ -182,15 +205,17 @@ class SurfacePositionModel (ModelFitter):
                 % (params[2],
                    self.info['max position']-0.02*self.info['position range']))
         if (self.info['guessed contact slope'] != None
-            and params[3] < 0.5 * self.info['guessed contact slope']):
-            raise PoorFit('Too far')
+            and abs(params[1]) < 0.5 * abs(self.info['guessed contact slope'])):
+            raise PoorFit('Too far (contact slope %g, but expected ~%g'
+                          % (params[3], self.info['guessed contact slope']))
         return params
 
-class VelocityClampPlugin (Builtin):
+class VelocityClampPlugin (Plugin):
     def __init__(self):
         super(VelocityClampPlugin, self).__init__(name='vclamp')
         self._commands = [
             SurfaceContactCommand(self), ForceCommand(self),
+            CantileverAdjustedExtensionCommand(self),
             ]
 
     def default_settings(self):
@@ -206,13 +231,13 @@ class VelocityClampPlugin (Builtin):
 class SurfaceContactCommand (Command):
     """Automatically determine a block's surface contact point.
 
-    Uses the block's `z piezo (m)` and `deflection (m)` arrays.
-    Stores the contact parameters in `block.info`'s `surface distance
-    offset (m)` and `surface deflection offset (m)`.  Model-specific
-    fitting information is stored in `surface detection parameters`.
+    You can select the contact point algorithm with the creatively
+    named `surface contact point algorithm` configuration setting.
+    Currently available options are:
 
-    The adjusted data columns `surface distance (m)` and `surface
-    adjusted deflection (m)` are also added to the block.
+    * fmms (:meth:`find_contact_point_fmms`)
+    * ms (:meth:`find_contact_point_ms`)
+    * wtk (:meth:`find_contact_point_wtk`)
     """
     def __init__(self, plugin):
         super(SurfaceContactCommand, self).__init__(
@@ -224,12 +249,47 @@ class SurfaceContactCommand (Command):
 Data block for which the force should be calculated.  For an
 approach/retract force curve, `0` selects the approaching curve and `1`
 selects the retracting curve.
+""".strip()),
+                Argument(name='input distance column', type='string',
+                         default='z piezo (m)',
+                         help="""
+Name of the column to use as the surface positioning input.
+""".strip()),
+                Argument(name='input deflection column', type='string',
+                         default='deflection (m)',
+                         help="""
+Name of the column to use as the deflection input.
+""".strip()),
+                Argument(name='output distance column', type='string',
+                         default='surface distance',
+                         help="""
+Name of the column (without units) to use as the surface positioning output.
+""".strip()),
+                Argument(name='output deflection column', type='string',
+                         default='surface deflection',
+                         help="""
+Name of the column (without units) to use as the deflection output.
+""".strip()),
+                Argument(name='distance info name', type='string',
+                         default='surface distance offset',
+                         help="""
+Name (without units) for storing the distance offset in the `.info` dictionary.
+""".strip()),
+                Argument(name='deflection info name', type='string',
+                         default='surface deflection offset',
+                         help="""
+Name (without units) for storing the deflection offset in the `.info` dictionary.
+""".strip()),
+                Argument(name='fit parameters info name', type='string',
+                         default='surface deflection offset',
+                         help="""
+Name (without units) for storing the deflection offset in the `.info` dictionary.
 """.strip()),
                 ],
             help=self.__doc__, plugin=plugin)
 
     def _run(self, hooke, inqueue, outqueue, params):
-        data = params['curve'].data[int(params['block'])] # HACK, int() should be handled by ui
+        data = params['curve'].data[params['block']]
         # HACK? rely on params['curve'] being bound to the local hooke
         # playlist (i.e. not a copy, as you would get by passing a
         # curve through the queue).  Ugh.  Stupid queues.  As an
@@ -238,103 +298,100 @@ selects the retracting curve.
         new = Data((data.shape[0], data.shape[1]+2), dtype=data.dtype)
         new.info = copy.deepcopy(data.info)
         new[:,:-2] = data
-        new.info['columns'].extend(
-            ['surface distance (m)', 'surface adjusted deflection (m)'])
-        z_data = data[:,data.info['columns'].index('z piezo (m)')]
-        d_data = data[:,data.info['columns'].index('deflection (m)')]
-        i,deflection_offset,ps = self._find_contact_point(
-            z_data, d_data, outqueue)
-        surface_offset = z_data[i]
-        new.info['surface distance offset (m)'] = surface_offset
-        new.info['surface deflection offset (m)'] = deflection_offset
-        new.info['surface detection parameters'] = ps
-        new[:,-2] = z_data - surface_offset
-        new[:,-1] = d_data - deflection_offset
-        data = params['curve'].data[int(params['block'])] # HACK, int() should be handled by ui
-        params['curve'].data[int(params['block'])] = new # HACK, int() should be handled by ui
-
-    def _find_contact_point(self, z_data, d_data, outqueue=None):
-        fn = getattr(self, '_find_contact_point_%s'
+        name,dist_units = split_data_label(params['input distance column'])
+        name,def_units = split_data_label(params['input deflection column'])
+        new.info['columns'].extend([
+                join_data_label(params['output distance column'], dist_units),
+                join_data_label(params['output deflection column'], def_units),
+                ])
+        dist_data = data[:,data.info['columns'].index(
+                params['input distance column'])]
+        def_data = data[:,data.info['columns'].index(
+                params['input deflection column'])]
+        i,def_offset,ps = self.find_contact_point(
+            params['curve'], dist_data, def_data, outqueue)
+        dist_offset = dist_data[i]
+        new.info[join_data_label(params['distance info name'], dist_units
+                                 )] = dist_offset
+        new.info[join_data_label(params['deflection info name'], def_units
+                                 )] = def_offset
+        new.info[params['fit parameters info name']] = ps
+        new[:,-2] = dist_data - dist_offset
+        new[:,-1] = def_data - def_offset
+        params['curve'].data[params['block']] = new
+
+    def find_contact_point(self, curve, z_data, d_data, outqueue=None):
+        """Railyard for the `find_contact_point_*` family.
+
+        Uses the `surface contact point algorithm` configuration
+        setting to call the appropriate backend algorithm.
+        """
+        fn = getattr(self, 'find_contact_point_%s'
                      % self.plugin.config['surface contact point algorithm'])
-        return fn(z_data, d_data, outqueue)
+        return fn(curve, z_data, d_data, outqueue)
 
-    def _find_contact_point_fmms(self, z_data, d_data, outqueue=None):
+    def find_contact_point_fmms(self, curve, z_data, d_data, outqueue=None):
         """Algorithm by Francesco Musiani and Massimo Sandal.
 
         Notes
         -----
         Algorithm:
 
-        * take care of the PicoForce trigger bug - exclude retraction portions with too high standard deviation
-        * fit the second half of the retraction curve to a line
-        * if the fit is not almost horizontal, take a smaller chunk and repeat
-        * otherwise, we have something horizontal
-        * so take the average of horizontal points and use it as a baseline
-
-        Then, start from the rise of the retraction curve and look at
-        the first point below the baseline.
+        0) Driver-specific workarounds, e.g. deal with the PicoForce
+          trigger bug by excluding retraction portions with excessive
+          deviation.
+        1) Select the second half (non-contact side) of the retraction
+          curve.
+        2) Fit the selection to a line.
+        3) If the fit is not almost horizontal, halve the selection
+          and retrun to (2).
+        4) Average the selection and use it as a baseline.
+        5) Slide in from the start (contact side) of the retraction
+        curve, until you find a point with greater than baseline
+        deflection.  That point is the contact point.
         """
-        if not plot:
-            plot=self.plots[0]
-
-        outplot=self.subtract_curves(1)
-        xret=outplot.vectors[1][0]
-        ydiff=outplot.vectors[1][1]
-
-        xext=plot.vectors[0][0]
-        yext=plot.vectors[0][1]
-        xret2=plot.vectors[1][0]
-        yret=plot.vectors[1][1]
-
-        #taking care of the picoforce trigger bug: we exclude portions of the curve that have too much
-        #standard deviation. yes, a lot of magic is here.
-        monster=True
-        monlength=len(xret)-int(len(xret)/20)
-        finalength=len(xret)
-        while monster:
-            monchunk=scipy.array(ydiff[monlength:finalength])
-            if abs(np.std(monchunk)) < 2e-10:
-                monster=False
-            else: #move away from the monster
-                monlength-=int(len(xret)/50)
-                finalength-=int(len(xret)/50)
-
-
-        #take half of the thing
-        endlength=int(len(xret)/2)
-
-        ok=False
-
-        while not ok:
-            xchunk=yext[endlength:monlength]
-            ychunk=yext[endlength:monlength]
-            regr=scipy.stats.linregress(xchunk,ychunk)[0:2]
-            #we stop if we found an almost-horizontal fit or if we're going too short...
-            #FIXME: 0.1 and 6 here are "magic numbers" (although reasonable)
-            if (abs(regr[1]) > 0.1) and ( endlength < len(xret)-int(len(xret)/6) ) :
-                endlength+=10
-            else:
-                ok=True
-
-
-        ymean=np.mean(ychunk) #baseline
-
-        index=0
-        point = ymean+1
-
-        #find the first point below the calculated baseline
-        while point > ymean:
-            try:
-                point=yret[index]
-                index+=1
-            except IndexError:
-                #The algorithm didn't find anything below the baseline! It should NEVER happen
-                index=0
-                return index
-
-        return index
-
-    def _find_contact_point_ms(self, z_data, d_data, outqueue=None):
+        if curve.info['filetype'] == 'picoforce':
+            # Take care of the picoforce trigger bug (TODO: example
+            # data file demonstrating the bug).  We exclude portions
+            # of the curve that have too much standard deviation.
+            # Yes, a lot of magic is here.
+            check_start = len(d_data)-len(d_data)/20
+            monster_start = len(d_data)
+            while True:
+                # look at the non-contact tail
+                non_monster = d_data[check_start:monster_start]
+                if non_monster.std() < 2e-10: # HACK: hardcoded cutoff
+                    break
+                else: # move further away from the monster
+                    check_start -= len(d_data)/50
+                    monster_start -= len(d_data)/50
+            z_data = z_data[:monster_start]
+            d_data = d_data[:monster_start]
+
+        # take half of the thing to start
+        selection_start = len(d_data)/2
+        while True:
+            z_chunk = z_data[selection_start:]
+            d_chunk = d_data[selection_start:]
+            slope,intercept,r,two_tailed_prob,stderr_of_the_estimate = \
+                scipy.stats.linregress(z_chunk, d_chunk)
+            # We stop if we found an almost-horizontal fit or if we're
+            # getting to small a selection.  FIXME: 0.1 and 5./6 here
+            # are "magic numbers" (although reasonable)
+            if (abs(slope) < 0.1  # deflection (m) / surface (m)
+                or selection_start > 5./6*len(d_data)):
+                break
+            selection_start += 10
+
+        d_baseline = d_chunk.mean()
+
+        # find the first point above the calculated baseline
+        i = 0
+        while i < len(d_data) and d_data[i] < ymean:
+            i += 1
+        return (i, d_baseline, {})
+
+    def find_contact_point_ms(self, curve, z_data, d_data, outqueue=None):
         """Algorithm by Massimo Sandal.
 
         Notes
@@ -346,13 +403,6 @@ selects the retracting curve.
 
         * ?
         """
-        #raw_plot=self.current.curve.default_plots()[0]
-        raw_plot=self.plots[0]
-        '''xext=self.plots[0].vectors[0][0]
-        yext=self.plots[0].vectors[0][1]
-        xret2=self.plots[0].vectors[1][0]
-        yret=self.plots[0].vectors[1][1]
-        '''
         xext=raw_plot.vectors[0][0]
         yext=raw_plot.vectors[0][1]
         xret2=raw_plot.vectors[1][0]
@@ -386,7 +436,7 @@ selects the retracting curve.
         else:
             return dummy.index
 
-    def _find_contact_point_wtk(self, z_data, d_data, outqueue=None):
+    def find_contact_point_wtk(self, curve, z_data, d_data, outqueue=None):
         """Algorithm by W. Trevor King.
 
         Notes
@@ -412,11 +462,9 @@ selects the retracting curve.
             surface_index = len(d_data)-1-surface_index
         return (numpy.round(surface_index), deflection_offset, info)
 
-class ForceCommand (Command):
-    """Calculate a block's `deflection (N)` array.
 
-    Uses the block's `deflection (m)` array and `spring constant
-    (N/m)`.
+class ForceCommand (Command):
+    """Convert a deflection column from meters to newtons.
     """
     def __init__(self, plugin):
         super(ForceCommand, self).__init__(
@@ -428,12 +476,27 @@ class ForceCommand (Command):
 Data block for which the force should be calculated.  For an
 approach/retract force curve, `0` selects the approaching curve and `1`
 selects the retracting curve.
+""".strip()),
+                Argument(name='input deflection column', type='string',
+                         default='surface deflection (m)',
+                         help="""
+Name of the column to use as the deflection input.
+""".strip()),
+                Argument(name='output deflection column', type='string',
+                         default='deflection',
+                         help="""
+Name of the column (without units) to use as the deflection output.
+""".strip()),
+                Argument(name='spring constant info name', type='string',
+                         default='spring constant (N/m)',
+                         help="""
+Name of the spring constant in the `.info` dictionary.
 """.strip()),
                 ],
             help=self.__doc__, plugin=plugin)
 
     def _run(self, hooke, inqueue, outqueue, params):
-        data = params['curve'].data[int(params['block'])] # HACK, int() should be handled by ui
+        data = params['curve'].data[params['block']]
         # HACK? rely on params['curve'] being bound to the local hooke
         # playlist (i.e. not a copy, as you would get by passing a
         # curve through the queue).  Ugh.  Stupid queues.  As an
@@ -442,159 +505,82 @@ selects the retracting curve.
         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('deflection (N)')
-        d_data = data[:,data.info['columns'].index('surface adjusted deflection (m)')]
-        new[:,-1] = d_data * data.info['spring constant (N/m)']
-        params['curve'].data[int(params['block'])] = new # HACK, int() should be handled by ui
-
-
-class generalvclampCommands(object):
-
-    def do_subtplot(self, args):
-        '''
-        SUBTPLOT
-        (procplots.py plugin)
-        Plots the difference between ret and ext current curve
-        -------
-        Syntax: subtplot
-        '''
-        #FIXME: sub_filter and sub_order must be args
-
-        if len(self.plots[0].vectors) != 2:
-            print 'This command only works on a curve with two different plots.'
-            pass
+        new.info['columns'].append(
+            join_data_label(params['output deflection column'], 'N'))
+        d_data = data[:,data.info['columns'].index(
+                params['input deflection column'])]
+        new[:,-1] = d_data * data.info[params['spring constant info name']]
+        params['curve'].data[params['block']] = new
 
-        outplot=self.subtract_curves(sub_order=1)
 
-        plot_graph=self.list_of_events['plot_graph']
-        wx.PostEvent(self.frame,plot_graph(plots=[outplot]))
-
-    def _plug_init(self):
-        self.basecurrent=None
-        self.basepoints=None
-        self.autofile=''
-
-    def do_distance(self,args):
-        '''
-        DISTANCE
-        (generalvclamp.py)
-        Measure the distance (in nm) between two points.
-        For a standard experiment this is the delta X distance.
-        For a force clamp experiment this is the delta Y distance (actually becomes
-        an alias of zpiezo)
-        -----------------
-        Syntax: distance
-        '''
-        if self.current.curve.experiment == 'clamp':
-            print 'You wanted to use zpiezo perhaps?'
-            return
-        else:
-            dx,unitx,dy,unity=self._delta(set=1)
-            print str(dx*(10**9))+' nm'
-            to_dump='distance '+self.current.path+' '+str(dx*(10**9))+' nm'
-            self.outlet.push(to_dump)
-
-
-    def do_force(self,args):
-        '''
-        FORCE
-        (generalvclamp.py)
-        Measure the force difference (in pN) between two points
-        ---------------
-        Syntax: force
-        '''
-        if self.current.curve.experiment == 'clamp':
-            print 'This command makes no sense for a force clamp experiment.'
-            return
-        dx,unitx,dy,unity=self._delta(set=1)
-        print str(dy*(10**12))+' pN'
-        to_dump='force '+self.current.path+' '+str(dy*(10**12))+' pN'
-        self.outlet.push(to_dump)
-
-
-    def do_forcebase(self,args):
-        '''
-        FORCEBASE
-        (generalvclamp.py)
-        Measures the difference in force (in pN) between a point and a baseline
-        took as the average between two points.
-
-        The baseline is fixed once for a given curve and different force measurements,
-        unless the user wants it to be recalculated
-        ------------
-        Syntax: forcebase [rebase]
-                rebase: Forces forcebase to ask again the baseline
-                max: Instead of asking for a point to measure, asks for two points and use
-                     the maximum peak in between
-        '''
-        rebase=False #if true=we select rebase
-        maxpoint=False #if true=we measure the maximum peak
-
-        plot=self._get_displayed_plot()
-        whatset=1 #fixme: for all sets
-        if 'rebase' in args or (self.basecurrent != self.current.path):
-            rebase=True
-        if 'max' in args:
-            maxpoint=True
-
-        if rebase:
-            print 'Select baseline'
-            self.basepoints=self._measure_N_points(N=2, whatset=whatset)
-            self.basecurrent=self.current.path
-
-        if maxpoint:
-            print 'Select two points'
-            points=self._measure_N_points(N=2, whatset=whatset)
-            boundpoints=[points[0].index, points[1].index]
-            boundpoints.sort()
-            try:
-                y=min(plot.vectors[whatset][1][boundpoints[0]:boundpoints[1]])
-            except ValueError:
-                print 'Chosen interval not valid. Try picking it again. Did you pick the same point as begin and end of interval?'
-        else:
-            print 'Select point to measure'
-            points=self._measure_N_points(N=1, whatset=whatset)
-            #whatplot=points[0].dest
-            y=points[0].graph_coords[1]
-
-        #fixme: code duplication
-        boundaries=[self.basepoints[0].index, self.basepoints[1].index]
-        boundaries.sort()
-        to_average=plot.vectors[whatset][1][boundaries[0]:boundaries[1]] #y points to average
-
-        avg=np.mean(to_average)
-        forcebase=abs(y-avg)
-        print str(forcebase*(10**12))+' pN'
-        to_dump='forcebase '+self.current.path+' '+str(forcebase*(10**12))+' pN'
-        self.outlet.push(to_dump)
-
-    def plotmanip_multiplier(self, plot, current):
-        '''
-        Multiplies all the Y values of an SMFS curve by a value stored in the 'force_multiplier'
-        configuration variable. Useful for calibrations and other stuff.
-        '''
+class CantileverAdjustedExtensionCommand (Command):
+    """Remove cantilever extension from a total extension column.
+    """
+    def __init__(self, plugin):
+        super(CantileverAdjustedExtensionCommand, self).__init__(
+            name='add block cantilever adjusted extension array',
+            arguments=[
+                CurveArgument,
+                Argument(name='block', aliases=['set'], type='int', default=0,
+                         help="""
+Data block for which the adjusted extension should be calculated.  For
+an approach/retract force curve, `0` selects the approaching curve and
+`1` selects the retracting curve.
+""".strip()),
+                Argument(name='input distance column', type='string',
+                         default='surface distance (m)',
+                         help="""
+Name of the column to use as the distance input.
+""".strip()),
+                Argument(name='input deflection column', type='string',
+                         default='deflection (N)',
+                         help="""
+Name of the column to use as the deflection input.
+""".strip()),
+                Argument(name='output distance column', type='string',
+                         default='cantilever adjusted extension',
+                         help="""
+Name of the column (without units) to use as the deflection output.
+""".strip()),
+                Argument(name='spring constant info name', type='string',
+                         default='spring constant (N/m)',
+                         help="""
+Name of the spring constant in the `.info` dictionary.
+""".strip()),
+                ],
+            help=self.__doc__, plugin=plugin)
 
-        #not a smfs curve...
-        if current.curve.experiment != 'smfs':
-            return plot
+    def _run(self, hooke, inqueue, outqueue, params):
+        data = params['curve'].data[params['block']]
+        # HACK? rely on params['curve'] being bound to the local hooke
+        # playlist (i.e. not a copy, as you would get by passing a
+        # curve through the queue).  Ugh.  Stupid queues.  As an
+        # alternative, we could pass lookup information through the
+        # queue.
+        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 distance column'], 'm'))
+        z_data = data[:,data.info['columns'].index(
+                params['input distance column'])]
+        d_data = data[:,data.info['columns'].index(
+                params['input deflection column'])]
+        k = data.info[params['spring constant info name']]
 
-        #only one set is present...
-        if len(self.plots[0].vectors) != 2:
-            return plot
+        z_name,z_unit = split_data_label(params['input distance column'])
+        assert z_unit == 'm', params['input distance column']
+        d_name,d_unit = split_data_label(params['input deflection column'])
+        assert d_unit == 'N', params['input deflection column']
+        k_name,k_unit = split_data_label(params['spring constant info name'])
+        assert k_unit == 'N/m', params['spring constant info name']
 
-        #multiplier is 1...
-        if (self.config['force_multiplier']==1):
-            return plot
+        new[:,-1] = z_data - d_data / k
+        params['curve'].data[params['block']] = new
 
-        for i in range(len(plot.vectors[0][1])):
-            plot.vectors[0][1][i]=plot.vectors[0][1][i]*self.config['force_multiplier']        
 
-        for i in range(len(plot.vectors[1][1])):
-            plot.vectors[1][1][i]=plot.vectors[1][1][i]*self.config['force_multiplier']
+class generalvclampCommands(object):
 
-        return plot            
-   
-    
     def plotmanip_flatten(self, plot, current, customvalue=False):
         '''
         Subtracts a polynomial fit to the non-contact part of the curve, as to flatten it.
@@ -673,192 +659,3 @@ class generalvclampCommands(object):
 
         return plot
 
-    #---SLOPE---
-    def do_slope(self,args):
-        '''
-        SLOPE
-        (generalvclamp.py)
-        Measures the slope of a delimited chunk on the return trace.
-        The chunk can be delimited either by two manual clicks, or have
-        a fixed width, given as an argument.
-        ---------------
-        Syntax: slope [width]
-                The facultative [width] parameter specifies how many
-                points will be considered for the fit. If [width] is
-                specified, only one click will be required.
-        (c) Marco Brucale, Massimo Sandal 2008
-        '''
-
-        # Reads the facultative width argument
-        try:
-            fitspan=int(args)
-        except:
-            fitspan=0
-
-        # Decides between the two forms of user input, as per (args)
-        if fitspan == 0:
-            # Gets the Xs of two clicked points as indexes on the current curve vector
-            print 'Click twice to delimit chunk'
-            points=self._measure_N_points(N=2,whatset=1)
-        else:
-            print 'Click once on the leftmost point of the chunk (i.e.usually the peak)'
-            points=self._measure_N_points(N=1,whatset=1)
-            
-        slope=self._slope(points,fitspan)
-
-        # Outputs the relevant slope parameter
-        print 'Slope:'
-        print str(slope)
-        to_dump='slope '+self.current.path+' '+str(slope)
-        self.outlet.push(to_dump)
-
-    def _slope(self,points,fitspan):
-        # Calls the function linefit_between
-        parameters=[0,0,[],[]]
-        try:
-            clickedpoints=[points[0].index,points[1].index]
-            clickedpoints.sort()
-        except:
-            clickedpoints=[points[0].index-fitspan,points[0].index]        
-
-        try:
-            parameters=self.linefit_between(clickedpoints[0],clickedpoints[1])
-        except:
-            print 'Cannot fit. Did you click twice the same point?'
-            return
-             
-        # Outputs the relevant slope parameter
-        print 'Slope:'
-        print str(parameters[0])
-        to_dump='slope '+self.curve.path+' '+str(parameters[0])
-        self.outlet.push(to_dump)
-
-        # Makes a vector with the fitted parameters and sends it to the GUI
-        xtoplot=parameters[2]
-        ytoplot=[]
-        x=0
-        for x in xtoplot:
-            ytoplot.append((x*parameters[0])+parameters[1])
-
-        clickvector_x, clickvector_y=[], []
-        for item in points:
-            clickvector_x.append(item.graph_coords[0])
-            clickvector_y.append(item.graph_coords[1])
-
-        lineplot=self._get_displayed_plot(0) #get topmost displayed plot
-
-        lineplot.add_set(xtoplot,ytoplot)
-        lineplot.add_set(clickvector_x, clickvector_y)
-
-
-        if lineplot.styles==[]:
-            lineplot.styles=[None,None,None,'scatter']
-        else:
-            lineplot.styles+=[None,'scatter']
-        if lineplot.colors==[]:
-            lineplot.colors=[None,None,'black',None]
-        else:
-            lineplot.colors+=['black',None]
-        
-        
-        self._send_plot([lineplot])
-
-        return parameters[0]
-
-
-    def linefit_between(self,index1,index2,whatset=1):
-        '''
-        Creates two vectors (xtofit,ytofit) slicing out from the
-        current return trace a portion delimited by the two indexes
-        given as arguments.
-        Then does a least squares linear fit on that slice.
-        Finally returns [0]=the slope, [1]=the intercept of the
-        fitted 1st grade polynomial, and [2,3]=the actual (x,y) vectors
-        used for the fit.
-        (c) Marco Brucale, Massimo Sandal 2008
-        '''
-        # Translates the indexes into two vectors containing the x,y data to fit
-        xtofit=self.plots[0].vectors[whatset][0][index1:index2]
-        ytofit=self.plots[0].vectors[whatset][1][index1:index2]
-
-        # Does the actual linear fitting (simple least squares with numpy.polyfit)
-        linefit=[]
-        linefit=np.polyfit(xtofit,ytofit,1)
-
-        return (linefit[0],linefit[1],xtofit,ytofit)
-
-
-    def fit_interval_nm(self,start_index,plot,nm,backwards):
-          '''
-          Calculates the number of points to fit, given a fit interval in nm
-          start_index: index of point
-          plot: plot to use
-          backwards: if true, finds a point backwards.
-          '''
-          whatset=1 #FIXME: should be decidable
-          x_vect=plot.vectors[1][0]
-          
-          c=0
-          i=start_index
-          start=x_vect[start_index]
-          maxlen=len(x_vect)
-          while abs(x_vect[i]-x_vect[start_index])*(10**9) < nm:
-              if i==0 or i==maxlen-1: #we reached boundaries of vector!
-                  return c
-              
-              if backwards:
-                  i-=1
-              else:
-                  i+=1
-              c+=1
-          return c
-
-
-
-    def find_current_peaks(self,noflatten, a=True, maxpeak=True):
-            #Find peaks.
-            if a==True:
-                  a=self.convfilt_config['mindeviation']
-            try:
-                  abs_devs=float(a)
-            except:
-                  print "Bad input, using default."
-                  abs_devs=self.convfilt_config['mindeviation']
-
-            defplot=self.current.curve.default_plots()[0]
-            if not noflatten:
-                flatten=self._find_plotmanip('flatten') #Extract flatten plotmanip
-                defplot=flatten(defplot, self.current, customvalue=1) #Flatten curve before feeding it to has_peaks
-            pk_location,peak_size=self.has_peaks(defplot, abs_devs, maxpeak)
-            return pk_location, peak_size
-
-
-    def pickup_contact_point(self,N=1,whatset=1):
-        '''macro to pick up the contact point by clicking'''
-        contact_point=self._measure_N_points(N=1, whatset=1)[0]
-        contact_point_index=contact_point.index
-        self.wlccontact_point=contact_point
-        self.wlccontact_index=contact_point.index
-        self.wlccurrent=self.current.path
-        return contact_point, contact_point_index
-
-
-    def baseline_points(self,peak_location, displayed_plot):
-        clicks=self.config['baseline_clicks']
-        if clicks==0:
-            self.basepoints=[]
-            base_index_0=peak_location[-1]+self.fit_interval_nm(peak_location[-1], displayed_plot, self.config['auto_right_baseline'],False)
-            self.basepoints.append(self._clickize(displayed_plot.vectors[1][0],displayed_plot.vectors[1][1],base_index_0))
-            base_index_1=self.basepoints[0].index+self.fit_interval_nm(self.basepoints[0].index, displayed_plot, self.config['auto_left_baseline'],False)
-            self.basepoints.append(self._clickize(displayed_plot.vectors[1][0],displayed_plot.vectors[1][1],base_index_1))
-        elif clicks>0:
-            print 'Select baseline'
-            if clicks==1:
-                self.basepoints=self._measure_N_points(N=1, whatset=1)
-                base_index_1=self.basepoints[0].index+self.fit_interval_nm(self.basepoints[0].index, displayed_plot, self.config['auto_left_baseline'], False)
-                self.basepoints.append(self._clickize(displayed_plot.vectors[1][0],displayed_plot.vectors[1][1],base_index_1))
-            else:
-                self.basepoints=self._measure_N_points(N=2, whatset=1)
-            
-        self.basecurrent=self.current.path
-        return self.basepoints