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