generalvclamp.py

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