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 def model(self, params):
69 """A continuous, bilinear model.
76 p_0 + p_1 x & \text{if $x <= p_2$}, \\
77 p_0 + p_1 p_2 + p_3 (x-p_2) & \text{if $x >= p_2$}.
80 Where :math:`p_0` is a vertical offset, :math:`p_1` is the slope
81 of the first region, :math:`p_2` is the transition location, and
82 :math:`p_3` is the slope of the second region.
84 p = params # convenient alias
85 rNC_ignore = self.info['ignore non-contact before index']
86 if self.info['force zero non-contact slope'] == True:
88 p.append(0.) # restore the non-contact slope parameter
89 r2 = numpy.round(abs(p[2]))
91 self._model_data[:r2] = p[0] + p[1] * numpy.arange(r2)
92 if r2 < len(self._data)-1:
93 self._model_data[r2:] = \
94 p[0] + p[1]*p[2] + p[3] * numpy.arange(len(self._data)-r2)
96 self._model_data[r2:rNC_ignore] = self._data[r2:rNC_ignore]
97 return self._model_data
99 def set_data(self, data, info=None, *args, **kwargs):
100 super(SurfacePositionModel, self).set_data(data, info, *args, **kwargs)
101 if self.info == None:
104 ('force zero non-contact slope', False),
105 ('ignore non-contact before index', -1),
106 ('min position', 0), # Store postions etc. to avoid recalculating.
107 ('max position', len(data)),
108 ('max deflection', data.max()),
109 ('min deflection', data.min()),
111 if key not in self.info:
112 self.info[key] = value
115 self.info['max position'] - self.info['min position']),
117 self.info['max deflection'] - self.info['min deflection']),
119 if key not in self.info:
120 self.info[key] = value
122 def guess_initial_params(self, outqueue=None):
123 """Guess the initial parameters.
127 We guess initial parameters such that the offset (:math:`p_1`)
128 matches the minimum deflection, the kink (:math:`p_2`) occurs in
129 the middle of the data, the initial (contact) slope (:math:`p_0`)
130 produces the maximum deflection at the left-most point, and the
131 final (non-contact) slope (:math:`p_3`) is zero.
133 left_offset = self.info['min deflection']
134 left_slope = 2*(self.info['deflection range']
135 /self.info['position range'])
136 kink_position = (self.info['max position']
137 +self.info['min position'])/2.0
139 self.info['guessed contact slope'] = left_slope
140 params = [left_offset, left_slope, kink_position, right_slope]
141 if self.info['force zero non-contact slope'] == True:
145 def fit(self, *args, **kwargs):
146 """Fit the model to the data.
150 We change the `epsfcn` default from :func:`scipy.optimize.leastsq`'s
151 `0` to `1e-3`, so the initial Jacobian estimate takes larger steps,
152 which helps avoid being trapped in noise-generated local minima.
154 self.info['guessed contact slope'] = None
155 if 'epsfcn' not in kwargs:
156 kwargs['epsfcn'] = 1e-3 # take big steps to estimate the Jacobian
157 params = super(SurfacePositionModel, self).fit(*args, **kwargs)
158 params[2] = abs(params[2])
159 if self.info['force zero non-contact slope'] == True:
160 params = list(params)
161 params.append(0.) # restore the non-contact slope parameter
163 # check that the fit is reasonable, see the :meth:`model` docstring
164 # for parameter descriptions. HACK: hardcoded cutoffs.
165 if abs(params[3]*10) > abs(params[1]) :
166 raise PoorFit('Slope in non-contact region, or no slope in contact')
167 if params[2] < self.info['min position']+0.02*self.info['position range']:
169 'No kink (kink %g less than %g, need more space to left)'
171 self.info['min position']+0.02*self.info['position range']))
172 if params[2] > self.info['max position']-0.02*self.info['position range']:
174 'No kink (kink %g more than %g, need more space to right)'
176 self.info['max position']-0.02*self.info['position range']))
177 if (self.info['guessed contact slope'] != None
178 and abs(params[1]) < 0.5 * abs(self.info['guessed contact slope'])):
179 raise PoorFit('Too far (contact slope %g, but expected ~%g'
180 % (params[3], self.info['guessed contact slope']))
184 class VelocityClampPlugin (Plugin):
186 super(VelocityClampPlugin, self).__init__(name='vclamp')
188 SurfaceContactCommand(self), ForceCommand(self),
189 CantileverAdjustedExtensionCommand(self), FlattenCommand(self),
192 def default_settings(self):
194 Setting(section=self.setting_section, help=self.__doc__),
195 Setting(section=self.setting_section,
196 option='surface contact point algorithm',
198 help='Select the surface contact point algorithm. See the documentation for descriptions of available algorithms.')
202 class SurfaceContactCommand (ColumnAddingCommand):
203 """Automatically determine a block's surface contact point.
205 You can select the contact point algorithm with the creatively
206 named `surface contact point algorithm` configuration setting.
207 Currently available options are:
209 * fmms (:meth:`find_contact_point_fmms`)
210 * ms (:meth:`find_contact_point_ms`)
211 * wtk (:meth:`find_contact_point_wtk`)
213 def __init__(self, plugin):
214 super(SurfaceContactCommand, self).__init__(
215 name='zero surface contact point',
217 ('distance column', 'z piezo (m)', """
218 Name of the column to use as the surface position input.
220 ('deflection column', 'deflection (m)', """
221 Name of the column to use as the deflection input.
225 ('output distance column', 'surface distance', """
226 Name of the column (without units) to use as the surface position output.
228 ('output deflection column', 'surface deflection', """
229 Name of the column (without units) to use as the deflection output.
233 Argument(name='ignore index', type='int', default=None,
235 Ignore the residual from the non-contact region before the indexed
236 point (for the `wtk` algorithm).
238 Argument(name='ignore after last peak info name',
239 type='string', default=None,
241 As an alternative to 'ignore index', ignore after the last peak in the
242 peak list stored in the `.info` dictionary.
244 Argument(name='distance info name', type='string',
245 default='surface distance offset',
247 Name (without units) for storing the distance offset in the `.info` dictionary.
249 Argument(name='deflection info name', type='string',
250 default='surface deflection offset',
252 Name (without units) for storing the deflection offset in the `.info` dictionary.
254 Argument(name='fit parameters info name', type='string',
255 default='surface deflection offset',
257 Name (without units) for storing fit parameters in the `.info` dictionary.
260 help=self.__doc__, plugin=plugin)
262 def _run(self, hooke, inqueue, outqueue, params):
263 self._add_to_command_stack(params)
264 params = self.__setup_params(hooke=hooke, params=params)
265 block = self._block(hooke=hooke, params=params)
266 dist_data = self._get_column(hooke=hooke, params=params,
267 column_name='distance column')
268 def_data = self._get_column(hooke=hooke, params=params,
269 column_name='deflection column')
270 i,def_offset,ps = self.find_contact_point(
271 params, dist_data, def_data, outqueue)
272 dist_offset = dist_data[i]
273 block.info[params['distance info name']] = dist_offset
274 block.info[params['deflection info name']] = def_offset
275 block.info[params['fit parameters info name']] = ps
276 self._set_column(hooke=hooke, params=params,
277 column_name='output distance column',
278 values=dist_data - dist_offset)
279 self._set_column(hooke=hooke, params=params,
280 column_name='output deflection column',
281 values=def_data - def_offset)
283 def __setup_params(self, hooke, params):
284 name,dist_unit = split_data_label(params['distance column'])
285 name,def_unit = split_data_label(params['deflection column'])
286 params['output distance column'] = join_data_label(
287 params['output distance column'], dist_unit)
288 params['output deflection column'] = join_data_label(
289 params['output deflection column'], def_unit)
290 params['distance info name'] = join_data_label(
291 params['distance info name'], dist_unit)
292 params['deflection info name'] = join_data_label(
293 params['deflection info name'], def_unit)
296 def find_contact_point(self, params, z_data, d_data, outqueue=None):
297 """Railyard for the `find_contact_point_*` family.
299 Uses the `surface contact point algorithm` configuration
300 setting to call the appropriate backend algorithm.
302 fn = getattr(self, 'find_contact_point_%s'
303 % self.plugin.config['surface contact point algorithm'])
304 return fn(params, z_data, d_data, outqueue)
306 def find_contact_point_fmms(self, params, z_data, d_data, outqueue=None):
307 """Algorithm by Francesco Musiani and Massimo Sandal.
313 0) Driver-specific workarounds, e.g. deal with the PicoForce
314 trigger bug by excluding retraction portions with excessive
316 1) Select the second half (non-contact side) of the retraction
318 2) Fit the selection to a line.
319 3) If the fit is not almost horizontal, halve the selection
321 4) Average the selection and use it as a baseline.
322 5) Slide in from the start (contact side) of the retraction
323 curve, until you find a point with greater than baseline
324 deflection. That point is the contact point.
326 if params['curve'].info['filetype'] == 'picoforce':
327 # Take care of the picoforce trigger bug (TODO: example
328 # data file demonstrating the bug). We exclude portions
329 # of the curve that have too much standard deviation.
330 # Yes, a lot of magic is here.
331 check_start = len(d_data)-len(d_data)/20
332 monster_start = len(d_data)
334 # look at the non-contact tail
335 non_monster = d_data[check_start:monster_start]
336 if non_monster.std() < 2e-10: # HACK: hardcoded cutoff
338 else: # move further away from the monster
339 check_start -= len(d_data)/50
340 monster_start -= len(d_data)/50
341 z_data = z_data[:monster_start]
342 d_data = d_data[:monster_start]
344 # take half of the thing to start
345 selection_start = len(d_data)/2
347 z_chunk = z_data[selection_start:]
348 d_chunk = d_data[selection_start:]
349 slope,intercept,r,two_tailed_prob,stderr_of_the_estimate = \
350 scipy.stats.linregress(z_chunk, d_chunk)
351 # We stop if we found an almost-horizontal fit or if we're
352 # getting to small a selection. FIXME: 0.1 and 5./6 here
353 # are "magic numbers" (although reasonable)
354 if (abs(slope) < 0.1 # deflection (m) / surface (m)
355 or selection_start > 5./6*len(d_data)):
357 selection_start += 10
359 d_baseline = d_chunk.mean()
361 # find the first point above the calculated baseline
363 while i < len(d_data) and d_data[i] < ymean:
365 return (i, d_baseline, {})
367 def find_contact_point_ms(self, params, z_data, d_data, outqueue=None):
368 """Algorithm by Massimo Sandal.
372 WTK: At least the commits are by Massimo, and I see no notes
373 attributing the algorithm to anyone else.
379 xext=raw_plot.vectors[0][0]
380 yext=raw_plot.vectors[0][1]
381 xret2=raw_plot.vectors[1][0]
382 yret=raw_plot.vectors[1][1]
384 first_point=[xext[0], yext[0]]
385 last_point=[xext[-1], yext[-1]]
387 #regr=scipy.polyfit(first_point, last_point,1)[0:2]
388 diffx=abs(first_point[0]-last_point[0])
389 diffy=abs(first_point[1]-last_point[1])
391 #using polyfit results in numerical errors. good old algebra.
393 b=first_point[1]-(a*first_point[0])
394 baseline=scipy.polyval((a,b), xext)
396 ysub=[item-basitem for item,basitem in zip(yext,baseline)]
398 contact=ysub.index(min(ysub))
400 return xext,ysub,contact
402 #now, exploit a ClickedPoint instance to calculate index...
404 dummy.absolute_coords=(x_intercept,y_intercept)
405 dummy.find_graph_coords(xret2,yret)
408 return dummy.index, regr, regr_contact
412 def find_contact_point_wtk(self, params, z_data, d_data, outqueue=None):
413 """Algorithm by W. Trevor King.
417 Uses :class:`SurfacePositionModel` internally.
419 reverse = z_data[0] > z_data[-1]
420 if reverse == True: # approaching, contact region on the right
421 d_data = d_data[::-1]
422 s = SurfacePositionModel(d_data, info={
423 'force zero non-contact slope':True},
426 if params['ignore index'] != None:
427 ignore_index = params['ignore index']
428 elif params['ignore after last peak info name'] != None:
429 peaks = z_data.info[params['ignore after last peak info name']]
430 if not len(peaks) > 0:
431 raise Failure('Need at least one peak in %s, not %s'
432 % (params['ignore after last peak info name'],
434 ignore_index = peaks[-1].post_index()
435 if ignore_index != None:
436 s.info['ignore non-contact before index'] = ignore_index
437 offset,contact_slope,surface_index,non_contact_slope = s.fit(
441 'contact slope': contact_slope,
442 'surface index': surface_index,
443 'non-contact slope': non_contact_slope,
446 deflection_offset = offset + contact_slope*surface_index,
448 surface_index = len(d_data)-1-surface_index
449 return (numpy.round(surface_index), deflection_offset, info)
452 class ForceCommand (ColumnAddingCommand):
453 """Convert a deflection column from meters to newtons.
455 def __init__(self, plugin):
456 super(ForceCommand, self).__init__(
457 name='convert distance to force',
459 ('deflection column', 'surface deflection (m)', """
460 Name of the column to use as the deflection input.
464 ('output deflection column', 'deflection', """
465 Name of the column (without units) to use as the deflection output.
469 Argument(name='spring constant info name', type='string',
470 default='spring constant (N/m)',
472 Name of the spring constant in the `.info` dictionary.
475 help=self.__doc__, plugin=plugin)
477 def _run(self, hooke, inqueue, outqueue, params):
478 self._add_to_command_stack(params)
479 params = self.__setup_params(hooke=hooke, params=params)
480 def_data = self._get_column(hooke=hooke, params=params,
481 column_name='deflection column')
482 out = def_data * def_data.info[params['spring constant info name']]
483 self._set_column(hooke=hooke, params=params,
484 column_name='output deflection column',
487 def __setup_params(self, hooke, params):
488 name,in_unit = split_data_label(params['deflection column'])
489 out_unit = 'N' # HACK: extract target units from k_unit.
490 params['output deflection column'] = join_data_label(
491 params['output deflection column'], out_unit)
492 name,k_unit = split_data_label(params['spring constant info name'])
493 expected_k_unit = '%s/%s' % (out_unit, in_unit)
494 if k_unit != expected_k_unit:
495 raise Failure('Cannot convert from %s to %s with %s'
496 % (params['deflection column'],
497 params['output deflection column'],
498 params['spring constant info name']))
502 class CantileverAdjustedExtensionCommand (ColumnAddingCommand):
503 """Remove cantilever extension from a total extension column.
505 If `distance column` and `deflection column` have the same units
506 (e.g. `z piezo (m)` and `deflection (m)`), `spring constant info
507 name` is ignored and a deflection/distance conversion factor of
510 def __init__(self, plugin):
511 super(CantileverAdjustedExtensionCommand, self).__init__(
512 name='remove cantilever from extension',
514 ('distance column', 'surface distance (m)', """
515 Name of the column to use as the surface position input.
517 ('deflection column', 'deflection (N)', """
518 Name of the column to use as the deflection input.
522 ('output distance column', 'cantilever adjusted extension', """
523 Name of the column (without units) to use as the surface position 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 dist_data = self._get_column(hooke=hooke, params=params,
541 column_name='distance column')
542 if params['spring constant info name'] == None:
543 k = 1.0 # distance and deflection in the same units
545 k = def_data.info[params['spring constant info name']]
546 self._set_column(hooke=hooke, params=params,
547 column_name='output distance column',
548 values=dist_data - def_data / k)
550 def __setup_params(self, hooke, params):
551 name,dist_unit = split_data_label(params['distance column'])
552 name,def_unit = split_data_label(params['deflection column'])
553 params['output distance column'] = join_data_label(
554 params['output distance column'], dist_unit)
555 if dist_unit == def_unit:
556 params['spring constant info name'] == None
558 name,k_unit = split_data_label(params['spring constant info name'])
559 expected_k_unit = '%s/%s' % (def_unit, dist_unit)
560 if k_unit != expected_k_unit:
561 raise Failure('Cannot convert from %s to %s with %s'
562 % (params['deflection column'],
563 params['output distance column'],
564 params['spring constant info name']))
568 class FlattenCommand (ColumnAddingCommand):
569 """Flatten a deflection column.
571 Subtracts a polynomial fit from the non-contact part of the curve
572 to flatten it. The best polynomial fit is chosen among
573 polynomials of degree 1 to `max degree`.
575 .. todo: Why does flattening use a polynomial fit and not a sinusoid?
576 Isn't most of the oscillation due to laser interference?
577 See Jaschke 1995 ( 10.1063/1.1146018 )
578 and the figure 4 caption of Weisenhorn 1992 ( 10.1103/PhysRevB.45.11226 )
580 def __init__(self, plugin):
581 super(FlattenCommand, self).__init__(
582 name='polynomial flatten',
584 ('distance column', 'surface distance (m)', """
585 Name of the column to use as the surface position input.
587 ('deflection column', 'deflection (N)', """
588 Name of the column to use as the deflection input.
592 ('output deflection column', 'flattened deflection', """
593 Name of the column (without units) to use as the deflection output.
597 Argument(name='degree', type='int',
600 Order of the polynomial used for flattening. Using values greater
601 than one usually doesn't help and can give artifacts. However, it
602 could be useful too. (TODO: Back this up with some theory...)
604 Argument(name='fit info name', type='string',
605 default='flatten fit',
607 Name of the flattening information in the `.info` dictionary.
610 help=self.__doc__, plugin=plugin)
612 def _run(self, hooke, inqueue, outqueue, params):
613 self._add_to_command_stack(params)
614 params = self.__setup_params(hooke=hooke, params=params)
615 block = self._block(hooke=hooke, params=params)
616 dist_data = self._get_column(hooke=hooke, params=params,
617 column_name='distance column')
618 def_data = self._get_column(hooke=hooke, params=params,
619 column_name='deflection column')
620 degree = params['degree']
622 indices = numpy.argwhere(mask)
623 if len(indices) == 0:
624 raise Failure('no positive distance values in %s'
625 % params['distance column'])
626 dist_nc = dist_data[indices].flatten()
627 def_nc = def_data[indices].flatten()
630 poly_values = scipy.polyfit(dist_nc, def_nc, degree)
631 def_nc_fit = scipy.polyval(poly_values, dist_nc)
633 raise Failure('failed to flatten with a degree %d polynomial: %s'
635 error = numpy.sqrt((def_nc_fit-def_nc)**2).sum() / len(def_nc)
636 block.info[params['fit info name']] = {
639 'polynomial values':poly_values,
641 out = def_data - mask*scipy.polyval(poly_values, dist_data)
642 self._set_column(hooke=hooke, params=params,
643 column_name='output deflection column',
646 def __setup_params(self, hooke, params):
647 d_name,d_unit = split_data_label(params['deflection column'])
648 params['output deflection column'] = join_data_label(
649 params['output deflection column'], d_unit)