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