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