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=None):
 151         '''
 152         Finds peak position in a force curve.
 153         FIXME: should be moved in libpeakspot.py
 154         '''
 155         if abs_devs==None:
 156             abs_devs=self.convfilt_config['mindeviation']
 157
 158
 159         xret=plot.vectors[1][0]
 160         yret=plot.vectors[1][1]
 161         #Calculate convolution.
 162         convoluted=lps.conv_dx(yret, self.convfilt_config['convolution'])
 163
 164         #surely cut everything before the contact point
 165         cut_index=self.find_contact_point(plot)
 166         #cut even more, before the blind window
 167         start_x=xret[cut_index]
 168         blind_index=0
 169         for value in xret[cut_index:]:
 170             if abs((value) - (start_x)) > self.convfilt_config['blindwindow']*(10**-9):
 171                 break
 172             blind_index+=1
 173         cut_index+=blind_index
 174         #do the dirty convolution-peak finding stuff
 175         noise_level=lps.noise_absdev(convoluted[cut_index:], self.convfilt_config['positive'], self.convfilt_config['maxcut'], self.convfilt_config['stable'])
 176         above=lps.abovenoise(convoluted,noise_level,cut_index,abs_devs)
 177         peak_location,peak_size=lps.find_peaks(above,seedouble=self.convfilt_config['seedouble'])
 178
 179         #take the maximum
 180         for i in range(len(peak_location)):
 181             peak=peak_location[i]
 182             maxpk=min(yret[peak-10:peak+10])
 183             index_maxpk=yret[peak-10:peak+10].index(maxpk)+(peak-10)
 184             peak_location[i]=index_maxpk
 185
 186         return peak_location,peak_size
 187
 188
 189     def exec_has_peaks(self,item,abs_devs):
 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=item.curve.default_plots()[0]
 197
 198         if 'flatten' in self.config['plotmanips']:
 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
 203         peak_location,peak_size=self.has_peaks(plot,abs_devs)
 204         #close all open files
 205         item.curve.close_all()
 206         #needed to avoid *big* memory leaks!
 207         del item.curve
 208         del item
 209         return peak_location, peak_size
 210
 211     #------------------------
 212     #------commands----------
 213     #------------------------
 214     def do_peaks(self,args):
 215         '''
 216         PEAKS
 217         (flatfilts.py)
 218         Test command for convolution filter / test.
 219         ----
 220         Syntax: peaks [deviations]
 221         absolute deviation = number of times the convolution signal is above the noise absolute deviation.
 222         Default is 5.
 223         '''
 224         if len(args)==0:
 225             args=self.convfilt_config['mindeviation']
 226
 227         try:
 228             abs_devs=float(args)
 229         except:
 230             print 'Wrong argument, using config value'
 231             abs_devs=float(self.convfilt_config['mindeviation'])
 232
 233         defplots=self.current.curve.default_plots()[0] #we need the raw, uncorrected plots
 234
 235         if 'flatten' in self.config['plotmanips']:
 236             flatten=self._find_plotmanip('flatten') #extract flatten plot manipulator
 237             defplots=flatten(defplots, self.current)
 238         else:
 239             print 'You have the flatten plot manipulator not loaded. Enabling it could give you better results.'
 240
 241         peak_location,peak_size=self.has_peaks(defplots,abs_devs)
 242         print 'Found '+str(len(peak_location))+' peaks.'
 243         to_dump='peaks '+self.current.path+' '+str(len(peak_location))
 244         self.outlet.push(to_dump)
 245         #print peak_location
 246
 247         #if no peaks, we have nothing to plot. exit.
 248         if len(peak_location)==0:
 249             return
 250
 251         #otherwise, we plot the peak locations.
 252         xplotted_ret=self.plots[0].vectors[1][0]
 253         yplotted_ret=self.plots[0].vectors[1][1]
 254         xgood=[xplotted_ret[index] for index in peak_location]
 255         ygood=[yplotted_ret[index] for index in peak_location]
 256
 257         recplot=self._get_displayed_plot()
 258         recplot.vectors.append([xgood,ygood])
 259         if recplot.styles==[]:
 260             recplot.styles=[None,None,'scatter']
 261             recplot.colors=[None,None,None]
 262         else:
 263             recplot.styles+=['scatter']
 264             recplot.colors+=[None]
 265
 266         self._send_plot([recplot])
 267
 268     def do_convfilt(self,args):
 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         min_npks=self.convfilt_config['minpeaks']
 289         min_deviation=self.convfilt_config['mindeviation']
 290
 291         args=args.split(' ')
 292         if len(args) == 2:
 293             min_npks=int(args[0])
 294             min_deviation=int(args[1])
 295         else:
 296             pass
 297
 298         print 'Processing playlist...'
 299         print '(Please wait)'
 300         notflat_list=[]
 301
 302         c=0
 303         for item in self.current_list:
 304             c+=1
 305
 306             try:
 307                 peak_location,peak_size=self.exec_has_peaks(item,min_deviation)
 308                 if len(peak_location)>=min_npks:
 309                     isok='+'
 310                 else:
 311                     isok=''
 312                 print 'Curve',item.path, 'is',c,'of',len(self.current_list),': found '+str(len(peak_location))+' peaks.'+isok
 313             except:
 314                 peak_location,peak_size=[],[]
 315                 print 'Curve',item.path, 'is',c,'of',len(self.current_list),': cannot be filtered. Probably unable to retrieve force data from corrupt file.'
 316
 317             if len(peak_location)>=min_npks:
 318                 item.peak_location=peak_location
 319                 item.peak_size=peak_size
 320                 item.curve=None #empty the item object, to further avoid memory leak
 321                 notflat_list.append(item)
 322
 323         #Warn that no flattening had been done.
 324         if not ('flatten' in self.config['plotmanips']):
 325             print 'Flatten manipulator was not found. Processing was done without flattening.'
 326             print 'Try to enable it in your configuration file for better results.'
 327
 328         if len(notflat_list)==0:
 329             print 'Found nothing interesting. Check your playlist, could be a bug or criteria could be too much stringent'
 330             return
 331         else:
 332             print 'Found ',len(notflat_list),' potentially interesting curves'
 333             print 'Regenerating playlist...'
 334             self.pointer=0
 335             self.current_list=notflat_list
 336             self.current=self.current_list[self.pointer]
 337             self.do_plot(0)
 338
 339
 340     def do_setconv(self,args):
 341         '''
 342         SETCONV
 343         (flatfilts.py)
 344         Sets the convfilt configuration variables
 345         ------
 346         Syntax: setconv variable value
 347         '''
 348         args=args.split()
 349         #FIXME: a general "set dictionary" function has to be built
 350         if len(args)==0:
 351             print self.convfilt_config
 352         else:
 353             if not (args[0] in self.convfilt_config.keys()):
 354                 print 'This is not an internal convfilt variable!'
 355                 print 'Run "setconv" without arguments to see a list of defined variables.'
 356                 return
 357
 358             if len(args)==1:
 359                 print self.convfilt_config[args[0]]
 360             elif len(args)>1:
 361                 try:
 362                     self.convfilt_config[args[0]]=eval(args[1])
 363                 except NameError: #we have a string argument
 364                     self.convfilt_config[args[0]]=args[1]
 365
 366
 367 #########################
 368 #HANDLING OF CONFIGURATION FILE
 369 class ConvfiltConfig:
 370     '''
 371     Handling of convfilt configuration file
 372
 373     Mostly based on the simple-yet-useful examples of the Python Library Reference
 374     about xml.dom.minidom
 375
 376     FIXME: starting to look a mess, should require refactoring
 377     '''
 378
 379     def __init__(self):
 380         self.config={}
 381
 382
 383     def load_config(self, filename):
 384         myconfig=file(filename)
 385         #the following 3 lines are needed to strip newlines. otherwise, since newlines
 386         #are XML elements too, the parser would read them (and re-save them, multiplying
 387         #newlines...)
 388         #yes, I'm an XML n00b
 389         the_file=myconfig.read()
 390         the_file_lines=the_file.split('\n')
 391         the_file=''.join(the_file_lines)
 392
 393         self.config_tree=xml.dom.minidom.parseString(the_file)
 394
 395         def getText(nodelist):
 396             #take the text from a nodelist
 397             #from Python Library Reference 13.7.2
 398             rc = ''
 399             for node in nodelist:
 400                 if node.nodeType == node.TEXT_NODE:
 401                     rc += node.data
 402             return rc
 403
 404         def handleConfig(config):
 405             noiseabsdev_elements=config.getElementsByTagName("noise_absdev")
 406             convfilt_elements=config.getElementsByTagName("convfilt")
 407             handleAbsdev(noiseabsdev_elements)
 408             handleConvfilt(convfilt_elements)
 409
 410         def handleAbsdev(noiseabsdev_elements):
 411             for element in noiseabsdev_elements:
 412                 for attribute in element.attributes.keys():
 413                     self.config[attribute]=element.getAttribute(attribute)
 414
 415         def handleConvfilt(convfilt_elements):
 416             for element in convfilt_elements:
 417                 for attribute in element.attributes.keys():
 418                     self.config[attribute]=element.getAttribute(attribute)
 419
 420         handleConfig(self.config_tree)
 421         #making items in the dictionary machine-readable
 422         for item in self.config.keys():
 423             try:
 424                 self.config[item]=eval(self.config[item])
 425             except NameError: #if it's an unreadable string, keep it as a string
 426                 pass
 427
 428         return self.config