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