1 # Copyright (C) 2008-2010 Alberto Gomez-Casado
4 # Massimo Sandal <devicerandom@gmail.com>
5 # W. Trevor King <wking@drexel.edu>
7 # This file is part of Hooke.
9 # Hooke is free software: you can redistribute it and/or
10 # modify it under the terms of the GNU Lesser General Public
11 # License as published by the Free Software Foundation, either
12 # version 3 of the License, or (at your option) any later version.
14 # Hooke is distributed in the hope that it will be useful,
15 # but WITHOUT ANY WARRANTY; without even the implied warranty of
16 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 # GNU Lesser General Public License for more details.
19 # You should have received a copy of the GNU Lesser General Public
20 # License along with Hooke. If not, see
21 # <http://www.gnu.org/licenses/>.
23 """The autopeak module provides :class:`Autopeak`, a
24 :class:`hooke.plugin.Plugin` for automatically extracting force peaks
25 (unfolding events) from force curves.
28 from hooke.libhooke import WX_GOOD
31 wxversion.select(WX_GOOD)
32 from wx import PostEvent
40 warnings.simplefilter('ignore',np.RankWarning)
42 #from .. import ui.gui.results as results
45 class autopeakCommands(object):
47 TODO: autopeak docstring.
49 def do_autopeak(self, args):
53 Automatically performs a number of analyses on the peaks of the given curve.
54 Currently it automatically:
55 - fits peaks with WLC or FJC function (depending on how the fit_function variable is set)
56 - measures peak maximum forces with a baseline
57 - measures slope in proximity of peak maximum
58 Requires flatten plotmanipulator , fit.py plugin , flatfilts.py plugin with convfilt
61 autopeak [rebase] [pl=value] [t=value] [noauto] [reclick]
63 rebase : Re-asks baseline interval
65 pl=[value] : Use a fixed persistent length for the fit. If pl is not given,
66 the fit will be a 2-variable
67 fit. DO NOT put spaces between 'pl', '=' and the value.
68 The value must be in meters.
69 Scientific notation like 0.35e-9 is fine.
71 t=[value] : Use a user-defined temperature. The value must be in
72 kelvins; by default it is 293 K.
73 DO NOT put spaces between 't', '=' and the value.
75 noauto : allows for clicking the contact point by
76 hand (otherwise it is automatically estimated) the first time.
77 If subsequent measurements are made, the same contact point
78 clicked the first time is used
80 reclick : redefines by hand the contact point, if noauto has been used before
81 but the user is unsatisfied of the previously choosen contact point.
83 usepoints : fit interval by number of points instead than by nanometers
85 noflatten : does not use the "flatten" plot manipulator
87 When you first issue the command, it will ask for the filename. If you are giving the filename
88 of an existing file, autopeak will resume it and append measurements to it. If you are giving
89 a new filename, it will create the file and append to it until you close Hooke.
92 Useful variables (to set with SET command):
94 fit_function = type of function to use for elasticity. If "wlc" worm-like chain is used, if "fjc" freely jointed
97 temperature= temperature of the system for wlc/fjc fit (in K)
99 auto_slope_span = number of points on which measure the slope, for slope
101 auto_fit_nm = number of nm to fit before the peak maximum, for WLC/FJC (if usepoints false)
102 auto_fit_points = number of points to fit before the peak maximum, for WLC/FJC (if usepoints true)
104 baseline_clicks = -1: no baseline, f=0 at the contact point (whether hand-picked or automatically found)
105 0: automatic baseline
106 1: decide baseline with a single click and length defined in auto_left_baseline
107 2: let user click points of baseline
108 auto_left_baseline = length in nm to use as baseline from the right point (if baseline_clicks=0 , 1)
109 auto_right_baseline = distance in nm of peak-most baseline point from last peak (if baseline_clicks = 0)
111 auto_min_p ; auto_max_p = Minimum and maximum persistence length (if using WLC) or Kuhn length (if using FJC)
112 outside of which the peak is automatically discarded (in nm)
115 #default fit etc. variables
117 T=self.config['temperature']
119 slope_span=int(self.config['auto_slope_span'])
121 rebase=False #if true=we select rebase
122 noflatten=False #if true=we avoid flattening
124 #initialize output data vectors
133 displayed_plot=self._get_displayed_plot(0)
135 #COMMAND LINE PARSING
136 #--Using points instead of nm interval
137 if 'usepoints' in args.split():
138 fit_points=int(self.config['auto_fit_points'])
143 #--Recalculate baseline
144 if 'rebase' in args or (self.basecurrent != self.current.path):
147 if 'noflatten' in args:
150 #--Custom persistent length / custom temperature
151 for arg in args.split():
152 #look for a persistent length argument.
154 pl_expression=arg.split('=')
155 pl_value=float(pl_expression[1]) #actual value
156 #look for a T argument. FIXME: spaces are not allowed between 'pl' and value
157 if ('t=' in arg[0:2]) or ('T=' in arg[0:2]):
158 t_expression=arg.split('=')
159 T=float(t_expression[1])
160 #--Contact point arguments
161 if 'reclick' in args.split():
162 print 'Click contact point'
163 contact_point, contact_point_index = self.pickup_contact_point()
164 elif 'noauto' in args.split():
165 if self.wlccontact_index==None or self.wlccurrent != self.current.path:
166 print 'Click contact point'
167 contact_point , contact_point_index = self.pickup_contact_point()
169 contact_point=self.wlccontact_point
170 contact_point_index=self.wlccontact_index
172 #Automatically find contact point
173 cindex=self.find_contact_point()
174 contact_point=self._clickize(displayed_plot.vectors[1][0], displayed_plot.vectors[1][1], cindex)
175 #--END COMMAND LINE PARSING--
178 peak_location, peak_size = self.find_current_peaks(noflatten)
180 if len(peak_location) == 0:
181 print 'No peaks to fit.'
184 fitplot=copy.deepcopy(displayed_plot)
186 #Pick up force baseline
188 self.basepoints=self.baseline_points(peak_location, displayed_plot)
190 boundaries=[self.basepoints[0].index, self.basepoints[1].index]
192 to_average=displayed_plot.vectors[1][1][boundaries[0]:boundaries[1]] #y points to average
193 avg=np.mean(to_average)
195 clicks=self.config['baseline_clicks']
198 avg=displayed_plot.vectors[1][1][contact_point_index]
200 avg=displayed_plot.vectors[1][1][cindex]
202 for peak in peak_location:
206 fit_points=self.fit_interval_nm(peak, displayed_plot, self.config['auto_fit_nm'], True)
207 peak_point=self._clickize(displayed_plot.vectors[1][0],displayed_plot.vectors[1][1],peak)
208 other_fit_point=self._clickize(displayed_plot.vectors[1][0],displayed_plot.vectors[1][1],peak-fit_points)
211 points=[contact_point, peak_point, other_fit_point]
213 if abs(peak_point.index-other_fit_point.index) < 2:
216 if self.config['fit_function']=='wlc':
218 params, yfit, xfit, fit_errors = self.wlc_fit(points, displayed_plot.vectors[1][0], displayed_plot.vectors[1][1], pl_value, T, return_errors=True)
219 elif self.config['fit_function']=='fjc':
220 params, yfit, xfit, fit_errors = self.fjc_fit(points, displayed_plot.vectors[1][0], displayed_plot.vectors[1][1], pl_value, T, return_errors=True)
222 print 'Unknown fit function'
223 print 'Please set fit_function as wlc or fjc'
228 delta_to_measure=displayed_plot.vectors[1][1][peak-delta_force:peak+delta_force]
229 y=min(delta_to_measure)
230 #save force values (pN)
232 slope=self.linefit_between(peak-slope_span,peak)[0]
235 #check fitted data and, if right, add peak to the measurement
236 #FIXME: code duplication
237 if len(params)==1: #if we did choose 1-value fit
238 p_lengths.append(pl_value)
239 c_lengths.append(params[0]*(1.0e+9))
240 sigma_p_lengths.append(0)
241 sigma_c_lengths.append(fit_errors[0]*(1.0e+9))
242 forces.append(abs(y-avg)*(1.0e+12))
244 #Add WLC fit lines to plot
245 fitplot.add_set(xfit,yfit)
246 if len(fitplot.styles)==0:
250 fitplot.styles.append(None)
251 fitplot.colors.append(None)
253 p_leng=params[1]*(1.0e+9)
254 #check if persistent length makes sense. otherwise, discard peak.
255 if p_leng>self.config['auto_min_p'] and p_leng<self.config['auto_max_p']:
256 p_lengths.append(p_leng)
257 c_lengths.append(params[0]*(1.0e+9))
258 sigma_c_lengths.append(fit_errors[0]*(1.0e+9))
259 sigma_p_lengths.append(fit_errors[1]*(1.0e+9))
260 forces.append(abs(y-avg)*(1.0e+12))
263 #Add WLC fit lines to plot
264 fitplot.add_set(xfit,yfit)
265 if len(fitplot.styles)==0:
269 fitplot.styles.append(None)
270 fitplot.colors.append(None)
275 #add basepoints to fitplot
276 fitplot.add_set([self.basepoints[0].graph_coords[0],self.basepoints[1].graph_coords[0]],[self.basepoints[0].graph_coords[1],self.basepoints[1].graph_coords[1]])
277 fitplot.styles.append('scatter')
278 fitplot.colors.append(None)
280 #Show wlc fits and peak locations
281 self._send_plot([fitplot])
283 print 'Using fit function: ',self.config['fit_function']
284 print 'Measurements for all peaks detected:'
285 print 'contour (nm)', c_lengths
286 print 'sigma contour (nm)',sigma_c_lengths
287 print 'p (nm)',p_lengths
288 print 'sigma p (nm)',sigma_p_lengths
289 print 'forces (pN)',forces
290 print 'slopes (N/m)',slopes
293 while controller==False:
294 #Ask the user what peaks to ignore from analysis.
295 print 'Peaks to ignore (0,1...n from contact point,return to take all)'
296 print 'N to discard measurement'
297 exclude_raw=raw_input('Input:')
303 if not exclude_raw=='':
304 exclude=exclude_raw.split(',')
306 exclude=[int(item) for item in exclude]
312 sigma_c_lengths[i]=None
313 sigma_p_lengths[i]=None
319 #Clean data vectors from ignored peaks
320 #FIXME:code duplication
321 c_lengths=[item for item in c_lengths if item != None]
322 p_lengths=[item for item in p_lengths if item != None]
323 forces=[item for item in forces if item != None]
324 slopes=[item for item in slopes if item != None]
325 sigma_c_lengths=[item for item in sigma_c_lengths if item != None]
326 sigma_p_lengths=[item for item in sigma_p_lengths if item != None]
328 print 'Measurements for chosen peaks:'
329 print 'contour (nm)',c_lengths
330 print 'sigma contour (nm)',sigma_c_lengths
331 print 'p (nm)',p_lengths
332 print 'sigma p (nm)',sigma_p_lengths
333 print 'forces (pN)',forces
334 print 'slopes (N/m)',slopes
337 if self.autofile=='':
338 self.autofile=raw_input('Autopeak filename? (return to ignore) ')
339 if self.autofile=='':
343 if not os.path.exists(self.autofile):
344 f=open(self.autofile,'w+')
345 f.write('Analysis started '+time.asctime()+'\n')
346 f.write('----------------------------------------\n')
347 f.write('; Contour length (nm) ; Persistence length (nm) ; Max.Force (pN) ; Slope (N/m) ; Sigma contour (nm) ; Sigma persistence (nm)\n')
351 f=open(self.autofile,'a+')
353 f.write(self.current.path+'\n')
354 for i in range(len(c_lengths)):
355 f.write(' ; '+str(c_lengths[i])+' ; '+str(p_lengths[i])+' ; '+str(forces[i])+' ; '+str(slopes[i])+' ; '+str(sigma_c_lengths[i])+' ; '+str(sigma_p_lengths[i])+'\n')
358 self.do_note('autopeak')