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['max 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 params[3] < 0.5 * self.info['guessed contact slope']):
 186             raise PoorFit('Too far')
 187         return params
 188
 189 class VelocityClampPlugin (Builtin):
 190     def __init__(self):
 191         super(VelocityClampPlugin, self).__init__(name='vclamp')
 192         self._commands = [
 193             SurfaceContactCommand(self), ForceCommand(self),
 194             ]
 195
 196     def default_settings(self):
 197         return [
 198             Setting(section=self.setting_section, help=self.__doc__),
 199             Setting(section=self.setting_section,
 200                     option='surface contact point algorithm',
 201                     value='wtk',
 202                     help='Select the surface contact point algorithm.  See the documentation for descriptions of available algorithms.')
 203             ]
 204
 205
 206 class SurfaceContactCommand (Command):
 207     """Automatically determine a block's surface contact point.
 208
 209     Uses the block's `z piezo (m)` and `deflection (m)` arrays.
 210     Stores the contact parameters in `block.info`'s `surface distance
 211     offset (m)` and `surface deflection offset (m)`.  Model-specific
 212     fitting information is stored in `surface detection parameters`.
 213
 214     The adjusted data columns `surface distance (m)` and `surface
 215     adjusted deflection (m)` are also added to the block.
 216     """
 217     def __init__(self, plugin):
 218         super(SurfaceContactCommand, self).__init__(
 219             name='zero block surface contact point',
 220             arguments=[
 221                 CurveArgument,
 222                 Argument(name='block', aliases=['set'], type='int', default=0,
 223                          help="""
 224 Data block for which the force should be calculated.  For an
 225 approach/retract force curve, `0` selects the approaching curve and `1`
 226 selects the retracting curve.
 227 """.strip()),
 228                 ],
 229             help=self.__doc__, plugin=plugin)
 230
 231     def _run(self, hooke, inqueue, outqueue, params):
 232         data = params['curve'].data[int(params['block'])] # HACK, int() should be handled by ui
 233         # HACK? rely on params['curve'] being bound to the local hooke
 234         # playlist (i.e. not a copy, as you would get by passing a
 235         # curve through the queue).  Ugh.  Stupid queues.  As an
 236         # alternative, we could pass lookup information through the
 237         # queue...
 238         new = Data((data.shape[0], data.shape[1]+2), dtype=data.dtype)
 239         new.info = copy.deepcopy(data.info)
 240         new[:,:-2] = data
 241         new.info['columns'].extend(
 242             ['surface distance (m)', 'surface adjusted deflection (m)'])
 243         z_data = data[:,data.info['columns'].index('z piezo (m)')]
 244         d_data = data[:,data.info['columns'].index('deflection (m)')]
 245         i,deflection_offset,ps = self._find_contact_point(
 246             z_data, d_data, outqueue)
 247         surface_offset = z_data[i]
 248         new.info['surface distance offset (m)'] = surface_offset
 249         new.info['surface deflection offset (m)'] = deflection_offset
 250         new.info['surface detection parameters'] = ps
 251         new[:,-2] = z_data - surface_offset
 252         new[:,-1] = d_data - deflection_offset
 253         data = params['curve'].data[int(params['block'])] # HACK, int() should be handled by ui
 254         params['curve'].data[int(params['block'])] = new # HACK, int() should be handled by ui
 255
 256     def _find_contact_point(self, z_data, d_data, outqueue=None):
 257         fn = getattr(self, '_find_contact_point_%s'
 258                      % self.plugin.config['surface contact point algorithm'])
 259         return fn(z_data, d_data, outqueue)
 260
 261     def _find_contact_point_fmms(self, z_data, d_data, outqueue=None):
 262         """Algorithm by Francesco Musiani and Massimo Sandal.
 263
 264         Notes
 265         -----
 266         Algorithm:
 267
 268         * take care of the PicoForce trigger bug - exclude retraction portions with too high standard deviation
 269         * fit the second half of the retraction curve to a line
 270         * if the fit is not almost horizontal, take a smaller chunk and repeat
 271         * otherwise, we have something horizontal
 272         * so take the average of horizontal points and use it as a baseline
 273
 274         Then, start from the rise of the retraction curve and look at
 275         the first point below the baseline.
 276         """
 277         if not plot:
 278             plot=self.plots[0]
 279
 280         outplot=self.subtract_curves(1)
 281         xret=outplot.vectors[1][0]
 282         ydiff=outplot.vectors[1][1]
 283
 284         xext=plot.vectors[0][0]
 285         yext=plot.vectors[0][1]
 286         xret2=plot.vectors[1][0]
 287         yret=plot.vectors[1][1]
 288
 289         #taking care of the picoforce trigger bug: we exclude portions of the curve that have too much
 290         #standard deviation. yes, a lot of magic is here.
 291         monster=True
 292         monlength=len(xret)-int(len(xret)/20)
 293         finalength=len(xret)
 294         while monster:
 295             monchunk=scipy.array(ydiff[monlength:finalength])
 296             if abs(np.std(monchunk)) < 2e-10:
 297                 monster=False
 298             else: #move away from the monster
 299                 monlength-=int(len(xret)/50)
 300                 finalength-=int(len(xret)/50)
 301
 302
 303         #take half of the thing
 304         endlength=int(len(xret)/2)
 305
 306         ok=False
 307
 308         while not ok:
 309             xchunk=yext[endlength:monlength]
 310             ychunk=yext[endlength:monlength]
 311             regr=scipy.stats.linregress(xchunk,ychunk)[0:2]
 312             #we stop if we found an almost-horizontal fit or if we're going too short...
 313             #FIXME: 0.1 and 6 here are "magic numbers" (although reasonable)
 314             if (abs(regr[1]) > 0.1) and ( endlength < len(xret)-int(len(xret)/6) ) :
 315                 endlength+=10
 316             else:
 317                 ok=True
 318
 319
 320         ymean=np.mean(ychunk) #baseline
 321
 322         index=0
 323         point = ymean+1
 324
 325         #find the first point below the calculated baseline
 326         while point > ymean:
 327             try:
 328                 point=yret[index]
 329                 index+=1
 330             except IndexError:
 331                 #The algorithm didn't find anything below the baseline! It should NEVER happen
 332                 index=0
 333                 return index
 334
 335         return index
 336
 337     def _find_contact_point_ms(self, z_data, d_data, outqueue=None):
 338         """Algorithm by Massimo Sandal.
 339
 340         Notes
 341         -----
 342         WTK: At least the commits are by Massimo, and I see no notes
 343         attributing the algorithm to anyone else.
 344
 345         Algorithm:
 346
 347         * ?
 348         """
 349         #raw_plot=self.current.curve.default_plots()[0]
 350         raw_plot=self.plots[0]
 351         '''xext=self.plots[0].vectors[0][0]
 352         yext=self.plots[0].vectors[0][1]
 353         xret2=self.plots[0].vectors[1][0]
 354         yret=self.plots[0].vectors[1][1]
 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, 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