hooke/plugin/vclamp.py

   1 # Copyright (C) 2008-2011 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'] is 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         if kink_position == 0:
 170             # jump vibration at the start of the retraction?
 171             start = int(min(max(20, 0.01 * len(self._data)), 0.5*len(self._data)))
 172             std = self._data[:start].std()
 173             left_offset = self._data[start].mean()
 174             stop = start
 175             while abs(self._data[stop] - left_offset) < 3*std:
 176                 stop += 1
 177             left_slope = (self._data[stop-start:stop].mean()
 178                           - left_offset) / (stop-start)
 179             left_offset -= left_slope * start/2
 180             kink_position = (self._data[-1] - left_offset)/left_slope
 181         else:
 182             left_slope = (self._data[-1] - self.info['min deflection']
 183                           )/kink_position
 184         right_slope = 0
 185         self.info['guessed contact slope'] = left_slope
 186         params = [left_offset, left_slope, kink_position, right_slope]
 187         if self.info['force zero non-contact slope'] == True:
 188             params = params[:-1]
 189         return params
 190
 191     def fit(self, *args, **kwargs):
 192         """Fit the model to the data.
 193
 194         Notes
 195         -----
 196         We change the `epsfcn` default from :func:`scipy.optimize.leastsq`'s
 197         `0` to `1e-3`, so the initial Jacobian estimate takes larger steps,
 198         which helps avoid being trapped in noise-generated local minima.
 199         """
 200         self.info['guessed contact slope'] = None
 201         if 'epsfcn' not in kwargs:
 202             kwargs['epsfcn'] = 1e-3  # take big steps to estimate the Jacobian
 203         params = super(SurfacePositionModel, self).fit(*args, **kwargs)
 204         params[2] = abs(params[2])
 205         if self.info['force zero non-contact slope'] == True:
 206             params = list(params)
 207             params.append(0.)  # restore the non-contact slope parameter
 208
 209         # check that the fit is reasonable, see the :meth:`model` docstring
 210         # for parameter descriptions.
 211         slope_ratio = abs(params[1]/params[3])
 212         if slope_ratio < self.info['min slope ratio']:
 213             raise PoorFit(
 214                'Slope in non-contact region, or no slope in contact (slope ratio %g less than %g)'
 215                % (slope_ratio, self.info['min slope ratio']))
 216         contact_fraction = ((params[2]-self.info['min position'])
 217                             /self.info['position range'])
 218         if contact_fraction < self.info['min contact fraction']:
 219             raise PoorFit(
 220                 'No kink (contact fraction %g less than %g)'
 221                 % (contact_fraction, self.info['min contact fraction']))
 222         if contact_fraction > self.info['max contact fraction']:
 223             raise PoorFit(
 224                 'No kink (contact fraction %g greater than %g)'
 225                 % (contact_fraction, self.info['max contact fraction']))
 226         slope_guess_ratio = abs(params[1]/self.info['guessed contact slope'])
 227         if (self.info['guessed contact slope'] != None
 228             and slope_guess_ratio < self.info['min slope guess ratio']):
 229             raise PoorFit(
 230                 'Too far (contact slope off guess by %g less than %g)'
 231                 % (slope_guess_ratio, self.info['min slope guess ratio']))
 232         return params
 233
 234
 235 class VelocityClampPlugin (Plugin):
 236     def __init__(self):
 237         super(VelocityClampPlugin, self).__init__(name='vclamp')
 238         self._commands = [
 239             SurfaceContactCommand(self), ForceCommand(self),
 240             CantileverAdjustedExtensionCommand(self), FlattenCommand(self),
 241             ]
 242
 243     def default_settings(self):
 244         return [
 245             Setting(section=self.setting_section, help=self.__doc__),
 246             Setting(section=self.setting_section,
 247                     option='surface contact point algorithm',
 248                     value='wtk',
 249                     help='Select the surface contact point algorithm.  See the documentation for descriptions of available algorithms.')
 250             ]
 251
 252
 253 class SurfaceContactCommand (ColumnAddingCommand):
 254     """Automatically determine a block's surface contact point.
 255
 256     You can select the contact point algorithm with the creatively
 257     named `surface contact point algorithm` configuration setting.
 258     Currently available options are:
 259
 260     * fmms (:meth:`find_contact_point_fmms`)
 261     * ms (:meth:`find_contact_point_ms`)
 262     * wtk (:meth:`find_contact_point_wtk`)
 263     """
 264     def __init__(self, plugin):
 265         self._wtk_fit_check_arguments = []
 266         for argument in SurfacePositionModel._fit_check_arguments:
 267             arg = copy.deepcopy(argument)
 268             arg._help += '  (wtk model)'
 269             self._wtk_fit_check_arguments.append(arg)
 270         super(SurfaceContactCommand, self).__init__(
 271             name='zero surface contact point',
 272             columns=[
 273                 ('distance column', 'z piezo (m)', """
 274 Name of the column to use as the surface position input.
 275 """.strip()),
 276                 ('deflection column', 'deflection (m)', """
 277 Name of the column to use as the deflection input.
 278 """.strip()),
 279                 ],
 280             new_columns=[
 281                 ('output distance column', 'surface distance', """
 282 Name of the column (without units) to use as the surface position output.
 283 """.strip()),
 284                 ('output deflection column', 'surface deflection', """
 285 Name of the column (without units) to use as the deflection output.
 286 """.strip()),
 287                 ],
 288             arguments=[
 289                 Argument(name='ignore index', type='int', default=None,
 290                          help="""
 291 Ignore the residual from the non-contact region before the indexed
 292 point (for the `wtk` algorithm).
 293 """.strip()),
 294                 Argument(name='ignore after last peak info name',
 295                          type='string', default=None,
 296                          help="""
 297 As an alternative to 'ignore index', ignore after the last peak in the
 298 peak list stored in the `.info` dictionary.
 299 """.strip()),
 300                 Argument(name='force zero non-contact slope', type='bool',
 301                          default=False, count=1,
 302                          help="""
 303 Fix the fitted non-contact slope at zero.
 304 """.strip()),
 305                 Argument(name='distance info name', type='string',
 306                          default='surface distance offset',
 307                          help="""
 308 Name (without units) for storing the distance offset in the `.info` dictionary.
 309 """.strip()),
 310                 Argument(name='deflection info name', type='string',
 311                          default='surface deflection offset',
 312                          help="""
 313 Name (without units) for storing the deflection offset in the `.info` dictionary.
 314 """.strip()),
 315                 Argument(name='fit parameters info name', type='string',
 316                          default='surface deflection offset',
 317                          help="""
 318 Name (without units) for storing fit parameters in the `.info` dictionary.
 319 """.strip()),
 320                 ] + self._wtk_fit_check_arguments,
 321             help=self.__doc__, plugin=plugin)
 322
 323     def _run(self, hooke, inqueue, outqueue, params):
 324         self._add_to_command_stack(params)
 325         params = self._setup_params(hooke=hooke, params=params)
 326         block = self._block(hooke=hooke, params=params)
 327         dist_data = self._get_column(hooke=hooke, params=params,
 328                                      column_name='distance column')
 329         def_data = self._get_column(hooke=hooke, params=params,
 330                                     column_name='deflection column')
 331         i,def_offset,ps = self.find_contact_point(
 332             params, dist_data, def_data, outqueue)
 333         dist_offset = dist_data[i]
 334         block.info[params['distance info name']] = dist_offset
 335         block.info[params['deflection info name']] = def_offset
 336         block.info[params['fit parameters info name']] = ps
 337         self._set_column(hooke=hooke, params=params,
 338                          column_name='output distance column',
 339                          values=dist_data - dist_offset)
 340         self._set_column(hooke=hooke, params=params,
 341                          column_name='output deflection column',
 342                          values=def_data - def_offset)
 343
 344     def _setup_params(self, hooke, params):
 345         name,dist_unit = split_data_label(params['distance column'])
 346         name,def_unit = split_data_label(params['deflection column'])
 347         params['output distance column'] = join_data_label(
 348             params['output distance column'], dist_unit)
 349         params['output deflection column'] = join_data_label(
 350             params['output deflection column'], def_unit)
 351         params['distance info name'] = join_data_label(
 352             params['distance info name'], dist_unit)
 353         params['deflection info name'] = join_data_label(
 354             params['deflection info name'], def_unit)
 355         return params
 356
 357     def find_contact_point(self, params, z_data, d_data, outqueue=None):
 358         """Railyard for the `find_contact_point_*` family.
 359
 360         Uses the `surface contact point algorithm` configuration
 361         setting to call the appropriate backend algorithm.
 362         """
 363         fn = getattr(self, 'find_contact_point_%s'
 364                      % self.plugin.config['surface contact point algorithm'])
 365         return fn(params, z_data, d_data, outqueue)
 366
 367     def find_contact_point_fmms(self, params, z_data, d_data, outqueue=None):
 368         """Algorithm by Francesco Musiani and Massimo Sandal.
 369
 370         Notes
 371         -----
 372         Algorithm:
 373
 374         0) Driver-specific workarounds, e.g. deal with the PicoForce
 375           trigger bug by excluding retraction portions with excessive
 376           deviation.
 377         1) Select the second half (non-contact side) of the retraction
 378           curve.
 379         2) Fit the selection to a line.
 380         3) If the fit is not almost horizontal, halve the selection
 381           and retrun to (2).
 382         4) Average the selection and use it as a baseline.
 383         5) Slide in from the start (contact side) of the retraction
 384         curve, until you find a point with greater than baseline
 385         deflection.  That point is the contact point.
 386         """
 387         if params['curve'].info['filetype'] == 'picoforce':
 388             # Take care of the picoforce trigger bug (TODO: example
 389             # data file demonstrating the bug).  We exclude portions
 390             # of the curve that have too much standard deviation.
 391             # Yes, a lot of magic is here.
 392             check_start = len(d_data)-len(d_data)/20
 393             monster_start = len(d_data)
 394             while True:
 395                 # look at the non-contact tail
 396                 non_monster = d_data[check_start:monster_start]
 397                 if non_monster.std() < 2e-10: # HACK: hardcoded cutoff
 398                     break
 399                 else: # move further away from the monster
 400                     check_start -= len(d_data)/50
 401                     monster_start -= len(d_data)/50
 402             z_data = z_data[:monster_start]
 403             d_data = d_data[:monster_start]
 404
 405         # take half of the thing to start
 406         selection_start = len(d_data)/2
 407         while True:
 408             z_chunk = z_data[selection_start:]
 409             d_chunk = d_data[selection_start:]
 410             slope,intercept,r,two_tailed_prob,stderr_of_the_estimate = \
 411                 scipy.stats.linregress(z_chunk, d_chunk)
 412             # We stop if we found an almost-horizontal fit or if we're
 413             # getting to small a selection.  FIXME: 0.1 and 5./6 here
 414             # are "magic numbers" (although reasonable)
 415             if (abs(slope) < 0.1  # deflection (m) / surface (m)
 416                 or selection_start > 5./6*len(d_data)):
 417                 break
 418             selection_start += 10
 419
 420         d_baseline = d_chunk.mean()
 421
 422         # find the first point above the calculated baseline
 423         i = 0
 424         while i < len(d_data) and d_data[i] < ymean:
 425             i += 1
 426         return (i, d_baseline, {})
 427
 428     def find_contact_point_ms(self, params, z_data, d_data, outqueue=None):
 429         """Algorithm by Massimo Sandal.
 430
 431         Notes
 432         -----
 433         WTK: At least the commits are by Massimo, and I see no notes
 434         attributing the algorithm to anyone else.
 435
 436         Algorithm:
 437
 438         * ?
 439         """
 440         xext=raw_plot.vectors[0][0]
 441         yext=raw_plot.vectors[0][1]
 442         xret2=raw_plot.vectors[1][0]
 443         yret=raw_plot.vectors[1][1]
 444
 445         first_point=[xext[0], yext[0]]
 446         last_point=[xext[-1], yext[-1]]
 447
 448         #regr=scipy.polyfit(first_point, last_point,1)[0:2]
 449         diffx=abs(first_point[0]-last_point[0])
 450         diffy=abs(first_point[1]-last_point[1])
 451
 452         #using polyfit results in numerical errors. good old algebra.
 453         a=diffy/diffx
 454         b=first_point[1]-(a*first_point[0])
 455         baseline=scipy.polyval((a,b), xext)
 456
 457         ysub=[item-basitem for item,basitem in zip(yext,baseline)]
 458
 459         contact=ysub.index(min(ysub))
 460
 461         return xext,ysub,contact
 462
 463         #now, exploit a ClickedPoint instance to calculate index...
 464         dummy=ClickedPoint()
 465         dummy.absolute_coords=(x_intercept,y_intercept)
 466         dummy.find_graph_coords(xret2,yret)
 467
 468         if debug:
 469             return dummy.index, regr, regr_contact
 470         else:
 471             return dummy.index
 472
 473     def find_contact_point_wtk(self, params, z_data, d_data, outqueue=None):
 474         """Algorithm by W. Trevor King.
 475
 476         Notes
 477         -----
 478         Uses :class:`SurfacePositionModel` internally.
 479         """
 480         reverse = z_data[0] > z_data[-1]
 481         if reverse == True:    # approaching, contact region on the right
 482             d_data = d_data[::-1]
 483         s = SurfacePositionModel(d_data, info={
 484                 'force zero non-contact slope':
 485                     params['force zero non-contact slope']},
 486                                  rescale=True)
 487         for argument in self._wtk_fit_check_arguments:
 488             s.info[argument.name] = params[argument.name]
 489         ignore_index = None
 490         if params['ignore index'] != None:
 491             ignore_index = params['ignore index']
 492         elif params['ignore after last peak info name'] != None:
 493             peaks = z_data.info[params['ignore after last peak info name']]
 494             if not len(peaks) > 0:
 495                 raise Failure('Need at least one peak in %s, not %s'
 496                               % (params['ignore after last peak info name'],
 497                                  peaks))
 498             ignore_index = peaks[-1].post_index()
 499         if ignore_index != None:
 500             s.info['ignore non-contact before index'] = ignore_index
 501         offset,contact_slope,surface_index,non_contact_slope = s.fit(
 502             outqueue=outqueue)
 503         deflection_offset = offset + contact_slope*surface_index
 504         delta_pos_per_point = z_data[1] - z_data[0]
 505         contact_slope /= delta_pos_per_point  # ddef/point -> ddev/dpos
 506         non_contact_slope /= delta_pos_per_point
 507         info = {
 508             'offset': offset,
 509             'contact slope': contact_slope,
 510             'surface index': surface_index,
 511             'non-contact slope': non_contact_slope,
 512             'reversed': reverse,
 513             }
 514         if reverse == True:
 515             surface_index = len(d_data)-1-surface_index
 516         return (numpy.round(surface_index), deflection_offset, info)
 517
 518
 519 class ForceCommand (ColumnAddingCommand):
 520     """Convert a deflection column from meters to newtons.
 521     """
 522     def __init__(self, plugin):
 523         super(ForceCommand, self).__init__(
 524             name='convert distance to force',
 525             columns=[
 526                 ('deflection column', 'surface deflection (m)', """
 527 Name of the column to use as the deflection input.
 528 """.strip()),
 529                 ],
 530             new_columns=[
 531                 ('output deflection column', 'deflection', """
 532 Name of the column (without units) to use as the deflection output.
 533 """.strip()),
 534                 ],
 535             arguments=[
 536                 Argument(name='spring constant info name', type='string',
 537                          default='spring constant (N/m)',
 538                          help="""
 539 Name of the spring constant in the `.info` dictionary.
 540 """.strip()),
 541                 ],
 542             help=self.__doc__, plugin=plugin)
 543
 544     def _run(self, hooke, inqueue, outqueue, params):
 545         self._add_to_command_stack(params)
 546         params = self._setup_params(hooke=hooke, params=params)
 547         # TODO: call .curve.ScaledColumnAdditionCommand
 548         def_data = self._get_column(hooke=hooke, params=params,
 549                                     column_name='deflection column')
 550         out = def_data * def_data.info[params['spring constant info name']]
 551         self._set_column(hooke=hooke, params=params,
 552                          column_name='output deflection column',
 553                          values=out)
 554
 555     def _setup_params(self, hooke, params):
 556         name,in_unit = split_data_label(params['deflection column'])
 557         out_unit = 'N'  # HACK: extract target units from k_unit.
 558         params['output deflection column'] = join_data_label(
 559             params['output deflection column'], out_unit)
 560         name,k_unit = split_data_label(params['spring constant info name'])
 561         expected_k_unit = '%s/%s' % (out_unit, in_unit)
 562         if k_unit != expected_k_unit:
 563             raise Failure('Cannot convert from %s to %s with %s'
 564                           % (params['deflection column'],
 565                              params['output deflection column'],
 566                              params['spring constant info name']))
 567         return params
 568
 569
 570 class CantileverAdjustedExtensionCommand (ColumnAddingCommand):
 571     """Remove cantilever extension from a total extension column.
 572
 573     If `distance column` and `deflection column` have the same units
 574     (e.g. `z piezo (m)` and `deflection (m)`), `spring constant info
 575     name` is ignored and a deflection/distance conversion factor of
 576     one is used.
 577     """
 578     def __init__(self, plugin):
 579         super(CantileverAdjustedExtensionCommand, self).__init__(
 580             name='remove cantilever from extension',
 581             columns=[
 582                 ('distance column', 'surface distance (m)', """
 583 Name of the column to use as the surface position input.
 584 """.strip()),
 585                 ('deflection column', 'deflection (N)', """
 586 Name of the column to use as the deflection input.
 587 """.strip()),
 588                 ],
 589             new_columns=[
 590                 ('output distance column', 'cantilever adjusted extension', """
 591 Name of the column (without units) to use as the surface position output.
 592 """.strip()),
 593                 ],
 594             arguments=[
 595                 Argument(name='spring constant info name', type='string',
 596                          default='spring constant (N/m)',
 597                          help="""
 598 Name of the spring constant in the `.info` dictionary.
 599 """.strip()),
 600                 ],
 601             help=self.__doc__, plugin=plugin)
 602
 603     def _run(self, hooke, inqueue, outqueue, params):
 604         self._add_to_command_stack(params)
 605         params = self._setup_params(hooke=hooke, params=params)
 606         def_data = self._get_column(hooke=hooke, params=params,
 607                                     column_name='deflection column')
 608         dist_data = self._get_column(hooke=hooke, params=params,
 609                                      column_name='distance column')
 610         if params['spring constant info name'] == None:
 611             k = 1.0  # distance and deflection in the same units
 612         else:
 613             k = def_data.info[params['spring constant info name']]
 614         self._set_column(hooke=hooke, params=params,
 615                          column_name='output distance column',
 616                          values=dist_data - def_data / k)
 617
 618     def _setup_params(self, hooke, params):
 619         name,dist_unit = split_data_label(params['distance column'])
 620         name,def_unit = split_data_label(params['deflection column'])
 621         params['output distance column'] = join_data_label(
 622             params['output distance column'], dist_unit)
 623         if dist_unit == def_unit:
 624             params['spring constant info name'] == None
 625         else:
 626             name,k_unit = split_data_label(params['spring constant info name'])
 627             expected_k_unit = '%s/%s' % (def_unit, dist_unit)
 628             if k_unit != expected_k_unit:
 629                 raise Failure('Cannot convert from %s to %s with %s'
 630                               % (params['deflection column'],
 631                                  params['output distance column'],
 632                                  params['spring constant info name']))
 633         return params
 634
 635
 636 class FlattenCommand (ColumnAddingCommand):
 637     """Flatten a deflection column.
 638
 639     Subtracts a polynomial fit from the non-contact part of the curve
 640     to flatten it.  The best polynomial fit is chosen among
 641     polynomials of degree 1 to `max degree`.
 642
 643     .. todo:  Why does flattening use a polynomial fit and not a sinusoid?
 644       Isn't most of the oscillation due to laser interference?
 645       See Jaschke 1995 ( 10.1063/1.1146018 )
 646       and the figure 4 caption of Weisenhorn 1992 ( 10.1103/PhysRevB.45.11226 )
 647     """
 648     def __init__(self, plugin):
 649         super(FlattenCommand, self).__init__(
 650             name='polynomial flatten',
 651             columns=[
 652                 ('distance column', 'surface distance (m)', """
 653 Name of the column to use as the surface position input.
 654 """.strip()),
 655                 ('deflection column', 'deflection (N)', """
 656 Name of the column to use as the deflection input.
 657 """.strip()),
 658                 ],
 659             new_columns=[
 660                 ('output deflection column', 'flattened deflection', """
 661 Name of the column (without units) to use as the deflection output.
 662 """.strip()),
 663                 ],
 664             arguments=[
 665                 Argument(name='degree', type='int',
 666                          default=1,
 667                          help="""
 668 Order of the polynomial used for flattening.  Using values greater
 669 than one usually doesn't help and can give artifacts.  However, it
 670 could be useful too.  (TODO: Back this up with some theory...)
 671 """.strip()),
 672                 Argument(name='fit info name', type='string',
 673                          default='flatten fit',
 674                          help="""
 675 Name of the flattening information in the `.info` dictionary.
 676 """.strip()),
 677                 ],
 678             help=self.__doc__, plugin=plugin)
 679
 680     def _run(self, hooke, inqueue, outqueue, params):
 681         self._add_to_command_stack(params)
 682         params = self._setup_params(hooke=hooke, params=params)
 683         block = self._block(hooke=hooke, params=params)
 684         dist_data = self._get_column(hooke=hooke, params=params,
 685                                      column_name='distance column')
 686         def_data = self._get_column(hooke=hooke, params=params,
 687                                     column_name='deflection column')
 688         degree = params['degree']
 689         mask = dist_data > 0
 690         indices = numpy.argwhere(mask)
 691         if len(indices) == 0:
 692             raise Failure('no positive distance values in %s'
 693                           % params['distance column'])
 694         dist_nc = dist_data[indices].flatten()
 695         def_nc = def_data[indices].flatten()
 696
 697         try:
 698             poly_values = scipy.polyfit(dist_nc, def_nc, degree)
 699             def_nc_fit = scipy.polyval(poly_values, dist_nc)
 700         except Exception, e:
 701             raise Failure('failed to flatten with a degree %d polynomial: %s'
 702                           % (degree, e))
 703         error = numpy.sqrt((def_nc_fit-def_nc)**2).sum() / len(def_nc)
 704         block.info[params['fit info name']] = {
 705             'error':error,
 706             'degree':degree,
 707             'polynomial values':poly_values,
 708             }
 709         out = (def_data
 710                - numpy.invert(mask)*scipy.polyval(poly_values[-1:], dist_data)
 711                - mask*scipy.polyval(poly_values, dist_data))
 712         self._set_column(hooke=hooke, params=params,
 713                          column_name='output deflection column',
 714                          values=out)
 715
 716     def _setup_params(self, hooke, params):
 717         d_name,d_unit = split_data_label(params['deflection column'])
 718         params['output deflection column'] = join_data_label(
 719             params['output deflection column'], d_unit)
 720         return params