hooke/plugin/vclamp.py

   1 # Copyright (C) 2008-2012 Alberto Gomez-Casado <a.gomezcasado@tnw.utwente.nl>
   2 #                         Marco Brucale <marco.brucale@unibo.it>
   3 #                         Massimo Sandal <devicerandom@gmail.com>
   4 #                         W. Trevor King <wking@tremily.us>
   5 #
   6 # This file is part of Hooke.
   7 #
   8 # Hooke is free software: you can redistribute it and/or modify it under the
   9 # terms of the GNU Lesser General Public License as published by the Free
  10 # Software Foundation, either version 3 of the License, or (at your option) any
  11 # later version.
  12 #
  13 # Hooke is distributed in the hope that it will be useful, but WITHOUT ANY
  14 # WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
  15 # A PARTICULAR PURPOSE.  See the GNU Lesser General Public License for more
  16 # details.
  17 #
  18 # You should have received a copy of the GNU Lesser General Public License
  19 # along with Hooke.  If not, see <http://www.gnu.org/licenses/>.
  20
  21 """The ``vclamp`` module provides :class:`VelocityClampPlugin` and
  22 several associated :class:`hooke.command.Command`\s for handling
  23 common velocity clamp analysis tasks.
  24 """
  25
  26 import copy
  27 import logging
  28
  29 import numpy
  30 import scipy
  31
  32 from ..command import Argument, Failure, NullQueue
  33 from ..config import Setting
  34 from ..curve import Data
  35 from ..util.fit import PoorFit, ModelFitter
  36 from ..util.si import join_data_label, split_data_label
  37 from . import Plugin
  38 from .curve import ColumnAddingCommand
  39
  40
  41 class SurfacePositionModel (ModelFitter):
  42     """Bilinear surface position model.
  43
  44     The bilinear model is symmetric, but the parameter guessing and
  45     sanity checks assume the contact region occurs for lower indicies
  46     ("left of") the non-contact region.  We also assume that
  47     tip-surface attractions produce positive deflections.
  48
  49     Notes
  50     -----
  51     Algorithm borrowed from WTK's `piezo package`_, specifically
  52     from :func:`piezo.z_piezo_utils.analyzeSurfPosData`.
  53
  54     .. _piezo package:
  55       http://www.physics.drexel.edu/~wking/code/git/git.php?p=piezo.git
  56
  57     Fits the data to the bilinear :method:`model`.
  58
  59     In order for this model to produce a satisfactory fit, there
  60     should be enough data in the off-surface region that interactions
  61     due to proteins, etc. will not seriously skew the fit in the
  62     off-surface region.  If you don't have much of a tail, you can set
  63     the `info` dict's `ignore non-contact before index` parameter to
  64     the index of the last surface- or protein-related feature.
  65     """
  66     _fit_check_arguments = [
  67         Argument('min slope ratio', type='float', default=10.,
  68                  help="""
  69 Minimum `contact-slope/non-contact-slope` ratio for a "good" fit.
  70 """.strip()),
  71         Argument('min contact fraction', type='float', default=0.02,
  72                  help="""
  73 Minimum fraction of points in the contact region for a "good" fit.
  74 """.strip()),
  75         Argument('max contact fraction', type='float', default=0.98,
  76                  help="""
  77 Maximum fraction of points in the contact region for a "good" fit.
  78 """.strip()),
  79         Argument('min slope guess ratio', type='float', default=0.5,
  80                  help="""
  81 Minimum `fit-contact-slope/guessed-contact-slope` ratio for a "good" fit.
  82 """.strip()),
  83         ]
  84
  85     def model(self, params):
  86         """A continuous, bilinear model.
  87
  88         Notes
  89         -----
  90         .. math::
  91
  92           y = \begin{cases}
  93             p_0 + p_1 x                 & \text{if $x <= p_2$}, \\
  94             p_0 + p_1 p_2 + p_3 (x-p_2) & \text{if $x >= p_2$}.
  95               \end{cases}
  96
  97         Where :math:`p_0` is a vertical offset, :math:`p_1` is the slope
  98         of the first region, :math:`p_2` is the transition location, and
  99         :math:`p_3` is the slope of the second region.
 100         """
 101         p = params  # convenient alias
 102         rNC_ignore = self.info['ignore non-contact before index']
 103         if self.info['force zero non-contact slope'] is True:
 104             p = list(p)
 105             p.append(0.)  # restore the non-contact slope parameter
 106         r2 = numpy.round(abs(p[2]))
 107         if r2 >= 1:
 108             r2 = min(r2, len(self._model_data))
 109             self._model_data[:r2] = p[0] + p[1] * numpy.arange(r2)
 110         if r2 < len(self._data)-1:
 111             self._model_data[r2:] = \
 112                 p[0] + p[1]*p[2] + p[3] * numpy.arange(len(self._data)-r2)
 113             if r2 < rNC_ignore:
 114                 self._model_data[r2:rNC_ignore] = self._data[r2:rNC_ignore]
 115         return self._model_data
 116
 117     def set_data(self, data, info=None, *args, **kwargs):
 118         super(SurfacePositionModel, self).set_data(data, info, *args, **kwargs)
 119         if info == None:
 120             info = {}
 121         if getattr(self, 'info', None) == None:
 122             self.info = {}
 123         self.info.update(info)
 124         for key,value in [
 125             ('force zero non-contact slope', False),
 126             ('ignore non-contact before index', -1),
 127             ('min position', 0),  # Store postions etc. to avoid recalculating.
 128             ('max position', len(data)),
 129             ('max deflection', data.max()),
 130             ('min deflection', data.min()),
 131             ]:
 132             if key not in self.info:
 133                 self.info[key] = value
 134         for key,value in [
 135             ('position range',
 136              self.info['max position'] - self.info['min position']),
 137             ('deflection range',
 138              self.info['max deflection'] - self.info['min deflection']),
 139             ]:
 140             if key not in self.info:
 141                 self.info[key] = value
 142         for argument in self._fit_check_arguments:
 143             if argument.name not in self.info:
 144                 self.info[argument.name] = argument.default
 145
 146     def guess_initial_params(self, outqueue=None):
 147         """Guess the initial parameters.
 148
 149         Notes
 150         -----
 151         We guess initial parameters such that the offset (:math:`p_1`)
 152         matches the minimum deflection, the kink (:math:`p_2`) occurs
 153         at the first point that the deflection crosses the middle of
 154         its range, the initial (contact) slope (:math:`p_0`) produces
 155         the right-most deflection at the kink point, and the final
 156         (non-contact) slope (:math:`p_3`) is zero.
 157
 158         In the event of a tie, :meth:`argmax` returns the lowest index
 159         to the maximum value.
 160         >>> (numpy.arange(10) >= 5).argmax()
 161         5
 162         """
 163         left_offset = self.info['min deflection']
 164         middle_deflection = (self.info['min deflection']
 165                              + self.info['deflection range']/2.)
 166         kink_position = 2*(self._data > middle_deflection).argmax()
 167         if kink_position == 0:
 168             # jump vibration at the start of the retraction?
 169             start = int(min(max(20, 0.01 * len(self._data)), 0.5*len(self._data)))
 170             std = self._data[:start].std()
 171             left_offset = self._data[start].mean()
 172             stop = start
 173             while abs(self._data[stop] - left_offset) < 3*std:
 174                 stop += 1
 175             left_slope = (self._data[stop-start:stop].mean()
 176                           - left_offset) / (stop-start)
 177             left_offset -= left_slope * start/2
 178             kink_position = (self._data[-1] - left_offset)/left_slope
 179         else:
 180             left_slope = (self._data[-1] - self.info['min deflection']
 181                           )/kink_position
 182         right_slope = 0
 183         self.info['guessed contact slope'] = left_slope
 184         params = [left_offset, left_slope, kink_position, right_slope]
 185         if self.info['force zero non-contact slope'] == True:
 186             params = params[:-1]
 187         return params
 188
 189     def fit(self, *args, **kwargs):
 190         """Fit the model to the data.
 191
 192         Notes
 193         -----
 194         We change the `epsfcn` default from :func:`scipy.optimize.leastsq`'s
 195         `0` to `1e-3`, so the initial Jacobian estimate takes larger steps,
 196         which helps avoid being trapped in noise-generated local minima.
 197         """
 198         self.info['guessed contact slope'] = None
 199         if 'epsfcn' not in kwargs:
 200             kwargs['epsfcn'] = 1e-3  # take big steps to estimate the Jacobian
 201         params = super(SurfacePositionModel, self).fit(*args, **kwargs)
 202         params[2] = abs(params[2])
 203         if self.info['force zero non-contact slope'] == True:
 204             params = list(params)
 205             params.append(0.)  # restore the non-contact slope parameter
 206
 207         # check that the fit is reasonable, see the :meth:`model` docstring
 208         # for parameter descriptions.
 209         slope_ratio = abs(params[1]/params[3])
 210         if slope_ratio < self.info['min slope ratio']:
 211             raise PoorFit(
 212                'Slope in non-contact region, or no slope in contact (slope ratio %g less than %g)'
 213                % (slope_ratio, self.info['min slope ratio']))
 214         contact_fraction = ((params[2]-self.info['min position'])
 215                             /self.info['position range'])
 216         if contact_fraction < self.info['min contact fraction']:
 217             raise PoorFit(
 218                 'No kink (contact fraction %g less than %g)'
 219                 % (contact_fraction, self.info['min contact fraction']))
 220         if contact_fraction > self.info['max contact fraction']:
 221             raise PoorFit(
 222                 'No kink (contact fraction %g greater than %g)'
 223                 % (contact_fraction, self.info['max contact fraction']))
 224         slope_guess_ratio = abs(params[1]/self.info['guessed contact slope'])
 225         if (self.info['guessed contact slope'] != None
 226             and slope_guess_ratio < self.info['min slope guess ratio']):
 227             raise PoorFit(
 228                 'Too far (contact slope off guess by %g less than %g)'
 229                 % (slope_guess_ratio, self.info['min slope guess ratio']))
 230         return params
 231
 232
 233 class VelocityClampPlugin (Plugin):
 234     def __init__(self):
 235         super(VelocityClampPlugin, self).__init__(name='vclamp')
 236         self._commands = [
 237             SurfaceContactCommand(self), ForceCommand(self),
 238             CantileverAdjustedExtensionCommand(self), FlattenCommand(self),
 239             ]
 240
 241     def default_settings(self):
 242         return [
 243             Setting(section=self.setting_section, help=self.__doc__),
 244             Setting(section=self.setting_section,
 245                     option='surface contact point algorithm',
 246                     value='wtk',
 247                     help='Select the surface contact point algorithm.  See the documentation for descriptions of available algorithms.')
 248             ]
 249
 250
 251 class SurfaceContactCommand (ColumnAddingCommand):
 252     """Automatically determine a block's surface contact point.
 253
 254     You can select the contact point algorithm with the creatively
 255     named `surface contact point algorithm` configuration setting.
 256     Currently available options are:
 257
 258     * fmms (:meth:`find_contact_point_fmms`)
 259     * ms (:meth:`find_contact_point_ms`)
 260     * wtk (:meth:`find_contact_point_wtk`)
 261     """
 262     def __init__(self, plugin):
 263         self._wtk_fit_check_arguments = []
 264         for argument in SurfacePositionModel._fit_check_arguments:
 265             arg = copy.deepcopy(argument)
 266             arg._help += '  (wtk model)'
 267             self._wtk_fit_check_arguments.append(arg)
 268         super(SurfaceContactCommand, self).__init__(
 269             name='zero surface contact point',
 270             columns=[
 271                 ('distance column', 'z piezo (m)', """
 272 Name of the column to use as the surface position input.
 273 """.strip()),
 274                 ('deflection column', 'deflection (m)', """
 275 Name of the column to use as the deflection input.
 276 """.strip()),
 277                 ],
 278             new_columns=[
 279                 ('output distance column', 'surface distance', """
 280 Name of the column (without units) to use as the surface position output.
 281 """.strip()),
 282                 ('output deflection column', 'surface deflection', """
 283 Name of the column (without units) to use as the deflection output.
 284 """.strip()),
 285                 ],
 286             arguments=[
 287                 Argument(name='ignore index', type='int', default=None,
 288                          help="""
 289 Ignore the residual from the non-contact region before the indexed
 290 point (for the `wtk` algorithm).
 291 """.strip()),
 292                 Argument(name='ignore after last peak info name',
 293                          type='string', default=None,
 294                          help="""
 295 As an alternative to 'ignore index', ignore after the last peak in the
 296 peak list stored in the `.info` dictionary.
 297 """.strip()),
 298                 Argument(name='force zero non-contact slope', type='bool',
 299                          default=False, count=1,
 300                          help="""
 301 Fix the fitted non-contact slope at zero.
 302 """.strip()),
 303                 Argument(name='distance info name', type='string',
 304                          default='surface distance offset',
 305                          help="""
 306 Name (without units) for storing the distance offset in the `.info` dictionary.
 307 """.strip()),
 308                 Argument(name='deflection info name', type='string',
 309                          default='surface deflection offset',
 310                          help="""
 311 Name (without units) for storing the deflection offset in the `.info` dictionary.
 312 """.strip()),
 313                 Argument(name='fit parameters info name', type='string',
 314                          default='surface deflection offset',
 315                          help="""
 316 Name (without units) for storing fit parameters in the `.info` dictionary.
 317 """.strip()),
 318                 ] + self._wtk_fit_check_arguments,
 319             help=self.__doc__, plugin=plugin)
 320
 321     def _run(self, hooke, inqueue, outqueue, params):
 322         self._add_to_command_stack(params)
 323         params = self._setup_params(hooke=hooke, params=params)
 324         block = self._block(hooke=hooke, params=params)
 325         dist_data = self._get_column(hooke=hooke, params=params,
 326                                      column_name='distance column')
 327         def_data = self._get_column(hooke=hooke, params=params,
 328                                     column_name='deflection column')
 329         i,def_offset,ps = self.find_contact_point(
 330             params, dist_data, def_data, outqueue)
 331         dist_offset = dist_data[i]
 332         block.info[params['distance info name']] = dist_offset
 333         block.info[params['deflection info name']] = def_offset
 334         block.info[params['fit parameters info name']] = ps
 335         self._set_column(hooke=hooke, params=params,
 336                          column_name='output distance column',
 337                          values=dist_data - dist_offset)
 338         self._set_column(hooke=hooke, params=params,
 339                          column_name='output deflection column',
 340                          values=def_data - def_offset)
 341
 342     def _setup_params(self, hooke, params):
 343         name,dist_unit = split_data_label(params['distance column'])
 344         name,def_unit = split_data_label(params['deflection column'])
 345         params['output distance column'] = join_data_label(
 346             params['output distance column'], dist_unit)
 347         params['output deflection column'] = join_data_label(
 348             params['output deflection column'], def_unit)
 349         params['distance info name'] = join_data_label(
 350             params['distance info name'], dist_unit)
 351         params['deflection info name'] = join_data_label(
 352             params['deflection info name'], def_unit)
 353         return params
 354
 355     def find_contact_point(self, params, z_data, d_data, outqueue=None):
 356         """Railyard for the `find_contact_point_*` family.
 357
 358         Uses the `surface contact point algorithm` configuration
 359         setting to call the appropriate backend algorithm.
 360         """
 361         fn = getattr(self, 'find_contact_point_%s'
 362                      % self.plugin.config['surface contact point algorithm'])
 363         return fn(params, z_data, d_data, outqueue)
 364
 365     def find_contact_point_fmms(self, params, z_data, d_data, outqueue=None):
 366         """Algorithm by Francesco Musiani and Massimo Sandal.
 367
 368         Notes
 369         -----
 370         Algorithm:
 371
 372         0) Driver-specific workarounds, e.g. deal with the PicoForce
 373           trigger bug by excluding retraction portions with excessive
 374           deviation.
 375         1) Select the second half (non-contact side) of the retraction
 376           curve.
 377         2) Fit the selection to a line.
 378         3) If the fit is not almost horizontal, halve the selection
 379           and retrun to (2).
 380         4) Average the selection and use it as a baseline.
 381         5) Slide in from the start (contact side) of the retraction
 382         curve, until you find a point with greater than baseline
 383         deflection.  That point is the contact point.
 384         """
 385         if params['curve'].info['filetype'] == 'picoforce':
 386             # Take care of the picoforce trigger bug (TODO: example
 387             # data file demonstrating the bug).  We exclude portions
 388             # of the curve that have too much standard deviation.
 389             # Yes, a lot of magic is here.
 390             check_start = len(d_data)-len(d_data)/20
 391             monster_start = len(d_data)
 392             while True:
 393                 # look at the non-contact tail
 394                 non_monster = d_data[check_start:monster_start]
 395                 if non_monster.std() < 2e-10: # HACK: hardcoded cutoff
 396                     break
 397                 else: # move further away from the monster
 398                     check_start -= len(d_data)/50
 399                     monster_start -= len(d_data)/50
 400             z_data = z_data[:monster_start]
 401             d_data = d_data[:monster_start]
 402
 403         # take half of the thing to start
 404         selection_start = len(d_data)/2
 405         while True:
 406             z_chunk = z_data[selection_start:]
 407             d_chunk = d_data[selection_start:]
 408             slope,intercept,r,two_tailed_prob,stderr_of_the_estimate = \
 409                 scipy.stats.linregress(z_chunk, d_chunk)
 410             # We stop if we found an almost-horizontal fit or if we're
 411             # getting to small a selection.  FIXME: 0.1 and 5./6 here
 412             # are "magic numbers" (although reasonable)
 413             if (abs(slope) < 0.1  # deflection (m) / surface (m)
 414                 or selection_start > 5./6*len(d_data)):
 415                 break
 416             selection_start += 10
 417
 418         d_baseline = d_chunk.mean()
 419
 420         # find the first point above the calculated baseline
 421         i = 0
 422         while i < len(d_data) and d_data[i] < ymean:
 423             i += 1
 424         return (i, d_baseline, {})
 425
 426     def find_contact_point_ms(self, params, z_data, d_data, outqueue=None):
 427         """Algorithm by Massimo Sandal.
 428
 429         Notes
 430         -----
 431         WTK: At least the commits are by Massimo, and I see no notes
 432         attributing the algorithm to anyone else.
 433
 434         Algorithm:
 435
 436         * ?
 437         """
 438         xext=raw_plot.vectors[0][0]
 439         yext=raw_plot.vectors[0][1]
 440         xret2=raw_plot.vectors[1][0]
 441         yret=raw_plot.vectors[1][1]
 442
 443         first_point=[xext[0], yext[0]]
 444         last_point=[xext[-1], yext[-1]]
 445
 446         #regr=scipy.polyfit(first_point, last_point,1)[0:2]
 447         diffx=abs(first_point[0]-last_point[0])
 448         diffy=abs(first_point[1]-last_point[1])
 449
 450         #using polyfit results in numerical errors. good old algebra.
 451         a=diffy/diffx
 452         b=first_point[1]-(a*first_point[0])
 453         baseline=scipy.polyval((a,b), xext)
 454
 455         ysub=[item-basitem for item,basitem in zip(yext,baseline)]
 456
 457         contact=ysub.index(min(ysub))
 458
 459         return xext,ysub,contact
 460
 461         #now, exploit a ClickedPoint instance to calculate index...
 462         dummy=ClickedPoint()
 463         dummy.absolute_coords=(x_intercept,y_intercept)
 464         dummy.find_graph_coords(xret2,yret)
 465
 466         if debug:
 467             return dummy.index, regr, regr_contact
 468         else:
 469             return dummy.index
 470
 471     def find_contact_point_wtk(self, params, z_data, d_data, outqueue=None):
 472         """Algorithm by W. Trevor King.
 473
 474         Notes
 475         -----
 476         Uses :class:`SurfacePositionModel` internally.
 477         """
 478         reverse = z_data[0] > z_data[-1]
 479         if reverse == True:    # approaching, contact region on the right
 480             d_data = d_data[::-1]
 481         s = SurfacePositionModel(d_data, info={
 482                 'force zero non-contact slope':
 483                     params['force zero non-contact slope']},
 484                                  rescale=True)
 485         for argument in self._wtk_fit_check_arguments:
 486             s.info[argument.name] = params[argument.name]
 487         ignore_index = None
 488         if params['ignore index'] != None:
 489             ignore_index = params['ignore index']
 490         elif params['ignore after last peak info name'] != None:
 491             peaks = z_data.info[params['ignore after last peak info name']]
 492             if not len(peaks) > 0:
 493                 raise Failure('Need at least one peak in %s, not %s'
 494                               % (params['ignore after last peak info name'],
 495                                  peaks))
 496             ignore_index = peaks[-1].post_index()
 497         if ignore_index != None:
 498             s.info['ignore non-contact before index'] = ignore_index
 499         offset,contact_slope,surface_index,non_contact_slope = s.fit(
 500             outqueue=outqueue)
 501         deflection_offset = offset + contact_slope*surface_index
 502         delta_pos_per_point = z_data[1] - z_data[0]
 503         contact_slope /= delta_pos_per_point  # ddef/point -> ddev/dpos
 504         non_contact_slope /= delta_pos_per_point
 505         info = {
 506             'offset': offset,
 507             'contact slope': contact_slope,
 508             'surface index': surface_index,
 509             'non-contact slope': non_contact_slope,
 510             'reversed': reverse,
 511             }
 512         if reverse == True:
 513             surface_index = len(d_data)-1-surface_index
 514         return (numpy.round(surface_index), deflection_offset, info)
 515
 516
 517 class ForceCommand (ColumnAddingCommand):
 518     """Convert a deflection column from meters to newtons.
 519     """
 520     def __init__(self, plugin):
 521         super(ForceCommand, self).__init__(
 522             name='convert distance to force',
 523             columns=[
 524                 ('deflection column', 'surface deflection (m)', """
 525 Name of the column to use as the deflection input.
 526 """.strip()),
 527                 ],
 528             new_columns=[
 529                 ('output deflection column', 'deflection', """
 530 Name of the column (without units) to use as the deflection output.
 531 """.strip()),
 532                 ],
 533             arguments=[
 534                 Argument(name='spring constant info name', type='string',
 535                          default='spring constant (N/m)',
 536                          help="""
 537 Name of the spring constant in the `.info` dictionary.
 538 """.strip()),
 539                 ],
 540             help=self.__doc__, plugin=plugin)
 541
 542     def _run(self, hooke, inqueue, outqueue, params):
 543         self._add_to_command_stack(params)
 544         params = self._setup_params(hooke=hooke, params=params)
 545         # TODO: call .curve.ScaledColumnAdditionCommand
 546         def_data = self._get_column(hooke=hooke, params=params,
 547                                     column_name='deflection column')
 548         out = def_data * def_data.info[params['spring constant info name']]
 549         self._set_column(hooke=hooke, params=params,
 550                          column_name='output deflection column',
 551                          values=out)
 552
 553     def _setup_params(self, hooke, params):
 554         name,in_unit = split_data_label(params['deflection column'])
 555         out_unit = 'N'  # HACK: extract target units from k_unit.
 556         params['output deflection column'] = join_data_label(
 557             params['output deflection column'], out_unit)
 558         name,k_unit = split_data_label(params['spring constant info name'])
 559         expected_k_unit = '%s/%s' % (out_unit, in_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 deflection column'],
 564                              params['spring constant info name']))
 565         return params
 566
 567
 568 class CantileverAdjustedExtensionCommand (ColumnAddingCommand):
 569     """Remove cantilever extension from a total extension column.
 570
 571     If `distance column` and `deflection column` have the same units
 572     (e.g. `z piezo (m)` and `deflection (m)`), `spring constant info
 573     name` is ignored and a deflection/distance conversion factor of
 574     one is used.
 575     """
 576     def __init__(self, plugin):
 577         super(CantileverAdjustedExtensionCommand, self).__init__(
 578             name='remove cantilever from extension',
 579             columns=[
 580                 ('distance column', 'surface distance (m)', """
 581 Name of the column to use as the surface position input.
 582 """.strip()),
 583                 ('deflection column', 'deflection (N)', """
 584 Name of the column to use as the deflection input.
 585 """.strip()),
 586                 ],
 587             new_columns=[
 588                 ('output distance column', 'cantilever adjusted extension', """
 589 Name of the column (without units) to use as the surface position output.
 590 """.strip()),
 591                 ],
 592             arguments=[
 593                 Argument(name='spring constant info name', type='string',
 594                          default='spring constant (N/m)',
 595                          help="""
 596 Name of the spring constant in the `.info` dictionary.
 597 """.strip()),
 598                 ],
 599             help=self.__doc__, plugin=plugin)
 600
 601     def _run(self, hooke, inqueue, outqueue, params):
 602         self._add_to_command_stack(params)
 603         params = self._setup_params(hooke=hooke, params=params)
 604         def_data = self._get_column(hooke=hooke, params=params,
 605                                     column_name='deflection column')
 606         dist_data = self._get_column(hooke=hooke, params=params,
 607                                      column_name='distance column')
 608         if params['spring constant info name'] == None:
 609             k = 1.0  # distance and deflection in the same units
 610         else:
 611             k = def_data.info[params['spring constant info name']]
 612         self._set_column(hooke=hooke, params=params,
 613                          column_name='output distance column',
 614                          values=dist_data - def_data / k)
 615
 616     def _setup_params(self, hooke, params):
 617         name,dist_unit = split_data_label(params['distance column'])
 618         name,def_unit = split_data_label(params['deflection column'])
 619         params['output distance column'] = join_data_label(
 620             params['output distance column'], dist_unit)
 621         if dist_unit == def_unit:
 622             params['spring constant info name'] == None
 623         else:
 624             name,k_unit = split_data_label(params['spring constant info name'])
 625             expected_k_unit = '%s/%s' % (def_unit, dist_unit)
 626             if k_unit != expected_k_unit:
 627                 raise Failure('Cannot convert from %s to %s with %s'
 628                               % (params['deflection column'],
 629                                  params['output distance column'],
 630                                  params['spring constant info name']))
 631         return params
 632
 633
 634 class FlattenCommand (ColumnAddingCommand):
 635     """Flatten a deflection column.
 636
 637     Subtracts a polynomial fit from the non-contact part of the curve
 638     to flatten it.  The best polynomial fit is chosen among
 639     polynomials of degree 1 to `max degree`.
 640
 641     .. todo:  Why does flattening use a polynomial fit and not a sinusoid?
 642       Isn't most of the oscillation due to laser interference?
 643       See Jaschke 1995 ( 10.1063/1.1146018 )
 644       and the figure 4 caption of Weisenhorn 1992 ( 10.1103/PhysRevB.45.11226 )
 645     """
 646     def __init__(self, plugin):
 647         super(FlattenCommand, self).__init__(
 648             name='polynomial flatten',
 649             columns=[
 650                 ('distance column', 'surface distance (m)', """
 651 Name of the column to use as the surface position input.
 652 """.strip()),
 653                 ('deflection column', 'deflection (N)', """
 654 Name of the column to use as the deflection input.
 655 """.strip()),
 656                 ],
 657             new_columns=[
 658                 ('output deflection column', 'flattened deflection', """
 659 Name of the column (without units) to use as the deflection output.
 660 """.strip()),
 661                 ],
 662             arguments=[
 663                 Argument(name='degree', type='int',
 664                          default=1,
 665                          help="""
 666 Order of the polynomial used for flattening.  Using values greater
 667 than one usually doesn't help and can give artifacts.  However, it
 668 could be useful too.  (TODO: Back this up with some theory...)
 669 """.strip()),
 670                 Argument(name='fit info name', type='string',
 671                          default='flatten fit',
 672                          help="""
 673 Name of the flattening information in the `.info` dictionary.
 674 """.strip()),
 675                 ],
 676             help=self.__doc__, plugin=plugin)
 677
 678     def _run(self, hooke, inqueue, outqueue, params):
 679         self._add_to_command_stack(params)
 680         params = self._setup_params(hooke=hooke, params=params)
 681         block = self._block(hooke=hooke, params=params)
 682         dist_data = self._get_column(hooke=hooke, params=params,
 683                                      column_name='distance column')
 684         def_data = self._get_column(hooke=hooke, params=params,
 685                                     column_name='deflection column')
 686         degree = params['degree']
 687         mask = dist_data > 0
 688         indices = numpy.argwhere(mask)
 689         if len(indices) == 0:
 690             raise Failure('no positive distance values in %s'
 691                           % params['distance column'])
 692         dist_nc = dist_data[indices].flatten()
 693         def_nc = def_data[indices].flatten()
 694
 695         try:
 696             poly_values = scipy.polyfit(dist_nc, def_nc, degree)
 697             def_nc_fit = scipy.polyval(poly_values, dist_nc)
 698         except Exception, e:
 699             raise Failure('failed to flatten with a degree %d polynomial: %s'
 700                           % (degree, e))
 701         error = numpy.sqrt((def_nc_fit-def_nc)**2).sum() / len(def_nc)
 702         block.info[params['fit info name']] = {
 703             'error':error,
 704             'degree':degree,
 705             'polynomial values':poly_values,
 706             }
 707         out = (def_data
 708                - numpy.invert(mask)*scipy.polyval(poly_values[-1:], dist_data)
 709                - mask*scipy.polyval(poly_values, dist_data))
 710         self._set_column(hooke=hooke, params=params,
 711                          column_name='output deflection column',
 712                          values=out)
 713
 714     def _setup_params(self, hooke, params):
 715         d_name,d_unit = split_data_label(params['deflection column'])
 716         params['output deflection column'] = join_data_label(
 717             params['output deflection column'], d_unit)
 718         return params