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