plugins/generalvclamp.py

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