hooke/plugin/generalvclamp.py

   1 #!/usr/bin/env python
   2
   3 '''
   4 generalvclamp.py
   5
   6 Plugin regarding general velocity clamp measurements
   7
   8 Copyright 2008 by Massimo Sandal, Fabrizio Benedetti, Marco Brucale, Bruno Samori (University of Bologna, Italy),
   9 and Alberto Gomez-Casado (University of Twente)
  10 with modifications by Dr. Rolf Schmidt (Concordia University, Canada)
  11
  12 This program is released under the GNU General Public License version 2.
  13 '''
  14
  15 import lib.libhooke as lh
  16 import wxversion
  17 wxversion.select(lh.WX_GOOD)
  18
  19 import numpy as np
  20 import scipy as sp
  21
  22 import warnings
  23 warnings.simplefilter('ignore', np.RankWarning)
  24
  25 import lib.curve
  26 import lib.prettyformat
  27
  28 class generalvclampCommands:
  29
  30     def _plug_init(self):
  31         self.basecurrent = ''
  32         self.basepoints = []
  33         #TODO: what is self.autofile for?
  34         #self.autofile = ''
  35
  36     def do_distance(self):
  37         '''
  38         DISTANCE
  39         (generalvclamp.py)
  40         Measure the distance (in nm) between two points.
  41         For a standard experiment this is the delta X distance.
  42         For a force clamp experiment this is the delta Y distance (actually becomes
  43         an alias of zpiezo)
  44         -----------------
  45         Syntax: distance
  46         '''
  47         color = self.GetColorFromConfig('generalvclamp', 'distance', 'color')
  48         decimals = self.GetIntFromConfig('generalvclamp', 'distance', 'decimals')
  49         prefix = self.GetStringFromConfig('generalvclamp', 'distance', 'prefix')
  50         multiplier = 10 ** lib.prettyformat.get_exponent(prefix)
  51         show =  self.GetBoolFromConfig('generalvclamp', 'distance', 'show')
  52         show_in_legend = self.GetBoolFromConfig('generalvclamp', 'distance', 'show_in_legend')
  53         size = self.GetIntFromConfig('generalvclamp', 'distance', 'size')
  54         whatset_str = self.GetStringFromConfig('generalvclamp', 'distance', 'whatset')
  55         whatset = 'retraction'
  56         if whatset_str == 'extension':
  57             whatset = lh.EXTENSION
  58         if whatset_str == 'retraction':
  59             whatset = lh.RETRACTION
  60
  61         active_file = self.GetActiveFile()
  62         if active_file.driver.experiment == 'clamp':
  63             self.AppendToOutput('You wanted to use zpiezo perhaps?')
  64             return
  65         plugin = lib.plugin.Plugin()
  66         plugin.name = 'generalvclamp'
  67         plugin.section = 'distance'
  68         delta = self._delta(message='Click 2 points to measure the distance.', whatset=whatset)
  69
  70         plot = self.GetDisplayedPlotCorrected()
  71         if show:
  72             #add the points to the plot
  73             points = lib.curve.Curve()
  74             points.color = color
  75             if show_in_legend:
  76                 points.label = 'distance'
  77             else:
  78                 points.label = '_nolegend_'
  79             points.size = size
  80             points.style = 'scatter'
  81             points.units.x = delta.units.x
  82             points.units.y = delta.units.y
  83             points.x = [delta.point1.x, delta.point2.x]
  84             points.y = [delta.point1.y, delta.point2.y]
  85             plot.curves.append(points)
  86
  87         self.UpdatePlot(plot)
  88
  89         output_str = lib.prettyformat.pretty_format(abs(delta.get_delta_x()), delta.units.x, decimals, multiplier)
  90         self.AppendToOutput(''.join(['Distance: ', output_str]))
  91
  92     def do_force(self):
  93         '''
  94         FORCE
  95         (generalvclamp.py)
  96         Measure the force difference (in pN) between two points
  97         ---------------
  98         Syntax: force
  99         '''
 100         color = self.GetColorFromConfig('generalvclamp', 'force', 'color')
 101         decimals = self.GetIntFromConfig('generalvclamp', 'force', 'decimals')
 102         prefix = self.GetStringFromConfig('generalvclamp', 'force', 'prefix')
 103         multiplier = 10 ** lib.prettyformat.get_exponent(prefix)
 104         show = self.GetBoolFromConfig('generalvclamp', 'force', 'show')
 105         show_in_legend = self.GetBoolFromConfig('generalvclamp', 'force', 'show_in_legend')
 106         size = self.GetIntFromConfig('generalvclamp', 'force', 'size')
 107         whatset_str = self.GetStringFromConfig('generalvclamp', 'force', 'whatset')
 108         whatset = 'retraction'
 109         if whatset_str == 'extension':
 110             whatset = lh.EXTENSION
 111         if whatset_str == 'retraction':
 112             whatset = lh.RETRACTION
 113
 114         active_file = self.GetActiveFile()
 115         if active_file.driver.experiment == 'clamp':
 116             self.AppendToOutput('This command makes no sense for a force clamp experiment.')
 117             return
 118         plugin = lib.plugin.Plugin()
 119         plugin.name = 'generalvclamp'
 120         plugin.section = 'force'
 121         delta = self._delta(message='Click 2 points to measure the force.', whatset=whatset)
 122
 123         plot = self.GetDisplayedPlotCorrected()
 124         if show:
 125             #add the points to the plot
 126             points = lib.curve.Curve()
 127             points.color = color
 128             if show_in_legend:
 129                 points.label = 'force'
 130             else:
 131                 points.label = '_nolegend_'
 132             points.size = size
 133             points.style = 'scatter'
 134             points.units.x = delta.units.x
 135             points.units.y = delta.units.y
 136             points.x = [delta.point1.x, delta.point2.x]
 137             points.y = [delta.point1.y, delta.point2.y]
 138             plot.curves.append(points)
 139
 140         self.UpdatePlot(plot)
 141
 142         output_str = lib.prettyformat.pretty_format(abs(delta.get_delta_y()), delta.units.y, decimals, multiplier)
 143         self.AppendToOutput(''.join(['Force: ', output_str]))
 144
 145     def do_forcebase(self):
 146         '''
 147         FORCEBASE
 148         (generalvclamp.py)
 149         Measures the difference in force (in pN) between a point and a baseline
 150         taken as the average between two points.
 151
 152         The baseline is fixed once for a given curve and different force measurements,
 153         unless the user wants it to be recalculated
 154         ------------
 155         Syntax: forcebase [rebase]
 156                 rebase: Forces forcebase to ask again the baseline
 157                 max: Instead of asking for a point to measure, asks for two points and use
 158                      the maximum peak in between
 159         '''
 160         baseline_color =  self.GetColorFromConfig('generalvclamp', 'forcebase', 'baseline_color')
 161         baseline_show = self.GetBoolFromConfig('generalvclamp', 'forcebase', 'baseline_show')
 162         baseline_show_in_legend = self.GetBoolFromConfig('generalvclamp', 'forcebase', 'baseline_show_in_legend')
 163         baseline_size = self.GetIntFromConfig('generalvclamp', 'forcebase', 'baseline_size')
 164         decimals = self.GetIntFromConfig('generalvclamp', 'forcebase', 'decimals')
 165         maximum_color =  self.GetColorFromConfig('generalvclamp', 'forcebase', 'maximum_color')
 166         maximum_show = self.GetBoolFromConfig('generalvclamp', 'forcebase', 'maximum_show')
 167         maximum_show_in_legend = self.GetBoolFromConfig('generalvclamp', 'forcebase', 'maximum_show_in_legend')
 168         maximum_size = self.GetIntFromConfig('generalvclamp', 'forcebase', 'maximum_size')
 169         maximumrange_color =  self.GetColorFromConfig('generalvclamp', 'forcebase', 'maximumrange_color')
 170         maximumrange_show = self.GetBoolFromConfig('generalvclamp', 'forcebase', 'maximumrange_show')
 171         maximumrange_show_in_legend = self.GetBoolFromConfig('generalvclamp', 'forcebase', 'maximumrange_show_in_legend')
 172         maximumrange_size = self.GetIntFromConfig('generalvclamp', 'forcebase', 'maximumrange_size')
 173         maxpoint = self.GetBoolFromConfig('generalvclamp', 'forcebase', 'max')
 174         prefix = self.GetStringFromConfig('generalvclamp', 'forcebase', 'prefix')
 175         multiplier = 10 ** lib.prettyformat.get_exponent(prefix)
 176         rebase = self.GetBoolFromConfig('generalvclamp', 'forcebase', 'rebase')
 177         whatset_str = self.GetStringFromConfig('generalvclamp', 'forcebase', 'whatset')
 178         whatset = 'retraction'
 179         if whatset_str == 'extension':
 180             whatset = lh.EXTENSION
 181         if whatset_str == 'retraction':
 182             whatset = lh.RETRACTION
 183
 184         plot = self.GetDisplayedPlotCorrected()
 185
 186         filename = self.GetActiveFile().name
 187         if rebase or (self.basecurrent != filename):
 188             self.basepoints = self._measure_N_points(N=2, message='Click on 2 points to select the baseline.', whatset=whatset)
 189             self.basecurrent = filename
 190
 191         #TODO: maxpoint does not seem to be picking up the 'real' maximum (at least not with test.hkp/default.000)
 192         maximumrange_points = []
 193         maximum_point = []
 194         if maxpoint:
 195             maximumrange_points = self._measure_N_points(N=2, message='Click 2 points to select the range for maximum detection.', whatset=whatset)
 196             boundpoints = [maximumrange_points[0].index, maximumrange_points[1].index]
 197             boundpoints.sort()
 198             try:
 199                 vector_x = plot.curves[whatset].x[boundpoints[0]:boundpoints[1]]
 200                 vector_y = plot.curves[whatset].y[boundpoints[0]:boundpoints[1]]
 201                 y = min(vector_y)
 202                 index = vector_y.index(y)
 203                 maximum_point = [self._clickize(vector_x, vector_y, index)]
 204             except ValueError:
 205                 self.AppendToOutput('Chosen interval not valid. Try picking it again. Did you pick the same point as begin and end of the interval?')
 206                 return
 207         else:
 208             maximum_point = self._measure_N_points(N=1, message='Click on the point to measure.', whatset=whatset)
 209             y = maximum_point[0].graph_coords[1]
 210
 211         boundaries = [self.basepoints[0].index, self.basepoints[1].index]
 212         boundaries.sort()
 213         to_average = plot.curves[whatset].y[boundaries[0]:boundaries[1]] #y points to average
 214
 215         avg = np.mean(to_average)
 216         forcebase = abs(y - avg)
 217
 218         curve = plot.curves[whatset]
 219         if self.basepoints and baseline_show:
 220             #add the baseline points to the plot
 221             baseline = lib.curve.Curve()
 222             baseline.color = baseline_color
 223             if baseline_show_in_legend:
 224                 baseline.label = 'baseline'
 225             else:
 226                 baseline.label = '_nolegend_'
 227             baseline.size = baseline_size
 228             baseline.style = 'scatter'
 229             baseline.units.x = curve.units.x
 230             baseline.units.y = curve.units.y
 231             for point in self.basepoints:
 232                 baseline.x += [point.graph_coords[0]]
 233                 baseline.y += [point.graph_coords[1]]
 234             plot.curves.append(baseline)
 235
 236         if maximumrange_points and maximumrange_show:
 237             #add the range points to the plot
 238             maximumrange = lib.curve.Curve()
 239             maximumrange.color = maximumrange_color
 240             if maximumrange_show_in_legend:
 241                 maximumrange.label = 'maximumrange'
 242             else:
 243                 maximumrange.label = '_nolegend_'
 244             maximumrange.size = maximumrange_size
 245             maximumrange.style = 'scatter'
 246             maximumrange.units.x = curve.units.x
 247             maximumrange.units.y = curve.units.y
 248             for point in maximumrange_points:
 249                 maximumrange.x += [point.graph_coords[0]]
 250                 maximumrange.y += [point.graph_coords[1]]
 251             plot.curves.append(maximumrange)
 252
 253         if maximum_show:
 254             #add the maximum to the plot
 255             maximum = lib.curve.Curve()
 256             maximum.color = maximum_color
 257             if maximum_show_in_legend:
 258                 maximum.label = 'maximum'
 259             else:
 260                 maximum.label = '_nolegend_'
 261             maximum.size = maximum_size
 262             maximum.style = 'scatter'
 263             maximum.units.x = curve.units.x
 264             maximum.units.y = curve.units.y
 265             maximum.x = [maximum_point[0].graph_coords[0]]
 266             maximum.y = [maximum_point[0].graph_coords[1]]
 267             plot.curves.append(maximum)
 268
 269         self.UpdatePlot(plot)
 270
 271         unit_str = plot.curves[whatset].units.y
 272         output_str = lib.prettyformat.pretty_format(forcebase, unit_str, decimals, multiplier)
 273         self.AppendToOutput(''.join(['Force: ', output_str]))
 274
 275     def plotmanip_multiplier(self, plot, current, customvalue=False):
 276         '''
 277         Multiplies all the Y values of an SMFS curve by a value stored in the 'force_multiplier'
 278         configuration variable. Useful for calibrations and other stuff.
 279         '''
 280
 281         #not a smfs curve...
 282         if current.driver.experiment != 'smfs':
 283             return plot
 284
 285         force_multiplier = self.GetFloatFromConfig('generalvclamp', 'force_multiplier')
 286         if force_multiplier == 1:
 287             return plot
 288
 289         plot.curves[lh.EXTENSION].y = [element * force_multiplier for element in plot.curves[lh.EXTENSION].y]
 290         plot.curves[lh.RETRACTION].y = [element * force_multiplier for element in plot.curves[lh.RETRACTION].y]
 291
 292         return plot
 293
 294     def plotmanip_flatten(self, plot, current, customvalue=0):
 295         '''
 296         Subtracts a polynomial fit to the non-contact part of the curve, as to flatten it.
 297         the best polynomial fit is chosen among polynomials of degree 1 to n, where n is
 298         given by the configuration file or by the customvalue.
 299
 300         customvalue = int (>0) --> starts the function even if config says no (default=0)
 301         '''
 302
 303         #not a smfs curve...
 304         if current.driver.experiment != 'smfs':
 305             return current
 306
 307         #config is not flatten, and customvalue flag is false too
 308         #if (not self.config['generalvclamp']['flatten'].as_bool('value')) and (customvalue == 0):
 309         ##TODO: do we need this?
 310         #if (not self.GetBoolFromConfig('generalvclamp', 'flatten')) and (customvalue == 0):
 311             #return plot
 312
 313         max_exponent = 12
 314         delta_contact = 0
 315
 316         if customvalue > 0:
 317             max_cycles = customvalue
 318         else:
 319             #Using > 1 usually doesn't help and can give artefacts. However, it could be useful too.
 320             max_cycles = self.GetIntFromConfig('generalvclamp', 'max_cycles')
 321
 322         contact_index = self.find_contact_point(plot)
 323
 324         valn = [[] for item in range(max_exponent)]
 325         yrn = [0.0 for item in range(max_exponent)]
 326         errn = [0.0 for item in range(max_exponent)]
 327
 328         for i in range(int(max_cycles)):
 329             x_ext = plot.curves[lh.EXTENSION].x[contact_index + delta_contact:]
 330             y_ext = plot.curves[lh.EXTENSION].y[contact_index + delta_contact:]
 331             x_ret = plot.curves[lh.RETRACTION].x[contact_index + delta_contact:]
 332             y_ret = plot.curves[lh.RETRACTION].y[contact_index + delta_contact:]
 333             for exponent in range(max_exponent):
 334                 try:
 335                     valn[exponent] = sp.polyfit(x_ext, y_ext, exponent)
 336                     yrn[exponent] = sp.polyval(valn[exponent], x_ret)
 337                     errn[exponent] = sp.sqrt(sum((yrn[exponent] - y_ext) ** 2) / float(len(y_ext)))
 338                 except Exception, e:
 339                     print 'Cannot flatten!'
 340                     print e
 341                     return current
 342
 343             best_exponent = errn.index(min(errn))
 344
 345             #extension
 346             ycorr_ext = y_ext - yrn[best_exponent] + y_ext[0] #noncontact part
 347             yjoin_ext = np.array(plot.curves[lh.EXTENSION].y[0:contact_index + delta_contact]) #contact part
 348             #retraction
 349             ycorr_ret = y_ret - yrn[best_exponent] + y_ext[0] #noncontact part
 350             yjoin_ret = np.array(plot.curves[lh.RETRACTION].y[0:contact_index + delta_contact]) #contact part
 351
 352             ycorr_ext = np.concatenate((yjoin_ext, ycorr_ext))
 353             ycorr_ret = np.concatenate((yjoin_ret, ycorr_ret))
 354
 355             plot.curves[lh.EXTENSION].y = list(ycorr_ext)
 356             plot.curves[lh.RETRACTION].y = list(ycorr_ret)
 357
 358         return plot
 359
 360     #---SLOPE---
 361     def do_slope(self):
 362         '''
 363         SLOPE
 364         (generalvclamp.py)
 365         Measures the slope of a delimited chunk on the return trace.
 366         The chunk can be delimited either by two manual clicks, or have
 367         a fixed width, given as an argument.
 368         ---------------
 369         Syntax: slope [width]
 370                 The facultative [width] parameter specifies how many
 371                 points will be considered for the fit. If [width] is
 372                 specified, only one click will be required.
 373         Copyright 2008 by Marco Brucale, Massimo Sandal
 374         '''
 375
 376         decimals = self.GetIntFromConfig('generalvclamp', 'slope', 'decimals')
 377         fitspan = self.GetIntFromConfig('generalvclamp', 'slope', 'fitspan')
 378         point_color = self.GetColorFromConfig('generalvclamp', 'slope', 'point_color')
 379         point_show = self.GetBoolFromConfig('generalvclamp', 'slope', 'point_show')
 380         point_show_in_legend = self.GetBoolFromConfig('generalvclamp', 'slope', 'point_show_in_legend')
 381         point_size = self.GetIntFromConfig('generalvclamp', 'slope', 'point_size')
 382         slope_color = self.GetColorFromConfig('generalvclamp', 'slope', 'slope_color')
 383         slope_linewidth = self.GetIntFromConfig('generalvclamp', 'slope', 'slope_linewidth')
 384         slope_show = self.GetBoolFromConfig('generalvclamp', 'slope', 'slope_show')
 385         slope_show_in_legend = self.GetBoolFromConfig('generalvclamp', 'slope', 'slope_show_in_legend')
 386         whatset_str = self.GetStringFromConfig('generalvclamp', 'slope', 'whatset')
 387         whatset = 'retraction'
 388         if whatset_str == 'extension':
 389             whatset = lh.EXTENSION
 390         if whatset_str == 'retraction':
 391             whatset = lh.RETRACTION
 392
 393         # Decides between the two forms of user input
 394         #TODO: add an option 'mode' with options 'chunk' and 'point'
 395         if fitspan == 0:
 396             # Gets the Xs of two clicked points as indexes on the curve curve vector
 397             clicked_points = []
 398             points = self._measure_N_points(N=2, message='Click 2 points to select the chunk.', whatset=whatset)
 399             clicked_points = [points[0].index, points[1].index]
 400             clicked_points.sort()
 401         else:
 402             clicked_points = []
 403             points = self._measure_N_points(N=1, message='Click on the leftmost point of the chunk (i.e.usually the peak).', whatset=whatset)
 404             clicked_points = [points[0].index - fitspan, points[0].index]
 405
 406         # Calls the function linefit_between
 407         parameters = [0, 0, [], []]
 408         try:
 409             parameters = self.linefit_between(clicked_points[0], clicked_points[1], whatset=whatset)
 410         except:
 411             self.AppendToOutput('Cannot fit. Did you click the same point twice?')
 412             return
 413
 414         plot = self.GetDisplayedPlotCorrected()
 415         # Makes a vector with the fitted parameters and sends it to the GUI
 416         xtoplot=parameters[2]
 417         ytoplot=[]
 418         x = 0
 419         for x in xtoplot:
 420             ytoplot.append((x * parameters[0]) + parameters[1])
 421
 422         clickvector_x = []
 423         clickvector_y = []
 424         for item in points:
 425             clickvector_x.append(item.graph_coords[0])
 426             clickvector_y.append(item.graph_coords[1])
 427
 428         if point_show:
 429             #add the clicked point to the plot
 430             point = lib.curve.Curve()
 431             point.color = point_color
 432             if point_show_in_legend:
 433                 point.label = 'clicked point'
 434             else:
 435                 point.label = '_nolegend_'
 436             point.size = point_size
 437             point.style = 'scatter'
 438             point.x = clickvector_x
 439             point.y = clickvector_y
 440             plot.curves.append(point)
 441
 442         if slope_show:
 443             #add the slope to the plot
 444             slope = lib.curve.Curve()
 445             slope.color = slope_color
 446             if slope_show_in_legend:
 447                 slope.label = 'slope'
 448             else:
 449                 slope.label = '_nolegend_'
 450             slope.linewidth = slope_linewidth
 451             slope.style = 'plot'
 452             slope.units.x = plot.curves[whatset].units.x
 453             slope.units.y = plot.curves[whatset].units.y
 454             slope.x = xtoplot
 455             slope.y = ytoplot
 456             plot.curves.append(slope)
 457
 458         self.UpdatePlot(plot)
 459
 460         # Outputs the relevant slope parameter
 461         unit_str = plot.curves[whatset].units.x + '/' + plot.curves[whatset].units.y
 462         output_str = lib.prettyformat.pretty_format(parameters[0], unit_str, decimals, 1)
 463         self.AppendToOutput(''.join(['Slope: ', output_str]))
 464
 465     def linefit_between(self, index1, index2, whatset=lh.RETRACTION):
 466         '''
 467         Creates two vectors (xtofit, ytofit) slicing out from the
 468         curve return trace a portion delimited by the two indeces
 469         given as arguments.
 470         Then does a least squares linear fit on that slice.
 471         Finally returns [0]=the slope, [1]=the intercept of the
 472         fitted 1st grade polynomial, and [2,3]=the actual (x,y) vectors
 473         used for the fit.
 474         Copyright 2008 by Marco Brucale, Massimo Sandal
 475         '''
 476         # Translates the indeces into two vectors containing the x, y data to fit
 477         plot = self.displayed_plot
 478         xtofit = plot.corrected_curves[whatset].x[index1:index2]
 479         ytofit = plot.corrected_curves[whatset].y[index1:index2]
 480
 481         # Does the actual linear fitting (simple least squares with numpy.polyfit)
 482         linefit = np.polyfit(xtofit, ytofit, 1)
 483
 484         return (linefit[0], linefit[1], xtofit, ytofit)
 485
 486
 487