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