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