- 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
-
- #---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]
+ Notes
+ -----
+ .. math::
+
+ y = \begin{cases}
+ p_0 + p_1 x & \text{if $x <= p_2$}, \\
+ p_0 + p_1 p_2 + p_3 (x-p_2) & \text{if $x >= p_2$}.
+ \end{cases}
+
+ Where :math:`p_0` is a vertical offset, :math:`p_1` is the slope
+ of the first region, :math:`p_2` is the transition location, and
+ :math:`p_3` is the slope of the second region.
+ """
+ p = params # convenient alias
+ rNC_ignore = self.info['ignore non-contact before index']
+ if self.info['force zero non-contact slope'] == True:
+ p = list(p)
+ p.append(0.) # restore the non-contact slope parameter
+ 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)
+ if r2 < rNC_ignore:
+ self._model_data[r2:rNC_ignore] = self._data[r2:rNC_ignore]
+ return self._model_data
+
+ def set_data(self, data, info=None, *args, **kwargs):
+ super(SurfacePositionModel, self).set_data(data, info, *args, **kwargs)
+ if self.info == None:
+ self.info = {}
+ for key,value in [
+ ('force zero non-contact slope', False),
+ ('ignore non-contact before index', 6158),
+ ('min position', 0), # Store postions etc. to avoid recalculating.
+ ('max position', len(data)),
+ ('max deflection', data.max()),
+ ('min deflection', data.min()),
+ ]:
+ if key not in self.info:
+ self.info[key] = value
+ for key,value in [
+ ('position range',
+ self.info['max position'] - self.info['min position']),
+ ('deflection range',
+ self.info['max deflection'] - self.info['min deflection']),
+ ]:
+ if key not in self.info:
+ self.info[key] = value
+
+ def guess_initial_params(self, outqueue=None):
+ """Guess the initial parameters.
+
+ Notes
+ -----
+ We guess initial parameters such that the offset (:math:`p_1`)
+ matches the minimum deflection, the kink (:math:`p_2`) occurs in
+ the middle of the data, the initial (contact) slope (:math:`p_0`)
+ produces the maximum deflection at the left-most point, and the
+ final (non-contact) slope (:math:`p_3`) is zero.
+ """
+ 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
+ self.info['guessed contact slope'] = left_slope
+ params = [left_offset, left_slope, kink_position, right_slope]
+ if self.info['force zero non-contact slope'] == True:
+ params = params[:-1]
+ return params
+
+ def fit(self, *args, **kwargs):
+ """Fit the model to the data.
+
+ Notes
+ -----
+ We change the `epsfcn` default from :func:`scipy.optimize.leastsq`'s
+ `0` to `1e-3`, so the initial Jacobian estimate takes larger steps,
+ which helps avoid being trapped in noise-generated local minima.
+ """
+ self.info['guessed contact slope'] = None
+ if 'epsfcn' not in kwargs:
+ kwargs['epsfcn'] = 1e-3 # take big steps to estimate the Jacobian
+ params = super(SurfacePositionModel, self).fit(*args, **kwargs)
+ params[2] = abs(params[2])
+ if self.info['force zero non-contact slope'] == 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. 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(
+ 'No kink (kink %g less than %g, need more space to left)'
+ % (params[2],
+ self.info['min position']+0.02*self.info['position range']))
+ if params[2] > self.info['max position']-0.02*self.info['position range']:
+ raise poorFit(
+ 'No kink (kink %g more than %g, need more space to right)'
+ % (params[2],
+ self.info['max position']-0.02*self.info['position range']))
+ if (self.info['guessed contact slope'] != None
+ 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 (Plugin):
+ def __init__(self):
+ super(VelocityClampPlugin, self).__init__(name='vclamp')
+ self._commands = [
+ SurfaceContactCommand(self), ForceCommand(self),
+ CantileverAdjustedExtensionCommand(self), FlattenCommand(self),
+ ]
+
+ def default_settings(self):
+ return [
+ Setting(section=self.setting_section, help=self.__doc__),
+ Setting(section=self.setting_section,
+ option='surface contact point algorithm',
+ value='wtk',
+ help='Select the surface contact point algorithm. See the documentation for descriptions of available algorithms.')
+ ]
+
+
+class SurfaceContactCommand (Command):
+ """Automatically determine a block's surface contact point.
+
+ You can select the contact point algorithm with the creatively
+ named `surface contact point algorithm` configuration setting.
+ Currently available options are:
+
+ * 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__(
+ name='zero block surface contact point',
+ arguments=[
+ CurveArgument,
+ Argument(name='block', aliases=['set'], type='int', default=0,
+ help="""
+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 position 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 position 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 fit parameters in the `.info` dictionary.
+""".strip()),
+ ],
+ help=self.__doc__, plugin=plugin)
+
+ 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]+2), dtype=data.dtype)
+ new.info = copy.deepcopy(data.info)
+ new[:,:-2] = data
+ 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(curve, z_data, d_data, outqueue)
+
+ def find_contact_point_fmms(self, curve, z_data, d_data, outqueue=None):
+ """Algorithm by Francesco Musiani and Massimo Sandal.
+
+ Notes
+ -----
+ Algorithm:
+
+ 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 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
+ -----
+ WTK: At least the commits are by Massimo, and I see no notes
+ attributing the algorithm to anyone else.
+
+ Algorithm:
+
+ * ?
+ """
+ xext=raw_plot.vectors[0][0]
+ yext=raw_plot.vectors[0][1]
+ xret2=raw_plot.vectors[1][0]
+ yret=raw_plot.vectors[1][1]
+
+ first_point=[xext[0], yext[0]]
+ last_point=[xext[-1], yext[-1]]
+
+ #regr=scipy.polyfit(first_point, last_point,1)[0:2]
+ diffx=abs(first_point[0]-last_point[0])
+ diffy=abs(first_point[1]-last_point[1])
+
+ #using polyfit results in numerical errors. good old algebra.
+ a=diffy/diffx
+ b=first_point[1]-(a*first_point[0])
+ baseline=scipy.polyval((a,b), xext)
+
+ ysub=[item-basitem for item,basitem in zip(yext,baseline)]
+
+ contact=ysub.index(min(ysub))
+
+ return xext,ysub,contact
+
+ #now, exploit a ClickedPoint instance to calculate index...
+ dummy=ClickedPoint()
+ dummy.absolute_coords=(x_intercept,y_intercept)
+ dummy.find_graph_coords(xret2,yret)
+
+ if debug:
+ return dummy.index, regr, regr_contact