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