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