From: W. Trevor King Date: Mon, 9 Aug 2010 01:59:16 +0000 (-0400) Subject: Convert the old flatten plotmanip into vclamp's FlattenCommand X-Git-Url: http://git.tremily.us/?a=commitdiff_plain;h=4b523539d61ee1f143e4d7c2d7ab9e150351e141;p=hooke.git Convert the old flatten plotmanip into vclamp's FlattenCommand --- diff --git a/hooke/plugin/vclamp.py b/hooke/plugin/vclamp.py index fbd6335..77887db 100644 --- a/hooke/plugin/vclamp.py +++ b/hooke/plugin/vclamp.py @@ -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