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