e0283b5e6e0e14b11b816c2849244ece6f7d1670
[hooke.git] / 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()
162         
163         #cut even more, before the blind window
164         start_x=xret[cut_index]
165         blind_index=0
166         for value in xret[cut_index:]:
167             if abs((value) - (start_x)) > self.convfilt_config['blindwindow']*(10**-9):
168                 break
169             blind_index+=1
170         cut_index+=blind_index
171         
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)
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            
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         else:
261             recplot.styles+=['scatter']
262         
263         self._send_plot([recplot])
264         
265     def do_convfilt(self,args):
266         '''
267         CONVFILT
268         (flatfilts.py)
269         Filters out flat (featureless) curves of the current playlist,
270         creating a playlist containing only the curves with potential
271         features.
272         ------------
273         Syntax:
274         convfilt [min_npks min_deviation]
275
276         min_npks = minmum number of peaks
277         (to set the default, see convfilt.conf file; CONVCONF and SETCONF commands)
278
279         min_deviation = minimum signal/noise ratio *in the convolution*
280         (to set the default, see convfilt.conf file; CONVCONF and SETCONF commands)
281
282         If called without arguments, it uses default values.
283         '''
284         
285         min_npks=self.convfilt_config['minpeaks']
286         min_deviation=self.convfilt_config['mindeviation']
287         
288         args=args.split(' ')
289         if len(args) == 2:
290             min_npks=int(args[0])
291             min_deviation=int(args[1])
292         else:
293             pass
294         
295         print 'Processing playlist...'
296         print '(Please wait)'
297         notflat_list=[]
298         
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     def do_setconv(self,args):
338         '''
339         SETCONV
340         (flatfilts.py)
341         Sets the convfilt configuration variables
342         ------
343         Syntax: setconv variable value
344         '''
345         args=args.split()
346         #FIXME: a general "set dictionary" function has to be built
347         if len(args)==0:
348             print self.convfilt_config
349         else:
350             if not (args[0] in self.convfilt_config.keys()):
351                 print 'This is not an internal convfilt variable!'
352                 print 'Run "setconv" without arguments to see a list of defined variables.'
353                 return
354             
355             if len(args)==1:
356                 print self.convfilt_config[args[0]]
357             elif len(args)>1:
358                 try:
359                     self.convfilt_config[args[0]]=eval(args[1])
360                 except NameError: #we have a string argument
361                     self.convfilt_config[args[0]]=args[1]
362
363
364 #########################
365 #HANDLING OF CONFIGURATION FILE
366 class ConvfiltConfig:
367     '''
368     Handling of convfilt configuration file
369     
370     Mostly based on the simple-yet-useful examples of the Python Library Reference
371     about xml.dom.minidom
372     
373     FIXME: starting to look a mess, should require refactoring
374     '''
375     
376     def __init__(self):
377         self.config={}
378         
379                 
380     def load_config(self, filename):
381         myconfig=file(filename)                    
382         #the following 3 lines are needed to strip newlines. otherwise, since newlines
383         #are XML elements too, the parser would read them (and re-save them, multiplying
384         #newlines...)
385         #yes, I'm an XML n00b
386         the_file=myconfig.read()
387         the_file_lines=the_file.split('\n')
388         the_file=''.join(the_file_lines)
389                        
390         self.config_tree=xml.dom.minidom.parseString(the_file)  
391         
392         def getText(nodelist):
393             #take the text from a nodelist
394             #from Python Library Reference 13.7.2
395             rc = ''
396             for node in nodelist:
397                 if node.nodeType == node.TEXT_NODE:
398                     rc += node.data
399             return rc
400         
401         def handleConfig(config):
402             noiseabsdev_elements=config.getElementsByTagName("noise_absdev")
403             convfilt_elements=config.getElementsByTagName("convfilt")
404             handleAbsdev(noiseabsdev_elements)
405             handleConvfilt(convfilt_elements)
406                         
407         def handleAbsdev(noiseabsdev_elements):
408             for element in noiseabsdev_elements:
409                 for attribute in element.attributes.keys():
410                     self.config[attribute]=element.getAttribute(attribute)
411                     
412         def handleConvfilt(convfilt_elements):
413             for element in convfilt_elements:
414                 for attribute in element.attributes.keys():
415                     self.config[attribute]=element.getAttribute(attribute)
416             
417         handleConfig(self.config_tree)
418         #making items in the dictionary machine-readable
419         for item in self.config.keys():
420             try:
421                 self.config[item]=eval(self.config[item])
422             except NameError: #if it's an unreadable string, keep it as a string
423                 pass
424             
425         return self.config