hooke/plugin/flatfilts-rolf.py

   1 #!/usr/bin/env python
   2
   3 '''
   4 FLATFILTS
   5
   6 Force spectroscopy curves filtering of flat curves
   7 Licensed under the GNU LGPL version 2
   8
   9 Other plugin dependencies:
  10 procplots.py (plot processing plugin)
  11 '''
  12
  13 import xml.dom.minidom
  14
  15 import wx
  16 #import scipy
  17 import numpy
  18 from numpy import diff
  19 import os.path
  20
  21 #import pickle
  22
  23 import libpeakspot as lps
  24 #import libhookecurve as lhc
  25 import hookecurve as lhc
  26 import libhooke as lh
  27 import wxversion
  28 wxversion.select(lh.WX_GOOD)
  29
  30
  31 class flatfiltsCommands:
  32
  33     def do_flatfilt(self):
  34         '''
  35         FLATFILT
  36         (flatfilts.py)
  37         Filters out flat (featureless) curves of the current playlist,
  38         creating a playlist containing only the curves with potential
  39         features.
  40         ------------
  41         Syntax:
  42         flatfilt [min_npks min_deviation]
  43
  44         min_npks = minmum number of points over the deviation
  45         (default=4)
  46
  47         min_deviation = minimum signal/noise ratio
  48         (default=9)
  49
  50         If called without arguments, it uses default values, that
  51         should work most of the times.
  52         '''
  53         #TODO: should this be optional?
  54         medianfilter = 7
  55
  56         self.AppendToOutput('Processing playlist...')
  57         self.AppendToOutput('(Please wait)')
  58         features = []
  59         playlist = self.GetActivePlaylist()
  60         curves = playlist.curves
  61         curve_index = 0
  62         for curve in curves:
  63             curve_index += 1
  64             try:
  65                 notflat = self.has_features(curve)
  66                 feature_string = ''
  67                 if notflat != 1:
  68                     if notflat > 0:
  69                         feature_string = str(notflat) + ' features'
  70                     else:
  71                         feature_string = 'no features'
  72                 else:
  73                     feature_string = '1 feature'
  74                 output_string = ''.join(['Curve ', curve.name, '(', str(curve_index), '/', str(len(curves)), '): ', feature_string])
  75             except:
  76                 notflat = False
  77                 output_string = ''.join(['Curve ', curve.name, '(', str(curve_index), '/', str(len(curves)), '): cannot be filtered. Probably unable to retrieve force data from corrupt file.'])
  78             self.AppendToOutput(output_string)
  79             if notflat:
  80                 curve.features = notflat
  81                 features.append(curve_index - 1)
  82         if not features:
  83             self.AppendToOutput('Found nothing interesting. Check the playlist, could be a bug or criteria could be too stringent.')
  84         else:
  85             if len(features) < playlist.count:
  86                 self.AppendToOutput(''.join(['Found ', str(len(features)), ' potentially interesting curves.']))
  87                 self.AppendToOutput('Regenerating playlist...')
  88                 playlist_filtered = playlist.filter_curves(features)
  89                 self.AddPlaylist(playlist_filtered, name='flatfilt')
  90             else:
  91                 self.AppendToOutput('No curves filtered. Try different filtering criteria.')
  92
  93     def has_features(self, curve):
  94         '''
  95         decides if a curve is flat enough to be rejected from analysis: it sees if there
  96         are at least min_npks points that are higher than min_deviation times the absolute value
  97         of noise.
  98
  99         Algorithm original idea by Francesco Musiani, with my tweaks and corrections.
 100         '''
 101         medianfilter = 7
 102         mindeviation = self.GetIntFromConfig('flatfilts', 'flatfilt', 'min_deviation')
 103         minpeaks = self.GetIntFromConfig('flatfilts', 'flatfilt', 'min_npks')
 104         #medianfilter = self.GetIntFromConfig('flatfilt', 'median_filter')
 105         #mindeviation = self.GetIntFromConfig('convfilt', 'mindeviation')
 106         #minpeaks = self.GetIntFromConfig('convfilt', 'minpeaks')
 107
 108         retvalue = 0
 109
 110         #item.identify(self.drivers)
 111         #we assume the first is the plot with the force curve
 112         #do the median to better resolve features from noise
 113         flat_plot = self.plotmanip_median(curve.driver.default_plots()[0], curve, customvalue=medianfilter)
 114         flat_vects = flat_plot.vectors
 115         curve.driver.close_all()
 116         #needed to avoid *big* memory leaks!
 117         #del item.driver
 118         #del item
 119
 120         #absolute value of derivate
 121         yretdiff=diff(flat_vects[1][1])
 122         yretdiff=[abs(value) for value in yretdiff]
 123         #average of derivate values
 124         diffmean=numpy.mean(yretdiff)
 125         yretdiff.sort()
 126         yretdiff.reverse()
 127         c_pks=0
 128         for value in yretdiff:
 129             if value/diffmean > mindeviation:
 130                 c_pks += 1
 131             else:
 132                 break
 133
 134         if c_pks >= minpeaks:
 135             retvalue = c_pks
 136
 137         del flat_plot, flat_vects, yretdiff
 138
 139         return retvalue
 140
 141     ################################################################
 142     #-----CONVFILT-------------------------------------------------
 143     #-----Convolution-based peak recognition and filtering.
 144     #Requires the libpeakspot.py library
 145
 146     def has_peaks(self, plot, curve=None):
 147         '''
 148         Finds peak position in a force curve.
 149         FIXME: should be moved to libpeakspot.py
 150         '''
 151
 152         blindwindow = self.GetFloatFromConfig('flatfilts', 'convfilt', 'blindwindow')
 153         #need to convert the string that contains the list into a list
 154         convolution = eval(self.GetStringFromConfig('flatfilts', 'convfilt', 'convolution'))
 155         maxcut = self.GetFloatFromConfig('flatfilts', 'convfilt', 'maxcut')
 156         mindeviation = self.GetFloatFromConfig('flatfilts', 'convfilt', 'mindeviation')
 157         positive = self.GetBoolFromConfig('flatfilts', 'convfilt', 'positive')
 158         seedouble = self.GetIntFromConfig('flatfilts', 'convfilt', 'seedouble')
 159         stable = self.GetFloatFromConfig('flatfilts', 'convfilt', 'stable')
 160
 161         xret = plot.vectors[1][0]
 162         yret = plot.vectors[1][1]
 163         #Calculate convolution
 164         convoluted = lps.conv_dx(yret, convolution)
 165
 166         #surely cut everything before the contact point
 167         cut_index = self.find_contact_point(plot, curve)
 168         #cut even more, before the blind window
 169         start_x = xret[cut_index]
 170         blind_index = 0
 171         for value in xret[cut_index:]:
 172             if abs((value) - (start_x)) > blindwindow * (10 ** -9):
 173                 break
 174             blind_index += 1
 175         cut_index += blind_index
 176         #do the dirty convolution-peak finding stuff
 177         noise_level = lps.noise_absdev(convoluted[cut_index:], positive, maxcut, stable)
 178         above = lps.abovenoise(convoluted, noise_level, cut_index, mindeviation)
 179         peak_location, peak_size = lps.find_peaks(above, seedouble=seedouble)
 180         #take the maximum
 181         for i in range(len(peak_location)):
 182             peak = peak_location[i]
 183             maxpk = min(yret[peak - 10:peak + 10])
 184             index_maxpk = yret[peak - 10:peak + 10].index(maxpk) + (peak - 10)
 185             peak_location[i] = index_maxpk
 186
 187         return peak_location, peak_size
 188
 189     def exec_has_peaks(self, curve):
 190         '''
 191         encapsulates has_peaks for the purpose of correctly treating the curve objects in the convfilt loop,
 192         to avoid memory leaks
 193         '''
 194         #item.identify(self.drivers)
 195         #we assume the first is the plot with the force curve
 196         plot = curve.driver.default_plots()[0]
 197
 198         if self.HasPlotmanipulator('plotmanip_flatten'):
 199             #If flatten is present, use it for better recognition of peaks...
 200             #flatten = self._find_plotmanip('flatten') #extract flatten plot manipulator
 201             #plot = flatten(plot, item, customvalue=1)
 202             plot = self.plotmanip_flatten(plot, curve, customvalue=1)
 203
 204         peak_location, peak_size = self.has_peaks(plot, curve)
 205         #close all open files
 206         curve.driver.close_all()
 207         #needed to avoid *big* memory leaks!
 208         #del item.driver
 209         #del item
 210         return peak_location, peak_size
 211
 212     #------------------------
 213     #------commands----------
 214     #------------------------
 215     def do_peaks(self):
 216         '''
 217         PEAKS
 218         (flatfilts.py)
 219         Test command for convolution filter / test.
 220         ----
 221         Syntax: peaks [deviations]
 222         absolute deviation = number of times the convolution signal is above the noise absolute deviation.
 223         Default is 5.
 224         '''
 225
 226         #TODO: check if the following line gives us what we need
 227         curve = self.GetActiveCurve()
 228         defplots = curve.driver.default_plots()[0] #we need the raw, uncorrected plots
 229
 230         #if 'flatten' in self.config['plotmanips']:
 231         if self.HasPlotmanipulator('plotmanip_flatten'):
 232             #flatten=self._find_plotmanip('flatten') #extract flatten plot manipulator
 233             #defplots=flatten(defplots, self.current)
 234             defplots = self.plotmanip_flatten(defplots, curve, customvalue=0)
 235         else:
 236             self.AppendToOutput('The flatten plot manipulator is not loaded. Enabling it could give better results.')
 237
 238         peak_location, peak_size = self.has_peaks(defplots, curve)
 239         self.AppendToOutput('Found ' + str(len(peak_location)) + ' peaks.')
 240         self.AppendToOutput('peaks ' + curve.filename + ' ' + str(len(peak_location)))
 241         #to_dump = 'peaks ' + current_curve.filename + ' ' + str(len(peak_location))
 242         #self.outlet.push(to_dump)
 243
 244         #if no peaks, we have nothing to plot. exit.
 245         if peak_location:
 246             #otherwise, we plot the peak locations.
 247             xplotted_ret = curve.plots[0].vectors[1][0]
 248             yplotted_ret = curve.plots[0].vectors[1][1]
 249             xgood = [xplotted_ret[index] for index in peak_location]
 250             ygood = [yplotted_ret[index] for index in peak_location]
 251
 252             curve.plots[0].add_set(xgood, ygood)
 253             curve.plots[0].styles.append('scatter')
 254             curve.plots[0].colors.append('indigo')
 255
 256             #recplot = self._get_displayed_plot()
 257             #recplot.vectors.append([xgood,ygood])
 258             #if recplot.styles == []:
 259                 #recplot.styles = [None, None, 'scatter']
 260                 #recplot.colors = [None, None, None]
 261             #else:
 262                 #recplot.styles += ['scatter']
 263                 #recplot.colors += [None]
 264
 265             #self._send_plot([recplot])
 266             self.UpdatePlot()
 267
 268     def do_convfilt(self):
 269         '''
 270         CONVFILT
 271         (flatfilts.py)
 272         Filters out flat (featureless) curves of the current playlist,
 273         creating a playlist containing only the curves with potential
 274         features.
 275         ------------
 276         Syntax:
 277         convfilt [min_npks min_deviation]
 278
 279         min_npks = minmum number of peaks
 280         (to set the default, see convfilt.conf file; CONVCONF and SETCONF commands)
 281
 282         min_deviation = minimum signal/noise ratio *in the convolution*
 283         (to set the default, see convfilt.conf file; CONVCONF and SETCONF commands)
 284
 285         If called without arguments, it uses default values.
 286         '''
 287
 288         self.AppendToOutput('Processing playlist...')
 289         self.AppendToOutput('(Please wait)')
 290         minpeaks = self.GetIntFromConfig('flatfilts', 'convfilt', 'minpeaks')
 291         features = []
 292         playlist = self.GetActivePlaylist()
 293
 294         curves = self.GetActivePlaylist().curves
 295         curve_index = 0
 296         for curve in curves:
 297             curve_index += 1
 298             try:
 299                 peak_location, peak_size = self.exec_has_peaks(curve)
 300                 number_of_peaks = len(peak_location)
 301                 if number_of_peaks != 1:
 302                     if number_of_peaks > 0:
 303                         feature_string = str(number_of_peaks) + ' features'
 304                     else:
 305                         feature_string = 'no features'
 306                 else:
 307                     feature_string = '1 feature'
 308                 if number_of_peaks >= minpeaks:
 309                     feature_string += '+'
 310                 output_string = ''.join(['Curve ', curve.name, '(', str(curve_index), '/', str(len(curves)), '): ', feature_string])
 311             except:
 312                 peak_location = []
 313                 peak_size = []
 314                 output_string = ''.join(['Curve ', curve.name, '(', str(curve_index), '/', str(len(curves)), '): cannot be filtered. Probably unable to retrieve force data from corrupt file.'])
 315             self.AppendToOutput(output_string)
 316             if number_of_peaks >= minpeaks:
 317                 curve.peak_location = peak_location
 318                 curve.peak_size = peak_size
 319                 features.append(curve_index - 1)
 320
 321         #TODO: do we need this? Flattening might not be necessary/desired
 322         #Warn that no flattening had been done.
 323         if not self.HasPlotmanipulator('plotmanip_flatten'):
 324             self.AppendToOutput('Flatten manipulator was not found. Processing was done without flattening.')
 325             self.AppendToOutput('Try to enable it in the configuration file for better results.')
 326         if not features:
 327             self.AppendToOutput('Found nothing interesting. Check the playlist, could be a bug or criteria could be too stringent.')
 328         else:
 329             if len(features) < playlist.count:
 330                 self.AppendToOutput(''.join(['Found ', str(len(features)), ' potentially interesting curves.']))
 331                 self.AppendToOutput('Regenerating playlist...')
 332                 playlist_filtered = playlist.filter_curves(features)
 333                 self.AddPlaylist(playlist_filtered, name='convfilt')
 334             else:
 335                 self.AppendToOutput('No curves filtered. Try different filtering criteria.')
 336