hooke/plugin/vclamp.py

   1 # Copyright (C) 2008-2010 Alberto Gomez-Casado
   2 #                         Fabrizio Benedetti
   3 #                         Marco Brucale
   4 #                         Massimo Sandal <devicerandom@gmail.com>
   5 #                         W. Trevor King <wking@drexel.edu>
   6 #
   7 # This file is part of Hooke.
   8 #
   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.
  13 #
  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.
  18 #
  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/>.
  22
  23 """The ``vclamp`` module provides :class:`VelocityClampPlugin` and
  24 several associated :class:`hooke.command.Command`\s for handling
  25 common velocity clamp analysis tasks.
  26 """
  27
  28 import copy
  29 import logging
  30
  31 import numpy
  32 import scipy
  33
  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
  39 from . import Plugin
  40 from .curve import ColumnAddingCommand
  41
  42
  43 class SurfacePositionModel (ModelFitter):
  44     """Bilinear surface position model.
  45
  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.
  50
  51     Notes
  52     -----
  53     Algorithm borrowed from WTK's `piezo package`_, specifically
  54     from :func:`piezo.z_piezo_utils.analyzeSurfPosData`.
  55
  56     .. _piezo package:
  57       http://www.physics.drexel.edu/~wking/code/git/git.php?p=piezo.git
  58
  59     Fits the data to the bilinear :method:`model`.
  60
  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.
  67     """
  68     _fit_check_arguments = [
  69         Argument('min slope ratio', type='float', default=10.,
  70                  help="""
  71 Minimum `contact-slope/non-contact-slope` ratio for a "good" fit.
  72 """.strip()),
  73         Argument('min contact fraction', type='float', default=0.02,
  74                  help="""
  75 Minimum fraction of points in the contact region for a "good" fit.
  76 """.strip()),
  77         Argument('max contact fraction', type='float', default=0.98,
  78                  help="""
  79 Maximum fraction of points in the contact region for a "good" fit.
  80 """.strip()),
  81         Argument('min slope guess ratio', type='float', default=0.5,
  82                  help="""
  83 Minimum `fit-contact-slope/guessed-contact-slope` ratio for a "good" fit.
  84 """.strip()),
  85         ]
  86
  87     def model(self, params):
  88         """A continuous, bilinear model.
  89
  90         Notes
  91         -----
  92         .. math::
  93
  94           y = \begin{cases}
  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$}.
  97               \end{cases}
  98
  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.
 102         """
 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:
 106             p = list(p)
 107             p.append(0.)  # restore the non-contact slope parameter
 108         r2 = numpy.round(abs(p[2]))
 109         if r2 >= 1:
 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)
 115             if r2 < rNC_ignore:
 116                 self._model_data[r2:rNC_ignore] = self._data[r2:rNC_ignore]
 117         return self._model_data
 118
 119     def set_data(self, data, info=None, *args, **kwargs):
 120         super(SurfacePositionModel, self).set_data(data, info, *args, **kwargs)
 121         if info == None:
 122             info = {}
 123         if getattr(self, 'info', None) == None:
 124             self.info = {}
 125         self.info.update(info)
 126         for key,value in [
 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()),
 133             ]:
 134             if key not in self.info:
 135                 self.info[key] = value
 136         for key,value in [
 137             ('position range',
 138              self.info['max position'] - self.info['min position']),
 139             ('deflection range',
 140              self.info['max deflection'] - self.info['min deflection']),
 141             ]:
 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
 147
 148     def guess_initial_params(self, outqueue=None):
 149         """Guess the initial parameters.
 150
 151         Notes
 152         -----
 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.
 159
 160         In the event of a tie, :meth:`argmax` returns the lowest index
 161         to the maximum value.
 162         >>> (numpy.arange(10) >= 5).argmax()
 163         5
 164         """
 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'])
 170                       /kink_position)
 171         right_slope = 0
 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:
 175             params = params[:-1]
 176         return params
 177
 178     def fit(self, *args, **kwargs):
 179         """Fit the model to the data.
 180
 181         Notes
 182         -----
 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.
 186         """
 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
 195
 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']:
 200             raise PoorFit(
 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']:
 206             raise PoorFit(
 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']:
 210             raise PoorFit(
 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']):
 216             raise PoorFit(
 217                 'Too far (contact slope off guess by %g less than %g)'
 218                 % (slope_guess_ratio, self.info['min slope guess ratio']))
 219         return params
 220
 221
 222 class VelocityClampPlugin (Plugin):
 223     def __init__(self):
 224         super(VelocityClampPlugin, self).__init__(name='vclamp')
 225         self._commands = [
 226             SurfaceContactCommand(self), ForceCommand(self),
 227             CantileverAdjustedExtensionCommand(self), FlattenCommand(self),
 228             ]
 229
 230     def default_settings(self):
 231         return [
 232             Setting(section=self.setting_section, help=self.__doc__),
 233             Setting(section=self.setting_section,
 234                     option='surface contact point algorithm',
 235                     value='wtk',
 236                     help='Select the surface contact point algorithm.  See the documentation for descriptions of available algorithms.')
 237             ]
 238
 239
 240 class SurfaceContactCommand (ColumnAddingCommand):
 241     """Automatically determine a block's surface contact point.
 242
 243     You can select the contact point algorithm with the creatively
 244     named `surface contact point algorithm` configuration setting.
 245     Currently available options are:
 246
 247     * fmms (:meth:`find_contact_point_fmms`)
 248     * ms (:meth:`find_contact_point_ms`)
 249     * wtk (:meth:`find_contact_point_wtk`)
 250     """
 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',
 259             columns=[
 260                 ('distance column', 'z piezo (m)', """
 261 Name of the column to use as the surface position input.
 262 """.strip()),
 263                 ('deflection column', 'deflection (m)', """
 264 Name of the column to use as the deflection input.
 265 """.strip()),
 266                 ],
 267             new_columns=[
 268                 ('output distance column', 'surface distance', """
 269 Name of the column (without units) to use as the surface position output.
 270 """.strip()),
 271                 ('output deflection column', 'surface deflection', """
 272 Name of the column (without units) to use as the deflection output.
 273 """.strip()),
 274                 ],
 275             arguments=[
 276                 Argument(name='ignore index', type='int', default=None,
 277                          help="""
 278 Ignore the residual from the non-contact region before the indexed
 279 point (for the `wtk` algorithm).
 280 """.strip()),
 281                 Argument(name='ignore after last peak info name',
 282                          type='string', default=None,
 283                          help="""
 284 As an alternative to 'ignore index', ignore after the last peak in the
 285 peak list stored in the `.info` dictionary.
 286 """.strip()),
 287                 Argument(name='distance info name', type='string',
 288                          default='surface distance offset',
 289                          help="""
 290 Name (without units) for storing the distance offset in the `.info` dictionary.
 291 """.strip()),
 292                 Argument(name='deflection info name', type='string',
 293                          default='surface deflection offset',
 294                          help="""
 295 Name (without units) for storing the deflection offset in the `.info` dictionary.
 296 """.strip()),
 297                 Argument(name='fit parameters info name', type='string',
 298                          default='surface deflection offset',
 299                          help="""
 300 Name (without units) for storing fit parameters in the `.info` dictionary.
 301 """.strip()),
 302                 ] + self._wtk_fit_check_arguments,
 303             help=self.__doc__, plugin=plugin)
 304
 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)
 325
 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)
 337         return params
 338
 339     def find_contact_point(self, params, z_data, d_data, outqueue=None):
 340         """Railyard for the `find_contact_point_*` family.
 341
 342         Uses the `surface contact point algorithm` configuration
 343         setting to call the appropriate backend algorithm.
 344         """
 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)
 348
 349     def find_contact_point_fmms(self, params, z_data, d_data, outqueue=None):
 350         """Algorithm by Francesco Musiani and Massimo Sandal.
 351
 352         Notes
 353         -----
 354         Algorithm:
 355
 356         0) Driver-specific workarounds, e.g. deal with the PicoForce
 357           trigger bug by excluding retraction portions with excessive
 358           deviation.
 359         1) Select the second half (non-contact side) of the retraction
 360           curve.
 361         2) Fit the selection to a line.
 362         3) If the fit is not almost horizontal, halve the selection
 363           and retrun to (2).
 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.
 368         """
 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)
 376             while True:
 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
 380                     break
 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]
 386
 387         # take half of the thing to start
 388         selection_start = len(d_data)/2
 389         while True:
 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)):
 399                 break
 400             selection_start += 10
 401
 402         d_baseline = d_chunk.mean()
 403
 404         # find the first point above the calculated baseline
 405         i = 0
 406         while i < len(d_data) and d_data[i] < ymean:
 407             i += 1
 408         return (i, d_baseline, {})
 409
 410     def find_contact_point_ms(self, params, z_data, d_data, outqueue=None):
 411         """Algorithm by Massimo Sandal.
 412
 413         Notes
 414         -----
 415         WTK: At least the commits are by Massimo, and I see no notes
 416         attributing the algorithm to anyone else.
 417
 418         Algorithm:
 419
 420         * ?
 421         """
 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]
 426
 427         first_point=[xext[0], yext[0]]
 428         last_point=[xext[-1], yext[-1]]
 429
 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])
 433
 434         #using polyfit results in numerical errors. good old algebra.
 435         a=diffy/diffx
 436         b=first_point[1]-(a*first_point[0])
 437         baseline=scipy.polyval((a,b), xext)
 438
 439         ysub=[item-basitem for item,basitem in zip(yext,baseline)]
 440
 441         contact=ysub.index(min(ysub))
 442
 443         return xext,ysub,contact
 444
 445         #now, exploit a ClickedPoint instance to calculate index...
 446         dummy=ClickedPoint()
 447         dummy.absolute_coords=(x_intercept,y_intercept)
 448         dummy.find_graph_coords(xret2,yret)
 449
 450         if debug:
 451             return dummy.index, regr, regr_contact
 452         else:
 453             return dummy.index
 454
 455     def find_contact_point_wtk(self, params, z_data, d_data, outqueue=None):
 456         """Algorithm by W. Trevor King.
 457
 458         Notes
 459         -----
 460         Uses :class:`SurfacePositionModel` internally.
 461         """
 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},
 467                                  rescale=True)
 468         for argument in self._wtk_fit_check_arguments:
 469             s.info[argument.name] = params[argument.name]
 470         ignore_index = None
 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'],
 478                                  peaks))
 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(
 483             outqueue=outqueue)
 484         info = {
 485             'offset': offset,
 486             'contact slope': contact_slope,
 487             'surface index': surface_index,
 488             'non-contact slope': non_contact_slope,
 489             'reversed': reverse,
 490             }
 491         deflection_offset = offset + contact_slope*surface_index,
 492         if reverse == True:
 493             surface_index = len(d_data)-1-surface_index
 494         return (numpy.round(surface_index), deflection_offset, info)
 495
 496
 497 class ForceCommand (ColumnAddingCommand):
 498     """Convert a deflection column from meters to newtons.
 499     """
 500     def __init__(self, plugin):
 501         super(ForceCommand, self).__init__(
 502             name='convert distance to force',
 503             columns=[
 504                 ('deflection column', 'surface deflection (m)', """
 505 Name of the column to use as the deflection input.
 506 """.strip()),
 507                 ],
 508             new_columns=[
 509                 ('output deflection column', 'deflection', """
 510 Name of the column (without units) to use as the deflection output.
 511 """.strip()),
 512                 ],
 513             arguments=[
 514                 Argument(name='spring constant info name', type='string',
 515                          default='spring constant (N/m)',
 516                          help="""
 517 Name of the spring constant in the `.info` dictionary.
 518 """.strip()),
 519                 ],
 520             help=self.__doc__, plugin=plugin)
 521
 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',
 530                          values=out)
 531
 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']))
 544         return params
 545
 546
 547 class CantileverAdjustedExtensionCommand (ColumnAddingCommand):
 548     """Remove cantilever extension from a total extension column.
 549
 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
 553     one is used.
 554     """
 555     def __init__(self, plugin):
 556         super(CantileverAdjustedExtensionCommand, self).__init__(
 557             name='remove cantilever from extension',
 558             columns=[
 559                 ('distance column', 'surface distance (m)', """
 560 Name of the column to use as the surface position input.
 561 """.strip()),
 562                 ('deflection column', 'deflection (N)', """
 563 Name of the column to use as the deflection input.
 564 """.strip()),
 565                 ],
 566             new_columns=[
 567                 ('output distance column', 'cantilever adjusted extension', """
 568 Name of the column (without units) to use as the surface position output.
 569 """.strip()),
 570                 ],
 571             arguments=[
 572                 Argument(name='spring constant info name', type='string',
 573                          default='spring constant (N/m)',
 574                          help="""
 575 Name of the spring constant in the `.info` dictionary.
 576 """.strip()),
 577                 ],
 578             help=self.__doc__, plugin=plugin)
 579
 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
 589         else:
 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)
 594
 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
 602         else:
 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']))
 610         return params
 611
 612
 613 class FlattenCommand (ColumnAddingCommand):
 614     """Flatten a deflection column.
 615
 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`.
 619
 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 )
 624     """
 625     def __init__(self, plugin):
 626         super(FlattenCommand, self).__init__(
 627             name='polynomial flatten',
 628             columns=[
 629                 ('distance column', 'surface distance (m)', """
 630 Name of the column to use as the surface position input.
 631 """.strip()),
 632                 ('deflection column', 'deflection (N)', """
 633 Name of the column to use as the deflection input.
 634 """.strip()),
 635                 ],
 636             new_columns=[
 637                 ('output deflection column', 'flattened deflection', """
 638 Name of the column (without units) to use as the deflection output.
 639 """.strip()),
 640                 ],
 641             arguments=[
 642                 Argument(name='degree', type='int',
 643                          default=1,
 644                          help="""
 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...)
 648 """.strip()),
 649                 Argument(name='fit info name', type='string',
 650                          default='flatten fit',
 651                          help="""
 652 Name of the flattening information in the `.info` dictionary.
 653 """.strip()),
 654                 ],
 655             help=self.__doc__, plugin=plugin)
 656
 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']
 666         mask = dist_data > 0
 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()
 673
 674         try:
 675             poly_values = scipy.polyfit(dist_nc, def_nc, degree)
 676             def_nc_fit = scipy.polyval(poly_values, dist_nc)
 677         except Exception, e:
 678             raise Failure('failed to flatten with a degree %d polynomial: %s'
 679                           % (degree, e))
 680         error = numpy.sqrt((def_nc_fit-def_nc)**2).sum() / len(def_nc)
 681         block.info[params['fit info name']] = {
 682             'error':error,
 683             'degree':degree,
 684             'polynomial values':poly_values,
 685             }
 686         out = (def_data
 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',
 691                          values=out)
 692
 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)
 697         return params