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 ???? by ?
  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):
 133         '''
 134         Finds peak position in a force curve.
 135         FIXME: should be moved to peakspot.py
 136         '''
 137
 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
 147         if plot is None:
 148             plot = self.GetActivePlot()
 149
 150         xret = plot.curves[lh.RETRACTION].x
 151         yret = plot.curves[lh.RETRACTION].y
 152         #Calculate convolution
 153         convoluted = lps.conv_dx(yret, convolution)
 154
 155         #surely cut everything before the contact point
 156         cut_index = self.find_contact_point(plot)
 157         #cut even more, before the blind window
 158         start_x = xret[cut_index]
 159         blind_index = 0
 160         for value in xret[cut_index:]:
 161             if abs((value) - (start_x)) > blindwindow * (10 ** -9):
 162                 break
 163             blind_index += 1
 164         cut_index += blind_index
 165         #do the dirty convolution-peak finding stuff
 166         noise_level = lps.noise_absdev(convoluted[cut_index:], positive, maxcut, stable)
 167         above = lps.abovenoise(convoluted, noise_level, cut_index, mindeviation)
 168         peak_location, peak_size = lps.find_peaks(above, seedouble=seedouble)
 169         #take the maximum
 170         for i in range(len(peak_location)):
 171             peak = peak_location[i]
 172             maxpk = min(yret[peak - 10:peak + 10])
 173             index_maxpk = yret[peak - 10:peak + 10].index(maxpk) + (peak - 10)
 174             peak_location[i] = index_maxpk
 175
 176         return peak_location, peak_size
 177
 178     def exec_has_peaks(self, item):
 179         '''
 180         encapsulates has_peaks for the purpose of correctly treating the curve objects in the convfilt loop,
 181         to avoid memory leaks
 182         '''
 183         #TODO: things have changed, check for the memory leak again
 184
 185         if self.HasPlotmanipulator('plotmanip_flatten'):
 186             #If flatten is present, use it for better recognition of peaks...
 187             flat_plot = self.plotmanip_flatten(item.plot, item, customvalue=1)
 188
 189         peak_location, peak_size = self.has_peaks(flat_plot)
 190
 191         return peak_location, peak_size
 192
 193     def do_peaks(self, plugin=None, peak_location=None, peak_size=None):
 194         '''
 195         PEAKS
 196         (flatfilts.py)
 197         Test command for convolution filter / test.
 198         ----
 199         Syntax: peaks [deviations]
 200         absolute deviation = number of times the convolution signal is above the noise absolute deviation.
 201         Default is 5.
 202         '''
 203
 204         if plugin is None:
 205             color = self.GetColorFromConfig('flatfilts', 'peaks', 'color')
 206             size = self.GetIntFromConfig('flatfilts', 'peaks', 'size')
 207         else:
 208             color = self.GetColorFromConfig(plugin.name, plugin.section, plugin.prefix + 'color')
 209             size = self.GetIntFromConfig(plugin.name, plugin.section, plugin.prefix + 'size')
 210
 211         plot = self.GetDisplayedPlotCorrected()
 212
 213         if peak_location is None and peak_size is None:
 214             if not self.AppliesPlotmanipulator('flatten'):
 215                 self.AppendToOutput('The flatten plot manipulator is not loaded. Enabling it could give better results.')
 216
 217             peak_location, peak_size = self.has_peaks(plot)
 218             if len(peak_location) != 1:
 219                 peak_str = ' peaks.'
 220             else:
 221                 peak_str = ' peak.'
 222             self.AppendToOutput('Found ' + str(len(peak_location)) + peak_str)
 223
 224         if peak_location:
 225             xplotted_ret = plot.curves[lh.RETRACTION].x
 226             yplotted_ret = plot.curves[lh.RETRACTION].y
 227             xgood = [xplotted_ret[index] for index in peak_location]
 228             ygood = [yplotted_ret[index] for index in peak_location]
 229
 230             peaks = lib.curve.Curve()
 231             peaks.color = color
 232             peaks.size = size
 233             peaks.style = 'scatter'
 234             peaks.x = xgood
 235             peaks.y = ygood
 236
 237             plot.curves.append(peaks)
 238
 239             self.UpdatePlot(plot)
 240
 241     def do_convfilt(self):
 242         '''
 243         CONVFILT
 244         (flatfilts.py)
 245         Filters out flat (featureless) files of the current playlist,
 246         creating a playlist containing only the files with potential
 247         features.
 248         ------------
 249         Syntax:
 250         convfilt [min_npks min_deviation]
 251
 252         min_npks = minmum number of peaks
 253         (to set the default, see convfilt.conf file; CONVCONF and SETCONF commands)
 254
 255         min_deviation = minimum signal/noise ratio *in the convolution*
 256         (to set the default, see convfilt.conf file; CONVCONF and SETCONF commands)
 257
 258         If called without arguments, it uses default values.
 259         '''
 260
 261         self.AppendToOutput('Processing playlist...')
 262         self.AppendToOutput('(Please wait)')
 263         minpeaks = self.GetIntFromConfig('flatfilts', 'convfilt', 'minpeaks')
 264         features = []
 265         playlist = self.GetActivePlaylist()
 266
 267         files = self.GetActivePlaylist().files
 268         file_index = 0
 269         for current_file in files:
 270             number_of_peaks = 0
 271             file_index += 1
 272             try:
 273                 current_file.identify(self.drivers)
 274                 peak_location, peak_size = self.exec_has_peaks(copy.deepcopy(current_file))
 275                 number_of_peaks = len(peak_location)
 276                 if number_of_peaks != 1:
 277                     if number_of_peaks > 0:
 278                         feature_string = str(number_of_peaks) + ' features'
 279                     else:
 280                         feature_string = 'no features'
 281                 else:
 282                     feature_string = '1 feature'
 283                 if number_of_peaks >= minpeaks:
 284                     feature_string += '+'
 285                 output_string = ''.join(['Curve ', current_file.name, '(', str(file_index), '/', str(len(files)), '): ', feature_string])
 286             except:
 287                 peak_location = []
 288                 peak_size = []
 289                 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.'])
 290             self.AppendToOutput(output_string)
 291             if number_of_peaks >= minpeaks:
 292                 current_file.peak_location = peak_location
 293                 current_file.peak_size = peak_size
 294                 features.append(file_index - 1)
 295
 296         #Warn that no flattening had been done.
 297         if not self.HasPlotmanipulator('plotmanip_flatten'):
 298             self.AppendToOutput('Flatten manipulator was not found. Processing was done without flattening.')
 299             self.AppendToOutput('Try to enable it in the configuration file for better results.')
 300         if not features:
 301             self.AppendToOutput('Found nothing interesting. Check the playlist, could be a bug or criteria could be too stringent.')
 302         else:
 303             if len(features) < playlist.count:
 304                 self.AppendToOutput(''.join(['Found ', str(len(features)), ' potentially interesting files.']))
 305                 self.AppendToOutput('Regenerating playlist...')
 306                 playlist_filtered = playlist.filter_curves(features)
 307                 self.AddPlaylist(playlist_filtered, name='convfilt')
 308             else:
 309                 self.AppendToOutput('No files filtered. Try different filtering criteria.')
 310