hooke/plugin/generalvclamp.py

   1 '''
   2 generalvclamp.py
   3
   4 Plugin regarding general velocity clamp measurements
   5 '''
   6
   7 from libhooke import WX_GOOD, ClickedPoint
   8 import wxversion
   9 wxversion.select(WX_GOOD)
  10 from wx import PostEvent
  11 import numpy as np
  12 import scipy as sp
  13 import copy
  14 import os.path
  15 import time
  16
  17 import warnings
  18 warnings.simplefilter('ignore',np.RankWarning)
  19
  20
  21 class generalvclampCommands:
  22
  23     def _plug_init(self):
  24         self.basecurrent=None
  25         self.basepoints=None
  26         self.autofile=''
  27
  28     def do_distance(self,args):
  29         '''
  30         DISTANCE
  31         (generalvclamp.py)
  32         Measure the distance (in nm) between two points.
  33         For a standard experiment this is the delta X distance.
  34         For a force clamp experiment this is the delta Y distance (actually becomes
  35         an alias of zpiezo)
  36         -----------------
  37         Syntax: distance
  38         '''
  39         if self.current.curve.experiment == 'clamp':
  40             print 'You wanted to use zpiezo perhaps?'
  41             return
  42         else:
  43             dx,unitx,dy,unity=self._delta(set=1)
  44             print str(dx*(10**9))+' nm'
  45             to_dump='distance '+self.current.path+' '+str(dx*(10**9))+' nm'
  46             self.outlet.push(to_dump)
  47
  48
  49     def do_force(self,args):
  50         '''
  51         FORCE
  52         (generalvclamp.py)
  53         Measure the force difference (in pN) between two points
  54         ---------------
  55         Syntax: force
  56         '''
  57         if self.current.curve.experiment == 'clamp':
  58             print 'This command makes no sense for a force clamp experiment.'
  59             return
  60         dx,unitx,dy,unity=self._delta(set=1)
  61         print str(dy*(10**12))+' pN'
  62         to_dump='force '+self.current.path+' '+str(dy*(10**12))+' pN'
  63         self.outlet.push(to_dump)
  64
  65
  66     def do_forcebase(self,args):
  67         '''
  68         FORCEBASE
  69         (generalvclamp.py)
  70         Measures the difference in force (in pN) between a point and a baseline
  71         took as the average between two points.
  72
  73         The baseline is fixed once for a given curve and different force measurements,
  74         unless the user wants it to be recalculated
  75         ------------
  76         Syntax: forcebase [rebase]
  77                 rebase: Forces forcebase to ask again the baseline
  78                 max: Instead of asking for a point to measure, asks for two points and use
  79                      the maximum peak in between
  80         '''
  81         rebase=False #if true=we select rebase
  82         maxpoint=False #if true=we measure the maximum peak
  83
  84         plot=self._get_displayed_plot()
  85         whatset=1 #fixme: for all sets
  86         if 'rebase' in args or (self.basecurrent != self.current.path):
  87             rebase=True
  88         if 'max' in args:
  89             maxpoint=True
  90
  91         if rebase:
  92             print 'Select baseline'
  93             self.basepoints=self._measure_N_points(N=2, whatset=whatset)
  94             self.basecurrent=self.current.path
  95
  96         if maxpoint:
  97             print 'Select two points'
  98             points=self._measure_N_points(N=2, whatset=whatset)
  99             boundpoints=[points[0].index, points[1].index]
 100             boundpoints.sort()
 101             try:
 102                 y=min(plot.vectors[whatset][1][boundpoints[0]:boundpoints[1]])
 103             except ValueError:
 104                 print 'Chosen interval not valid. Try picking it again. Did you pick the same point as begin and end of interval?'
 105         else:
 106             print 'Select point to measure'
 107             points=self._measure_N_points(N=1, whatset=whatset)
 108             #whatplot=points[0].dest
 109             y=points[0].graph_coords[1]
 110
 111         #fixme: code duplication
 112         boundaries=[self.basepoints[0].index, self.basepoints[1].index]
 113         boundaries.sort()
 114         to_average=plot.vectors[whatset][1][boundaries[0]:boundaries[1]] #y points to average
 115
 116         avg=np.mean(to_average)
 117         forcebase=abs(y-avg)
 118         print str(forcebase*(10**12))+' pN'
 119         to_dump='forcebase '+self.current.path+' '+str(forcebase*(10**12))+' pN'
 120         self.outlet.push(to_dump)
 121 \r
 122     def plotmanip_multiplier(self, plot, current):
 123         '''
 124         Multiplies all the Y values of an SMFS curve by a value stored in the 'force_multiplier'\r
 125         configuration variable. Useful for calibrations and other stuff.
 126         '''
 127
 128         #not a smfs curve...
 129         if current.curve.experiment != 'smfs':
 130             return plot
 131
 132         #only one set is present...
 133         if len(self.plots[0].vectors) != 2:
 134             return plot
 135
 136         #multiplier is 1...
 137         if (self.config['force_multiplier']==1):
 138             return plot\r
 139 \r
 140         for i in range(len(plot.vectors[0][1])):
 141             plot.vectors[0][1][i]=plot.vectors[0][1][i]*self.config['force_multiplier']        \r
 142 \r
 143         for i in range(len(plot.vectors[1][1])):
 144             plot.vectors[1][1][i]=plot.vectors[1][1][i]*self.config['force_multiplier']\r
 145 \r
 146         return plot            \r
 147
 148
 149     def plotmanip_flatten(self, plot, current, customvalue=False):
 150         '''
 151         Subtracts a polynomial fit to the non-contact part of the curve, as to flatten it.
 152         the best polynomial fit is chosen among polynomials of degree 1 to n, where n is
 153         given by the configuration file or by the customvalue.
 154
 155         customvalue= int (>0) --> starts the function even if config says no (default=False)
 156         '''
 157
 158         #not a smfs curve...
 159         if current.curve.experiment != 'smfs':
 160             return plot
 161
 162         #only one set is present...
 163         if len(self.plots[0].vectors) != 2:
 164             return plot
 165
 166         #config is not flatten, and customvalue flag is false too
 167         if (not self.config['flatten']) and (not customvalue):
 168             return plot
 169
 170         max_exponent=12
 171         delta_contact=0
 172
 173         if customvalue:
 174             max_cycles=customvalue
 175         else:
 176             max_cycles=self.config['flatten'] #Using > 1 usually doesn't help and can give artefacts. However, it could be useful too.
 177
 178         contact_index=self.find_contact_point()
 179
 180         valn=[[] for item in range(max_exponent)]
 181         yrn=[0.0 for item in range(max_exponent)]
 182         errn=[0.0 for item in range(max_exponent)]
 183
 184         for i in range(int(max_cycles)):
 185
 186             x_ext=plot.vectors[0][0][contact_index+delta_contact:]
 187             y_ext=plot.vectors[0][1][contact_index+delta_contact:]
 188             x_ret=plot.vectors[1][0][contact_index+delta_contact:]
 189             y_ret=plot.vectors[1][1][contact_index+delta_contact:]
 190             for exponent in range(max_exponent):
 191                 try:
 192                     valn[exponent]=sp.polyfit(x_ext,y_ext,exponent)
 193                     yrn[exponent]=sp.polyval(valn[exponent],x_ret)
 194                     errn[exponent]=sp.sqrt(sum((yrn[exponent]-y_ext)**2)/float(len(y_ext)))
 195                 except Exception,e:
 196                     print 'Cannot flatten!'
 197                     print e
 198                     return plot
 199
 200             best_exponent=errn.index(min(errn))
 201
 202             #extension
 203             ycorr_ext=y_ext-yrn[best_exponent]+y_ext[0] #noncontact part
 204             yjoin_ext=np.array(plot.vectors[0][1][0:contact_index+delta_contact]) #contact part
 205             #retraction
 206             ycorr_ret=y_ret-yrn[best_exponent]+y_ext[0] #noncontact part
 207             yjoin_ret=np.array(plot.vectors[1][1][0:contact_index+delta_contact]) #contact part
 208
 209             ycorr_ext=np.concatenate((yjoin_ext, ycorr_ext))
 210             ycorr_ret=np.concatenate((yjoin_ret, ycorr_ret))
 211
 212             plot.vectors[0][1]=list(ycorr_ext)
 213             plot.vectors[1][1]=list(ycorr_ret)
 214
 215         return plot
 216
 217     #---SLOPE---
 218     def do_slope(self,args):
 219         '''
 220         SLOPE
 221         (generalvclamp.py)
 222         Measures the slope of a delimited chunk on the return trace.
 223         The chunk can be delimited either by two manual clicks, or have
 224         a fixed width, given as an argument.
 225         ---------------
 226         Syntax: slope [width]
 227                 The facultative [width] parameter specifies how many
 228                 points will be considered for the fit. If [width] is
 229                 specified, only one click will be required.
 230         (c) Marco Brucale, Massimo Sandal 2008
 231         '''
 232
 233         # Reads the facultative width argument
 234         try:
 235             fitspan=int(args)
 236         except:
 237             fitspan=0
 238
 239         # Decides between the two forms of user input, as per (args)
 240         if fitspan == 0:
 241             # Gets the Xs of two clicked points as indexes on the current curve vector
 242             print 'Click twice to delimit chunk'
 243             clickedpoints=[]
 244             points=self._measure_N_points(N=2,whatset=1)
 245             clickedpoints=[points[0].index,points[1].index]
 246             clickedpoints.sort()
 247         else:
 248             print 'Click once on the leftmost point of the chunk (i.e.usually the peak)'
 249             clickedpoints=[]
 250             points=self._measure_N_points(N=1,whatset=1)
 251             clickedpoints=[points[0].index-fitspan,points[0].index]
 252
 253         # Calls the function linefit_between
 254         parameters=[0,0,[],[]]
 255         try:
 256             parameters=self.linefit_between(clickedpoints[0],clickedpoints[1])
 257         except:
 258             print 'Cannot fit. Did you click twice the same point?'
 259             return
 260
 261         # Outputs the relevant slope parameter
 262         print 'Slope:'
 263         print str(parameters[0])
 264         to_dump='slope '+self.current.path+' '+str(parameters[0])
 265         self.outlet.push(to_dump)
 266
 267         # Makes a vector with the fitted parameters and sends it to the GUI
 268         xtoplot=parameters[2]
 269         ytoplot=[]
 270         x=0
 271         for x in xtoplot:
 272             ytoplot.append((x*parameters[0])+parameters[1])
 273
 274         clickvector_x, clickvector_y=[], []
 275         for item in points:
 276             clickvector_x.append(item.graph_coords[0])
 277             clickvector_y.append(item.graph_coords[1])
 278
 279         lineplot=self._get_displayed_plot(0) #get topmost displayed plot
 280
 281         lineplot.add_set(xtoplot,ytoplot)
 282         lineplot.add_set(clickvector_x, clickvector_y)
 283
 284
 285         if lineplot.styles==[]:
 286             lineplot.styles=[None,None,None,'scatter']
 287         else:
 288             lineplot.styles+=[None,'scatter']
 289         if lineplot.colors==[]:
 290             lineplot.styles=[None,None,None,None]
 291         else:
 292             lineplot.colors+=[None,None]
 293
 294
 295         self._send_plot([lineplot])
 296
 297     def linefit_between(self,index1,index2,whatset=1):
 298         '''
 299         Creates two vectors (xtofit,ytofit) slicing out from the
 300         current return trace a portion delimited by the two indexes
 301         given as arguments.
 302         Then does a least squares linear fit on that slice.
 303         Finally returns [0]=the slope, [1]=the intercept of the
 304         fitted 1st grade polynomial, and [2,3]=the actual (x,y) vectors
 305         used for the fit.
 306         (c) Marco Brucale, Massimo Sandal 2008
 307         '''
 308         # Translates the indexes into two vectors containing the x,y data to fit
 309         xtofit=self.plots[0].vectors[whatset][0][index1:index2]
 310         ytofit=self.plots[0].vectors[whatset][1][index1:index2]
 311
 312         # Does the actual linear fitting (simple least squares with numpy.polyfit)
 313         linefit=[]
 314         linefit=np.polyfit(xtofit,ytofit,1)
 315
 316         return (linefit[0],linefit[1],xtofit,ytofit)