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