"""
import copy
+import logging
import numpy
import scipy
super(VelocityClampPlugin, self).__init__(name='vclamp')
self._commands = [
SurfaceContactCommand(self), ForceCommand(self),
- CantileverAdjustedExtensionCommand(self),
+ CantileverAdjustedExtensionCommand(self), FlattenCommand(self),
]
def default_settings(self):
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
projects / hooke.git / commitdiff