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
  30 import numpy
  31 import scipy
  32
  33 from ..command import Command, Argument, Failure, NullQueue
  34 from ..config import Setting
  35 from ..curve import Data
  36 from ..plugin import Plugin
  37 from ..util.fit import PoorFit, ModelFitter
  38 from ..util.si import join_data_label, split_data_label
  39 from .curve import CurveArgument
  40
  41
  42 def scale(hooke, curve, block=None):
  43     """Run 'add block force array' on `block`.
  44
  45     Runs 'zero block surface contact point' first, if necessary.  Does
  46     not run either command if the columns they add to the block are
  47     already present.
  48
  49     If `block` is `None`, scale all blocks in `curve`.
  50     """
  51     commands = hooke.commands
  52     contact = [c for c in hooke.commands
  53                if c.name == 'zero block surface contact point'][0]
  54     force = [c for c in hooke.commands if c.name == 'add block force array'][0]
  55     cant_adjust = [c for c in hooke.commands
  56                if c.name == 'add block cantilever adjusted extension array'][0]
  57     inqueue = None
  58     outqueue = NullQueue()
  59     if block == None:
  60         for i in range(len(curve.data)):
  61             scale(hooke, curve, block=i)
  62     else:
  63         params = {'curve':curve, 'block':block}
  64         b = curve.data[block]
  65         if ('surface distance (m)' not in b.info['columns']
  66             or 'surface deflection (m)' not in b.info['columns']):
  67             try:
  68                 contact.run(hooke, inqueue, outqueue, **params)
  69             except PoorFit, e:
  70                 raise PoorFit('Could not fit %s %s: %s'
  71                               % (curve.path, block, str(e)))
  72         if ('deflection (N)' not in b.info['columns']):
  73             force.run(hooke, inqueue, outqueue, **params)
  74         if ('cantilever adjusted extension (m)' not in b.info['columns']):
  75             cant_adjust.run(hooke, inqueue, outqueue, **params)
  76     return curve
  77
  78 class SurfacePositionModel (ModelFitter):
  79     """Bilinear surface position model.
  80
  81     The bilinear model is symmetric, but the parameter guessing and
  82     sanity checks assume the contact region occurs for lower indicies
  83     ("left of") the non-contact region.  We also assume that
  84     tip-surface attractions produce positive deflections.
  85
  86     Notes
  87     -----
  88     Algorithm borrowed from WTK's `piezo package`_, specifically
  89     from :func:`piezo.z_piezo_utils.analyzeSurfPosData`.
  90
  91     .. _piezo package:
  92       http://www.physics.drexel.edu/~wking/code/git/git.php?p=piezo.git
  93
  94     Fits the data to the bilinear :method:`model`.
  95
  96     In order for this model to produce a satisfactory fit, there
  97     should be enough data in the off-surface region that interactions
  98     due to proteins, etc. will not seriously skew the fit in the
  99     off-surface region.
 100     """
 101     def model(self, params):
 102         """A continuous, bilinear model.
 103
 104         Notes
 105         -----
 106         .. math::
 107
 108           y = \begin{cases}
 109             p_0 + p_1 x                 & \text{if $x <= p_2$}, \\
 110             p_0 + p_1 p_2 + p_3 (x-p_2) & \text{if $x >= p_2$}.
 111               \end{cases}
 112
 113         Where :math:`p_0` is a vertical offset, :math:`p_1` is the slope
 114         of the first region, :math:`p_2` is the transition location, and
 115         :math:`p_3` is the slope of the second region.
 116         """
 117         p = params  # convenient alias
 118         if self.info.get('force zero non-contact slope', None) == True:
 119             p = list(p)
 120             p.append(0.)  # restore the non-contact slope parameter
 121         r2 = numpy.round(abs(p[2]))
 122         if r2 >= 1:
 123             self._model_data[:r2] = p[0] + p[1] * numpy.arange(r2)
 124         if r2 < len(self._data)-1:
 125             self._model_data[r2:] = \
 126                 p[0] + p[1]*p[2] + p[3] * numpy.arange(len(self._data)-r2)
 127         return self._model_data
 128
 129     def set_data(self, data, info=None):
 130         super(SurfacePositionModel, self).set_data(data, info)
 131         if info == None:
 132             info = {}
 133         self.info = info
 134         self.info['min position'] = 0
 135         self.info['max position'] = len(data)
 136         self.info['max deflection'] = data.max()
 137         self.info['min deflection'] = data.min()
 138         self.info['position range'] = self.info['max position'] - self.info['min position']
 139         self.info['deflection range'] = self.info['max deflection'] - self.info['min deflection']
 140
 141     def guess_initial_params(self, outqueue=None):
 142         """Guess the initial parameters.
 143
 144         Notes
 145         -----
 146         We guess initial parameters such that the offset (:math:`p_1`)
 147         matches the minimum deflection, the kink (:math:`p_2`) occurs in
 148         the middle of the data, the initial (contact) slope (:math:`p_0`)
 149         produces the maximum deflection at the left-most point, and the
 150         final (non-contact) slope (:math:`p_3`) is zero.
 151         """
 152         left_offset = self.info['min deflection']
 153         left_slope = 2*(self.info['deflection range']
 154                         /self.info['position range'])
 155         kink_position = (self.info['max position']
 156                          +self.info['min position'])/2.0
 157         right_slope = 0
 158         self.info['guessed contact slope'] = left_slope
 159         params = [left_offset, left_slope, kink_position, right_slope]
 160         if self.info.get('force zero non-contact slope', None) == True:
 161             params = params[:-1]
 162         return params
 163
 164     def guess_scale(self, params, outqueue=None):
 165         """Guess the parameter scales.
 166
 167         Notes
 168         -----
 169         We guess offset scale (:math:`p_0`) as one tenth of the total
 170         deflection range, the kink scale (:math:`p_2`) as one tenth of
 171         the total index range, the initial (contact) slope scale
 172         (:math:`p_1`) as one tenth of the contact slope estimation,
 173         and the final (non-contact) slope scale (:math:`p_3`) is as
 174         one tenth of the initial slope scale.
 175         """
 176         offset_scale = self.info['deflection range']/10.
 177         left_slope_scale = abs(params[1])/10.
 178         kink_scale = self.info['position range']/10.
 179         right_slope_scale = left_slope_scale/10.
 180         scale = [offset_scale, left_slope_scale, kink_scale, right_slope_scale]
 181         if self.info.get('force zero non-contact slope', None) == True:
 182             scale = scale[:-1]
 183         return scale
 184
 185     def fit(self, *args, **kwargs):
 186         self.info['guessed contact slope'] = None
 187         params = super(SurfacePositionModel, self).fit(*args, **kwargs)
 188         params[2] = abs(params[2])
 189         if self.info.get('force zero non-contact slope', None) == True:
 190             params = list(params)
 191             params.append(0.)  # restore the non-contact slope parameter
 192
 193         # check that the fit is reasonable, see the :meth:`model` docstring
 194         # for parameter descriptions.  HACK: hardcoded cutoffs.
 195         if abs(params[3]*10) > abs(params[1]) :
 196             raise PoorFit('Slope in non-contact region, or no slope in contact')
 197         if params[2] < self.info['min position']+0.02*self.info['position range']:
 198             raise PoorFit(
 199                 'No kink (kink %g less than %g, need more space to left)'
 200                 % (params[2],
 201                    self.info['min position']+0.02*self.info['position range']))
 202         if params[2] > self.info['max position']-0.02*self.info['position range']:
 203             raise poorFit(
 204                 'No kink (kink %g more than %g, need more space to right)'
 205                 % (params[2],
 206                    self.info['max position']-0.02*self.info['position range']))
 207         if (self.info['guessed contact slope'] != None
 208             and abs(params[1]) < 0.5 * abs(self.info['guessed contact slope'])):
 209             raise PoorFit('Too far (contact slope %g, but expected ~%g'
 210                           % (params[3], self.info['guessed contact slope']))
 211         return params
 212
 213 class VelocityClampPlugin (Plugin):
 214     def __init__(self):
 215         super(VelocityClampPlugin, self).__init__(name='vclamp')
 216         self._commands = [
 217             SurfaceContactCommand(self), ForceCommand(self),
 218             CantileverAdjustedExtensionCommand(self),
 219             ]
 220
 221     def default_settings(self):
 222         return [
 223             Setting(section=self.setting_section, help=self.__doc__),
 224             Setting(section=self.setting_section,
 225                     option='surface contact point algorithm',
 226                     value='wtk',
 227                     help='Select the surface contact point algorithm.  See the documentation for descriptions of available algorithms.')
 228             ]
 229
 230
 231 class SurfaceContactCommand (Command):
 232     """Automatically determine a block's surface contact point.
 233
 234     You can select the contact point algorithm with the creatively
 235     named `surface contact point algorithm` configuration setting.
 236     Currently available options are:
 237
 238     * fmms (:meth:`find_contact_point_fmms`)
 239     * ms (:meth:`find_contact_point_ms`)
 240     * wtk (:meth:`find_contact_point_wtk`)
 241     """
 242     def __init__(self, plugin):
 243         super(SurfaceContactCommand, self).__init__(
 244             name='zero block surface contact point',
 245             arguments=[
 246                 CurveArgument,
 247                 Argument(name='block', aliases=['set'], type='int', default=0,
 248                          help="""
 249 Data block for which the force should be calculated.  For an
 250 approach/retract force curve, `0` selects the approaching curve and `1`
 251 selects the retracting curve.
 252 """.strip()),
 253                 Argument(name='input distance column', type='string',
 254                          default='z piezo (m)',
 255                          help="""
 256 Name of the column to use as the surface positioning input.
 257 """.strip()),
 258                 Argument(name='input deflection column', type='string',
 259                          default='deflection (m)',
 260                          help="""
 261 Name of the column to use as the deflection input.
 262 """.strip()),
 263                 Argument(name='output distance column', type='string',
 264                          default='surface distance',
 265                          help="""
 266 Name of the column (without units) to use as the surface positioning output.
 267 """.strip()),
 268                 Argument(name='output deflection column', type='string',
 269                          default='surface deflection',
 270                          help="""
 271 Name of the column (without units) to use as the deflection output.
 272 """.strip()),
 273                 Argument(name='distance info name', type='string',
 274                          default='surface distance offset',
 275                          help="""
 276 Name (without units) for storing the distance offset in the `.info` dictionary.
 277 """.strip()),
 278                 Argument(name='deflection info name', type='string',
 279                          default='surface deflection offset',
 280                          help="""
 281 Name (without units) for storing the deflection offset in the `.info` dictionary.
 282 """.strip()),
 283                 Argument(name='fit parameters info name', type='string',
 284                          default='surface deflection offset',
 285                          help="""
 286 Name (without units) for storing the deflection offset in the `.info` dictionary.
 287 """.strip()),
 288                 ],
 289             help=self.__doc__, plugin=plugin)
 290
 291     def _run(self, hooke, inqueue, outqueue, params):
 292         data = params['curve'].data[params['block']]
 293         # HACK? rely on params['curve'] being bound to the local hooke
 294         # playlist (i.e. not a copy, as you would get by passing a
 295         # curve through the queue).  Ugh.  Stupid queues.  As an
 296         # alternative, we could pass lookup information through the
 297         # queue...
 298         new = Data((data.shape[0], data.shape[1]+2), dtype=data.dtype)
 299         new.info = copy.deepcopy(data.info)
 300         new[:,:-2] = data
 301         name,dist_units = split_data_label(params['input distance column'])
 302         name,def_units = split_data_label(params['input deflection column'])
 303         new.info['columns'].extend([
 304                 join_data_label(params['output distance column'], dist_units),
 305                 join_data_label(params['output deflection column'], def_units),
 306                 ])
 307         dist_data = data[:,data.info['columns'].index(
 308                 params['input distance column'])]
 309         def_data = data[:,data.info['columns'].index(
 310                 params['input deflection column'])]
 311         i,def_offset,ps = self.find_contact_point(
 312             params['curve'], dist_data, def_data, outqueue)
 313         dist_offset = dist_data[i]
 314         new.info[join_data_label(params['distance info name'], dist_units
 315                                  )] = dist_offset
 316         new.info[join_data_label(params['deflection info name'], def_units
 317                                  )] = def_offset
 318         new.info[params['fit parameters info name']] = ps
 319         new[:,-2] = dist_data - dist_offset
 320         new[:,-1] = def_data - def_offset
 321         params['curve'].data[params['block']] = new
 322
 323     def find_contact_point(self, curve, z_data, d_data, outqueue=None):
 324         """Railyard for the `find_contact_point_*` family.
 325
 326         Uses the `surface contact point algorithm` configuration
 327         setting to call the appropriate backend algorithm.
 328         """
 329         fn = getattr(self, 'find_contact_point_%s'
 330                      % self.plugin.config['surface contact point algorithm'])
 331         return fn(curve, z_data, d_data, outqueue)
 332
 333     def find_contact_point_fmms(self, curve, z_data, d_data, outqueue=None):
 334         """Algorithm by Francesco Musiani and Massimo Sandal.
 335
 336         Notes
 337         -----
 338         Algorithm:
 339
 340         0) Driver-specific workarounds, e.g. deal with the PicoForce
 341           trigger bug by excluding retraction portions with excessive
 342           deviation.
 343         1) Select the second half (non-contact side) of the retraction
 344           curve.
 345         2) Fit the selection to a line.
 346         3) If the fit is not almost horizontal, halve the selection
 347           and retrun to (2).
 348         4) Average the selection and use it as a baseline.
 349         5) Slide in from the start (contact side) of the retraction
 350         curve, until you find a point with greater than baseline
 351         deflection.  That point is the contact point.
 352         """
 353         if curve.info['filetype'] == 'picoforce':
 354             # Take care of the picoforce trigger bug (TODO: example
 355             # data file demonstrating the bug).  We exclude portions
 356             # of the curve that have too much standard deviation.
 357             # Yes, a lot of magic is here.
 358             check_start = len(d_data)-len(d_data)/20
 359             monster_start = len(d_data)
 360             while True:
 361                 # look at the non-contact tail
 362                 non_monster = d_data[check_start:monster_start]
 363                 if non_monster.std() < 2e-10: # HACK: hardcoded cutoff
 364                     break
 365                 else: # move further away from the monster
 366                     check_start -= len(d_data)/50
 367                     monster_start -= len(d_data)/50
 368             z_data = z_data[:monster_start]
 369             d_data = d_data[:monster_start]
 370
 371         # take half of the thing to start
 372         selection_start = len(d_data)/2
 373         while True:
 374             z_chunk = z_data[selection_start:]
 375             d_chunk = d_data[selection_start:]
 376             slope,intercept,r,two_tailed_prob,stderr_of_the_estimate = \
 377                 scipy.stats.linregress(z_chunk, d_chunk)
 378             # We stop if we found an almost-horizontal fit or if we're
 379             # getting to small a selection.  FIXME: 0.1 and 5./6 here
 380             # are "magic numbers" (although reasonable)
 381             if (abs(slope) < 0.1  # deflection (m) / surface (m)
 382                 or selection_start > 5./6*len(d_data)):
 383                 break
 384             selection_start += 10
 385
 386         d_baseline = d_chunk.mean()
 387
 388         # find the first point above the calculated baseline
 389         i = 0
 390         while i < len(d_data) and d_data[i] < ymean:
 391             i += 1
 392         return (i, d_baseline, {})
 393
 394     def find_contact_point_ms(self, curve, z_data, d_data, outqueue=None):
 395         """Algorithm by Massimo Sandal.
 396
 397         Notes
 398         -----
 399         WTK: At least the commits are by Massimo, and I see no notes
 400         attributing the algorithm to anyone else.
 401
 402         Algorithm:
 403
 404         * ?
 405         """
 406         xext=raw_plot.vectors[0][0]
 407         yext=raw_plot.vectors[0][1]
 408         xret2=raw_plot.vectors[1][0]
 409         yret=raw_plot.vectors[1][1]
 410
 411         first_point=[xext[0], yext[0]]
 412         last_point=[xext[-1], yext[-1]]
 413
 414         #regr=scipy.polyfit(first_point, last_point,1)[0:2]
 415         diffx=abs(first_point[0]-last_point[0])
 416         diffy=abs(first_point[1]-last_point[1])
 417
 418         #using polyfit results in numerical errors. good old algebra.
 419         a=diffy/diffx
 420         b=first_point[1]-(a*first_point[0])
 421         baseline=scipy.polyval((a,b), xext)
 422
 423         ysub=[item-basitem for item,basitem in zip(yext,baseline)]
 424
 425         contact=ysub.index(min(ysub))
 426
 427         return xext,ysub,contact
 428
 429         #now, exploit a ClickedPoint instance to calculate index...
 430         dummy=ClickedPoint()
 431         dummy.absolute_coords=(x_intercept,y_intercept)
 432         dummy.find_graph_coords(xret2,yret)
 433
 434         if debug:
 435             return dummy.index, regr, regr_contact
 436         else:
 437             return dummy.index
 438
 439     def find_contact_point_wtk(self, curve, z_data, d_data, outqueue=None):
 440         """Algorithm by W. Trevor King.
 441
 442         Notes
 443         -----
 444         Uses :func:`analyze_surf_pos_data_wtk` internally.
 445         """
 446         reverse = z_data[0] > z_data[-1]
 447         if reverse == True:    # approaching, contact region on the right
 448             d_data = d_data[::-1]
 449         s = SurfacePositionModel(d_data)
 450         s.info['force zero non-contact slope'] = True
 451         offset,contact_slope,surface_index,non_contact_slope = s.fit(
 452             outqueue=outqueue)
 453         info = {
 454             'offset': offset,
 455             'contact slope': contact_slope,
 456             'surface index': surface_index,
 457             'non-contact slope': non_contact_slope,
 458             'reversed': reverse,
 459             }
 460         deflection_offset = offset + contact_slope*surface_index,
 461         if reverse == True:
 462             surface_index = len(d_data)-1-surface_index
 463         return (numpy.round(surface_index), deflection_offset, info)
 464
 465
 466 class ForceCommand (Command):
 467     """Convert a deflection column from meters to newtons.
 468     """
 469     def __init__(self, plugin):
 470         super(ForceCommand, self).__init__(
 471             name='add block force array',
 472             arguments=[
 473                 CurveArgument,
 474                 Argument(name='block', aliases=['set'], type='int', default=0,
 475                          help="""
 476 Data block for which the force should be calculated.  For an
 477 approach/retract force curve, `0` selects the approaching curve and `1`
 478 selects the retracting curve.
 479 """.strip()),
 480                 Argument(name='input deflection column', type='string',
 481                          default='surface deflection (m)',
 482                          help="""
 483 Name of the column to use as the deflection input.
 484 """.strip()),
 485                 Argument(name='output deflection column', type='string',
 486                          default='deflection',
 487                          help="""
 488 Name of the column (without units) to use as the deflection output.
 489 """.strip()),
 490                 Argument(name='spring constant info name', type='string',
 491                          default='spring constant (N/m)',
 492                          help="""
 493 Name of the spring constant in the `.info` dictionary.
 494 """.strip()),
 495                 ],
 496             help=self.__doc__, plugin=plugin)
 497
 498     def _run(self, hooke, inqueue, outqueue, params):
 499         data = params['curve'].data[params['block']]
 500         # HACK? rely on params['curve'] being bound to the local hooke
 501         # playlist (i.e. not a copy, as you would get by passing a
 502         # curve through the queue).  Ugh.  Stupid queues.  As an
 503         # alternative, we could pass lookup information through the
 504         # queue.
 505         new = Data((data.shape[0], data.shape[1]+1), dtype=data.dtype)
 506         new.info = copy.deepcopy(data.info)
 507         new[:,:-1] = data
 508         new.info['columns'].append(
 509             join_data_label(params['output deflection column'], 'N'))
 510         d_data = data[:,data.info['columns'].index(
 511                 params['input deflection column'])]
 512         new[:,-1] = d_data * data.info[params['spring constant info name']]
 513         params['curve'].data[params['block']] = new
 514
 515
 516 class CantileverAdjustedExtensionCommand (Command):
 517     """Remove cantilever extension from a total extension column.
 518     """
 519     def __init__(self, plugin):
 520         super(CantileverAdjustedExtensionCommand, self).__init__(
 521             name='add block cantilever adjusted extension array',
 522             arguments=[
 523                 CurveArgument,
 524                 Argument(name='block', aliases=['set'], type='int', default=0,
 525                          help="""
 526 Data block for which the adjusted extension should be calculated.  For
 527 an approach/retract force curve, `0` selects the approaching curve and
 528 `1` selects the retracting curve.
 529 """.strip()),
 530                 Argument(name='input distance column', type='string',
 531                          default='surface distance (m)',
 532                          help="""
 533 Name of the column to use as the distance input.
 534 """.strip()),
 535                 Argument(name='input deflection column', type='string',
 536                          default='deflection (N)',
 537                          help="""
 538 Name of the column to use as the deflection input.
 539 """.strip()),
 540                 Argument(name='output distance column', type='string',
 541                          default='cantilever adjusted extension',
 542                          help="""
 543 Name of the column (without units) to use as the deflection output.
 544 """.strip()),
 545                 Argument(name='spring constant info name', type='string',
 546                          default='spring constant (N/m)',
 547                          help="""
 548 Name of the spring constant in the `.info` dictionary.
 549 """.strip()),
 550                 ],
 551             help=self.__doc__, plugin=plugin)
 552
 553     def _run(self, hooke, inqueue, outqueue, params):
 554         data = params['curve'].data[params['block']]
 555         # HACK? rely on params['curve'] being bound to the local hooke
 556         # playlist (i.e. not a copy, as you would get by passing a
 557         # curve through the queue).  Ugh.  Stupid queues.  As an
 558         # alternative, we could pass lookup information through the
 559         # queue.
 560         new = Data((data.shape[0], data.shape[1]+1), dtype=data.dtype)
 561         new.info = copy.deepcopy(data.info)
 562         new[:,:-1] = data
 563         new.info['columns'].append(
 564             join_data_label(params['output distance column'], 'm'))
 565         z_data = data[:,data.info['columns'].index(
 566                 params['input distance column'])]
 567         d_data = data[:,data.info['columns'].index(
 568                 params['input deflection column'])]
 569         k = data.info[params['spring constant info name']]
 570
 571         z_name,z_unit = split_data_label(params['input distance column'])
 572         assert z_unit == 'm', params['input distance column']
 573         d_name,d_unit = split_data_label(params['input deflection column'])
 574         assert d_unit == 'N', params['input deflection column']
 575         k_name,k_unit = split_data_label(params['spring constant info name'])
 576         assert k_unit == 'N/m', params['spring constant info name']
 577
 578         new[:,-1] = z_data - d_data / k
 579         params['curve'].data[params['block']] = new
 580
 581
 582 class generalvclampCommands(object):
 583
 584     def plotmanip_flatten(self, plot, current, customvalue=False):
 585         '''
 586         Subtracts a polynomial fit to the non-contact part of the curve, as to flatten it.
 587         the best polynomial fit is chosen among polynomials of degree 1 to n, where n is
 588         given by the configuration file or by the customvalue.
 589
 590         customvalue= int (>0) --> starts the function even if config says no (default=False)
 591         '''
 592
 593         #not a smfs curve...
 594         if current.curve.experiment != 'smfs':
 595             return plot
 596
 597         #only one set is present...
 598         if len(self.plots[0].vectors) != 2:
 599             return plot
 600
 601         #config is not flatten, and customvalue flag is false too
 602         if (not self.config['flatten']) and (not customvalue):
 603             return plot
 604
 605         max_exponent=12
 606         delta_contact=0
 607
 608         if customvalue:
 609             max_cycles=customvalue
 610         else:
 611             max_cycles=self.config['flatten'] #Using > 1 usually doesn't help and can give artefacts. However, it could be useful too.
 612
 613         contact_index=self.find_contact_point()
 614
 615         valn=[[] for item in range(max_exponent)]
 616         yrn=[0.0 for item in range(max_exponent)]
 617         errn=[0.0 for item in range(max_exponent)]
 618
 619         #Check if we have a proper numerical value
 620         try:
 621             zzz=int(max_cycles)
 622         except:
 623             #Loudly and annoyingly complain if it's not a number, then fallback to zero
 624             print '''Warning: flatten value is not a number!
 625             Use "set flatten" or edit hooke.conf to set it properly
 626             Using zero.'''
 627             max_cycles=0
 628
 629         for i in range(int(max_cycles)):
 630
 631             x_ext=plot.vectors[0][0][contact_index+delta_contact:]
 632             y_ext=plot.vectors[0][1][contact_index+delta_contact:]
 633             x_ret=plot.vectors[1][0][contact_index+delta_contact:]
 634             y_ret=plot.vectors[1][1][contact_index+delta_contact:]
 635             for exponent in range(max_exponent):
 636                 try:
 637                     valn[exponent]=sp.polyfit(x_ext,y_ext,exponent)
 638                     yrn[exponent]=sp.polyval(valn[exponent],x_ret)
 639                     errn[exponent]=sp.sqrt(sum((yrn[exponent]-y_ext)**2)/float(len(y_ext)))
 640                 except Exception,e:
 641                     print 'Cannot flatten!'
 642                     print e
 643                     return plot
 644
 645             best_exponent=errn.index(min(errn))
 646
 647             #extension
 648             ycorr_ext=y_ext-yrn[best_exponent]+y_ext[0] #noncontact part
 649             yjoin_ext=np.array(plot.vectors[0][1][0:contact_index+delta_contact]) #contact part
 650             #retraction
 651             ycorr_ret=y_ret-yrn[best_exponent]+y_ext[0] #noncontact part
 652             yjoin_ret=np.array(plot.vectors[1][1][0:contact_index+delta_contact]) #contact part
 653
 654             ycorr_ext=np.concatenate((yjoin_ext, ycorr_ext))
 655             ycorr_ret=np.concatenate((yjoin_ret, ycorr_ret))
 656
 657             plot.vectors[0][1]=list(ycorr_ext)
 658             plot.vectors[1][1]=list(ycorr_ret)
 659
 660         return plot
 661