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 ..plugin import Plugin
  38 from ..util.fit import PoorFit, ModelFitter
  39 from ..util.si import join_data_label, split_data_label
  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         params = self.__setup_params(hooke=hooke, params=params)
 264         block = self._block(hooke=hooke, params=params)
 265         dist_data = self._get_column(hooke=hooke, params=params,
 266                                      column_name='distance column')
 267         def_data = self._get_column(hooke=hooke, params=params,
 268                                     column_name='deflection column')
 269         i,def_offset,ps = self.find_contact_point(
 270             params, dist_data, def_data, outqueue)
 271         dist_offset = dist_data[i]
 272         block.info[params['distance info name']] = dist_offset
 273         block.info[params['deflection info name']] = def_offset
 274         block.info[params['fit parameters info name']] = ps
 275         self._set_column(hooke=hooke, params=params,
 276                          column_name='output distance column',
 277                          values=dist_data - dist_offset)
 278         self._set_column(hooke=hooke, params=params,
 279                          column_name='output deflection column',
 280                          values=def_data - def_offset)
 281
 282     def __setup_params(self, hooke, params):
 283         name,dist_unit = split_data_label(params['distance column'])
 284         name,def_unit = split_data_label(params['deflection column'])
 285         params['output distance column'] = join_data_label(
 286             params['output distance column'], dist_unit)
 287         params['output deflection column'] = join_data_label(
 288             params['output deflection column'], def_unit)
 289         params['distance info name'] = join_data_label(
 290             params['distance info name'], dist_unit)
 291         params['deflection info name'] = join_data_label(
 292             params['deflection info name'], def_unit)
 293         return params
 294
 295     def find_contact_point(self, params, z_data, d_data, outqueue=None):
 296         """Railyard for the `find_contact_point_*` family.
 297
 298         Uses the `surface contact point algorithm` configuration
 299         setting to call the appropriate backend algorithm.
 300         """
 301         fn = getattr(self, 'find_contact_point_%s'
 302                      % self.plugin.config['surface contact point algorithm'])
 303         return fn(params, z_data, d_data, outqueue)
 304
 305     def find_contact_point_fmms(self, params, z_data, d_data, outqueue=None):
 306         """Algorithm by Francesco Musiani and Massimo Sandal.
 307
 308         Notes
 309         -----
 310         Algorithm:
 311
 312         0) Driver-specific workarounds, e.g. deal with the PicoForce
 313           trigger bug by excluding retraction portions with excessive
 314           deviation.
 315         1) Select the second half (non-contact side) of the retraction
 316           curve.
 317         2) Fit the selection to a line.
 318         3) If the fit is not almost horizontal, halve the selection
 319           and retrun to (2).
 320         4) Average the selection and use it as a baseline.
 321         5) Slide in from the start (contact side) of the retraction
 322         curve, until you find a point with greater than baseline
 323         deflection.  That point is the contact point.
 324         """
 325         if params['curve'].info['filetype'] == 'picoforce':
 326             # Take care of the picoforce trigger bug (TODO: example
 327             # data file demonstrating the bug).  We exclude portions
 328             # of the curve that have too much standard deviation.
 329             # Yes, a lot of magic is here.
 330             check_start = len(d_data)-len(d_data)/20
 331             monster_start = len(d_data)
 332             while True:
 333                 # look at the non-contact tail
 334                 non_monster = d_data[check_start:monster_start]
 335                 if non_monster.std() < 2e-10: # HACK: hardcoded cutoff
 336                     break
 337                 else: # move further away from the monster
 338                     check_start -= len(d_data)/50
 339                     monster_start -= len(d_data)/50
 340             z_data = z_data[:monster_start]
 341             d_data = d_data[:monster_start]
 342
 343         # take half of the thing to start
 344         selection_start = len(d_data)/2
 345         while True:
 346             z_chunk = z_data[selection_start:]
 347             d_chunk = d_data[selection_start:]
 348             slope,intercept,r,two_tailed_prob,stderr_of_the_estimate = \
 349                 scipy.stats.linregress(z_chunk, d_chunk)
 350             # We stop if we found an almost-horizontal fit or if we're
 351             # getting to small a selection.  FIXME: 0.1 and 5./6 here
 352             # are "magic numbers" (although reasonable)
 353             if (abs(slope) < 0.1  # deflection (m) / surface (m)
 354                 or selection_start > 5./6*len(d_data)):
 355                 break
 356             selection_start += 10
 357
 358         d_baseline = d_chunk.mean()
 359
 360         # find the first point above the calculated baseline
 361         i = 0
 362         while i < len(d_data) and d_data[i] < ymean:
 363             i += 1
 364         return (i, d_baseline, {})
 365
 366     def find_contact_point_ms(self, params, z_data, d_data, outqueue=None):
 367         """Algorithm by Massimo Sandal.
 368
 369         Notes
 370         -----
 371         WTK: At least the commits are by Massimo, and I see no notes
 372         attributing the algorithm to anyone else.
 373
 374         Algorithm:
 375
 376         * ?
 377         """
 378         xext=raw_plot.vectors[0][0]
 379         yext=raw_plot.vectors[0][1]
 380         xret2=raw_plot.vectors[1][0]
 381         yret=raw_plot.vectors[1][1]
 382
 383         first_point=[xext[0], yext[0]]
 384         last_point=[xext[-1], yext[-1]]
 385
 386         #regr=scipy.polyfit(first_point, last_point,1)[0:2]
 387         diffx=abs(first_point[0]-last_point[0])
 388         diffy=abs(first_point[1]-last_point[1])
 389
 390         #using polyfit results in numerical errors. good old algebra.
 391         a=diffy/diffx
 392         b=first_point[1]-(a*first_point[0])
 393         baseline=scipy.polyval((a,b), xext)
 394
 395         ysub=[item-basitem for item,basitem in zip(yext,baseline)]
 396
 397         contact=ysub.index(min(ysub))
 398
 399         return xext,ysub,contact
 400
 401         #now, exploit a ClickedPoint instance to calculate index...
 402         dummy=ClickedPoint()
 403         dummy.absolute_coords=(x_intercept,y_intercept)
 404         dummy.find_graph_coords(xret2,yret)
 405
 406         if debug:
 407             return dummy.index, regr, regr_contact
 408         else:
 409             return dummy.index
 410
 411     def find_contact_point_wtk(self, params, z_data, d_data, outqueue=None):
 412         """Algorithm by W. Trevor King.
 413
 414         Notes
 415         -----
 416         Uses :class:`SurfacePositionModel` internally.
 417         """
 418         reverse = z_data[0] > z_data[-1]
 419         if reverse == True:    # approaching, contact region on the right
 420             d_data = d_data[::-1]
 421         s = SurfacePositionModel(d_data, info={
 422                 'force zero non-contact slope':True},
 423                                  rescale=True)
 424         ignore_index = None
 425         if params['ignore index'] != None:
 426             ignore_index = params['ignore index']
 427         elif params['ignore after last peak info name'] != None:
 428             peaks = z_data.info[params['ignore after last peak info name']]
 429             if not len(peaks) > 0:
 430                 raise Failure('Need at least one peak in %s, not %s'
 431                               % (params['ignore after last peak info name'],
 432                                  peaks))
 433             ignore_index = peaks[-1].post_index()
 434         if ignore_index != None:
 435             s.info['ignore non-contact before index'] = ignore_index
 436         offset,contact_slope,surface_index,non_contact_slope = s.fit(
 437             outqueue=outqueue)
 438         info = {
 439             'offset': offset,
 440             'contact slope': contact_slope,
 441             'surface index': surface_index,
 442             'non-contact slope': non_contact_slope,
 443             'reversed': reverse,
 444             }
 445         deflection_offset = offset + contact_slope*surface_index,
 446         if reverse == True:
 447             surface_index = len(d_data)-1-surface_index
 448         return (numpy.round(surface_index), deflection_offset, info)
 449
 450
 451 class ForceCommand (ColumnAddingCommand):
 452     """Convert a deflection column from meters to newtons.
 453     """
 454     def __init__(self, plugin):
 455         super(ForceCommand, self).__init__(
 456             name='convert distance to force',
 457             columns=[
 458                 ('deflection column', 'surface deflection (m)', """
 459 Name of the column to use as the deflection input.
 460 """.strip()),
 461                 ],
 462             new_columns=[
 463                 ('output deflection column', 'deflection', """
 464 Name of the column (without units) to use as the deflection output.
 465 """.strip()),
 466                 ],
 467             arguments=[
 468                 Argument(name='spring constant info name', type='string',
 469                          default='spring constant (N/m)',
 470                          help="""
 471 Name of the spring constant in the `.info` dictionary.
 472 """.strip()),
 473                 ],
 474             help=self.__doc__, plugin=plugin)
 475
 476     def _run(self, hooke, inqueue, outqueue, params):
 477         params = self.__setup_params(hooke=hooke, params=params)
 478         def_data = self._get_column(hooke=hooke, params=params,
 479                                     column_name='deflection column')
 480         out = def_data * def_data.info[params['spring constant info name']]
 481         self._set_column(hooke=hooke, params=params,
 482                          column_name='output deflection column',
 483                          values=out)
 484
 485     def __setup_params(self, hooke, params):
 486         name,in_unit = split_data_label(params['deflection column'])
 487         out_unit = 'N'  # HACK: extract target units from k_unit.
 488         params['output deflection column'] = join_data_label(
 489             params['output deflection column'], out_unit)
 490         name,k_unit = split_data_label(params['spring constant info name'])
 491         expected_k_unit = '%s/%s' % (out_unit, in_unit)
 492         if k_unit != expected_k_unit:
 493             raise Failure('Cannot convert from %s to %s with %s'
 494                           % (params['deflection column'],
 495                              params['output deflection column'],
 496                              params['spring constant info name']))
 497         return params
 498
 499
 500 class CantileverAdjustedExtensionCommand (ColumnAddingCommand):
 501     """Remove cantilever extension from a total extension column.
 502
 503     If `distance column` and `deflection column` have the same units
 504     (e.g. `z piezo (m)` and `deflection (m)`), `spring constant info
 505     name` is ignored and a deflection/distance conversion factor of
 506     one is used.
 507     """
 508     def __init__(self, plugin):
 509         super(CantileverAdjustedExtensionCommand, self).__init__(
 510             name='remove cantilever from extension',
 511             columns=[
 512                 ('distance column', 'surface distance (m)', """
 513 Name of the column to use as the surface position input.
 514 """.strip()),
 515                 ('deflection column', 'deflection (N)', """
 516 Name of the column to use as the deflection input.
 517 """.strip()),
 518                 ],
 519             new_columns=[
 520                 ('output distance column', 'cantilever adjusted extension', """
 521 Name of the column (without units) to use as the surface position output.
 522 """.strip()),
 523                 ],
 524             arguments=[
 525                 Argument(name='spring constant info name', type='string',
 526                          default='spring constant (N/m)',
 527                          help="""
 528 Name of the spring constant in the `.info` dictionary.
 529 """.strip()),
 530                 ],
 531             help=self.__doc__, plugin=plugin)
 532
 533     def _run(self, hooke, inqueue, outqueue, params):
 534         params = self.__setup_params(hooke=hooke, params=params)
 535         def_data = self._get_column(hooke=hooke, params=params,
 536                                     column_name='deflection column')
 537         dist_data = self._get_column(hooke=hooke, params=params,
 538                                      column_name='distance column')
 539         if params['spring constant info name'] == None:
 540             k = 1.0  # distance and deflection in the same units
 541         else:
 542             k = def_data.info[params['spring constant info name']]
 543         self._set_column(hooke=hooke, params=params,
 544                          column_name='output distance column',
 545                          values=dist_data - def_data / k)
 546
 547     def __setup_params(self, hooke, params):
 548         name,dist_unit = split_data_label(params['distance column'])
 549         name,def_unit = split_data_label(params['deflection column'])
 550         params['output distance column'] = join_data_label(
 551             params['output distance column'], dist_unit)
 552         if dist_unit == def_unit:
 553             params['spring constant info name'] == None
 554         else:
 555             name,k_unit = split_data_label(params['spring constant info name'])
 556             expected_k_unit = '%s/%s' % (def_unit, dist_unit)
 557             if k_unit != expected_k_unit:
 558                 raise Failure('Cannot convert from %s to %s with %s'
 559                               % (params['deflection column'],
 560                                  params['output distance column'],
 561                                  params['spring constant info name']))
 562         return params
 563
 564
 565 class FlattenCommand (ColumnAddingCommand):
 566     """Flatten a deflection column.
 567
 568     Subtracts a polynomial fit from the non-contact part of the curve
 569     to flatten it.  The best polynomial fit is chosen among
 570     polynomials of degree 1 to `max degree`.
 571
 572     .. todo:  Why does flattening use a polynomial fit and not a sinusoid?
 573       Isn't most of the oscillation due to laser interference?
 574       See Jaschke 1995 ( 10.1063/1.1146018 )
 575       and the figure 4 caption of Weisenhorn 1992 ( 10.1103/PhysRevB.45.11226 )
 576     """
 577     def __init__(self, plugin):
 578         super(FlattenCommand, self).__init__(
 579             name='polynomial flatten',
 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 deflection column', 'flattened deflection', """
 590 Name of the column (without units) to use as the deflection output.
 591 """.strip()),
 592                 ],
 593             arguments=[
 594                 Argument(name='degree', type='int',
 595                          default=1,
 596                          help="""
 597 Order of the polynomial used for flattening.  Using values greater
 598 than one usually doesn't help and can give artifacts.  However, it
 599 could be useful too.  (TODO: Back this up with some theory...)
 600 """.strip()),
 601                 Argument(name='fit info name', type='string',
 602                          default='flatten fit',
 603                          help="""
 604 Name of the flattening information in the `.info` dictionary.
 605 """.strip()),
 606                 ],
 607             help=self.__doc__, plugin=plugin)
 608
 609     def _run(self, hooke, inqueue, outqueue, params):
 610         params = self.__setup_params(hooke=hooke, params=params)
 611         block = self._block(hooke=hooke, params=params)
 612         dist_data = self._get_column(hooke=hooke, params=params,
 613                                      column_name='distance column')
 614         def_data = self._get_column(hooke=hooke, params=params,
 615                                     column_name='deflection column')
 616         degree = params['degree']
 617         mask = dist_data > 0
 618         indices = numpy.argwhere(mask)
 619         if len(indices) == 0:
 620             raise Failure('no positive distance values in %s'
 621                           % params['distance column'])
 622         dist_nc = dist_data[indices].flatten()
 623         def_nc = def_data[indices].flatten()
 624
 625         try:
 626             poly_values = scipy.polyfit(dist_nc, def_nc, degree)
 627             def_nc_fit = scipy.polyval(poly_values, dist_nc)
 628         except Exception, e:
 629             raise Failure('failed to flatten with a degree %d polynomial: %s'
 630                           % (degree, e))
 631         error = numpy.sqrt((def_nc_fit-def_nc)**2).sum() / len(def_nc)
 632         block.info[params['fit info name']] = {
 633             'error':error,
 634             'degree':degree,
 635             'polynomial values':poly_values,
 636             }
 637         out = def_data - mask*scipy.polyval(poly_values, dist_data)
 638         self._set_column(hooke=hooke, params=params,
 639                          column_name='output deflection column',
 640                          values=out)
 641
 642     def __setup_params(self, hooke, params):
 643         d_name,d_unit = split_data_label(params['deflection column'])
 644         params['output deflection column'] = join_data_label(
 645             params['output deflection column'], d_unit)
 646         return params