hooke/plugin/vclamp.py

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