hooke/plugin/vclamp.py

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