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 do_subtplot(self, args):
 509         '''
 510         SUBTPLOT
 511         (procplots.py plugin)
 512         Plots the difference between ret and ext current curve
 513         -------
 514         Syntax: subtplot
 515         '''
 516         #FIXME: sub_filter and sub_order must be args
 517
 518         if len(self.plots[0].vectors) != 2:
 519             print 'This command only works on a curve with two different plots.'
 520             pass
 521
 522         outplot=self.subtract_curves(sub_order=1)
 523
 524         plot_graph=self.list_of_events['plot_graph']
 525         wx.PostEvent(self.frame,plot_graph(plots=[outplot]))
 526
 527     def _plug_init(self):
 528         self.basecurrent=None
 529         self.basepoints=None
 530         self.autofile=''
 531
 532     def do_distance(self,args):
 533         '''
 534         DISTANCE
 535         (generalvclamp.py)
 536         Measure the distance (in nm) between two points.
 537         For a standard experiment this is the delta X distance.
 538         For a force clamp experiment this is the delta Y distance (actually becomes
 539         an alias of zpiezo)
 540         -----------------
 541         Syntax: distance
 542         '''
 543         if self.current.curve.experiment == 'clamp':
 544             print 'You wanted to use zpiezo perhaps?'
 545             return
 546         else:
 547             dx,unitx,dy,unity=self._delta(set=1)
 548             print str(dx*(10**9))+' nm'
 549             to_dump='distance '+self.current.path+' '+str(dx*(10**9))+' nm'
 550             self.outlet.push(to_dump)
 551
 552
 553     def do_force(self,args):
 554         '''
 555         FORCE
 556         (generalvclamp.py)
 557         Measure the force difference (in pN) between two points
 558         ---------------
 559         Syntax: force
 560         '''
 561         if self.current.curve.experiment == 'clamp':
 562             print 'This command makes no sense for a force clamp experiment.'
 563             return
 564         dx,unitx,dy,unity=self._delta(set=1)
 565         print str(dy*(10**12))+' pN'
 566         to_dump='force '+self.current.path+' '+str(dy*(10**12))+' pN'
 567         self.outlet.push(to_dump)
 568
 569
 570     def do_forcebase(self,args):
 571         '''
 572         FORCEBASE
 573         (generalvclamp.py)
 574         Measures the difference in force (in pN) between a point and a baseline
 575         took as the average between two points.
 576
 577         The baseline is fixed once for a given curve and different force measurements,
 578         unless the user wants it to be recalculated
 579         ------------
 580         Syntax: forcebase [rebase]
 581                 rebase: Forces forcebase to ask again the baseline
 582                 max: Instead of asking for a point to measure, asks for two points and use
 583                      the maximum peak in between
 584         '''
 585         rebase=False #if true=we select rebase
 586         maxpoint=False #if true=we measure the maximum peak
 587
 588         plot=self._get_displayed_plot()
 589         whatset=1 #fixme: for all sets
 590         if 'rebase' in args or (self.basecurrent != self.current.path):
 591             rebase=True
 592         if 'max' in args:
 593             maxpoint=True
 594
 595         if rebase:
 596             print 'Select baseline'
 597             self.basepoints=self._measure_N_points(N=2, whatset=whatset)
 598             self.basecurrent=self.current.path
 599
 600         if maxpoint:
 601             print 'Select two points'
 602             points=self._measure_N_points(N=2, whatset=whatset)
 603             boundpoints=[points[0].index, points[1].index]
 604             boundpoints.sort()
 605             try:
 606                 y=min(plot.vectors[whatset][1][boundpoints[0]:boundpoints[1]])
 607             except ValueError:
 608                 print 'Chosen interval not valid. Try picking it again. Did you pick the same point as begin and end of interval?'
 609         else:
 610             print 'Select point to measure'
 611             points=self._measure_N_points(N=1, whatset=whatset)
 612             #whatplot=points[0].dest
 613             y=points[0].graph_coords[1]
 614
 615         #fixme: code duplication
 616         boundaries=[self.basepoints[0].index, self.basepoints[1].index]
 617         boundaries.sort()
 618         to_average=plot.vectors[whatset][1][boundaries[0]:boundaries[1]] #y points to average
 619
 620         avg=np.mean(to_average)
 621         forcebase=abs(y-avg)
 622         print str(forcebase*(10**12))+' pN'
 623         to_dump='forcebase '+self.current.path+' '+str(forcebase*(10**12))+' pN'
 624         self.outlet.push(to_dump)
 625
 626     def plotmanip_multiplier(self, plot, current):
 627         '''
 628         Multiplies all the Y values of an SMFS curve by a value stored in the 'force_multiplier'
 629         configuration variable. Useful for calibrations and other stuff.
 630         '''
 631
 632         #not a smfs curve...
 633         if current.curve.experiment != 'smfs':
 634             return plot
 635
 636         #only one set is present...
 637         if len(self.plots[0].vectors) != 2:
 638             return plot
 639
 640         #multiplier is 1...
 641         if (self.config['force_multiplier']==1):
 642             return plot
 643
 644         for i in range(len(plot.vectors[0][1])):
 645             plot.vectors[0][1][i]=plot.vectors[0][1][i]*self.config['force_multiplier']
 646
 647         for i in range(len(plot.vectors[1][1])):
 648             plot.vectors[1][1][i]=plot.vectors[1][1][i]*self.config['force_multiplier']
 649
 650         return plot
 651
 652
 653     def plotmanip_flatten(self, plot, current, customvalue=False):
 654         '''
 655         Subtracts a polynomial fit to the non-contact part of the curve, as to flatten it.
 656         the best polynomial fit is chosen among polynomials of degree 1 to n, where n is
 657         given by the configuration file or by the customvalue.
 658
 659         customvalue= int (>0) --> starts the function even if config says no (default=False)
 660         '''
 661
 662         #not a smfs curve...
 663         if current.curve.experiment != 'smfs':
 664             return plot
 665
 666         #only one set is present...
 667         if len(self.plots[0].vectors) != 2:
 668             return plot
 669
 670         #config is not flatten, and customvalue flag is false too
 671         if (not self.config['flatten']) and (not customvalue):
 672             return plot
 673
 674         max_exponent=12
 675         delta_contact=0
 676
 677         if customvalue:
 678             max_cycles=customvalue
 679         else:
 680             max_cycles=self.config['flatten'] #Using > 1 usually doesn't help and can give artefacts. However, it could be useful too.
 681
 682         contact_index=self.find_contact_point()
 683
 684         valn=[[] for item in range(max_exponent)]
 685         yrn=[0.0 for item in range(max_exponent)]
 686         errn=[0.0 for item in range(max_exponent)]
 687
 688         #Check if we have a proper numerical value
 689         try:
 690             zzz=int(max_cycles)
 691         except:
 692             #Loudly and annoyingly complain if it's not a number, then fallback to zero
 693             print '''Warning: flatten value is not a number!
 694             Use "set flatten" or edit hooke.conf to set it properly
 695             Using zero.'''
 696             max_cycles=0
 697
 698         for i in range(int(max_cycles)):
 699
 700             x_ext=plot.vectors[0][0][contact_index+delta_contact:]
 701             y_ext=plot.vectors[0][1][contact_index+delta_contact:]
 702             x_ret=plot.vectors[1][0][contact_index+delta_contact:]
 703             y_ret=plot.vectors[1][1][contact_index+delta_contact:]
 704             for exponent in range(max_exponent):
 705                 try:
 706                     valn[exponent]=sp.polyfit(x_ext,y_ext,exponent)
 707                     yrn[exponent]=sp.polyval(valn[exponent],x_ret)
 708                     errn[exponent]=sp.sqrt(sum((yrn[exponent]-y_ext)**2)/float(len(y_ext)))
 709                 except Exception,e:
 710                     print 'Cannot flatten!'
 711                     print e
 712                     return plot
 713
 714             best_exponent=errn.index(min(errn))
 715
 716             #extension
 717             ycorr_ext=y_ext-yrn[best_exponent]+y_ext[0] #noncontact part
 718             yjoin_ext=np.array(plot.vectors[0][1][0:contact_index+delta_contact]) #contact part
 719             #retraction
 720             ycorr_ret=y_ret-yrn[best_exponent]+y_ext[0] #noncontact part
 721             yjoin_ret=np.array(plot.vectors[1][1][0:contact_index+delta_contact]) #contact part
 722
 723             ycorr_ext=np.concatenate((yjoin_ext, ycorr_ext))
 724             ycorr_ret=np.concatenate((yjoin_ret, ycorr_ret))
 725
 726             plot.vectors[0][1]=list(ycorr_ext)
 727             plot.vectors[1][1]=list(ycorr_ret)
 728
 729         return plot
 730
 731     #---SLOPE---
 732     def do_slope(self,args):
 733         '''
 734         SLOPE
 735         (generalvclamp.py)
 736         Measures the slope of a delimited chunk on the return trace.
 737         The chunk can be delimited either by two manual clicks, or have
 738         a fixed width, given as an argument.
 739         ---------------
 740         Syntax: slope [width]
 741                 The facultative [width] parameter specifies how many
 742                 points will be considered for the fit. If [width] is
 743                 specified, only one click will be required.
 744         (c) Marco Brucale, Massimo Sandal 2008
 745         '''
 746
 747         # Reads the facultative width argument
 748         try:
 749             fitspan=int(args)
 750         except:
 751             fitspan=0
 752
 753         # Decides between the two forms of user input, as per (args)
 754         if fitspan == 0:
 755             # Gets the Xs of two clicked points as indexes on the current curve vector
 756             print 'Click twice to delimit chunk'
 757             points=self._measure_N_points(N=2,whatset=1)
 758         else:
 759             print 'Click once on the leftmost point of the chunk (i.e.usually the peak)'
 760             points=self._measure_N_points(N=1,whatset=1)
 761
 762         slope=self._slope(points,fitspan)
 763
 764         # Outputs the relevant slope parameter
 765         print 'Slope:'
 766         print str(slope)
 767         to_dump='slope '+self.current.path+' '+str(slope)
 768         self.outlet.push(to_dump)
 769
 770     def _slope(self,points,fitspan):
 771         # Calls the function linefit_between
 772         parameters=[0,0,[],[]]
 773         try:
 774             clickedpoints=[points[0].index,points[1].index]
 775             clickedpoints.sort()
 776         except:
 777             clickedpoints=[points[0].index-fitspan,points[0].index]
 778
 779         try:
 780             parameters=self.linefit_between(clickedpoints[0],clickedpoints[1])
 781         except:
 782             print 'Cannot fit. Did you click twice the same point?'
 783             return
 784
 785         # Outputs the relevant slope parameter
 786         print 'Slope:'
 787         print str(parameters[0])
 788         to_dump='slope '+self.curve.path+' '+str(parameters[0])
 789         self.outlet.push(to_dump)
 790
 791         # Makes a vector with the fitted parameters and sends it to the GUI
 792         xtoplot=parameters[2]
 793         ytoplot=[]
 794         x=0
 795         for x in xtoplot:
 796             ytoplot.append((x*parameters[0])+parameters[1])
 797
 798         clickvector_x, clickvector_y=[], []
 799         for item in points:
 800             clickvector_x.append(item.graph_coords[0])
 801             clickvector_y.append(item.graph_coords[1])
 802
 803         lineplot=self._get_displayed_plot(0) #get topmost displayed plot
 804
 805         lineplot.add_set(xtoplot,ytoplot)
 806         lineplot.add_set(clickvector_x, clickvector_y)
 807
 808
 809         if lineplot.styles==[]:
 810             lineplot.styles=[None,None,None,'scatter']
 811         else:
 812             lineplot.styles+=[None,'scatter']
 813         if lineplot.colors==[]:
 814             lineplot.colors=[None,None,'black',None]
 815         else:
 816             lineplot.colors+=['black',None]
 817
 818
 819         self._send_plot([lineplot])
 820
 821         return parameters[0]
 822
 823
 824     def linefit_between(self,index1,index2,whatset=1):
 825         '''
 826         Creates two vectors (xtofit,ytofit) slicing out from the
 827         current return trace a portion delimited by the two indexes
 828         given as arguments.
 829         Then does a least squares linear fit on that slice.
 830         Finally returns [0]=the slope, [1]=the intercept of the
 831         fitted 1st grade polynomial, and [2,3]=the actual (x,y) vectors
 832         used for the fit.
 833         (c) Marco Brucale, Massimo Sandal 2008
 834         '''
 835         # Translates the indexes into two vectors containing the x,y data to fit
 836         xtofit=self.plots[0].vectors[whatset][0][index1:index2]
 837         ytofit=self.plots[0].vectors[whatset][1][index1:index2]
 838
 839         # Does the actual linear fitting (simple least squares with numpy.polyfit)
 840         linefit=[]
 841         linefit=np.polyfit(xtofit,ytofit,1)
 842
 843         return (linefit[0],linefit[1],xtofit,ytofit)
 844
 845
 846     def fit_interval_nm(self,start_index,plot,nm,backwards):
 847           '''
 848           Calculates the number of points to fit, given a fit interval in nm
 849           start_index: index of point
 850           plot: plot to use
 851           backwards: if true, finds a point backwards.
 852           '''
 853           whatset=1 #FIXME: should be decidable
 854           x_vect=plot.vectors[1][0]
 855
 856           c=0
 857           i=start_index
 858           start=x_vect[start_index]
 859           maxlen=len(x_vect)
 860           while abs(x_vect[i]-x_vect[start_index])*(10**9) < nm:
 861               if i==0 or i==maxlen-1: #we reached boundaries of vector!
 862                   return c
 863
 864               if backwards:
 865                   i-=1
 866               else:
 867                   i+=1
 868               c+=1
 869           return c
 870
 871
 872
 873     def find_current_peaks(self,noflatten, a=True, maxpeak=True):
 874             #Find peaks.
 875             if a==True:
 876                   a=self.convfilt_config['mindeviation']
 877             try:
 878                   abs_devs=float(a)
 879             except:
 880                   print "Bad input, using default."
 881                   abs_devs=self.convfilt_config['mindeviation']
 882
 883             defplot=self.current.curve.default_plots()[0]
 884             if not noflatten:
 885                 flatten=self._find_plotmanip('flatten') #Extract flatten plotmanip
 886                 defplot=flatten(defplot, self.current, customvalue=1) #Flatten curve before feeding it to has_peaks
 887             pk_location,peak_size=self.has_peaks(defplot, abs_devs, maxpeak)
 888             return pk_location, peak_size
 889
 890
 891     def pickup_contact_point(self,N=1,whatset=1):
 892         '''macro to pick up the contact point by clicking'''
 893         contact_point=self._measure_N_points(N=1, whatset=1)[0]
 894         contact_point_index=contact_point.index
 895         self.wlccontact_point=contact_point
 896         self.wlccontact_index=contact_point.index
 897         self.wlccurrent=self.current.path
 898         return contact_point, contact_point_index
 899
 900
 901     def baseline_points(self,peak_location, displayed_plot):
 902         clicks=self.config['baseline_clicks']
 903         if clicks==0:
 904             self.basepoints=[]
 905             base_index_0=peak_location[-1]+self.fit_interval_nm(peak_location[-1], displayed_plot, self.config['auto_right_baseline'],False)
 906             self.basepoints.append(self._clickize(displayed_plot.vectors[1][0],displayed_plot.vectors[1][1],base_index_0))
 907             base_index_1=self.basepoints[0].index+self.fit_interval_nm(self.basepoints[0].index, displayed_plot, self.config['auto_left_baseline'],False)
 908             self.basepoints.append(self._clickize(displayed_plot.vectors[1][0],displayed_plot.vectors[1][1],base_index_1))
 909         elif clicks>0:
 910             print 'Select baseline'
 911             if clicks==1:
 912                 self.basepoints=self._measure_N_points(N=1, whatset=1)
 913                 base_index_1=self.basepoints[0].index+self.fit_interval_nm(self.basepoints[0].index, displayed_plot, self.config['auto_left_baseline'], False)
 914                 self.basepoints.append(self._clickize(displayed_plot.vectors[1][0],displayed_plot.vectors[1][1],base_index_1))
 915             else:
 916                 self.basepoints=self._measure_N_points(N=2, whatset=1)
 917
 918         self.basecurrent=self.current.path
 919         return self.basepoints