hooke/plugin/flatfilts-rolf.py

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