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
  10 # modify it under the terms of the GNU Lesser General Public
  11 # License as published by the Free Software Foundation, either
  12 # version 3 of the License, or (at your option) any later version.
  13 #
  14 # Hooke is distributed in the hope that it will be useful,
  15 # but WITHOUT ANY WARRANTY; without even the implied warranty of
  16 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  17 # GNU Lesser General 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 from ..command import Command, Argument, Failure
  31 from ..curve import Data
  32 from ..plugin import Builtin
  33
  34
  35 # Utility functions
  36
  37 def scale(curve):
  38     for i,block in enumerate(curve.data):
  39         data = Data((block.shape[0], block.shape[1]+2), dtype=block.dtype)
  40         data.info = copy.deepcopy(block.info)
  41         data[:,:-2] = block
  42         data.info['columns'].extend(['surface z piezo (m)', 'deflection (N)'])
  43         z_data = data[:,data.info['columns'].index('z piezo (m)')]
  44         d_data = data[:,data.info['columns'].index('deflection (m)')]
  45         surface_offset = min(z_data) # TODO
  46         data.info['surface offset (m)'] = surface_offset
  47         data[:,-2] = z_data - surface_offset
  48         data[:,-1] = d_data * data.info['spring constant (N/m)']
  49         curve.data[i] = data
  50     return curve
  51
  52
  53 class VelocityClampPlugin (Builtin):
  54     def __init__(self):
  55         super(VelocityClampPlugin, self).__init__(name='vclamp')
  56         self._commands = []
  57 #            NextCommand(self), PreviousCommand(self), JumpCommand(self),
  58 #            ]
  59
  60 class generalvclampCommands(object):
  61
  62     def do_subtplot(self, args):
  63         '''
  64         SUBTPLOT
  65         (procplots.py plugin)
  66         Plots the difference between ret and ext current curve
  67         -------
  68         Syntax: subtplot
  69         '''
  70         #FIXME: sub_filter and sub_order must be args
  71
  72         if len(self.plots[0].vectors) != 2:
  73             print 'This command only works on a curve with two different plots.'
  74             pass
  75
  76         outplot=self.subtract_curves(sub_order=1)
  77
  78         plot_graph=self.list_of_events['plot_graph']
  79         wx.PostEvent(self.frame,plot_graph(plots=[outplot]))
  80
  81     def _plug_init(self):
  82         self.basecurrent=None
  83         self.basepoints=None
  84         self.autofile=''
  85
  86     def do_distance(self,args):
  87         '''
  88         DISTANCE
  89         (generalvclamp.py)
  90         Measure the distance (in nm) between two points.
  91         For a standard experiment this is the delta X distance.
  92         For a force clamp experiment this is the delta Y distance (actually becomes
  93         an alias of zpiezo)
  94         -----------------
  95         Syntax: distance
  96         '''
  97         if self.current.curve.experiment == 'clamp':
  98             print 'You wanted to use zpiezo perhaps?'
  99             return
 100         else:
 101             dx,unitx,dy,unity=self._delta(set=1)
 102             print str(dx*(10**9))+' nm'
 103             to_dump='distance '+self.current.path+' '+str(dx*(10**9))+' nm'
 104             self.outlet.push(to_dump)
 105
 106
 107     def do_force(self,args):
 108         '''
 109         FORCE
 110         (generalvclamp.py)
 111         Measure the force difference (in pN) between two points
 112         ---------------
 113         Syntax: force
 114         '''
 115         if self.current.curve.experiment == 'clamp':
 116             print 'This command makes no sense for a force clamp experiment.'
 117             return
 118         dx,unitx,dy,unity=self._delta(set=1)
 119         print str(dy*(10**12))+' pN'
 120         to_dump='force '+self.current.path+' '+str(dy*(10**12))+' pN'
 121         self.outlet.push(to_dump)
 122
 123
 124     def do_forcebase(self,args):
 125         '''
 126         FORCEBASE
 127         (generalvclamp.py)
 128         Measures the difference in force (in pN) between a point and a baseline
 129         took as the average between two points.
 130
 131         The baseline is fixed once for a given curve and different force measurements,
 132         unless the user wants it to be recalculated
 133         ------------
 134         Syntax: forcebase [rebase]
 135                 rebase: Forces forcebase to ask again the baseline
 136                 max: Instead of asking for a point to measure, asks for two points and use
 137                      the maximum peak in between
 138         '''
 139         rebase=False #if true=we select rebase
 140         maxpoint=False #if true=we measure the maximum peak
 141
 142         plot=self._get_displayed_plot()
 143         whatset=1 #fixme: for all sets
 144         if 'rebase' in args or (self.basecurrent != self.current.path):
 145             rebase=True
 146         if 'max' in args:
 147             maxpoint=True
 148
 149         if rebase:
 150             print 'Select baseline'
 151             self.basepoints=self._measure_N_points(N=2, whatset=whatset)
 152             self.basecurrent=self.current.path
 153
 154         if maxpoint:
 155             print 'Select two points'
 156             points=self._measure_N_points(N=2, whatset=whatset)
 157             boundpoints=[points[0].index, points[1].index]
 158             boundpoints.sort()
 159             try:
 160                 y=min(plot.vectors[whatset][1][boundpoints[0]:boundpoints[1]])
 161             except ValueError:
 162                 print 'Chosen interval not valid. Try picking it again. Did you pick the same point as begin and end of interval?'
 163         else:
 164             print 'Select point to measure'
 165             points=self._measure_N_points(N=1, whatset=whatset)
 166             #whatplot=points[0].dest
 167             y=points[0].graph_coords[1]
 168
 169         #fixme: code duplication
 170         boundaries=[self.basepoints[0].index, self.basepoints[1].index]
 171         boundaries.sort()
 172         to_average=plot.vectors[whatset][1][boundaries[0]:boundaries[1]] #y points to average
 173
 174         avg=np.mean(to_average)
 175         forcebase=abs(y-avg)
 176         print str(forcebase*(10**12))+' pN'
 177         to_dump='forcebase '+self.current.path+' '+str(forcebase*(10**12))+' pN'
 178         self.outlet.push(to_dump)
 179
 180     def plotmanip_multiplier(self, plot, current):
 181         '''
 182         Multiplies all the Y values of an SMFS curve by a value stored in the 'force_multiplier'
 183         configuration variable. Useful for calibrations and other stuff.
 184         '''
 185
 186         #not a smfs curve...
 187         if current.curve.experiment != 'smfs':
 188             return plot
 189
 190         #only one set is present...
 191         if len(self.plots[0].vectors) != 2:
 192             return plot
 193
 194         #multiplier is 1...
 195         if (self.config['force_multiplier']==1):
 196             return plot
 197
 198         for i in range(len(plot.vectors[0][1])):
 199             plot.vectors[0][1][i]=plot.vectors[0][1][i]*self.config['force_multiplier']
 200
 201         for i in range(len(plot.vectors[1][1])):
 202             plot.vectors[1][1][i]=plot.vectors[1][1][i]*self.config['force_multiplier']
 203
 204         return plot
 205
 206
 207     def plotmanip_flatten(self, plot, current, customvalue=False):
 208         '''
 209         Subtracts a polynomial fit to the non-contact part of the curve, as to flatten it.
 210         the best polynomial fit is chosen among polynomials of degree 1 to n, where n is
 211         given by the configuration file or by the customvalue.
 212
 213         customvalue= int (>0) --> starts the function even if config says no (default=False)
 214         '''
 215
 216         #not a smfs curve...
 217         if current.curve.experiment != 'smfs':
 218             return plot
 219
 220         #only one set is present...
 221         if len(self.plots[0].vectors) != 2:
 222             return plot
 223
 224         #config is not flatten, and customvalue flag is false too
 225         if (not self.config['flatten']) and (not customvalue):
 226             return plot
 227
 228         max_exponent=12
 229         delta_contact=0
 230
 231         if customvalue:
 232             max_cycles=customvalue
 233         else:
 234             max_cycles=self.config['flatten'] #Using > 1 usually doesn't help and can give artefacts. However, it could be useful too.
 235
 236         contact_index=self.find_contact_point()
 237
 238         valn=[[] for item in range(max_exponent)]
 239         yrn=[0.0 for item in range(max_exponent)]
 240         errn=[0.0 for item in range(max_exponent)]
 241
 242         #Check if we have a proper numerical value
 243         try:
 244             zzz=int(max_cycles)
 245         except:
 246             #Loudly and annoyingly complain if it's not a number, then fallback to zero
 247             print '''Warning: flatten value is not a number!
 248             Use "set flatten" or edit hooke.conf to set it properly
 249             Using zero.'''
 250             max_cycles=0
 251
 252         for i in range(int(max_cycles)):
 253
 254             x_ext=plot.vectors[0][0][contact_index+delta_contact:]
 255             y_ext=plot.vectors[0][1][contact_index+delta_contact:]
 256             x_ret=plot.vectors[1][0][contact_index+delta_contact:]
 257             y_ret=plot.vectors[1][1][contact_index+delta_contact:]
 258             for exponent in range(max_exponent):
 259                 try:
 260                     valn[exponent]=sp.polyfit(x_ext,y_ext,exponent)
 261                     yrn[exponent]=sp.polyval(valn[exponent],x_ret)
 262                     errn[exponent]=sp.sqrt(sum((yrn[exponent]-y_ext)**2)/float(len(y_ext)))
 263                 except Exception,e:
 264                     print 'Cannot flatten!'
 265                     print e
 266                     return plot
 267
 268             best_exponent=errn.index(min(errn))
 269
 270             #extension
 271             ycorr_ext=y_ext-yrn[best_exponent]+y_ext[0] #noncontact part
 272             yjoin_ext=np.array(plot.vectors[0][1][0:contact_index+delta_contact]) #contact part
 273             #retraction
 274             ycorr_ret=y_ret-yrn[best_exponent]+y_ext[0] #noncontact part
 275             yjoin_ret=np.array(plot.vectors[1][1][0:contact_index+delta_contact]) #contact part
 276
 277             ycorr_ext=np.concatenate((yjoin_ext, ycorr_ext))
 278             ycorr_ret=np.concatenate((yjoin_ret, ycorr_ret))
 279
 280             plot.vectors[0][1]=list(ycorr_ext)
 281             plot.vectors[1][1]=list(ycorr_ret)
 282
 283         return plot
 284
 285     #---SLOPE---
 286     def do_slope(self,args):
 287         '''
 288         SLOPE
 289         (generalvclamp.py)
 290         Measures the slope of a delimited chunk on the return trace.
 291         The chunk can be delimited either by two manual clicks, or have
 292         a fixed width, given as an argument.
 293         ---------------
 294         Syntax: slope [width]
 295                 The facultative [width] parameter specifies how many
 296                 points will be considered for the fit. If [width] is
 297                 specified, only one click will be required.
 298         (c) Marco Brucale, Massimo Sandal 2008
 299         '''
 300
 301         # Reads the facultative width argument
 302         try:
 303             fitspan=int(args)
 304         except:
 305             fitspan=0
 306
 307         # Decides between the two forms of user input, as per (args)
 308         if fitspan == 0:
 309             # Gets the Xs of two clicked points as indexes on the current curve vector
 310             print 'Click twice to delimit chunk'
 311             points=self._measure_N_points(N=2,whatset=1)
 312         else:
 313             print 'Click once on the leftmost point of the chunk (i.e.usually the peak)'
 314             points=self._measure_N_points(N=1,whatset=1)
 315
 316         slope=self._slope(points,fitspan)
 317
 318         # Outputs the relevant slope parameter
 319         print 'Slope:'
 320         print str(slope)
 321         to_dump='slope '+self.current.path+' '+str(slope)
 322         self.outlet.push(to_dump)
 323
 324     def _slope(self,points,fitspan):
 325         # Calls the function linefit_between
 326         parameters=[0,0,[],[]]
 327         try:
 328             clickedpoints=[points[0].index,points[1].index]
 329             clickedpoints.sort()
 330         except:
 331             clickedpoints=[points[0].index-fitspan,points[0].index]
 332
 333         try:
 334             parameters=self.linefit_between(clickedpoints[0],clickedpoints[1])
 335         except:
 336             print 'Cannot fit. Did you click twice the same point?'
 337             return
 338
 339         # Outputs the relevant slope parameter
 340         print 'Slope:'
 341         print str(parameters[0])
 342         to_dump='slope '+self.curve.path+' '+str(parameters[0])
 343         self.outlet.push(to_dump)
 344
 345         # Makes a vector with the fitted parameters and sends it to the GUI
 346         xtoplot=parameters[2]
 347         ytoplot=[]
 348         x=0
 349         for x in xtoplot:
 350             ytoplot.append((x*parameters[0])+parameters[1])
 351
 352         clickvector_x, clickvector_y=[], []
 353         for item in points:
 354             clickvector_x.append(item.graph_coords[0])
 355             clickvector_y.append(item.graph_coords[1])
 356
 357         lineplot=self._get_displayed_plot(0) #get topmost displayed plot
 358
 359         lineplot.add_set(xtoplot,ytoplot)
 360         lineplot.add_set(clickvector_x, clickvector_y)
 361
 362
 363         if lineplot.styles==[]:
 364             lineplot.styles=[None,None,None,'scatter']
 365         else:
 366             lineplot.styles+=[None,'scatter']
 367         if lineplot.colors==[]:
 368             lineplot.colors=[None,None,'black',None]
 369         else:
 370             lineplot.colors+=['black',None]
 371
 372
 373         self._send_plot([lineplot])
 374
 375         return parameters[0]
 376
 377
 378     def linefit_between(self,index1,index2,whatset=1):
 379         '''
 380         Creates two vectors (xtofit,ytofit) slicing out from the
 381         current return trace a portion delimited by the two indexes
 382         given as arguments.
 383         Then does a least squares linear fit on that slice.
 384         Finally returns [0]=the slope, [1]=the intercept of the
 385         fitted 1st grade polynomial, and [2,3]=the actual (x,y) vectors
 386         used for the fit.
 387         (c) Marco Brucale, Massimo Sandal 2008
 388         '''
 389         # Translates the indexes into two vectors containing the x,y data to fit
 390         xtofit=self.plots[0].vectors[whatset][0][index1:index2]
 391         ytofit=self.plots[0].vectors[whatset][1][index1:index2]
 392
 393         # Does the actual linear fitting (simple least squares with numpy.polyfit)
 394         linefit=[]
 395         linefit=np.polyfit(xtofit,ytofit,1)
 396
 397         return (linefit[0],linefit[1],xtofit,ytofit)
 398
 399
 400     def fit_interval_nm(self,start_index,plot,nm,backwards):
 401           '''
 402           Calculates the number of points to fit, given a fit interval in nm
 403           start_index: index of point
 404           plot: plot to use
 405           backwards: if true, finds a point backwards.
 406           '''
 407           whatset=1 #FIXME: should be decidable
 408           x_vect=plot.vectors[1][0]
 409
 410           c=0
 411           i=start_index
 412           start=x_vect[start_index]
 413           maxlen=len(x_vect)
 414           while abs(x_vect[i]-x_vect[start_index])*(10**9) < nm:
 415               if i==0 or i==maxlen-1: #we reached boundaries of vector!
 416                   return c
 417
 418               if backwards:
 419                   i-=1
 420               else:
 421                   i+=1
 422               c+=1
 423           return c
 424
 425
 426
 427     def find_current_peaks(self,noflatten, a=True, maxpeak=True):
 428             #Find peaks.
 429             if a==True:
 430                   a=self.convfilt_config['mindeviation']
 431             try:
 432                   abs_devs=float(a)
 433             except:
 434                   print "Bad input, using default."
 435                   abs_devs=self.convfilt_config['mindeviation']
 436
 437             defplot=self.current.curve.default_plots()[0]
 438             if not noflatten:
 439                 flatten=self._find_plotmanip('flatten') #Extract flatten plotmanip
 440                 defplot=flatten(defplot, self.current, customvalue=1) #Flatten curve before feeding it to has_peaks
 441             pk_location,peak_size=self.has_peaks(defplot, abs_devs, maxpeak)
 442             return pk_location, peak_size
 443
 444
 445     def pickup_contact_point(self,N=1,whatset=1):
 446         '''macro to pick up the contact point by clicking'''
 447         contact_point=self._measure_N_points(N=1, whatset=1)[0]
 448         contact_point_index=contact_point.index
 449         self.wlccontact_point=contact_point
 450         self.wlccontact_index=contact_point.index
 451         self.wlccurrent=self.current.path
 452         return contact_point, contact_point_index
 453
 454
 455
 456     def baseline_points(self,peak_location, displayed_plot):
 457         clicks=self.config['baseline_clicks']
 458         if clicks==0:
 459             self.basepoints=[]
 460             base_index_0=peak_location[-1]+self.fit_interval_nm(peak_location[-1], displayed_plot, self.config['auto_right_baseline'],False)
 461             self.basepoints.append(self._clickize(displayed_plot.vectors[1][0],displayed_plot.vectors[1][1],base_index_0))
 462             base_index_1=self.basepoints[0].index+self.fit_interval_nm(self.basepoints[0].index, displayed_plot, self.config['auto_left_baseline'],False)
 463             self.basepoints.append(self._clickize(displayed_plot.vectors[1][0],displayed_plot.vectors[1][1],base_index_1))
 464         elif clicks>0:
 465             print 'Select baseline'
 466             if clicks==1:
 467                 self.basepoints=self._measure_N_points(N=1, whatset=1)
 468                 base_index_1=self.basepoints[0].index+self.fit_interval_nm(self.basepoints[0].index, displayed_plot, self.config['auto_left_baseline'], False)
 469                 self.basepoints.append(self._clickize(displayed_plot.vectors[1][0],displayed_plot.vectors[1][1],base_index_1))
 470             else:
 471                 self.basepoints=self._measure_N_points(N=2, whatset=1)
 472
 473         self.basecurrent=self.current.path
 474         return self.basepoints