hooke/plugin/flatfilts.py

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