flatfilts.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 from libhooke import WX_GOOD
  13 import wxversion
  14 wxversion.select(WX_GOOD)
  15
  16 import xml.dom.minidom
  17
  18 import wx
  19 import scipy
  20 import numpy
  21 from numpy import diff
  22
  23 #import pickle
  24
  25 import libpeakspot as lps
  26 import libhookecurve as lhc
  27
  28
  29 class flatfiltsCommands:
  30
  31     def _plug_init(self):
  32         #configurate convfilt variables
  33         convfilt_configurator=ConvfiltConfig()
  34
  35         #different OSes have different path conventions
  36         if self.config['hookedir'][0]=='/':
  37             slash='/' #a Unix or Unix-like system
  38         else:
  39             slash='\\' #it's a drive letter, we assume it's Windows
  40
  41         self.convfilt_config=convfilt_configurator.load_config(self.config['hookedir']+slash+'convfilt.conf')
  42
  43     def do_flatfilt(self,args):
  44         '''
  45         FLATFILT
  46         (flatfilts.py)
  47         Filters out flat (featureless) curves of the current playlist,
  48         creating a playlist containing only the curves with potential
  49         features.
  50         ------------
  51         Syntax:
  52         flatfilt [min_npks min_deviation]
  53
  54         min_npks = minmum number of points over the deviation
  55         (default=4)
  56
  57         min_deviation = minimum signal/noise ratio
  58         (default=9)
  59
  60         If called without arguments, it uses default values, that
  61         should work most of the times.
  62         '''
  63         median_filter=7
  64         min_npks=4
  65         min_deviation=9
  66
  67         args=args.split(' ')
  68         if len(args) == 2:
  69             min_npks=int(args[0])
  70             min_deviation=int(args[1])
  71         else:
  72             pass
  73
  74         print 'Processing playlist...'
  75         notflat_list=[]
  76
  77         c=0
  78         for item in self.current_list:
  79             c+=1
  80
  81             try:
  82                 notflat=self.has_features(item, median_filter, min_npks, min_deviation)
  83                 print 'Curve',item.path, 'is',c,'of',len(self.current_list),': features are ',notflat
  84             except:
  85                 notflat=False
  86                 print 'Curve',item.path, 'is',c,'of',len(self.current_list),': cannot be filtered. Probably unable to retrieve force data from corrupt file.'
  87
  88             if notflat:
  89                 item.features=notflat
  90                 item.curve=None #empty the item object, to further avoid memory leak
  91                 notflat_list.append(item)
  92
  93         if len(notflat_list)==0:
  94             print 'Found nothing interesting. Check your playlist, could be a bug or criteria could be too much stringent'
  95             return
  96         else:
  97             print 'Found ',len(notflat_list),' potentially interesting curves'
  98             print 'Regenerating playlist...'
  99             self.pointer=0
 100             self.current_list=notflat_list
 101             self.current=self.current_list[self.pointer]
 102             self.do_plot(0)
 103
 104     def has_features(self,item,median_filter,min_npks,min_deviation):
 105         '''
 106         decides if a curve is flat enough to be rejected from analysis: it sees if there
 107         are at least min_npks points that are higher than min_deviation times the absolute value
 108         of noise.
 109
 110         Algorithm original idea by Francesco Musiani, with my tweaks and corrections.
 111         '''
 112         retvalue=False
 113
 114         item.identify(self.drivers)
 115         #we assume the first is the plot with the force curve
 116         #do the median to better resolve features from noise
 117         flat_plot=self.plotmanip_median(item.curve.default_plots()[0], item, customvalue=median_filter)
 118         flat_vects=flat_plot.vectors
 119         item.curve.close_all()
 120         #needed to avoid *big* memory leaks!
 121         del item.curve
 122         del item
 123
 124         #absolute value of derivate
 125         yretdiff=diff(flat_vects[1][1])
 126         yretdiff=[abs(value) for value in yretdiff]
 127         #average of derivate values
 128         diffmean=numpy.mean(yretdiff)
 129         yretdiff.sort()
 130         yretdiff.reverse()
 131         c_pks=0
 132         for value in yretdiff:
 133             if value/diffmean > min_deviation:
 134                 c_pks+=1
 135             else:
 136                 break
 137
 138         if c_pks>=min_npks:
 139             retvalue = c_pks
 140
 141         del flat_plot, flat_vects, yretdiff
 142
 143         return retvalue
 144
 145     ################################################################
 146     #-----CONVFILT-------------------------------------------------
 147     #-----Convolution-based peak recognition and filtering.
 148     #Requires the libpeakspot.py library
 149
 150     def has_peaks(self, plot, abs_devs):
 151         '''
 152         Finds peak position in a force curve.
 153         FIXME: should be moved in libpeakspot.py
 154         '''
 155         xret=plot.vectors[1][0]
 156         yret=plot.vectors[1][1]
 157         #Calculate convolution.
 158         convoluted=lps.conv_dx(yret, self.convfilt_config['convolution'])
 159
 160         #surely cut everything before the contact point
 161         cut_index=self.find_contact_point(plot)
 162         #cut even more, before the blind window
 163         start_x=xret[cut_index]
 164         blind_index=0
 165         for value in xret[cut_index:]:
 166             if abs((value) - (start_x)) > self.convfilt_config['blindwindow']*(10**-9):
 167                 break
 168             blind_index+=1
 169         cut_index+=blind_index
 170         #do the dirty convolution-peak finding stuff
 171         noise_level=lps.noise_absdev(convoluted[cut_index:], self.convfilt_config['positive'], self.convfilt_config['maxcut'], self.convfilt_config['stable'])
 172         above=lps.abovenoise(convoluted,noise_level,cut_index,abs_devs)
 173         peak_location,peak_size=lps.find_peaks(above,seedouble=self.convfilt_config['seedouble'])
 174
 175         #take the maximum
 176         for i in range(len(peak_location)):
 177             peak=peak_location[i]
 178             maxpk=min(yret[peak-10:peak+10])
 179             index_maxpk=yret[peak-10:peak+10].index(maxpk)+(peak-10)
 180             peak_location[i]=index_maxpk
 181
 182         return peak_location,peak_size
 183
 184
 185     def exec_has_peaks(self,item,abs_devs):
 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=item.curve.default_plots()[0]
 193
 194         if 'flatten' in self.config['plotmanips']:
 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
 199         peak_location,peak_size=self.has_peaks(plot,abs_devs)
 200         #close all open files
 201         item.curve.close_all()
 202         #needed to avoid *big* memory leaks!
 203         del item.curve
 204         del item
 205         return peak_location, peak_size
 206
 207     #------------------------
 208     #------commands----------
 209     #------------------------
 210     def do_peaks(self,args):
 211         '''
 212         PEAKS
 213         (flatfilts.py)
 214         Test command for convolution filter / test.
 215         ----
 216         Syntax: peaks [deviations]
 217         absolute deviation = number of times the convolution signal is above the noise absolute deviation.
 218         Default is 5.
 219         '''
 220         if len(args)==0:
 221             args=self.convfilt_config['mindeviation']
 222
 223         try:
 224             abs_devs=float(args)
 225         except:
 226             pass
 227
 228         defplots=self.current.curve.default_plots()[0] #we need the raw, uncorrected plots
 229
 230         if 'flatten' in self.config['plotmanips']:
 231             flatten=self._find_plotmanip('flatten') #extract flatten plot manipulator
 232             defplots=flatten(defplots, self.current)
 233         else:
 234             print 'You have the flatten plot manipulator not loaded. Enabling it could give you better results.'
 235
 236         peak_location,peak_size=self.has_peaks(defplots,abs_devs)
 237         print 'Found '+str(len(peak_location))+' peaks.'
 238         to_dump='peaks '+self.current.path+' '+str(len(peak_location))
 239         self.outlet.push(to_dump)
 240         #print peak_location
 241
 242         #if no peaks, we have nothing to plot. exit.
 243         if len(peak_location)==0:
 244             return
 245
 246         #otherwise, we plot the peak locations.
 247         xplotted_ret=self.plots[0].vectors[1][0]
 248         yplotted_ret=self.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         recplot=self._get_displayed_plot()
 253         recplot.vectors.append([xgood,ygood])
 254         if recplot.styles==[]:
 255             recplot.styles=[None,None,'scatter']
 256         else:
 257             recplot.styles+=['scatter']
 258
 259         self._send_plot([recplot])
 260
 261     def do_convfilt(self,args):
 262         '''
 263         CONVFILT
 264         (flatfilts.py)
 265         Filters out flat (featureless) curves of the current playlist,
 266         creating a playlist containing only the curves with potential
 267         features.
 268         ------------
 269         Syntax:
 270         convfilt [min_npks min_deviation]
 271
 272         min_npks = minmum number of peaks
 273         (to set the default, see convfilt.conf file; CONVCONF and SETCONF commands)
 274
 275         min_deviation = minimum signal/noise ratio *in the convolution*
 276         (to set the default, see convfilt.conf file; CONVCONF and SETCONF commands)
 277
 278         If called without arguments, it uses default values.
 279         '''
 280
 281         min_npks=self.convfilt_config['minpeaks']
 282         min_deviation=self.convfilt_config['mindeviation']
 283
 284         args=args.split(' ')
 285         if len(args) == 2:
 286             min_npks=int(args[0])
 287             min_deviation=int(args[1])
 288         else:
 289             pass
 290
 291         print 'Processing playlist...'
 292         print '(Please wait)'
 293         notflat_list=[]
 294
 295         c=0
 296         for item in self.current_list:
 297             c+=1
 298
 299             try:
 300                 peak_location,peak_size=self.exec_has_peaks(item,min_deviation)
 301                 if len(peak_location)>=min_npks:
 302                     isok='+'
 303                 else:
 304                     isok=''
 305                 print 'Curve',item.path, 'is',c,'of',len(self.current_list),': found '+str(len(peak_location))+' peaks.'+isok
 306             except:
 307                 peak_location,peak_size=[],[]
 308                 print 'Curve',item.path, 'is',c,'of',len(self.current_list),': cannot be filtered. Probably unable to retrieve force data from corrupt file.'
 309
 310             if len(peak_location)>=min_npks:
 311                 item.peak_location=peak_location
 312                 item.peak_size=peak_size
 313                 item.curve=None #empty the item object, to further avoid memory leak
 314                 notflat_list.append(item)
 315
 316         #Warn that no flattening had been done.
 317         if not ('flatten' in self.config['plotmanips']):
 318             print 'Flatten manipulator was not found. Processing was done without flattening.'
 319             print 'Try to enable it in your configuration file for better results.'
 320
 321         if len(notflat_list)==0:
 322             print 'Found nothing interesting. Check your playlist, could be a bug or criteria could be too much stringent'
 323             return
 324         else:
 325             print 'Found ',len(notflat_list),' potentially interesting curves'
 326             print 'Regenerating playlist...'
 327             self.pointer=0
 328             self.current_list=notflat_list
 329             self.current=self.current_list[self.pointer]
 330             self.do_plot(0)
 331
 332
 333     def do_setconv(self,args):
 334         '''
 335         SETCONV
 336         (flatfilts.py)
 337         Sets the convfilt configuration variables
 338         ------
 339         Syntax: setconv variable value
 340         '''
 341         args=args.split()
 342         #FIXME: a general "set dictionary" function has to be built
 343         if len(args)==0:
 344             print self.convfilt_config
 345         else:
 346             if not (args[0] in self.convfilt_config.keys()):
 347                 print 'This is not an internal convfilt variable!'
 348                 print 'Run "setconv" without arguments to see a list of defined variables.'
 349                 return
 350
 351             if len(args)==1:
 352                 print self.convfilt_config[args[0]]
 353             elif len(args)>1:
 354                 try:
 355                     self.convfilt_config[args[0]]=eval(args[1])
 356                 except NameError: #we have a string argument
 357                     self.convfilt_config[args[0]]=args[1]
 358
 359
 360 #########################
 361 #HANDLING OF CONFIGURATION FILE
 362 class ConvfiltConfig:
 363     '''
 364     Handling of convfilt configuration file
 365
 366     Mostly based on the simple-yet-useful examples of the Python Library Reference
 367     about xml.dom.minidom
 368
 369     FIXME: starting to look a mess, should require refactoring
 370     '''
 371
 372     def __init__(self):
 373         self.config={}
 374
 375
 376     def load_config(self, filename):
 377         myconfig=file(filename)
 378         #the following 3 lines are needed to strip newlines. otherwise, since newlines
 379         #are XML elements too, the parser would read them (and re-save them, multiplying
 380         #newlines...)
 381         #yes, I'm an XML n00b
 382         the_file=myconfig.read()
 383         the_file_lines=the_file.split('\n')
 384         the_file=''.join(the_file_lines)
 385
 386         self.config_tree=xml.dom.minidom.parseString(the_file)
 387
 388         def getText(nodelist):
 389             #take the text from a nodelist
 390             #from Python Library Reference 13.7.2
 391             rc = ''
 392             for node in nodelist:
 393                 if node.nodeType == node.TEXT_NODE:
 394                     rc += node.data
 395             return rc
 396
 397         def handleConfig(config):
 398             noiseabsdev_elements=config.getElementsByTagName("noise_absdev")
 399             convfilt_elements=config.getElementsByTagName("convfilt")
 400             handleAbsdev(noiseabsdev_elements)
 401             handleConvfilt(convfilt_elements)
 402
 403         def handleAbsdev(noiseabsdev_elements):
 404             for element in noiseabsdev_elements:
 405                 for attribute in element.attributes.keys():
 406                     self.config[attribute]=element.getAttribute(attribute)
 407
 408         def handleConvfilt(convfilt_elements):
 409             for element in convfilt_elements:
 410                 for attribute in element.attributes.keys():
 411                     self.config[attribute]=element.getAttribute(attribute)
 412
 413         handleConfig(self.config_tree)
 414         #making items in the dictionary machine-readable
 415         for item in self.config.keys():
 416             try:
 417                 self.config[item]=eval(self.config[item])
 418             except NameError: #if it's an unreadable string, keep it as a string
 419                 pass
 420
 421         return self.config