1 # Copyright (C) 2008-2010 Alberto Gomez-Casado
4 # Massimo Sandal <devicerandom@gmail.com>
5 # W. Trevor King <wking@drexel.edu>
7 # This file is part of Hooke.
9 # Hooke is free software: you can redistribute it and/or modify it
10 # under the terms of the GNU Lesser General Public License as
11 # published by the Free Software Foundation, either version 3 of the
12 # License, or (at your option) any later version.
14 # Hooke is distributed in the hope that it will be useful, but WITHOUT
15 # ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
16 # or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General
17 # Public License for more details.
19 # You should have received a copy of the GNU Lesser General Public
20 # License along with Hooke. If not, see
21 # <http://www.gnu.org/licenses/>.
23 """The ``vclamp`` module provides :class:`VelocityClampPlugin` and
24 several associated :class:`hooke.command.Command`\s for handling
25 common velocity clamp analysis tasks.
34 from ..command import Argument, Failure, NullQueue
35 from ..config import Setting
36 from ..curve import Data
37 from ..util.fit import PoorFit, ModelFitter
38 from ..util.si import join_data_label, split_data_label
40 from .curve import ColumnAddingCommand
43 class SurfacePositionModel (ModelFitter):
44 """Bilinear surface position model.
46 The bilinear model is symmetric, but the parameter guessing and
47 sanity checks assume the contact region occurs for lower indicies
48 ("left of") the non-contact region. We also assume that
49 tip-surface attractions produce positive deflections.
53 Algorithm borrowed from WTK's `piezo package`_, specifically
54 from :func:`piezo.z_piezo_utils.analyzeSurfPosData`.
57 http://www.physics.drexel.edu/~wking/code/git/git.php?p=piezo.git
59 Fits the data to the bilinear :method:`model`.
61 In order for this model to produce a satisfactory fit, there
62 should be enough data in the off-surface region that interactions
63 due to proteins, etc. will not seriously skew the fit in the
64 off-surface region. If you don't have much of a tail, you can set
65 the `info` dict's `ignore non-contact before index` parameter to
66 the index of the last surface- or protein-related feature.
68 _fit_check_arguments = [
69 Argument('min slope ratio', type='float', default=10.,
71 Minimum `contact-slope/non-contact-slope` ratio for a "good" fit.
73 Argument('min contact fraction', type='float', default=0.02,
75 Minimum fraction of points in the contact region for a "good" fit.
77 Argument('max contact fraction', type='float', default=0.98,
79 Maximum fraction of points in the contact region for a "good" fit.
81 Argument('min slope guess ratio', type='float', default=0.5,
83 Minimum `fit-contact-slope/guessed-contact-slope` ratio for a "good" fit.
87 def model(self, params):
88 """A continuous, bilinear model.
95 p_0 + p_1 x & \text{if $x <= p_2$}, \\
96 p_0 + p_1 p_2 + p_3 (x-p_2) & \text{if $x >= p_2$}.
99 Where :math:`p_0` is a vertical offset, :math:`p_1` is the slope
100 of the first region, :math:`p_2` is the transition location, and
101 :math:`p_3` is the slope of the second region.
103 p = params # convenient alias
104 rNC_ignore = self.info['ignore non-contact before index']
105 if self.info['force zero non-contact slope'] == True:
107 p.append(0.) # restore the non-contact slope parameter
108 r2 = numpy.round(abs(p[2]))
110 r2 = min(r2, len(self._model_data))
111 self._model_data[:r2] = p[0] + p[1] * numpy.arange(r2)
112 if r2 < len(self._data)-1:
113 self._model_data[r2:] = \
114 p[0] + p[1]*p[2] + p[3] * numpy.arange(len(self._data)-r2)
116 self._model_data[r2:rNC_ignore] = self._data[r2:rNC_ignore]
117 return self._model_data
119 def set_data(self, data, info=None, *args, **kwargs):
120 super(SurfacePositionModel, self).set_data(data, info, *args, **kwargs)
123 if getattr(self, 'info', None) == None:
125 self.info.update(info)
127 ('force zero non-contact slope', False),
128 ('ignore non-contact before index', -1),
129 ('min position', 0), # Store postions etc. to avoid recalculating.
130 ('max position', len(data)),
131 ('max deflection', data.max()),
132 ('min deflection', data.min()),
134 if key not in self.info:
135 self.info[key] = value
138 self.info['max position'] - self.info['min position']),
140 self.info['max deflection'] - self.info['min deflection']),
142 if key not in self.info:
143 self.info[key] = value
144 for argument in self._fit_check_arguments:
145 if argument.name not in self.info:
146 self.info[argument.name] = argument.default
148 def guess_initial_params(self, outqueue=None):
149 """Guess the initial parameters.
153 We guess initial parameters such that the offset (:math:`p_1`)
154 matches the minimum deflection, the kink (:math:`p_2`) occurs
155 at the first point that the deflection crosses the middle of
156 its range, the initial (contact) slope (:math:`p_0`) produces
157 the right-most deflection at the kink point, and the final
158 (non-contact) slope (:math:`p_3`) is zero.
160 In the event of a tie, :meth:`argmax` returns the lowest index
161 to the maximum value.
162 >>> (numpy.arange(10) >= 5).argmax()
165 left_offset = self.info['min deflection']
166 middle_deflection = (self.info['min deflection']
167 + self.info['deflection range']/2.)
168 kink_position = 2*(self._data > middle_deflection).argmax()
169 if kink_position == 0:
170 # jump vibration at the start of the retraction?
171 start = int(min(max(20, 0.01 * len(self._data)), 0.5*len(self._data)))
172 std = self._data[:start].std()
173 left_offset = self._data[start].mean()
175 while abs(self._data[stop] - left_offset) < 3*std:
177 left_slope = (self._data[stop-start:stop].mean()
178 - left_offset) / (stop-start)
179 left_offset -= left_slope * start/2
180 kink_position = (self._data[-1] - left_offset)/left_slope
182 left_slope = (self._data[-1] - self.info['min deflection']
185 self.info['guessed contact slope'] = left_slope
186 params = [left_offset, left_slope, kink_position, right_slope]
187 if self.info['force zero non-contact slope'] == True:
191 def fit(self, *args, **kwargs):
192 """Fit the model to the data.
196 We change the `epsfcn` default from :func:`scipy.optimize.leastsq`'s
197 `0` to `1e-3`, so the initial Jacobian estimate takes larger steps,
198 which helps avoid being trapped in noise-generated local minima.
200 self.info['guessed contact slope'] = None
201 if 'epsfcn' not in kwargs:
202 kwargs['epsfcn'] = 1e-3 # take big steps to estimate the Jacobian
203 params = super(SurfacePositionModel, self).fit(*args, **kwargs)
204 params[2] = abs(params[2])
205 if self.info['force zero non-contact slope'] == True:
206 params = list(params)
207 params.append(0.) # restore the non-contact slope parameter
209 # check that the fit is reasonable, see the :meth:`model` docstring
210 # for parameter descriptions.
211 slope_ratio = abs(params[1]/params[3])
212 if slope_ratio < self.info['min slope ratio']:
214 'Slope in non-contact region, or no slope in contact (slope ratio %g less than %g)'
215 % (slope_ratio, self.info['min slope ratio']))
216 contact_fraction = ((params[2]-self.info['min position'])
217 /self.info['position range'])
218 if contact_fraction < self.info['min contact fraction']:
220 'No kink (contact fraction %g less than %g)'
221 % (contact_fraction, self.info['min contact fraction']))
222 if contact_fraction > self.info['max contact fraction']:
224 'No kink (contact fraction %g greater than %g)'
225 % (contact_fraction, self.info['max contact fraction']))
226 slope_guess_ratio = abs(params[1]/self.info['guessed contact slope'])
227 if (self.info['guessed contact slope'] != None
228 and slope_guess_ratio < self.info['min slope guess ratio']):
230 'Too far (contact slope off guess by %g less than %g)'
231 % (slope_guess_ratio, self.info['min slope guess ratio']))
235 class VelocityClampPlugin (Plugin):
237 super(VelocityClampPlugin, self).__init__(name='vclamp')
239 SurfaceContactCommand(self), ForceCommand(self),
240 CantileverAdjustedExtensionCommand(self), FlattenCommand(self),
243 def default_settings(self):
245 Setting(section=self.setting_section, help=self.__doc__),
246 Setting(section=self.setting_section,
247 option='surface contact point algorithm',
249 help='Select the surface contact point algorithm. See the documentation for descriptions of available algorithms.')
253 class SurfaceContactCommand (ColumnAddingCommand):
254 """Automatically determine a block's surface contact point.
256 You can select the contact point algorithm with the creatively
257 named `surface contact point algorithm` configuration setting.
258 Currently available options are:
260 * fmms (:meth:`find_contact_point_fmms`)
261 * ms (:meth:`find_contact_point_ms`)
262 * wtk (:meth:`find_contact_point_wtk`)
264 def __init__(self, plugin):
265 self._wtk_fit_check_arguments = []
266 for argument in SurfacePositionModel._fit_check_arguments:
267 arg = copy.deepcopy(argument)
268 arg._help += ' (wtk model)'
269 self._wtk_fit_check_arguments.append(arg)
270 super(SurfaceContactCommand, self).__init__(
271 name='zero surface contact point',
273 ('distance column', 'z piezo (m)', """
274 Name of the column to use as the surface position input.
276 ('deflection column', 'deflection (m)', """
277 Name of the column to use as the deflection input.
281 ('output distance column', 'surface distance', """
282 Name of the column (without units) to use as the surface position output.
284 ('output deflection column', 'surface deflection', """
285 Name of the column (without units) to use as the deflection output.
289 Argument(name='ignore index', type='int', default=None,
291 Ignore the residual from the non-contact region before the indexed
292 point (for the `wtk` algorithm).
294 Argument(name='ignore after last peak info name',
295 type='string', default=None,
297 As an alternative to 'ignore index', ignore after the last peak in the
298 peak list stored in the `.info` dictionary.
300 Argument(name='distance info name', type='string',
301 default='surface distance offset',
303 Name (without units) for storing the distance offset in the `.info` dictionary.
305 Argument(name='deflection info name', type='string',
306 default='surface deflection offset',
308 Name (without units) for storing the deflection offset in the `.info` dictionary.
310 Argument(name='fit parameters info name', type='string',
311 default='surface deflection offset',
313 Name (without units) for storing fit parameters in the `.info` dictionary.
315 ] + self._wtk_fit_check_arguments,
316 help=self.__doc__, plugin=plugin)
318 def _run(self, hooke, inqueue, outqueue, params):
319 self._add_to_command_stack(params)
320 params = self._setup_params(hooke=hooke, params=params)
321 block = self._block(hooke=hooke, params=params)
322 dist_data = self._get_column(hooke=hooke, params=params,
323 column_name='distance column')
324 def_data = self._get_column(hooke=hooke, params=params,
325 column_name='deflection column')
326 i,def_offset,ps = self.find_contact_point(
327 params, dist_data, def_data, outqueue)
328 dist_offset = dist_data[i]
329 block.info[params['distance info name']] = dist_offset
330 block.info[params['deflection info name']] = def_offset
331 block.info[params['fit parameters info name']] = ps
332 self._set_column(hooke=hooke, params=params,
333 column_name='output distance column',
334 values=dist_data - dist_offset)
335 self._set_column(hooke=hooke, params=params,
336 column_name='output deflection column',
337 values=def_data - def_offset)
339 def _setup_params(self, hooke, params):
340 name,dist_unit = split_data_label(params['distance column'])
341 name,def_unit = split_data_label(params['deflection column'])
342 params['output distance column'] = join_data_label(
343 params['output distance column'], dist_unit)
344 params['output deflection column'] = join_data_label(
345 params['output deflection column'], def_unit)
346 params['distance info name'] = join_data_label(
347 params['distance info name'], dist_unit)
348 params['deflection info name'] = join_data_label(
349 params['deflection info name'], def_unit)
352 def find_contact_point(self, params, z_data, d_data, outqueue=None):
353 """Railyard for the `find_contact_point_*` family.
355 Uses the `surface contact point algorithm` configuration
356 setting to call the appropriate backend algorithm.
358 fn = getattr(self, 'find_contact_point_%s'
359 % self.plugin.config['surface contact point algorithm'])
360 return fn(params, z_data, d_data, outqueue)
362 def find_contact_point_fmms(self, params, z_data, d_data, outqueue=None):
363 """Algorithm by Francesco Musiani and Massimo Sandal.
369 0) Driver-specific workarounds, e.g. deal with the PicoForce
370 trigger bug by excluding retraction portions with excessive
372 1) Select the second half (non-contact side) of the retraction
374 2) Fit the selection to a line.
375 3) If the fit is not almost horizontal, halve the selection
377 4) Average the selection and use it as a baseline.
378 5) Slide in from the start (contact side) of the retraction
379 curve, until you find a point with greater than baseline
380 deflection. That point is the contact point.
382 if params['curve'].info['filetype'] == 'picoforce':
383 # Take care of the picoforce trigger bug (TODO: example
384 # data file demonstrating the bug). We exclude portions
385 # of the curve that have too much standard deviation.
386 # Yes, a lot of magic is here.
387 check_start = len(d_data)-len(d_data)/20
388 monster_start = len(d_data)
390 # look at the non-contact tail
391 non_monster = d_data[check_start:monster_start]
392 if non_monster.std() < 2e-10: # HACK: hardcoded cutoff
394 else: # move further away from the monster
395 check_start -= len(d_data)/50
396 monster_start -= len(d_data)/50
397 z_data = z_data[:monster_start]
398 d_data = d_data[:monster_start]
400 # take half of the thing to start
401 selection_start = len(d_data)/2
403 z_chunk = z_data[selection_start:]
404 d_chunk = d_data[selection_start:]
405 slope,intercept,r,two_tailed_prob,stderr_of_the_estimate = \
406 scipy.stats.linregress(z_chunk, d_chunk)
407 # We stop if we found an almost-horizontal fit or if we're
408 # getting to small a selection. FIXME: 0.1 and 5./6 here
409 # are "magic numbers" (although reasonable)
410 if (abs(slope) < 0.1 # deflection (m) / surface (m)
411 or selection_start > 5./6*len(d_data)):
413 selection_start += 10
415 d_baseline = d_chunk.mean()
417 # find the first point above the calculated baseline
419 while i < len(d_data) and d_data[i] < ymean:
421 return (i, d_baseline, {})
423 def find_contact_point_ms(self, params, z_data, d_data, outqueue=None):
424 """Algorithm by Massimo Sandal.
428 WTK: At least the commits are by Massimo, and I see no notes
429 attributing the algorithm to anyone else.
435 xext=raw_plot.vectors[0][0]
436 yext=raw_plot.vectors[0][1]
437 xret2=raw_plot.vectors[1][0]
438 yret=raw_plot.vectors[1][1]
440 first_point=[xext[0], yext[0]]
441 last_point=[xext[-1], yext[-1]]
443 #regr=scipy.polyfit(first_point, last_point,1)[0:2]
444 diffx=abs(first_point[0]-last_point[0])
445 diffy=abs(first_point[1]-last_point[1])
447 #using polyfit results in numerical errors. good old algebra.
449 b=first_point[1]-(a*first_point[0])
450 baseline=scipy.polyval((a,b), xext)
452 ysub=[item-basitem for item,basitem in zip(yext,baseline)]
454 contact=ysub.index(min(ysub))
456 return xext,ysub,contact
458 #now, exploit a ClickedPoint instance to calculate index...
460 dummy.absolute_coords=(x_intercept,y_intercept)
461 dummy.find_graph_coords(xret2,yret)
464 return dummy.index, regr, regr_contact
468 def find_contact_point_wtk(self, params, z_data, d_data, outqueue=None):
469 """Algorithm by W. Trevor King.
473 Uses :class:`SurfacePositionModel` internally.
475 reverse = z_data[0] > z_data[-1]
476 if reverse == True: # approaching, contact region on the right
477 d_data = d_data[::-1]
478 s = SurfacePositionModel(d_data, info={
479 'force zero non-contact slope':True},
481 for argument in self._wtk_fit_check_arguments:
482 s.info[argument.name] = params[argument.name]
484 if params['ignore index'] != None:
485 ignore_index = params['ignore index']
486 elif params['ignore after last peak info name'] != None:
487 peaks = z_data.info[params['ignore after last peak info name']]
488 if not len(peaks) > 0:
489 raise Failure('Need at least one peak in %s, not %s'
490 % (params['ignore after last peak info name'],
492 ignore_index = peaks[-1].post_index()
493 if ignore_index != None:
494 s.info['ignore non-contact before index'] = ignore_index
495 offset,contact_slope,surface_index,non_contact_slope = s.fit(
499 'contact slope': contact_slope,
500 'surface index': surface_index,
501 'non-contact slope': non_contact_slope,
504 deflection_offset = offset + contact_slope*surface_index,
506 surface_index = len(d_data)-1-surface_index
507 return (numpy.round(surface_index), deflection_offset, info)
510 class ForceCommand (ColumnAddingCommand):
511 """Convert a deflection column from meters to newtons.
513 def __init__(self, plugin):
514 super(ForceCommand, self).__init__(
515 name='convert distance to force',
517 ('deflection column', 'surface deflection (m)', """
518 Name of the column to use as the deflection input.
522 ('output deflection column', 'deflection', """
523 Name of the column (without units) to use as the deflection output.
527 Argument(name='spring constant info name', type='string',
528 default='spring constant (N/m)',
530 Name of the spring constant in the `.info` dictionary.
533 help=self.__doc__, plugin=plugin)
535 def _run(self, hooke, inqueue, outqueue, params):
536 self._add_to_command_stack(params)
537 params = self._setup_params(hooke=hooke, params=params)
538 def_data = self._get_column(hooke=hooke, params=params,
539 column_name='deflection column')
540 out = def_data * def_data.info[params['spring constant info name']]
541 self._set_column(hooke=hooke, params=params,
542 column_name='output deflection column',
545 def _setup_params(self, hooke, params):
546 name,in_unit = split_data_label(params['deflection column'])
547 out_unit = 'N' # HACK: extract target units from k_unit.
548 params['output deflection column'] = join_data_label(
549 params['output deflection column'], out_unit)
550 name,k_unit = split_data_label(params['spring constant info name'])
551 expected_k_unit = '%s/%s' % (out_unit, in_unit)
552 if k_unit != expected_k_unit:
553 raise Failure('Cannot convert from %s to %s with %s'
554 % (params['deflection column'],
555 params['output deflection column'],
556 params['spring constant info name']))
560 class CantileverAdjustedExtensionCommand (ColumnAddingCommand):
561 """Remove cantilever extension from a total extension column.
563 If `distance column` and `deflection column` have the same units
564 (e.g. `z piezo (m)` and `deflection (m)`), `spring constant info
565 name` is ignored and a deflection/distance conversion factor of
568 def __init__(self, plugin):
569 super(CantileverAdjustedExtensionCommand, self).__init__(
570 name='remove cantilever from extension',
572 ('distance column', 'surface distance (m)', """
573 Name of the column to use as the surface position input.
575 ('deflection column', 'deflection (N)', """
576 Name of the column to use as the deflection input.
580 ('output distance column', 'cantilever adjusted extension', """
581 Name of the column (without units) to use as the surface position output.
585 Argument(name='spring constant info name', type='string',
586 default='spring constant (N/m)',
588 Name of the spring constant in the `.info` dictionary.
591 help=self.__doc__, plugin=plugin)
593 def _run(self, hooke, inqueue, outqueue, params):
594 self._add_to_command_stack(params)
595 params = self._setup_params(hooke=hooke, params=params)
596 def_data = self._get_column(hooke=hooke, params=params,
597 column_name='deflection column')
598 dist_data = self._get_column(hooke=hooke, params=params,
599 column_name='distance column')
600 if params['spring constant info name'] == None:
601 k = 1.0 # distance and deflection in the same units
603 k = def_data.info[params['spring constant info name']]
604 self._set_column(hooke=hooke, params=params,
605 column_name='output distance column',
606 values=dist_data - def_data / k)
608 def _setup_params(self, hooke, params):
609 name,dist_unit = split_data_label(params['distance column'])
610 name,def_unit = split_data_label(params['deflection column'])
611 params['output distance column'] = join_data_label(
612 params['output distance column'], dist_unit)
613 if dist_unit == def_unit:
614 params['spring constant info name'] == None
616 name,k_unit = split_data_label(params['spring constant info name'])
617 expected_k_unit = '%s/%s' % (def_unit, dist_unit)
618 if k_unit != expected_k_unit:
619 raise Failure('Cannot convert from %s to %s with %s'
620 % (params['deflection column'],
621 params['output distance column'],
622 params['spring constant info name']))
626 class FlattenCommand (ColumnAddingCommand):
627 """Flatten a deflection column.
629 Subtracts a polynomial fit from the non-contact part of the curve
630 to flatten it. The best polynomial fit is chosen among
631 polynomials of degree 1 to `max degree`.
633 .. todo: Why does flattening use a polynomial fit and not a sinusoid?
634 Isn't most of the oscillation due to laser interference?
635 See Jaschke 1995 ( 10.1063/1.1146018 )
636 and the figure 4 caption of Weisenhorn 1992 ( 10.1103/PhysRevB.45.11226 )
638 def __init__(self, plugin):
639 super(FlattenCommand, self).__init__(
640 name='polynomial flatten',
642 ('distance column', 'surface distance (m)', """
643 Name of the column to use as the surface position input.
645 ('deflection column', 'deflection (N)', """
646 Name of the column to use as the deflection input.
650 ('output deflection column', 'flattened deflection', """
651 Name of the column (without units) to use as the deflection output.
655 Argument(name='degree', type='int',
658 Order of the polynomial used for flattening. Using values greater
659 than one usually doesn't help and can give artifacts. However, it
660 could be useful too. (TODO: Back this up with some theory...)
662 Argument(name='fit info name', type='string',
663 default='flatten fit',
665 Name of the flattening information in the `.info` dictionary.
668 help=self.__doc__, plugin=plugin)
670 def _run(self, hooke, inqueue, outqueue, params):
671 self._add_to_command_stack(params)
672 params = self._setup_params(hooke=hooke, params=params)
673 block = self._block(hooke=hooke, params=params)
674 dist_data = self._get_column(hooke=hooke, params=params,
675 column_name='distance column')
676 def_data = self._get_column(hooke=hooke, params=params,
677 column_name='deflection column')
678 degree = params['degree']
680 indices = numpy.argwhere(mask)
681 if len(indices) == 0:
682 raise Failure('no positive distance values in %s'
683 % params['distance column'])
684 dist_nc = dist_data[indices].flatten()
685 def_nc = def_data[indices].flatten()
688 poly_values = scipy.polyfit(dist_nc, def_nc, degree)
689 def_nc_fit = scipy.polyval(poly_values, dist_nc)
691 raise Failure('failed to flatten with a degree %d polynomial: %s'
693 error = numpy.sqrt((def_nc_fit-def_nc)**2).sum() / len(def_nc)
694 block.info[params['fit info name']] = {
697 'polynomial values':poly_values,
700 - numpy.invert(mask)*scipy.polyval(poly_values[-1:], dist_data)
701 - mask*scipy.polyval(poly_values, dist_data))
702 self._set_column(hooke=hooke, params=params,
703 column_name='output deflection column',
706 def _setup_params(self, hooke, params):
707 d_name,d_unit = split_data_label(params['deflection column'])
708 params['output deflection column'] = join_data_label(
709 params['output deflection column'], d_unit)