1 # Copyright (C) 2008-2010 Alberto Gomez-Casado
3 # Massimo Sandal <devicerandom@gmail.com>
4 # W. Trevor King <wking@drexel.edu>
6 # This file is part of Hooke.
8 # Hooke is free software: you can redistribute it and/or
9 # modify it under the terms of the GNU Lesser General Public
10 # License as published by the Free Software Foundation, either
11 # version 3 of the License, or (at your option) any later version.
13 # Hooke is distributed in the hope that it will be useful,
14 # but WITHOUT ANY WARRANTY; without even the implied warranty of
15 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 # GNU Lesser General Public License for more details.
18 # You should have received a copy of the GNU Lesser General Public
19 # License along with Hooke. If not, see
20 # <http://www.gnu.org/licenses/>.
22 """Force spectroscopy curves filtering of flat curves
24 Other plugin dependencies:
25 procplots.py (plot processing plugin)
28 from hooke.libhooke import WX_GOOD
31 wxversion.select(WX_GOOD)
32 import xml.dom.minidom
36 from numpy import diff
39 from .. import libpeakspot as lps
40 from .. import curve as lhc
43 class flatfiltsCommands(object):
46 #configurate convfilt variables
47 convfilt_configurator=ConvfiltConfig()
48 self.convfilt_config=convfilt_configurator.load_config('convfilt.conf')
50 def do_flatfilt(self,args):
54 Filters out flat (featureless) curves of the current playlist,
55 creating a playlist containing only the curves with potential
59 flatfilt [min_npks min_deviation]
61 min_npks = minmum number of points over the deviation
64 min_deviation = minimum signal/noise ratio
67 If called without arguments, it uses default values, that
68 should work most of the times.
77 min_deviation=int(args[1])
81 print 'Processing playlist...'
85 for item in self.current_list:
89 notflat=self.has_features(item, median_filter, min_npks, min_deviation)
90 print 'Curve',item.path, 'is',c,'of',len(self.current_list),': features are ',notflat
93 print 'Curve',item.path, 'is',c,'of',len(self.current_list),': cannot be filtered. Probably unable to retrieve force data from corrupt file.'
97 item.curve=None #empty the item object, to further avoid memory leak
98 notflat_list.append(item)
100 if len(notflat_list)==0:
101 print 'Found nothing interesting. Check your playlist, could be a bug or criteria could be too much stringent'
104 print 'Found ',len(notflat_list),' potentially interesting curves'
105 print 'Regenerating playlist...'
107 self.current_list=notflat_list
108 self.current=self.current_list[self.pointer]
111 def has_features(self,item,median_filter,min_npks,min_deviation):
113 decides if a curve is flat enough to be rejected from analysis: it sees if there
114 are at least min_npks points that are higher than min_deviation times the absolute value
117 Algorithm original idea by Francesco Musiani, with my tweaks and corrections.
121 item.identify(self.drivers)
122 #we assume the first is the plot with the force curve
123 #do the median to better resolve features from noise
124 flat_plot=self.plotmanip_median(item.curve.default_plots()[0], item, customvalue=median_filter)
125 flat_vects=flat_plot.vectors
126 item.curve.close_all()
127 #needed to avoid *big* memory leaks!
131 #absolute value of derivate
132 yretdiff=diff(flat_vects[1][1])
133 yretdiff=[abs(value) for value in yretdiff]
134 #average of derivate values
135 diffmean=numpy.mean(yretdiff)
139 for value in yretdiff:
140 if value/diffmean > min_deviation:
148 del flat_plot, flat_vects, yretdiff
152 ################################################################
153 #-----CONVFILT-------------------------------------------------
154 #-----Convolution-based peak recognition and filtering.
155 #Requires the libpeakspot.py library
157 def has_peaks(self, plot, abs_devs=None, maxpeak=True, window=10):
159 Finds peak position in a force curve.
160 FIXME: should be moved in libpeakspot.py
163 abs_devs=self.convfilt_config['mindeviation']
166 xret=plot.vectors[1][0]
167 yret=plot.vectors[1][1]
168 #Calculate convolution.
169 convoluted=lps.conv_dx(yret, self.convfilt_config['convolution'])
171 #surely cut everything before the contact point
172 cut_index=self.find_contact_point(plot)
173 #cut even more, before the blind window
174 start_x=xret[cut_index]
176 for value in xret[cut_index:]:
177 if abs((value) - (start_x)) > self.convfilt_config['blindwindow']*(10**-9):
180 cut_index+=blind_index
181 #do the dirty convolution-peak finding stuff
182 noise_level=lps.noise_absdev(convoluted[cut_index:], self.convfilt_config['positive'], self.convfilt_config['maxcut'], self.convfilt_config['stable'])
183 above=lps.abovenoise(convoluted,noise_level,cut_index,abs_devs)
184 peak_location,peak_size=lps.find_peaks(above,seedouble=self.convfilt_config['seedouble'])
186 #take the minimum or the maximum of a peak
187 for i in range(len(peak_location)):
188 peak=peak_location[i]
189 valpk=min(yret[peak-window:peak+window]) #maximum in force (near the unfolding point)
190 index_pk=yret[peak-window:peak+window].index(valpk)+(peak-window)
193 valpk=max(yret[peak:peak+window]) #minimum in force, near the baseline
194 index_pk=yret[peak:peak+window].index(valpk)+(peak)
196 # Let's explain that for the minimum. Immaging that we know that there is a peak at position/region 100 and you have found its y-value,
197 # Now you look in the array, from 100-10 to 100+10 (if the window is 10).
198 # This "100-10 to 100+10" is substancially a new array with its index. In this array you have 20
199 # elements, so the index of your y-value will be 10.
200 # Now to find the index in the TOTAL array you have to add the "position" of the "region" (that in this case
201 # correspond to 100) and also substract the window size ---> (+100-10)
203 peak_location[i]=index_pk
205 return peak_location,peak_size
208 def exec_has_peaks(self,item,abs_devs):
210 encapsulates has_peaks for the purpose of correctly treating the curve objects in the convfilt loop,
211 to avoid memory leaks
213 item.identify(self.drivers)
214 #we assume the first is the plot with the force curve
215 plot=item.curve.default_plots()[0]
217 if 'flatten' in self.config['plotmanips']:
218 #If flatten is present, use it for better recognition of peaks...
219 flatten=self._find_plotmanip('flatten') #extract flatten plot manipulator
220 plot=flatten(plot, item, customvalue=1)
222 peak_location,peak_size=self.has_peaks(plot,abs_devs)
223 #close all open files
224 item.curve.close_all()
225 #needed to avoid *big* memory leaks!
228 return peak_location, peak_size
230 #------------------------
231 #------commands----------
232 #------------------------
233 def do_peaks(self,args):
237 Test command for convolution filter / test.
239 Syntax: peaks [deviations]
240 absolute deviation = number of times the convolution signal is above the noise absolute deviation.
244 args=self.convfilt_config['mindeviation']
249 print 'Wrong argument, using config value'
250 abs_devs=float(self.convfilt_config['mindeviation'])
252 defplots=self.current.curve.default_plots()[0] #we need the raw, uncorrected plots
254 if 'flatten' in self.config['plotmanips']:
255 flatten=self._find_plotmanip('flatten') #extract flatten plot manipulator
256 defplots=flatten(defplots, self.current)
258 print 'You have the flatten plot manipulator not loaded. Enabling it could give you better results.'
260 peak_location,peak_size=self.has_peaks(defplots,abs_devs)
261 print 'Found '+str(len(peak_location))+' peaks.'
262 to_dump='peaks '+self.current.path+' '+str(len(peak_location))
263 self.outlet.push(to_dump)
266 #if no peaks, we have nothing to plot. exit.
267 if len(peak_location)==0:
270 #otherwise, we plot the peak locations.
271 xplotted_ret=self.plots[0].vectors[1][0]
272 yplotted_ret=self.plots[0].vectors[1][1]
273 xgood=[xplotted_ret[index] for index in peak_location]
274 ygood=[yplotted_ret[index] for index in peak_location]
276 recplot=self._get_displayed_plot()
277 recplot.vectors.append([xgood,ygood])
278 if recplot.styles==[]:
279 recplot.styles=[None,None,'scatter']
280 recplot.colors=[None,None,None]
282 recplot.styles+=['scatter']
283 recplot.colors+=[None]
285 self._send_plot([recplot])
287 def do_convfilt(self,args):
291 Filters out flat (featureless) curves of the current playlist,
292 creating a playlist containing only the curves with potential
296 convfilt [min_npks min_deviation]
298 min_npks = minmum number of peaks
299 (to set the default, see convfilt.conf file; CONVCONF and SETCONF commands)
301 min_deviation = minimum signal/noise ratio *in the convolution*
302 (to set the default, see convfilt.conf file; CONVCONF and SETCONF commands)
304 If called without arguments, it uses default values.
307 min_npks=self.convfilt_config['minpeaks']
308 min_deviation=self.convfilt_config['mindeviation']
312 min_npks=int(args[0])
313 min_deviation=int(args[1])
317 print 'Processing playlist...'
318 print '(Please wait)'
322 for item in self.current_list:
326 peak_location,peak_size=self.exec_has_peaks(item,min_deviation)
327 if len(peak_location)>=min_npks:
331 print 'Curve',item.path, 'is',c,'of',len(self.current_list),': found '+str(len(peak_location))+' peaks.'+isok
333 peak_location,peak_size=[],[]
334 print 'Curve',item.path, 'is',c,'of',len(self.current_list),': cannot be filtered. Probably unable to retrieve force data from corrupt file.'
336 if len(peak_location)>=min_npks:
337 item.peak_location=peak_location
338 item.peak_size=peak_size
339 item.curve=None #empty the item object, to further avoid memory leak
340 notflat_list.append(item)
342 #Warn that no flattening had been done.
343 if not ('flatten' in self.config['plotmanips']):
344 print 'Flatten manipulator was not found. Processing was done without flattening.'
345 print 'Try to enable it in your configuration file for better results.'
347 if len(notflat_list)==0:
348 print 'Found nothing interesting. Check your playlist, could be a bug or criteria could be too much stringent'
351 print 'Found ',len(notflat_list),' potentially interesting curves'
352 print 'Regenerating playlist...'
354 self.current_list=notflat_list
355 self.current=self.current_list[self.pointer]
359 def do_setconv(self,args):
363 Sets the convfilt configuration variables
365 Syntax: setconv variable value
368 #FIXME: a general "set dictionary" function has to be built
370 print self.convfilt_config
372 if not (args[0] in self.convfilt_config.keys()):
373 print 'This is not an internal convfilt variable!'
374 print 'Run "setconv" without arguments to see a list of defined variables.'
378 print self.convfilt_config[args[0]]
381 self.convfilt_config[args[0]]=eval(args[1])
382 except NameError: #we have a string argument
383 self.convfilt_config[args[0]]=args[1]
386 #########################
387 #HANDLING OF CONFIGURATION FILE
388 class ConvfiltConfig(object):
390 Handling of convfilt configuration file
392 Mostly based on the simple-yet-useful examples of the Python Library Reference
393 about xml.dom.minidom
395 FIXME: starting to look a mess, should require refactoring
402 def load_config(self, filename):
403 myconfig=file(filename)
404 #the following 3 lines are needed to strip newlines. otherwise, since newlines
405 #are XML elements too, the parser would read them (and re-save them, multiplying
407 #yes, I'm an XML n00b
408 the_file=myconfig.read()
409 the_file_lines=the_file.split('\n')
410 the_file=''.join(the_file_lines)
412 self.config_tree=xml.dom.minidom.parseString(the_file)
414 def getText(nodelist):
415 #take the text from a nodelist
416 #from Python Library Reference 13.7.2
418 for node in nodelist:
419 if node.nodeType == node.TEXT_NODE:
423 def handleConfig(config):
424 noiseabsdev_elements=config.getElementsByTagName("noise_absdev")
425 convfilt_elements=config.getElementsByTagName("convfilt")
426 handleAbsdev(noiseabsdev_elements)
427 handleConvfilt(convfilt_elements)
429 def handleAbsdev(noiseabsdev_elements):
430 for element in noiseabsdev_elements:
431 for attribute in element.attributes.keys():
432 self.config[attribute]=element.getAttribute(attribute)
434 def handleConvfilt(convfilt_elements):
435 for element in convfilt_elements:
436 for attribute in element.attributes.keys():
437 self.config[attribute]=element.getAttribute(attribute)
439 handleConfig(self.config_tree)
440 #making items in the dictionary machine-readable
441 for item in self.config.keys():
443 self.config[item]=eval(self.config[item])
444 except NameError: #if it's an unreadable string, keep it as a string