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