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
  10 # modify it under the terms of the GNU Lesser General Public
  11 # License as published by the Free Software Foundation, either
  12 # version 3 of the License, or (at your option) any later version.
  13 #
  14 # Hooke is distributed in the hope that it will be useful,
  15 # but WITHOUT ANY WARRANTY; without even the implied warranty of
  16 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  17 # GNU Lesser General 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 Builtin
  37 from ..util.fit import PoorFit, ModelFitter
  38 from .curve import CurveArgument
  39
  40
  41 def scale(hooke, curve):
  42     commands = hooke.commands
  43     contact = [c for c in hooke.commands
  44                if c.name == 'zero block surface contact point'][0]
  45     force = [c for c in hooke.commands if c.name == 'add block force array'][0]
  46     inqueue = None
  47     outqueue = NullQueue()
  48     for i,block in enumerate(curve.data):
  49         numpy.savetxt(open('curve.dat', 'w'), block, delimiter='\t')
  50         params = {'curve':curve, 'block':i}
  51         try:
  52             contact._run(hooke, inqueue, outqueue, params)
  53         except PoorFit, e:
  54             raise PoorFit('Could not fit %s %s: %s'
  55                           % (curve.path, i, str(e)))
  56         force._run(hooke, inqueue, outqueue, params)
  57     return curve
  58
  59 class SurfacePositionModel (ModelFitter):
  60     """
  61
  62     The bilinear model is symmetric, but the parameter guessing and
  63     sanity checks assume the contact region occurs for lower indicies
  64     ("left of") the non-contact region.  We also assume that
  65     tip-surface attractions produce positive deflections.
  66
  67     Notes
  68     -----
  69     Algorithm borrowed from WTK's `piezo package`_, specifically
  70     from :func:`piezo.z_piezo_utils.analyzeSurfPosData`.
  71
  72     .. _piezo package:
  73       http://www.physics.drexel.edu/~wking/code/git/git.php?p=piezo.git
  74
  75     Fits the data to the bilinear :method:`model`.
  76
  77     In order for this model to produce a satisfactory fit, there
  78     should be enough data in the off-surface region that interactions
  79     due to proteins, etc. will not seriously skew the fit in the
  80     off-surface region.
  81
  82
  83     We guess
  84     """
  85     def model(self, params):
  86         """A continuous, bilinear model.
  87
  88         Notes
  89         -----
  90         .. math::
  91
  92           y = \begin{cases}
  93             p_0 + p_1 x                 & \text{if $x <= p_2$}, \\
  94             p_0 + p_1 p_2 + p_3 (x-p_2) & \text{if $x >= p_2$}.
  95               \end{cases}
  96
  97         Where :math:`p_0` is a vertical offset, :math:`p_1` is the slope
  98         of the first region, :math:`p_2` is the transition location, and
  99         :math:`p_3` is the slope of the second region.
 100         """
 101         p = params  # convenient alias
 102         if self.info.get('force zero non-contact slope', None) == True:
 103             p = list(p)
 104             p.append(0.)  # restore the non-contact slope parameter
 105         r2 = numpy.round(abs(p[2]))
 106         if r2 >= 1:
 107             self._model_data[:r2] = p[0] + p[1] * numpy.arange(r2)
 108         if r2 < len(self._data)-1:
 109             self._model_data[r2:] = \
 110                 p[0] + p[1]*p[2] + p[3] * numpy.arange(len(self._data)-r2)
 111         return self._model_data
 112
 113     def set_data(self, data, info=None):
 114         super(SurfacePositionModel, self).set_data(data, info)
 115         if info == None:
 116             info = {}
 117         self.info = info
 118         self.info['min position'] = 0
 119         self.info['max position'] = len(data)
 120         self.info['max deflection'] = data.max()
 121         self.info['min deflection'] = data.min()
 122         self.info['position range'] = self.info['max position'] - self.info['min position']
 123         self.info['deflection range'] = self.info['max deflection'] - self.info['min deflection']
 124
 125     def guess_initial_params(self, outqueue=None):
 126         """Guess the initial parameters.
 127
 128         Notes
 129         -----
 130         We guess initial parameters such that the offset (:math:`p_1`)
 131         matches the minimum deflection, the kink (:math:`p_2`) occurs in
 132         the middle of the data, the initial (contact) slope (:math:`p_0`)
 133         produces the maximum deflection at the left-most point, and the
 134         final (non-contact) slope (:math:`p_3`) is zero.
 135         """
 136         left_offset = self.info['min deflection']
 137         left_slope = 2*(self.info['deflection range']
 138                         /self.info['position range'])
 139         kink_position = (self.info['max position']
 140                          +self.info['min position'])/2.0
 141         right_slope = 0
 142         self.info['guessed contact slope'] = left_slope
 143         params = [left_offset, left_slope, kink_position, right_slope]
 144         if self.info.get('force zero non-contact slope', None) == True:
 145             params = params[:-1]
 146         return params
 147
 148     def guess_scale(self, params, outqueue=None):
 149         """Guess the parameter scales.
 150
 151         Notes
 152         -----
 153
 154         We guess offset scale (:math:`p_0`) as one tenth of the total
 155         deflection range, the kink scale (:math:`p_2`) as one tenth of
 156         the total index range, the initial (contact) slope scale
 157         (:math:`p_1`) as one tenth of the contact slope estimation,
 158         and the final (non-contact) slope scale (:math:`p_3`) is as
 159         one tenth of the initial slope scale.
 160         """
 161         offset_scale = self.info['deflection range']/10.
 162         left_slope_scale = abs(params[1])/10.
 163         kink_scale = self.info['position range']/10.
 164         right_slope_scale = left_slope_scale/10.
 165         scale = [offset_scale, left_slope_scale, kink_scale, right_slope_scale]
 166         if self.info.get('force zero non-contact slope', None) == True:
 167             scale = scale[:-1]
 168         return scale
 169
 170     def fit(self, *args, **kwargs):
 171         self.info['guessed contact slope'] = None
 172         params = super(SurfacePositionModel, self).fit(*args, **kwargs)
 173         params[2] = abs(params[2])
 174         if self.info.get('force zero non-contact slope', None) == True:
 175             params = list(params)
 176             params.append(0.)  # restore the non-contact slope parameter
 177
 178         # check that the fit is reasonable, see the :meth:`model` docstring
 179         # for parameter descriptions.  HACK: hardcoded cutoffs.
 180         if abs(params[3]*10) > abs(params[1]) :
 181             raise PoorFit('Slope in non-contact region, or no slope in contact')
 182         if params[2] < self.info['min position']+0.02*self.info['position range']:
 183             raise PoorFit(
 184                 'No kink (kink %g less than %g, need more space to left)'
 185                 % (params[2],
 186                    self.info['min position']+0.02*self.info['position range']))
 187         if params[2] > self.info['max position']-0.02*self.info['position range']:
 188             raise poorFit(
 189                 'No kink (kink %g more than %g, need more space to right)'
 190                 % (params[2],
 191                    self.info['max position']-0.02*self.info['position range']))
 192         if (self.info['guessed contact slope'] != None
 193             and abs(params[1]) < 0.5 * abs(self.info['guessed contact slope'])):
 194             raise PoorFit('Too far (contact slope %g, but expected ~%g'
 195                           % (params[3], self.info['guessed contact slope']))
 196         return params
 197
 198 class VelocityClampPlugin (Builtin):
 199     def __init__(self):
 200         super(VelocityClampPlugin, self).__init__(name='vclamp')
 201         self._commands = [
 202             SurfaceContactCommand(self), ForceCommand(self),
 203             ]
 204
 205     def default_settings(self):
 206         return [
 207             Setting(section=self.setting_section, help=self.__doc__),
 208             Setting(section=self.setting_section,
 209                     option='surface contact point algorithm',
 210                     value='wtk',
 211                     help='Select the surface contact point algorithm.  See the documentation for descriptions of available algorithms.')
 212             ]
 213
 214
 215 class SurfaceContactCommand (Command):
 216     """Automatically determine a block's surface contact point.
 217
 218     Uses the block's `z piezo (m)` and `deflection (m)` arrays.
 219     Stores the contact parameters in `block.info`'s `surface distance
 220     offset (m)` and `surface deflection offset (m)`.  Model-specific
 221     fitting information is stored in `surface detection parameters`.
 222
 223     The adjusted data columns `surface distance (m)` and `surface
 224     adjusted deflection (m)` are also added to the block.
 225
 226     You can select the contact point algorithm with the creatively
 227     named `surface contact point algorithm` configuration setting.
 228     Currently available options are:
 229
 230     * fmms (:meth:`find_contact_point_fmms`)
 231     * ms (:meth:`find_contact_point_ms`)
 232     * wtk (:meth:`find_contact_point_wtk`)
 233     """
 234     def __init__(self, plugin):
 235         super(SurfaceContactCommand, self).__init__(
 236             name='zero block surface contact point',
 237             arguments=[
 238                 CurveArgument,
 239                 Argument(name='block', aliases=['set'], type='int', default=0,
 240                          help="""
 241 Data block for which the force should be calculated.  For an
 242 approach/retract force curve, `0` selects the approaching curve and `1`
 243 selects the retracting curve.
 244 """.strip()),
 245                 ],
 246             help=self.__doc__, plugin=plugin)
 247
 248     def _run(self, hooke, inqueue, outqueue, params):
 249         data = params['curve'].data[int(params['block'])] # HACK, int() should be handled by ui
 250         # HACK? rely on params['curve'] being bound to the local hooke
 251         # playlist (i.e. not a copy, as you would get by passing a
 252         # curve through the queue).  Ugh.  Stupid queues.  As an
 253         # alternative, we could pass lookup information through the
 254         # queue...
 255         new = Data((data.shape[0], data.shape[1]+2), dtype=data.dtype)
 256         new.info = copy.deepcopy(data.info)
 257         new[:,:-2] = data
 258         new.info['columns'].extend(
 259             ['surface distance (m)', 'surface adjusted deflection (m)'])
 260         z_data = data[:,data.info['columns'].index('z piezo (m)')]
 261         d_data = data[:,data.info['columns'].index('deflection (m)')]
 262         i,deflection_offset,ps = self.find_contact_point(
 263             params['curve'], z_data, d_data, outqueue)
 264         surface_offset = z_data[i]
 265         new.info['surface distance offset (m)'] = surface_offset
 266         new.info['surface deflection offset (m)'] = deflection_offset
 267         new.info['surface detection parameters'] = ps
 268         new[:,-2] = z_data - surface_offset
 269         new[:,-1] = d_data - deflection_offset
 270         data = params['curve'].data[int(params['block'])] # HACK, int() should be handled by ui
 271         params['curve'].data[int(params['block'])] = new # HACK, int() should be handled by ui
 272
 273     def find_contact_point(self, curve, z_data, d_data, outqueue=None):
 274         """Railyard for the `find_contact_point_*` family.
 275
 276         Uses the `surface contact point algorithm` configuration
 277         setting to call the appropriate backend algorithm.
 278         """
 279         fn = getattr(self, 'find_contact_point_%s'
 280                      % self.plugin.config['surface contact point algorithm'])
 281         return fn(curve, z_data, d_data, outqueue)
 282
 283     def find_contact_point_fmms(self, curve, z_data, d_data, outqueue=None):
 284         """Algorithm by Francesco Musiani and Massimo Sandal.
 285
 286         Notes
 287         -----
 288         Algorithm:
 289
 290         0) Driver-specific workarounds, e.g. deal with the PicoForce
 291           trigger bug by excluding retraction portions with excessive
 292           deviation.
 293         1) Select the second half (non-contact side) of the retraction
 294           curve.
 295         2) Fit the selection to a line.
 296         3) If the fit is not almost horizontal, halve the selection
 297           and retrun to (2).
 298         4) Average the selection and use it as a baseline.
 299         5) Slide in from the start (contact side) of the retraction
 300         curve, until you find a point with greater than baseline
 301         deflection.  That point is the contact point.
 302         """
 303         if curve.info['filetype'] == 'picoforce':
 304             # Take care of the picoforce trigger bug (TODO: example
 305             # data file demonstrating the bug).  We exclude portions
 306             # of the curve that have too much standard deviation.
 307             # Yes, a lot of magic is here.
 308             check_start = len(d_data)-len(d_data)/20
 309             monster_start = len(d_data)
 310             while True:
 311                 # look at the non-contact tail
 312                 non_monster = d_data[check_start:monster_start]
 313                 if non_monster.std() < 2e-10: # HACK: hardcoded cutoff
 314                     break
 315                 else: # move further away from the monster
 316                     check_start -= len(d_data)/50
 317                     monster_start -= len(d_data)/50
 318             z_data = z_data[:monster_start]
 319             d_data = d_data[:monster_start]
 320
 321         # take half of the thing to start
 322         selection_start = len(d_data)/2
 323         while True:
 324             z_chunk = z_data[selection_start:]
 325             d_chunk = d_data[selection_start:]
 326             slope,intercept,r,two_tailed_prob,stderr_of_the_estimate = \
 327                 scipy.stats.linregress(z_chunk, d_chunk)
 328             # We stop if we found an almost-horizontal fit or if we're
 329             # getting to small a selection.  FIXME: 0.1 and 5./6 here
 330             # are "magic numbers" (although reasonable)
 331             if (abs(slope) < 0.1  # deflection (m) / surface (m)
 332                 or selection_start > 5./6*len(d_data)):
 333                 break
 334             selection_start += 10
 335
 336         d_baseline = d_chunk.mean()
 337
 338         # find the first point above the calculated baseline
 339         i = 0
 340         while i < len(d_data) and d_data[i] < ymean:
 341             i += 1
 342         return (i, d_baseline, {})
 343
 344     def find_contact_point_ms(self, curve, z_data, d_data, outqueue=None):
 345         """Algorithm by Massimo Sandal.
 346
 347         Notes
 348         -----
 349         WTK: At least the commits are by Massimo, and I see no notes
 350         attributing the algorithm to anyone else.
 351
 352         Algorithm:
 353
 354         * ?
 355         """
 356         xext=raw_plot.vectors[0][0]
 357         yext=raw_plot.vectors[0][1]
 358         xret2=raw_plot.vectors[1][0]
 359         yret=raw_plot.vectors[1][1]
 360
 361         first_point=[xext[0], yext[0]]
 362         last_point=[xext[-1], yext[-1]]
 363
 364         #regr=scipy.polyfit(first_point, last_point,1)[0:2]
 365         diffx=abs(first_point[0]-last_point[0])
 366         diffy=abs(first_point[1]-last_point[1])
 367
 368         #using polyfit results in numerical errors. good old algebra.
 369         a=diffy/diffx
 370         b=first_point[1]-(a*first_point[0])
 371         baseline=scipy.polyval((a,b), xext)
 372
 373         ysub=[item-basitem for item,basitem in zip(yext,baseline)]
 374
 375         contact=ysub.index(min(ysub))
 376
 377         return xext,ysub,contact
 378
 379         #now, exploit a ClickedPoint instance to calculate index...
 380         dummy=ClickedPoint()
 381         dummy.absolute_coords=(x_intercept,y_intercept)
 382         dummy.find_graph_coords(xret2,yret)
 383
 384         if debug:
 385             return dummy.index, regr, regr_contact
 386         else:
 387             return dummy.index
 388
 389     def find_contact_point_wtk(self, curve, z_data, d_data, outqueue=None):
 390         """Algorithm by W. Trevor King.
 391
 392         Notes
 393         -----
 394         Uses :func:`analyze_surf_pos_data_wtk` internally.
 395         """
 396         reverse = z_data[0] > z_data[-1]
 397         if reverse == True:    # approaching, contact region on the right
 398             d_data = d_data[::-1]
 399         s = SurfacePositionModel(d_data)
 400         s.info['force zero non-contact slope'] = True
 401         offset,contact_slope,surface_index,non_contact_slope = s.fit(
 402             outqueue=outqueue)
 403         info = {
 404             'offset': offset,
 405             'contact slope': contact_slope,
 406             'surface index': surface_index,
 407             'non-contact slope': non_contact_slope,
 408             'reversed': reverse,
 409             }
 410         deflection_offset = offset + contact_slope*surface_index,
 411         if reverse == True:
 412             surface_index = len(d_data)-1-surface_index
 413         return (numpy.round(surface_index), deflection_offset, info)
 414
 415 class ForceCommand (Command):
 416     """Calculate a block's `deflection (N)` array.
 417
 418     Uses the block's `deflection (m)` array and `spring constant
 419     (N/m)`.
 420     """
 421     def __init__(self, plugin):
 422         super(ForceCommand, self).__init__(
 423             name='add block force array',
 424             arguments=[
 425                 CurveArgument,
 426                 Argument(name='block', aliases=['set'], type='int', default=0,
 427                          help="""
 428 Data block for which the force should be calculated.  For an
 429 approach/retract force curve, `0` selects the approaching curve and `1`
 430 selects the retracting curve.
 431 """.strip()),
 432                 ],
 433             help=self.__doc__, plugin=plugin)
 434
 435     def _run(self, hooke, inqueue, outqueue, params):
 436         data = params['curve'].data[int(params['block'])] # HACK, int() should be handled by ui
 437         # HACK? rely on params['curve'] being bound to the local hooke
 438         # playlist (i.e. not a copy, as you would get by passing a
 439         # curve through the queue).  Ugh.  Stupid queues.  As an
 440         # alternative, we could pass lookup information through the
 441         # queue.
 442         new = Data((data.shape[0], data.shape[1]+1), dtype=data.dtype)
 443         new.info = copy.deepcopy(data.info)
 444         new[:,:-1] = data
 445         new.info['columns'].append('deflection (N)')
 446         d_data = data[:,data.info['columns'].index('surface adjusted deflection (m)')]
 447         new[:,-1] = d_data * data.info['spring constant (N/m)']
 448         params['curve'].data[int(params['block'])] = new # HACK, int() should be handled by ui
 449
 450
 451 class generalvclampCommands(object):
 452
 453     def do_subtplot(self, args):
 454         '''
 455         SUBTPLOT
 456         (procplots.py plugin)
 457         Plots the difference between ret and ext current curve
 458         -------
 459         Syntax: subtplot
 460         '''
 461         #FIXME: sub_filter and sub_order must be args
 462
 463         if len(self.plots[0].vectors) != 2:
 464             print 'This command only works on a curve with two different plots.'
 465             pass
 466
 467         outplot=self.subtract_curves(sub_order=1)
 468
 469         plot_graph=self.list_of_events['plot_graph']
 470         wx.PostEvent(self.frame,plot_graph(plots=[outplot]))
 471
 472     def _plug_init(self):
 473         self.basecurrent=None
 474         self.basepoints=None
 475         self.autofile=''
 476
 477     def do_distance(self,args):
 478         '''
 479         DISTANCE
 480         (generalvclamp.py)
 481         Measure the distance (in nm) between two points.
 482         For a standard experiment this is the delta X distance.
 483         For a force clamp experiment this is the delta Y distance (actually becomes
 484         an alias of zpiezo)
 485         -----------------
 486         Syntax: distance
 487         '''
 488         if self.current.curve.experiment == 'clamp':
 489             print 'You wanted to use zpiezo perhaps?'
 490             return
 491         else:
 492             dx,unitx,dy,unity=self._delta(set=1)
 493             print str(dx*(10**9))+' nm'
 494             to_dump='distance '+self.current.path+' '+str(dx*(10**9))+' nm'
 495             self.outlet.push(to_dump)
 496
 497
 498     def do_force(self,args):
 499         '''
 500         FORCE
 501         (generalvclamp.py)
 502         Measure the force difference (in pN) between two points
 503         ---------------
 504         Syntax: force
 505         '''
 506         if self.current.curve.experiment == 'clamp':
 507             print 'This command makes no sense for a force clamp experiment.'
 508             return
 509         dx,unitx,dy,unity=self._delta(set=1)
 510         print str(dy*(10**12))+' pN'
 511         to_dump='force '+self.current.path+' '+str(dy*(10**12))+' pN'
 512         self.outlet.push(to_dump)
 513
 514
 515     def do_forcebase(self,args):
 516         '''
 517         FORCEBASE
 518         (generalvclamp.py)
 519         Measures the difference in force (in pN) between a point and a baseline
 520         took as the average between two points.
 521
 522         The baseline is fixed once for a given curve and different force measurements,
 523         unless the user wants it to be recalculated
 524         ------------
 525         Syntax: forcebase [rebase]
 526                 rebase: Forces forcebase to ask again the baseline
 527                 max: Instead of asking for a point to measure, asks for two points and use
 528                      the maximum peak in between
 529         '''
 530         rebase=False #if true=we select rebase
 531         maxpoint=False #if true=we measure the maximum peak
 532
 533         plot=self._get_displayed_plot()
 534         whatset=1 #fixme: for all sets
 535         if 'rebase' in args or (self.basecurrent != self.current.path):
 536             rebase=True
 537         if 'max' in args:
 538             maxpoint=True
 539
 540         if rebase:
 541             print 'Select baseline'
 542             self.basepoints=self._measure_N_points(N=2, whatset=whatset)
 543             self.basecurrent=self.current.path
 544
 545         if maxpoint:
 546             print 'Select two points'
 547             points=self._measure_N_points(N=2, whatset=whatset)
 548             boundpoints=[points[0].index, points[1].index]
 549             boundpoints.sort()
 550             try:
 551                 y=min(plot.vectors[whatset][1][boundpoints[0]:boundpoints[1]])
 552             except ValueError:
 553                 print 'Chosen interval not valid. Try picking it again. Did you pick the same point as begin and end of interval?'
 554         else:
 555             print 'Select point to measure'
 556             points=self._measure_N_points(N=1, whatset=whatset)
 557             #whatplot=points[0].dest
 558             y=points[0].graph_coords[1]
 559
 560         #fixme: code duplication
 561         boundaries=[self.basepoints[0].index, self.basepoints[1].index]
 562         boundaries.sort()
 563         to_average=plot.vectors[whatset][1][boundaries[0]:boundaries[1]] #y points to average
 564
 565         avg=np.mean(to_average)
 566         forcebase=abs(y-avg)
 567         print str(forcebase*(10**12))+' pN'
 568         to_dump='forcebase '+self.current.path+' '+str(forcebase*(10**12))+' pN'
 569         self.outlet.push(to_dump)
 570
 571     def plotmanip_multiplier(self, plot, current):
 572         '''
 573         Multiplies all the Y values of an SMFS curve by a value stored in the 'force_multiplier'
 574         configuration variable. Useful for calibrations and other stuff.
 575         '''
 576
 577         #not a smfs curve...
 578         if current.curve.experiment != 'smfs':
 579             return plot
 580
 581         #only one set is present...
 582         if len(self.plots[0].vectors) != 2:
 583             return plot
 584
 585         #multiplier is 1...
 586         if (self.config['force_multiplier']==1):
 587             return plot
 588
 589         for i in range(len(plot.vectors[0][1])):
 590             plot.vectors[0][1][i]=plot.vectors[0][1][i]*self.config['force_multiplier']
 591
 592         for i in range(len(plot.vectors[1][1])):
 593             plot.vectors[1][1][i]=plot.vectors[1][1][i]*self.config['force_multiplier']
 594
 595         return plot
 596
 597
 598     def plotmanip_flatten(self, plot, current, customvalue=False):
 599         '''
 600         Subtracts a polynomial fit to the non-contact part of the curve, as to flatten it.
 601         the best polynomial fit is chosen among polynomials of degree 1 to n, where n is
 602         given by the configuration file or by the customvalue.
 603
 604         customvalue= int (>0) --> starts the function even if config says no (default=False)
 605         '''
 606
 607         #not a smfs curve...
 608         if current.curve.experiment != 'smfs':
 609             return plot
 610
 611         #only one set is present...
 612         if len(self.plots[0].vectors) != 2:
 613             return plot
 614
 615         #config is not flatten, and customvalue flag is false too
 616         if (not self.config['flatten']) and (not customvalue):
 617             return plot
 618
 619         max_exponent=12
 620         delta_contact=0
 621
 622         if customvalue:
 623             max_cycles=customvalue
 624         else:
 625             max_cycles=self.config['flatten'] #Using > 1 usually doesn't help and can give artefacts. However, it could be useful too.
 626
 627         contact_index=self.find_contact_point()
 628
 629         valn=[[] for item in range(max_exponent)]
 630         yrn=[0.0 for item in range(max_exponent)]
 631         errn=[0.0 for item in range(max_exponent)]
 632
 633         #Check if we have a proper numerical value
 634         try:
 635             zzz=int(max_cycles)
 636         except:
 637             #Loudly and annoyingly complain if it's not a number, then fallback to zero
 638             print '''Warning: flatten value is not a number!
 639             Use "set flatten" or edit hooke.conf to set it properly
 640             Using zero.'''
 641             max_cycles=0
 642
 643         for i in range(int(max_cycles)):
 644
 645             x_ext=plot.vectors[0][0][contact_index+delta_contact:]
 646             y_ext=plot.vectors[0][1][contact_index+delta_contact:]
 647             x_ret=plot.vectors[1][0][contact_index+delta_contact:]
 648             y_ret=plot.vectors[1][1][contact_index+delta_contact:]
 649             for exponent in range(max_exponent):
 650                 try:
 651                     valn[exponent]=sp.polyfit(x_ext,y_ext,exponent)
 652                     yrn[exponent]=sp.polyval(valn[exponent],x_ret)
 653                     errn[exponent]=sp.sqrt(sum((yrn[exponent]-y_ext)**2)/float(len(y_ext)))
 654                 except Exception,e:
 655                     print 'Cannot flatten!'
 656                     print e
 657                     return plot
 658
 659             best_exponent=errn.index(min(errn))
 660
 661             #extension
 662             ycorr_ext=y_ext-yrn[best_exponent]+y_ext[0] #noncontact part
 663             yjoin_ext=np.array(plot.vectors[0][1][0:contact_index+delta_contact]) #contact part
 664             #retraction
 665             ycorr_ret=y_ret-yrn[best_exponent]+y_ext[0] #noncontact part
 666             yjoin_ret=np.array(plot.vectors[1][1][0:contact_index+delta_contact]) #contact part
 667
 668             ycorr_ext=np.concatenate((yjoin_ext, ycorr_ext))
 669             ycorr_ret=np.concatenate((yjoin_ret, ycorr_ret))
 670
 671             plot.vectors[0][1]=list(ycorr_ext)
 672             plot.vectors[1][1]=list(ycorr_ret)
 673
 674         return plot
 675
 676     #---SLOPE---
 677     def do_slope(self,args):
 678         '''
 679         SLOPE
 680         (generalvclamp.py)
 681         Measures the slope of a delimited chunk on the return trace.
 682         The chunk can be delimited either by two manual clicks, or have
 683         a fixed width, given as an argument.
 684         ---------------
 685         Syntax: slope [width]
 686                 The facultative [width] parameter specifies how many
 687                 points will be considered for the fit. If [width] is
 688                 specified, only one click will be required.
 689         (c) Marco Brucale, Massimo Sandal 2008
 690         '''
 691
 692         # Reads the facultative width argument
 693         try:
 694             fitspan=int(args)
 695         except:
 696             fitspan=0
 697
 698         # Decides between the two forms of user input, as per (args)
 699         if fitspan == 0:
 700             # Gets the Xs of two clicked points as indexes on the current curve vector
 701             print 'Click twice to delimit chunk'
 702             points=self._measure_N_points(N=2,whatset=1)
 703         else:
 704             print 'Click once on the leftmost point of the chunk (i.e.usually the peak)'
 705             points=self._measure_N_points(N=1,whatset=1)
 706
 707         slope=self._slope(points,fitspan)
 708
 709         # Outputs the relevant slope parameter
 710         print 'Slope:'
 711         print str(slope)
 712         to_dump='slope '+self.current.path+' '+str(slope)
 713         self.outlet.push(to_dump)
 714
 715     def _slope(self,points,fitspan):
 716         # Calls the function linefit_between
 717         parameters=[0,0,[],[]]
 718         try:
 719             clickedpoints=[points[0].index,points[1].index]
 720             clickedpoints.sort()
 721         except:
 722             clickedpoints=[points[0].index-fitspan,points[0].index]
 723
 724         try:
 725             parameters=self.linefit_between(clickedpoints[0],clickedpoints[1])
 726         except:
 727             print 'Cannot fit. Did you click twice the same point?'
 728             return
 729
 730         # Outputs the relevant slope parameter
 731         print 'Slope:'
 732         print str(parameters[0])
 733         to_dump='slope '+self.curve.path+' '+str(parameters[0])
 734         self.outlet.push(to_dump)
 735
 736         # Makes a vector with the fitted parameters and sends it to the GUI
 737         xtoplot=parameters[2]
 738         ytoplot=[]
 739         x=0
 740         for x in xtoplot:
 741             ytoplot.append((x*parameters[0])+parameters[1])
 742
 743         clickvector_x, clickvector_y=[], []
 744         for item in points:
 745             clickvector_x.append(item.graph_coords[0])
 746             clickvector_y.append(item.graph_coords[1])
 747
 748         lineplot=self._get_displayed_plot(0) #get topmost displayed plot
 749
 750         lineplot.add_set(xtoplot,ytoplot)
 751         lineplot.add_set(clickvector_x, clickvector_y)
 752
 753
 754         if lineplot.styles==[]:
 755             lineplot.styles=[None,None,None,'scatter']
 756         else:
 757             lineplot.styles+=[None,'scatter']
 758         if lineplot.colors==[]:
 759             lineplot.colors=[None,None,'black',None]
 760         else:
 761             lineplot.colors+=['black',None]
 762
 763
 764         self._send_plot([lineplot])
 765
 766         return parameters[0]
 767
 768
 769     def linefit_between(self,index1,index2,whatset=1):
 770         '''
 771         Creates two vectors (xtofit,ytofit) slicing out from the
 772         current return trace a portion delimited by the two indexes
 773         given as arguments.
 774         Then does a least squares linear fit on that slice.
 775         Finally returns [0]=the slope, [1]=the intercept of the
 776         fitted 1st grade polynomial, and [2,3]=the actual (x,y) vectors
 777         used for the fit.
 778         (c) Marco Brucale, Massimo Sandal 2008
 779         '''
 780         # Translates the indexes into two vectors containing the x,y data to fit
 781         xtofit=self.plots[0].vectors[whatset][0][index1:index2]
 782         ytofit=self.plots[0].vectors[whatset][1][index1:index2]
 783
 784         # Does the actual linear fitting (simple least squares with numpy.polyfit)
 785         linefit=[]
 786         linefit=np.polyfit(xtofit,ytofit,1)
 787
 788         return (linefit[0],linefit[1],xtofit,ytofit)
 789
 790
 791     def fit_interval_nm(self,start_index,plot,nm,backwards):
 792           '''
 793           Calculates the number of points to fit, given a fit interval in nm
 794           start_index: index of point
 795           plot: plot to use
 796           backwards: if true, finds a point backwards.
 797           '''
 798           whatset=1 #FIXME: should be decidable
 799           x_vect=plot.vectors[1][0]
 800
 801           c=0
 802           i=start_index
 803           start=x_vect[start_index]
 804           maxlen=len(x_vect)
 805           while abs(x_vect[i]-x_vect[start_index])*(10**9) < nm:
 806               if i==0 or i==maxlen-1: #we reached boundaries of vector!
 807                   return c
 808
 809               if backwards:
 810                   i-=1
 811               else:
 812                   i+=1
 813               c+=1
 814           return c
 815
 816
 817
 818     def find_current_peaks(self,noflatten, a=True, maxpeak=True):
 819             #Find peaks.
 820             if a==True:
 821                   a=self.convfilt_config['mindeviation']
 822             try:
 823                   abs_devs=float(a)
 824             except:
 825                   print "Bad input, using default."
 826                   abs_devs=self.convfilt_config['mindeviation']
 827
 828             defplot=self.current.curve.default_plots()[0]
 829             if not noflatten:
 830                 flatten=self._find_plotmanip('flatten') #Extract flatten plotmanip
 831                 defplot=flatten(defplot, self.current, customvalue=1) #Flatten curve before feeding it to has_peaks
 832             pk_location,peak_size=self.has_peaks(defplot, abs_devs, maxpeak)
 833             return pk_location, peak_size
 834
 835
 836     def pickup_contact_point(self,N=1,whatset=1):
 837         '''macro to pick up the contact point by clicking'''
 838         contact_point=self._measure_N_points(N=1, whatset=1)[0]
 839         contact_point_index=contact_point.index
 840         self.wlccontact_point=contact_point
 841         self.wlccontact_index=contact_point.index
 842         self.wlccurrent=self.current.path
 843         return contact_point, contact_point_index
 844
 845
 846     def baseline_points(self,peak_location, displayed_plot):
 847         clicks=self.config['baseline_clicks']
 848         if clicks==0:
 849             self.basepoints=[]
 850             base_index_0=peak_location[-1]+self.fit_interval_nm(peak_location[-1], displayed_plot, self.config['auto_right_baseline'],False)
 851             self.basepoints.append(self._clickize(displayed_plot.vectors[1][0],displayed_plot.vectors[1][1],base_index_0))
 852             base_index_1=self.basepoints[0].index+self.fit_interval_nm(self.basepoints[0].index, displayed_plot, self.config['auto_left_baseline'],False)
 853             self.basepoints.append(self._clickize(displayed_plot.vectors[1][0],displayed_plot.vectors[1][1],base_index_1))
 854         elif clicks>0:
 855             print 'Select baseline'
 856             if clicks==1:
 857                 self.basepoints=self._measure_N_points(N=1, whatset=1)
 858                 base_index_1=self.basepoints[0].index+self.fit_interval_nm(self.basepoints[0].index, displayed_plot, self.config['auto_left_baseline'], False)
 859                 self.basepoints.append(self._clickize(displayed_plot.vectors[1][0],displayed_plot.vectors[1][1],base_index_1))
 860             else:
 861                 self.basepoints=self._measure_N_points(N=2, whatset=1)
 862
 863         self.basecurrent=self.current.path
 864         return self.basepoints