generalvclamp.py

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