hooke/plugin/vclamp.py

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