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 left_slope = ((self._data[-1] - self.info['min deflection'])
172 self.info['guessed contact slope'] = left_slope
173 params = [left_offset, left_slope, kink_position, right_slope]
174 if self.info['force zero non-contact slope'] == True:
178 def fit(self, *args, **kwargs):
179 """Fit the model to the data.
183 We change the `epsfcn` default from :func:`scipy.optimize.leastsq`'s
184 `0` to `1e-3`, so the initial Jacobian estimate takes larger steps,
185 which helps avoid being trapped in noise-generated local minima.
187 self.info['guessed contact slope'] = None
188 if 'epsfcn' not in kwargs:
189 kwargs['epsfcn'] = 1e-3 # take big steps to estimate the Jacobian
190 params = super(SurfacePositionModel, self).fit(*args, **kwargs)
191 params[2] = abs(params[2])
192 if self.info['force zero non-contact slope'] == True:
193 params = list(params)
194 params.append(0.) # restore the non-contact slope parameter
196 # check that the fit is reasonable, see the :meth:`model` docstring
197 # for parameter descriptions.
198 slope_ratio = abs(params[1]/params[3])
199 if slope_ratio < self.info['min slope ratio']:
201 'Slope in non-contact region, or no slope in contact (slope ratio %g less than %g)'
202 % (slope_ratio, self.info['min slope ratio']))
203 contact_fraction = ((params[2]-self.info['min position'])
204 /self.info['position range'])
205 if contact_fraction < self.info['min contact fraction']:
207 'No kink (contact fraction %g less than %g)'
208 % (contact_fraction, self.info['min contact fraction']))
209 if contact_fraction > self.info['max contact fraction']:
211 'No kink (contact fraction %g greater than %g)'
212 % (contact_fraction, self.info['max contact fraction']))
213 slope_guess_ratio = abs(params[1]/self.info['guessed contact slope'])
214 if (self.info['guessed contact slope'] != None
215 and slope_guess_ratio < self.info['min slope guess ratio']):
217 'Too far (contact slope off guess by %g less than %g)'
218 % (slope_guess_ratio, self.info['min slope guess ratio']))
222 class VelocityClampPlugin (Plugin):
224 super(VelocityClampPlugin, self).__init__(name='vclamp')
226 SurfaceContactCommand(self), ForceCommand(self),
227 CantileverAdjustedExtensionCommand(self), FlattenCommand(self),
230 def default_settings(self):
232 Setting(section=self.setting_section, help=self.__doc__),
233 Setting(section=self.setting_section,
234 option='surface contact point algorithm',
236 help='Select the surface contact point algorithm. See the documentation for descriptions of available algorithms.')
240 class SurfaceContactCommand (ColumnAddingCommand):
241 """Automatically determine a block's surface contact point.
243 You can select the contact point algorithm with the creatively
244 named `surface contact point algorithm` configuration setting.
245 Currently available options are:
247 * fmms (:meth:`find_contact_point_fmms`)
248 * ms (:meth:`find_contact_point_ms`)
249 * wtk (:meth:`find_contact_point_wtk`)
251 def __init__(self, plugin):
252 self._wtk_fit_check_arguments = []
253 for argument in SurfacePositionModel._fit_check_arguments:
254 arg = copy.deepcopy(argument)
255 arg._help += ' (wtk model)'
256 self._wtk_fit_check_arguments.append(arg)
257 super(SurfaceContactCommand, self).__init__(
258 name='zero surface contact point',
260 ('distance column', 'z piezo (m)', """
261 Name of the column to use as the surface position input.
263 ('deflection column', 'deflection (m)', """
264 Name of the column to use as the deflection input.
268 ('output distance column', 'surface distance', """
269 Name of the column (without units) to use as the surface position output.
271 ('output deflection column', 'surface deflection', """
272 Name of the column (without units) to use as the deflection output.
276 Argument(name='ignore index', type='int', default=None,
278 Ignore the residual from the non-contact region before the indexed
279 point (for the `wtk` algorithm).
281 Argument(name='ignore after last peak info name',
282 type='string', default=None,
284 As an alternative to 'ignore index', ignore after the last peak in the
285 peak list stored in the `.info` dictionary.
287 Argument(name='distance info name', type='string',
288 default='surface distance offset',
290 Name (without units) for storing the distance offset in the `.info` dictionary.
292 Argument(name='deflection info name', type='string',
293 default='surface deflection offset',
295 Name (without units) for storing the deflection offset in the `.info` dictionary.
297 Argument(name='fit parameters info name', type='string',
298 default='surface deflection offset',
300 Name (without units) for storing fit parameters in the `.info` dictionary.
302 ] + self._wtk_fit_check_arguments,
303 help=self.__doc__, plugin=plugin)
305 def _run(self, hooke, inqueue, outqueue, params):
306 self._add_to_command_stack(params)
307 params = self.__setup_params(hooke=hooke, params=params)
308 block = self._block(hooke=hooke, params=params)
309 dist_data = self._get_column(hooke=hooke, params=params,
310 column_name='distance column')
311 def_data = self._get_column(hooke=hooke, params=params,
312 column_name='deflection column')
313 i,def_offset,ps = self.find_contact_point(
314 params, dist_data, def_data, outqueue)
315 dist_offset = dist_data[i]
316 block.info[params['distance info name']] = dist_offset
317 block.info[params['deflection info name']] = def_offset
318 block.info[params['fit parameters info name']] = ps
319 self._set_column(hooke=hooke, params=params,
320 column_name='output distance column',
321 values=dist_data - dist_offset)
322 self._set_column(hooke=hooke, params=params,
323 column_name='output deflection column',
324 values=def_data - def_offset)
326 def __setup_params(self, hooke, params):
327 name,dist_unit = split_data_label(params['distance column'])
328 name,def_unit = split_data_label(params['deflection column'])
329 params['output distance column'] = join_data_label(
330 params['output distance column'], dist_unit)
331 params['output deflection column'] = join_data_label(
332 params['output deflection column'], def_unit)
333 params['distance info name'] = join_data_label(
334 params['distance info name'], dist_unit)
335 params['deflection info name'] = join_data_label(
336 params['deflection info name'], def_unit)
339 def find_contact_point(self, params, z_data, d_data, outqueue=None):
340 """Railyard for the `find_contact_point_*` family.
342 Uses the `surface contact point algorithm` configuration
343 setting to call the appropriate backend algorithm.
345 fn = getattr(self, 'find_contact_point_%s'
346 % self.plugin.config['surface contact point algorithm'])
347 return fn(params, z_data, d_data, outqueue)
349 def find_contact_point_fmms(self, params, z_data, d_data, outqueue=None):
350 """Algorithm by Francesco Musiani and Massimo Sandal.
356 0) Driver-specific workarounds, e.g. deal with the PicoForce
357 trigger bug by excluding retraction portions with excessive
359 1) Select the second half (non-contact side) of the retraction
361 2) Fit the selection to a line.
362 3) If the fit is not almost horizontal, halve the selection
364 4) Average the selection and use it as a baseline.
365 5) Slide in from the start (contact side) of the retraction
366 curve, until you find a point with greater than baseline
367 deflection. That point is the contact point.
369 if params['curve'].info['filetype'] == 'picoforce':
370 # Take care of the picoforce trigger bug (TODO: example
371 # data file demonstrating the bug). We exclude portions
372 # of the curve that have too much standard deviation.
373 # Yes, a lot of magic is here.
374 check_start = len(d_data)-len(d_data)/20
375 monster_start = len(d_data)
377 # look at the non-contact tail
378 non_monster = d_data[check_start:monster_start]
379 if non_monster.std() < 2e-10: # HACK: hardcoded cutoff
381 else: # move further away from the monster
382 check_start -= len(d_data)/50
383 monster_start -= len(d_data)/50
384 z_data = z_data[:monster_start]
385 d_data = d_data[:monster_start]
387 # take half of the thing to start
388 selection_start = len(d_data)/2
390 z_chunk = z_data[selection_start:]
391 d_chunk = d_data[selection_start:]
392 slope,intercept,r,two_tailed_prob,stderr_of_the_estimate = \
393 scipy.stats.linregress(z_chunk, d_chunk)
394 # We stop if we found an almost-horizontal fit or if we're
395 # getting to small a selection. FIXME: 0.1 and 5./6 here
396 # are "magic numbers" (although reasonable)
397 if (abs(slope) < 0.1 # deflection (m) / surface (m)
398 or selection_start > 5./6*len(d_data)):
400 selection_start += 10
402 d_baseline = d_chunk.mean()
404 # find the first point above the calculated baseline
406 while i < len(d_data) and d_data[i] < ymean:
408 return (i, d_baseline, {})
410 def find_contact_point_ms(self, params, z_data, d_data, outqueue=None):
411 """Algorithm by Massimo Sandal.
415 WTK: At least the commits are by Massimo, and I see no notes
416 attributing the algorithm to anyone else.
422 xext=raw_plot.vectors[0][0]
423 yext=raw_plot.vectors[0][1]
424 xret2=raw_plot.vectors[1][0]
425 yret=raw_plot.vectors[1][1]
427 first_point=[xext[0], yext[0]]
428 last_point=[xext[-1], yext[-1]]
430 #regr=scipy.polyfit(first_point, last_point,1)[0:2]
431 diffx=abs(first_point[0]-last_point[0])
432 diffy=abs(first_point[1]-last_point[1])
434 #using polyfit results in numerical errors. good old algebra.
436 b=first_point[1]-(a*first_point[0])
437 baseline=scipy.polyval((a,b), xext)
439 ysub=[item-basitem for item,basitem in zip(yext,baseline)]
441 contact=ysub.index(min(ysub))
443 return xext,ysub,contact
445 #now, exploit a ClickedPoint instance to calculate index...
447 dummy.absolute_coords=(x_intercept,y_intercept)
448 dummy.find_graph_coords(xret2,yret)
451 return dummy.index, regr, regr_contact
455 def find_contact_point_wtk(self, params, z_data, d_data, outqueue=None):
456 """Algorithm by W. Trevor King.
460 Uses :class:`SurfacePositionModel` internally.
462 reverse = z_data[0] > z_data[-1]
463 if reverse == True: # approaching, contact region on the right
464 d_data = d_data[::-1]
465 s = SurfacePositionModel(d_data, info={
466 'force zero non-contact slope':True},
468 for argument in self._wtk_fit_check_arguments:
469 s.info[argument.name] = params[argument.name]
471 if params['ignore index'] != None:
472 ignore_index = params['ignore index']
473 elif params['ignore after last peak info name'] != None:
474 peaks = z_data.info[params['ignore after last peak info name']]
475 if not len(peaks) > 0:
476 raise Failure('Need at least one peak in %s, not %s'
477 % (params['ignore after last peak info name'],
479 ignore_index = peaks[-1].post_index()
480 if ignore_index != None:
481 s.info['ignore non-contact before index'] = ignore_index
482 offset,contact_slope,surface_index,non_contact_slope = s.fit(
486 'contact slope': contact_slope,
487 'surface index': surface_index,
488 'non-contact slope': non_contact_slope,
491 deflection_offset = offset + contact_slope*surface_index,
493 surface_index = len(d_data)-1-surface_index
494 return (numpy.round(surface_index), deflection_offset, info)
497 class ForceCommand (ColumnAddingCommand):
498 """Convert a deflection column from meters to newtons.
500 def __init__(self, plugin):
501 super(ForceCommand, self).__init__(
502 name='convert distance to force',
504 ('deflection column', 'surface deflection (m)', """
505 Name of the column to use as the deflection input.
509 ('output deflection column', 'deflection', """
510 Name of the column (without units) to use as the deflection output.
514 Argument(name='spring constant info name', type='string',
515 default='spring constant (N/m)',
517 Name of the spring constant in the `.info` dictionary.
520 help=self.__doc__, plugin=plugin)
522 def _run(self, hooke, inqueue, outqueue, params):
523 self._add_to_command_stack(params)
524 params = self.__setup_params(hooke=hooke, params=params)
525 def_data = self._get_column(hooke=hooke, params=params,
526 column_name='deflection column')
527 out = def_data * def_data.info[params['spring constant info name']]
528 self._set_column(hooke=hooke, params=params,
529 column_name='output deflection column',
532 def __setup_params(self, hooke, params):
533 name,in_unit = split_data_label(params['deflection column'])
534 out_unit = 'N' # HACK: extract target units from k_unit.
535 params['output deflection column'] = join_data_label(
536 params['output deflection column'], out_unit)
537 name,k_unit = split_data_label(params['spring constant info name'])
538 expected_k_unit = '%s/%s' % (out_unit, in_unit)
539 if k_unit != expected_k_unit:
540 raise Failure('Cannot convert from %s to %s with %s'
541 % (params['deflection column'],
542 params['output deflection column'],
543 params['spring constant info name']))
547 class CantileverAdjustedExtensionCommand (ColumnAddingCommand):
548 """Remove cantilever extension from a total extension column.
550 If `distance column` and `deflection column` have the same units
551 (e.g. `z piezo (m)` and `deflection (m)`), `spring constant info
552 name` is ignored and a deflection/distance conversion factor of
555 def __init__(self, plugin):
556 super(CantileverAdjustedExtensionCommand, self).__init__(
557 name='remove cantilever from extension',
559 ('distance column', 'surface distance (m)', """
560 Name of the column to use as the surface position input.
562 ('deflection column', 'deflection (N)', """
563 Name of the column to use as the deflection input.
567 ('output distance column', 'cantilever adjusted extension', """
568 Name of the column (without units) to use as the surface position output.
572 Argument(name='spring constant info name', type='string',
573 default='spring constant (N/m)',
575 Name of the spring constant in the `.info` dictionary.
578 help=self.__doc__, plugin=plugin)
580 def _run(self, hooke, inqueue, outqueue, params):
581 self._add_to_command_stack(params)
582 params = self.__setup_params(hooke=hooke, params=params)
583 def_data = self._get_column(hooke=hooke, params=params,
584 column_name='deflection column')
585 dist_data = self._get_column(hooke=hooke, params=params,
586 column_name='distance column')
587 if params['spring constant info name'] == None:
588 k = 1.0 # distance and deflection in the same units
590 k = def_data.info[params['spring constant info name']]
591 self._set_column(hooke=hooke, params=params,
592 column_name='output distance column',
593 values=dist_data - def_data / k)
595 def __setup_params(self, hooke, params):
596 name,dist_unit = split_data_label(params['distance column'])
597 name,def_unit = split_data_label(params['deflection column'])
598 params['output distance column'] = join_data_label(
599 params['output distance column'], dist_unit)
600 if dist_unit == def_unit:
601 params['spring constant info name'] == None
603 name,k_unit = split_data_label(params['spring constant info name'])
604 expected_k_unit = '%s/%s' % (def_unit, dist_unit)
605 if k_unit != expected_k_unit:
606 raise Failure('Cannot convert from %s to %s with %s'
607 % (params['deflection column'],
608 params['output distance column'],
609 params['spring constant info name']))
613 class FlattenCommand (ColumnAddingCommand):
614 """Flatten a deflection column.
616 Subtracts a polynomial fit from the non-contact part of the curve
617 to flatten it. The best polynomial fit is chosen among
618 polynomials of degree 1 to `max degree`.
620 .. todo: Why does flattening use a polynomial fit and not a sinusoid?
621 Isn't most of the oscillation due to laser interference?
622 See Jaschke 1995 ( 10.1063/1.1146018 )
623 and the figure 4 caption of Weisenhorn 1992 ( 10.1103/PhysRevB.45.11226 )
625 def __init__(self, plugin):
626 super(FlattenCommand, self).__init__(
627 name='polynomial flatten',
629 ('distance column', 'surface distance (m)', """
630 Name of the column to use as the surface position input.
632 ('deflection column', 'deflection (N)', """
633 Name of the column to use as the deflection input.
637 ('output deflection column', 'flattened deflection', """
638 Name of the column (without units) to use as the deflection output.
642 Argument(name='degree', type='int',
645 Order of the polynomial used for flattening. Using values greater
646 than one usually doesn't help and can give artifacts. However, it
647 could be useful too. (TODO: Back this up with some theory...)
649 Argument(name='fit info name', type='string',
650 default='flatten fit',
652 Name of the flattening information in the `.info` dictionary.
655 help=self.__doc__, plugin=plugin)
657 def _run(self, hooke, inqueue, outqueue, params):
658 self._add_to_command_stack(params)
659 params = self.__setup_params(hooke=hooke, params=params)
660 block = self._block(hooke=hooke, params=params)
661 dist_data = self._get_column(hooke=hooke, params=params,
662 column_name='distance column')
663 def_data = self._get_column(hooke=hooke, params=params,
664 column_name='deflection column')
665 degree = params['degree']
667 indices = numpy.argwhere(mask)
668 if len(indices) == 0:
669 raise Failure('no positive distance values in %s'
670 % params['distance column'])
671 dist_nc = dist_data[indices].flatten()
672 def_nc = def_data[indices].flatten()
675 poly_values = scipy.polyfit(dist_nc, def_nc, degree)
676 def_nc_fit = scipy.polyval(poly_values, dist_nc)
678 raise Failure('failed to flatten with a degree %d polynomial: %s'
680 error = numpy.sqrt((def_nc_fit-def_nc)**2).sum() / len(def_nc)
681 block.info[params['fit info name']] = {
684 'polynomial values':poly_values,
687 - numpy.invert(mask)*scipy.polyval(poly_values[-1:], dist_data)
688 - mask*scipy.polyval(poly_values, dist_data))
689 self._set_column(hooke=hooke, params=params,
690 column_name='output deflection column',
693 def __setup_params(self, hooke, params):
694 d_name,d_unit = split_data_label(params['deflection column'])
695 params['output deflection column'] = join_data_label(
696 params['output deflection column'], d_unit)