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             pass
 231
 232         defplots=self.current.curve.default_plots()[0] #we need the raw, uncorrected plots
 233
 234         if 'flatten' in self.config['plotmanips']:
 235             flatten=self._find_plotmanip('flatten') #extract flatten plot manipulator
 236             defplots=flatten(defplots, self.current)
 237         else:
 238             print 'You have the flatten plot manipulator not loaded. Enabling it could give you better results.'
 239
 240         peak_location,peak_size=self.has_peaks(defplots,abs_devs)
 241         print 'Found '+str(len(peak_location))+' peaks.'
 242         to_dump='peaks '+self.current.path+' '+str(len(peak_location))
 243         self.outlet.push(to_dump)
 244         #print peak_location
 245
 246         #if no peaks, we have nothing to plot. exit.
 247         if len(peak_location)==0:
 248             return
 249
 250         #otherwise, we plot the peak locations.
 251         xplotted_ret=self.plots[0].vectors[1][0]
 252         yplotted_ret=self.plots[0].vectors[1][1]
 253         xgood=[xplotted_ret[index] for index in peak_location]
 254         ygood=[yplotted_ret[index] for index in peak_location]
 255
 256         recplot=self._get_displayed_plot()
 257         recplot.vectors.append([xgood,ygood])
 258         if recplot.styles==[]:
 259             recplot.styles=[None,None,'scatter']
 260             recplot.colors=[None,None,None]
 261         else:
 262             recplot.styles+=['scatter']
 263             recplot.colors+=[None]
 264
 265         self._send_plot([recplot])
 266
 267     def do_convfilt(self,args):
 268         '''
 269         CONVFILT
 270         (flatfilts.py)
 271         Filters out flat (featureless) curves of the current playlist,
 272         creating a playlist containing only the curves with potential
 273         features.
 274         ------------
 275         Syntax:
 276         convfilt [min_npks min_deviation]
 277
 278         min_npks = minmum number of peaks
 279         (to set the default, see convfilt.conf file; CONVCONF and SETCONF commands)
 280
 281         min_deviation = minimum signal/noise ratio *in the convolution*
 282         (to set the default, see convfilt.conf file; CONVCONF and SETCONF commands)
 283
 284         If called without arguments, it uses default values.
 285         '''
 286
 287         min_npks=self.convfilt_config['minpeaks']
 288         min_deviation=self.convfilt_config['mindeviation']
 289
 290         args=args.split(' ')
 291         if len(args) == 2:
 292             min_npks=int(args[0])
 293             min_deviation=int(args[1])
 294         else:
 295             pass
 296
 297         print 'Processing playlist...'
 298         print '(Please wait)'
 299         notflat_list=[]
 300
 301         c=0
 302         for item in self.current_list:
 303             c+=1
 304
 305             try:
 306                 peak_location,peak_size=self.exec_has_peaks(item,min_deviation)
 307                 if len(peak_location)>=min_npks:
 308                     isok='+'
 309                 else:
 310                     isok=''
 311                 print 'Curve',item.path, 'is',c,'of',len(self.current_list),': found '+str(len(peak_location))+' peaks.'+isok
 312             except:
 313                 peak_location,peak_size=[],[]
 314                 print 'Curve',item.path, 'is',c,'of',len(self.current_list),': cannot be filtered. Probably unable to retrieve force data from corrupt file.'
 315
 316             if len(peak_location)>=min_npks:
 317                 item.peak_location=peak_location
 318                 item.peak_size=peak_size
 319                 item.curve=None #empty the item object, to further avoid memory leak
 320                 notflat_list.append(item)
 321
 322         #Warn that no flattening had been done.
 323         if not ('flatten' in self.config['plotmanips']):
 324             print 'Flatten manipulator was not found. Processing was done without flattening.'
 325             print 'Try to enable it in your configuration file for better results.'
 326
 327         if len(notflat_list)==0:
 328             print 'Found nothing interesting. Check your playlist, could be a bug or criteria could be too much stringent'
 329             return
 330         else:
 331             print 'Found ',len(notflat_list),' potentially interesting curves'
 332             print 'Regenerating playlist...'
 333             self.pointer=0
 334             self.current_list=notflat_list
 335             self.current=self.current_list[self.pointer]
 336             self.do_plot(0)
 337
 338
 339     def do_setconv(self,args):
 340         '''
 341         SETCONV
 342         (flatfilts.py)
 343         Sets the convfilt configuration variables
 344         ------
 345         Syntax: setconv variable value
 346         '''
 347         args=args.split()
 348         #FIXME: a general "set dictionary" function has to be built
 349         if len(args)==0:
 350             print self.convfilt_config
 351         else:
 352             if not (args[0] in self.convfilt_config.keys()):
 353                 print 'This is not an internal convfilt variable!'
 354                 print 'Run "setconv" without arguments to see a list of defined variables.'
 355                 return
 356
 357             if len(args)==1:
 358                 print self.convfilt_config[args[0]]
 359             elif len(args)>1:
 360                 try:
 361                     self.convfilt_config[args[0]]=eval(args[1])
 362                 except NameError: #we have a string argument
 363                     self.convfilt_config[args[0]]=args[1]
 364
 365
 366 #########################
 367 #HANDLING OF CONFIGURATION FILE
 368 class ConvfiltConfig:
 369     '''
 370     Handling of convfilt configuration file
 371
 372     Mostly based on the simple-yet-useful examples of the Python Library Reference
 373     about xml.dom.minidom
 374
 375     FIXME: starting to look a mess, should require refactoring
 376     '''
 377
 378     def __init__(self):
 379         self.config={}
 380
 381
 382     def load_config(self, filename):
 383         myconfig=file(filename)
 384         #the following 3 lines are needed to strip newlines. otherwise, since newlines
 385         #are XML elements too, the parser would read them (and re-save them, multiplying
 386         #newlines...)
 387         #yes, I'm an XML n00b
 388         the_file=myconfig.read()
 389         the_file_lines=the_file.split('\n')
 390         the_file=''.join(the_file_lines)
 391
 392         self.config_tree=xml.dom.minidom.parseString(the_file)
 393
 394         def getText(nodelist):
 395             #take the text from a nodelist
 396             #from Python Library Reference 13.7.2
 397             rc = ''
 398             for node in nodelist:
 399                 if node.nodeType == node.TEXT_NODE:
 400                     rc += node.data
 401             return rc
 402
 403         def handleConfig(config):
 404             noiseabsdev_elements=config.getElementsByTagName("noise_absdev")
 405             convfilt_elements=config.getElementsByTagName("convfilt")
 406             handleAbsdev(noiseabsdev_elements)
 407             handleConvfilt(convfilt_elements)
 408
 409         def handleAbsdev(noiseabsdev_elements):
 410             for element in noiseabsdev_elements:
 411                 for attribute in element.attributes.keys():
 412                     self.config[attribute]=element.getAttribute(attribute)
 413
 414         def handleConvfilt(convfilt_elements):
 415             for element in convfilt_elements:
 416                 for attribute in element.attributes.keys():
 417                     self.config[attribute]=element.getAttribute(attribute)
 418
 419         handleConfig(self.config_tree)
 420         #making items in the dictionary machine-readable
 421         for item in self.config.keys():
 422             try:
 423                 self.config[item]=eval(self.config[item])
 424             except NameError: #if it's an unreadable string, keep it as a string
 425                 pass
 426
 427         return self.config