plugins/flatfilts.py

   1 #!/usr/bin/env python
   2
   3 '''
   4 flatfilts.py
   5
   6 Force spectroscopy files filtering of flat files.
   7
   8 Plugin dependencies:
   9 procplots.py (plot processing plugin)
  10
  11 Copyright 2008 Massimo Sandal, Fabrizio Benedetti
  12 with modifications by Dr. Rolf Schmidt (Concordia University, Canada)
  13
  14 This program is released under the GNU General Public License version 2.
  15 '''
  16
  17 import lib.libhooke as lh
  18 import wxversion
  19 wxversion.select(lh.WX_GOOD)
  20
  21 import copy
  22 from numpy import diff, mean
  23
  24 import lib.peakspot as lps
  25 import lib.curve
  26
  27 class flatfiltsCommands:
  28
  29     def do_flatfilt(self):
  30         '''
  31         FLATFILT
  32         (flatfilts.py)
  33         Filters out flat (featureless) files of the current playlist,
  34         creating a playlist containing only the files 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
  50         self.AppendToOutput('Processing playlist...')
  51         self.AppendToOutput('(Please wait)')
  52         features = []
  53         playlist = self.GetActivePlaylist()
  54         files = playlist.files
  55         file_index = 0
  56         for current_file in files:
  57             file_index += 1
  58             try:
  59                 current_file.identify(self.drivers)
  60                 notflat = self.has_features(copy.deepcopy(current_file))
  61                 feature_string = ''
  62                 if notflat != 1:
  63                     if notflat > 0:
  64                         feature_string = str(notflat) + ' features'
  65                     else:
  66                         feature_string = 'no features'
  67                 else:
  68                     feature_string = '1 feature'
  69                 output_string = ''.join(['Curve ', current_file.name, '(', str(file_index), '/', str(len(files)), '): ', feature_string])
  70             except:
  71                 notflat = False
  72                 output_string = ''.join(['Curve ', current_file.name, '(', str(file_index), '/', str(len(files)), '): cannot be filtered. Probably unable to retrieve force data from corrupt file.'])
  73             self.AppendToOutput(output_string)
  74             if notflat:
  75                 current_file.features = notflat
  76                 features.append(file_index - 1)
  77         if not features:
  78             self.AppendToOutput('Found nothing interesting. Check the playlist, could be a bug or criteria could be too stringent.')
  79         else:
  80             if len(features) < playlist.count:
  81                 self.AppendToOutput(''.join(['Found ', str(len(features)), ' potentially interesting files.']))
  82                 self.AppendToOutput('Regenerating playlist...')
  83                 playlist_filtered = playlist.filter_curves(features)
  84                 self.AddPlaylist(playlist_filtered, name='flatfilt')
  85             else:
  86                 self.AppendToOutput('No files filtered. Try different filtering criteria.')
  87
  88     def has_features(self, current_file):
  89         '''
  90         decides if a curve is flat enough to be rejected from analysis: it sees if there
  91         are at least min_npks points that are higher than min_deviation times the absolute value
  92         of noise.
  93
  94         Algorithm original idea by Francesco Musiani, with my tweaks and corrections.
  95         '''
  96         #TODO: shoudl medianfilter be variable?
  97         medianfilter = 7
  98         #medianfilter = self.GetIntFromConfig('flatfilts', 'flatfilt', 'median_filter')
  99         mindeviation = self.GetIntFromConfig('flatfilts', 'flatfilt', 'min_deviation')
 100         minpeaks = self.GetIntFromConfig('flatfilts', 'flatfilt', 'min_npks')
 101
 102         retvalue = 0
 103
 104         #we assume the first is the plot with the force curve
 105         #do the median to better resolve features from noise
 106         flat_curve = self.plotmanip_median(current_file.plot, current_file, customvalue=medianfilter)
 107
 108         #absolute value of derivate
 109         yretdiff = diff(flat_curve.curves[lh.RETRACTION].y)
 110         yretdiff = [abs(value) for value in yretdiff]
 111         #average of derivate values
 112         diffmean = mean(yretdiff)
 113         yretdiff.sort()
 114         yretdiff.reverse()
 115         c_pks = 0
 116         for value in yretdiff:
 117             if value / diffmean > mindeviation:
 118                 c_pks += 1
 119             else:
 120                 break
 121
 122         if c_pks >= minpeaks:
 123             retvalue = c_pks
 124
 125         return retvalue
 126
 127     ################################################################
 128     #-----CONVFILT-------------------------------------------------
 129     #-----Convolution-based peak recognition and filtering.
 130     #Requires the peakspot.py library
 131
 132     def has_peaks(self, plot=None, plugin=None):
 133         '''
 134         Finds peak position in a force curve.
 135         '''
 136
 137         if plugin is None:
 138             blindwindow = self.GetFloatFromConfig('flatfilts', 'convfilt', 'blindwindow')
 139             #need to convert the string that contains the list into a list
 140             convolution = eval(self.GetStringFromConfig('flatfilts', 'convfilt', 'convolution'))
 141             maxcut = self.GetFloatFromConfig('flatfilts', 'convfilt', 'maxcut')
 142             mindeviation = self.GetFloatFromConfig('flatfilts', 'convfilt', 'mindeviation')
 143             positive = self.GetBoolFromConfig('flatfilts', 'convfilt', 'positive')
 144             seedouble = self.GetIntFromConfig('flatfilts', 'convfilt', 'seedouble')
 145             stable = self.GetFloatFromConfig('flatfilts', 'convfilt', 'stable')
 146         else:
 147             blindwindow = self.GetFloatFromConfig(plugin.name, plugin.section, plugin.prefix + 'blindwindow')
 148             #need to convert the string that contains the list into a list
 149             convolution = eval(self.GetStringFromConfig(plugin.name, plugin.section, plugin.prefix + 'convolution'))
 150             maxcut = self.GetFloatFromConfig(plugin.name, plugin.section, plugin.prefix + 'maxcut')
 151             mindeviation = self.GetFloatFromConfig(plugin.name, plugin.section, plugin.prefix + 'mindeviation')
 152             positive = self.GetBoolFromConfig(plugin.name, plugin.section, plugin.prefix + 'positive')
 153             seedouble = self.GetIntFromConfig(plugin.name, plugin.section, plugin.prefix + 'seedouble')
 154             stable = self.GetFloatFromConfig(plugin.name, plugin.section, plugin.prefix + 'stable')
 155
 156         if plot is None:
 157             plot = self.GetDisplayedPlotCorrected()
 158
 159         retraction = plot.curves[lh.RETRACTION]
 160         #Calculate convolution
 161         convoluted = lps.conv_dx(retraction.y, convolution)
 162
 163         #surely cut everything before the contact point
 164         cut_index = self.find_contact_point(plot)
 165         #cut even more, before the blind window
 166         start_x = retraction.x[cut_index]
 167         blind_index = 0
 168         for value in retraction.x[cut_index:]:
 169             if abs((value) - (start_x)) > blindwindow * (10 ** -9):
 170                 break
 171             blind_index += 1
 172         cut_index += blind_index
 173         #do the dirty convolution-peak finding stuff
 174         noise_level = lps.noise_absdev(convoluted[cut_index:], positive, maxcut, stable)
 175         above = lps.abovenoise(convoluted, noise_level, cut_index, mindeviation)
 176         peak_location, peak_size = lps.find_peaks(above, seedouble=seedouble)
 177         #take the maximum
 178         for i in range(len(peak_location)):
 179             peak = peak_location[i]
 180             maxpk = min(retraction.y[peak - 10:peak + 10])
 181             index_maxpk = retraction.y[peak - 10:peak + 10].index(maxpk) + (peak - 10)
 182             peak_location[i] = index_maxpk
 183
 184         return peak_location, peak_size
 185
 186     def do_peaks(self, plugin=None, peak_location=None, peak_size=None):
 187         '''
 188         Test command for convolution filter.
 189         ----
 190         Syntax: peaks [deviations]
 191         absolute deviation = number of times the convolution signal is above the noise absolute deviation.
 192         Default is 5.
 193         '''
 194
 195         if plugin is None:
 196             color = self.GetColorFromConfig('flatfilts', 'peaks', 'color')
 197             size = self.GetIntFromConfig('flatfilts', 'peaks', 'size')
 198         else:
 199             color = self.GetColorFromConfig(plugin.name, plugin.section, plugin.prefix + 'color')
 200             size = self.GetIntFromConfig(plugin.name, plugin.section, plugin.prefix + 'size')
 201
 202         plot = self.GetDisplayedPlotCorrected()
 203
 204         if peak_location is None and peak_size is None:
 205             if not self.AppliesPlotmanipulator('flatten'):
 206                 self.AppendToOutput('The flatten plot manipulator is not loaded. Enabling it could give better results.')
 207
 208             peak_location, peak_size = self.has_peaks(plot)
 209             if len(peak_location) != 1:
 210                 peak_str = ' peaks.'
 211             else:
 212                 peak_str = ' peak.'
 213             self.AppendToOutput('Found ' + str(len(peak_location)) + peak_str)
 214
 215         if peak_location:
 216             retraction = plot.curves[lh.RETRACTION]
 217
 218             peaks = lib.curve.Curve()
 219             peaks.color = color
 220             peaks.size = size
 221             peaks.style = 'scatter'
 222             peaks.title = 'Peaks'
 223             peaks.x = [retraction.x[index] for index in peak_location]
 224             peaks.y = [retraction.y[index] for index in peak_location]
 225
 226             plot.curves.append(peaks)
 227
 228             self.UpdatePlot(plot)
 229
 230     def do_convfilt(self):
 231         '''
 232         Filters out flat (featureless) files of the current playlist,
 233         creating a playlist containing only the files with potential
 234         features.
 235         ------------
 236         min_npks: minmum number of peaks
 237         min_deviation: minimum signal/noise ratio *in the convolution*
 238         '''
 239
 240         self.AppendToOutput('Processing playlist...')
 241         self.AppendToOutput('(Please wait)')
 242         apply_plotmanipulators = self.GetStringFromConfig('flatfilts', 'convfilt', 'apply_plotmanipulators')
 243         minpeaks = self.GetIntFromConfig('flatfilts', 'convfilt', 'minpeaks')
 244         features = []
 245         playlist = self.GetActivePlaylist()
 246
 247         files = self.GetActivePlaylist().files
 248         file_index = 0
 249         for current_file in files:
 250             number_of_peaks = 0
 251             file_index += 1
 252             try:
 253                 current_file.identify(self.drivers)
 254                 if apply_plotmanipulators == 'all':
 255                     plot = self.ApplyPlotmanipulators(current_file.plot, current_file)
 256                 if apply_plotmanipulators == 'flatten':
 257                     plotmanipulator = self.GetPlotmanipulator('flatten')
 258                     plot = plotmanipulator.method(current_file.plot, current_file)
 259                 if apply_plotmanipulators == 'none':
 260                     plot = copy.deepcopy(current_file.plot)
 261
 262                 peak_location, peak_size = self.has_peaks(plot)
 263                 number_of_peaks = len(peak_location)
 264                 if number_of_peaks != 1:
 265                     if number_of_peaks > 0:
 266                         feature_string = str(number_of_peaks) + ' features'
 267                     else:
 268                         feature_string = 'no features'
 269                 else:
 270                     feature_string = '1 feature'
 271                 if number_of_peaks >= minpeaks:
 272                     feature_string += '+'
 273                 output_string = ''.join(['Curve ', current_file.name, '(', str(file_index), '/', str(len(files)), '): ', feature_string])
 274             except:
 275                 peak_location = []
 276                 peak_size = []
 277                 output_string = ''.join(['Curve ', current_file.name, '(', str(file_index), '/', str(len(files)), '): cannot be filtered. Probably unable to retrieve force data from corrupt file.'])
 278             self.AppendToOutput(output_string)
 279             if number_of_peaks >= minpeaks:
 280                 current_file.peak_location = peak_location
 281                 current_file.peak_size = peak_size
 282                 features.append(file_index - 1)
 283
 284         #Warn that no flattening had been done.
 285         if not self.HasPlotmanipulator('plotmanip_flatten'):
 286             self.AppendToOutput('Flatten manipulator was not found. Processing was done without flattening.')
 287         else:
 288             if not self.AppliesPlotmanipulator('flatten'):
 289                 self.AppendToOutput('Flatten manipulator was not applied.')
 290                 self.AppendToOutput('Try to enable the flatten plotmanipulator for better results.')
 291
 292         if not features:
 293             self.AppendToOutput('Found nothing interesting. Check the playlist, could be a bug or criteria could be too stringent.')
 294         else:
 295             if len(features) < playlist.count:
 296                 self.AppendToOutput(''.join(['Found ', str(len(features)), ' potentially interesting files.']))
 297                 self.AppendToOutput('Regenerating playlist...')
 298                 playlist_filtered = playlist.filter_curves(features)
 299                 self.AddPlaylist(playlist_filtered, name='convfilt')
 300             else:
 301                 self.AppendToOutput('No files filtered. Try different filtering criteria.')
 302