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