Convert the old flatten plotmanip into vclamp's FlattenCommand
authorW. Trevor King <wking@drexel.edu>
Mon, 9 Aug 2010 01:59:16 +0000 (21:59 -0400)
committerW. Trevor King <wking@drexel.edu>
Mon, 9 Aug 2010 01:59:16 +0000 (21:59 -0400)
hooke/plugin/vclamp.py

index fbd6335aec334c78950ad04a19361c1f32aad65b..77887db22a6b14894506a20cc15db038a34cb7f6 100644 (file)
@@ -26,6 +26,7 @@ common velocity clamp analysis tasks.
 """
 
 import copy
+import logging
 
 import numpy
 import scipy
@@ -215,7 +216,7 @@ class VelocityClampPlugin (Plugin):
         super(VelocityClampPlugin, self).__init__(name='vclamp')
         self._commands = [
             SurfaceContactCommand(self), ForceCommand(self),
-            CantileverAdjustedExtensionCommand(self),
+            CantileverAdjustedExtensionCommand(self), FlattenCommand(self),
             ]
 
     def default_settings(self):
@@ -579,83 +580,109 @@ Name of the spring constant in the `.info` dictionary.
         params['curve'].data[params['block']] = new
 
 
-class generalvclampCommands(object):
+class FlattenCommand (Command):
+    """Flatten a deflection column.
 
-    def plotmanip_flatten(self, plot, current, customvalue=False):
-        '''
-        Subtracts a polynomial fit to the non-contact part of the curve, as to flatten it.
-        the best polynomial fit is chosen among polynomials of degree 1 to n, where n is
-        given by the configuration file or by the customvalue.
+    Subtracts a polynomial fit from the non-contact part of the curve
+    to flatten it.  The best polynomial fit is chosen among
+    polynomials of degree 1 to `max degree`.
 
-        customvalue= int (>0) --> starts the function even if config says no (default=False)
-        '''
-
-        #not a smfs curve...
-        if current.curve.experiment != 'smfs':
-            return plot
-
-        #only one set is present...
-        if len(self.plots[0].vectors) != 2:
-            return plot
-
-        #config is not flatten, and customvalue flag is false too
-        if (not self.config['flatten']) and (not customvalue):
-            return plot
-
-        max_exponent=12
-        delta_contact=0
+    .. todo:  Why does flattening use a polynomial fit and not a sinusoid?
+      Isn't most of the oscillation due to laser interference?
+      See Jaschke 1995 ( 10.1063/1.1146018 )
+      and the figure 4 caption of Weisenhorn 1992 ( 10.1103/PhysRevB.45.11226 )
+    """
+    def __init__(self, plugin):
+        super(FlattenCommand, self).__init__(
+            name='add flattened 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='max degree', type='int',
+                         default=1,
+                         help="""
+Highest order polynomial to consider for flattening.  Using values
+greater than one usually doesn't help and can give artifacts.
+However, it could be useful too.  (TODO: Back this up with some
+theory...)
+""".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 deflection column', type='string',
+                         default='flattened deflection',
+                         help="""
+Name of the column (without units) to use as the deflection output.
+""".strip()),
+                Argument(name='fit info name', type='string',
+                         default='flatten fit',
+                         help="""
+Name of the flattening information in the `.info` dictionary.
+""".strip()),
+                ],
+            help=self.__doc__, plugin=plugin)
 
-        if customvalue:
-            max_cycles=customvalue
-        else:
-            max_cycles=self.config['flatten'] #Using > 1 usually doesn't help and can give artefacts. However, it could be useful too.
-
-        contact_index=self.find_contact_point()
-
-        valn=[[] for item in range(max_exponent)]
-        yrn=[0.0 for item in range(max_exponent)]
-        errn=[0.0 for item in range(max_exponent)]
-        
-        #Check if we have a proper numerical value
-        try:
-            zzz=int(max_cycles)
-        except:
-            #Loudly and annoyingly complain if it's not a number, then fallback to zero
-            print '''Warning: flatten value is not a number!
-            Use "set flatten" or edit hooke.conf to set it properly
-            Using zero.'''
-            max_cycles=0
-        
-        for i in range(int(max_cycles)):
-
-            x_ext=plot.vectors[0][0][contact_index+delta_contact:]
-            y_ext=plot.vectors[0][1][contact_index+delta_contact:]
-            x_ret=plot.vectors[1][0][contact_index+delta_contact:]
-            y_ret=plot.vectors[1][1][contact_index+delta_contact:]
-            for exponent in range(max_exponent):
-                try:
-                    valn[exponent]=sp.polyfit(x_ext,y_ext,exponent)
-                    yrn[exponent]=sp.polyval(valn[exponent],x_ret)
-                    errn[exponent]=sp.sqrt(sum((yrn[exponent]-y_ext)**2)/float(len(y_ext)))
-                except Exception,e:
-                    print 'Cannot flatten!'
-                    print e
-                    return plot
-
-            best_exponent=errn.index(min(errn))
-
-            #extension
-            ycorr_ext=y_ext-yrn[best_exponent]+y_ext[0] #noncontact part
-            yjoin_ext=np.array(plot.vectors[0][1][0:contact_index+delta_contact]) #contact part
-            #retraction
-            ycorr_ret=y_ret-yrn[best_exponent]+y_ext[0] #noncontact part
-            yjoin_ret=np.array(plot.vectors[1][1][0:contact_index+delta_contact]) #contact part
-
-            ycorr_ext=np.concatenate((yjoin_ext, ycorr_ext))
-            ycorr_ret=np.concatenate((yjoin_ret, ycorr_ret))
-
-            plot.vectors[0][1]=list(ycorr_ext)
-            plot.vectors[1][1]=list(ycorr_ret)
-
-        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
+        z_data = data[:,data.info['columns'].index(
+                params['input distance column'])]
+        d_data = data[:,data.info['columns'].index(
+                params['input deflection column'])]
 
+        d_name,d_unit = split_data_label(params['input deflection column'])
+        new.info['columns'].append(
+            join_data_label(params['output deflection column'], d_unit))
+
+        contact_index = numpy.absolute(z_data).argmin()
+        mask = z_data > 0
+        indices = numpy.argwhere(mask)
+        z_nc = z_data[indices].flatten()
+        d_nc = d_data[indices].flatten()
+
+        min_err = numpy.inf
+        degree = poly_values = None
+        log = logging.getLogger('hooke')
+        for deg in range(params['max degree']):
+            try:
+                pv = scipy.polyfit(z_nc, d_nc, deg)
+                df = scipy.polyval(pv, z_nc)
+                err = numpy.sqrt((df-d_nc)**2).sum()
+            except Exception,e:
+                log.warn('failed to flatten with a degree %d polynomial: %s'
+                         % (deg, e))
+                continue
+            if err < min_err:  # new best fit
+                min_err = err
+                degree = deg
+                poly_values = pv
+
+        if degree == None:
+            raise Failure('failed to flatten with all degrees')
+        new.info[params['fit info name']] = {
+            'error':min_err/len(z_nc),
+            'degree':degree,
+            'max degree':params['max degree'],
+            'polynomial values':poly_values,
+            }
+        new[:,-1] = d_data - mask*scipy.polyval(poly_values, z_data)
+        params['curve'].data[params['block']] = new