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