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 modify it
10 # under the terms of the GNU Lesser General Public License as
11 # published by the Free Software Foundation, either version 3 of the
12 # License, or (at your option) any later version.
14 # Hooke is distributed in the hope that it will be useful, but WITHOUT
15 # ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
16 # or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General
17 # 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 Autopeak carries out force curve fitting with a chosen model:
52 def do_autopeak(self, args):
56 Automatically performs a number of analyses on the peaks of the given curve.
57 Currently it automatically:
58 - fits peaks with WLC or FJC function (depending on how the fit_function variable is set)
59 - measures peak maximum forces with a baseline
60 - measures slope in proximity of peak maximum
61 Requires flatten plotmanipulator , fit.py plugin , flatfilts.py plugin with convfilt
64 autopeak [rebase] [pl=value] [t=value] [noauto] [reclick]
66 rebase : Re-asks baseline interval
68 pl=[value] : Use a fixed persistent length for the fit. If pl is not given,
69 the fit will be a 2-variable
70 fit. DO NOT put spaces between 'pl', '=' and the value.
71 The value must be in meters.
72 Scientific notation like 0.35e-9 is fine.
74 t=[value] : Use a user-defined temperature. The value must be in
75 kelvins; by default it is 293 K.
76 DO NOT put spaces between 't', '=' and the value.
78 noauto : allows for clicking the contact point by
79 hand (otherwise it is automatically estimated) the first time.
80 If subsequent measurements are made, the same contact point
81 clicked the first time is used
83 reclick : redefines by hand the contact point, if noauto has been used before
84 but the user is unsatisfied of the previously choosen contact point.
86 usepoints : fit interval by number of points instead than by nanometers
88 noflatten : does not use the "flatten" plot manipulator
90 When you first issue the command, it will ask for the filename. If you are giving the filename
91 of an existing file, autopeak will resume it and append measurements to it. If you are giving
92 a new filename, it will create the file and append to it until you close Hooke.
95 Useful variables (to set with SET command):
97 fit_function = type of function to use for elasticity. If "wlc" worm-like chain is used, if "fjc" freely jointed
100 temperature= temperature of the system for wlc/fjc fit (in K)
102 auto_slope_span = number of points on which measure the slope, for slope
104 auto_fit_nm = number of nm to fit before the peak maximum, for WLC/FJC (if usepoints false)
105 auto_fit_points = number of points to fit before the peak maximum, for WLC/FJC (if usepoints true)
107 baseline_clicks = -1: no baseline, f=0 at the contact point (whether hand-picked or automatically found)
108 0: automatic baseline
109 1: decide baseline with a single click and length defined in auto_left_baseline
110 2: let user click points of baseline
111 auto_left_baseline = length in nm to use as baseline from the right point (if baseline_clicks=0 , 1)
112 auto_right_baseline = distance in nm of peak-most baseline point from last peak (if baseline_clicks = 0)
114 auto_min_p ; auto_max_p = Minimum and maximum persistence length (if using WLC) or Kuhn length (if using FJC)
115 outside of which the peak is automatically discarded (in nm)
118 #default fit etc. variables
120 T=self.config['temperature']
122 slope_span=int(self.config['auto_slope_span'])
124 rebase=False #if true=we select rebase
125 noflatten=False #if true=we avoid flattening
127 #initialize output data vectors
136 displayed_plot=self._get_displayed_plot(0)
138 #COMMAND LINE PARSING
139 #--Using points instead of nm interval
140 if 'usepoints' in args.split():
141 fit_points=int(self.config['auto_fit_points'])
146 #--Recalculate baseline
147 if 'rebase' in args or (self.basecurrent != self.current.path):
150 if 'noflatten' in args:
153 #--Custom persistent length / custom temperature
154 for arg in args.split():
155 #look for a persistent length argument.
157 pl_expression=arg.split('=')
158 pl_value=float(pl_expression[1]) #actual value
159 #look for a T argument. FIXME: spaces are not allowed between 'pl' and value
160 if ('t=' in arg[0:2]) or ('T=' in arg[0:2]):
161 t_expression=arg.split('=')
162 T=float(t_expression[1])
163 #--Contact point arguments
164 if 'reclick' in args.split():
165 print 'Click contact point'
166 contact_point, contact_point_index = self.pickup_contact_point()
167 elif 'noauto' in args.split():
168 if self.wlccontact_index==None or self.wlccurrent != self.current.path:
169 print 'Click contact point'
170 contact_point , contact_point_index = self.pickup_contact_point()
172 contact_point=self.wlccontact_point
173 contact_point_index=self.wlccontact_index
175 #Automatically find contact point
176 cindex=self.find_contact_point()
177 contact_point=self._clickize(displayed_plot.vectors[1][0], displayed_plot.vectors[1][1], cindex)
178 #--END COMMAND LINE PARSING--
181 peak_location, peak_size = self.find_current_peaks(noflatten)
183 if len(peak_location) == 0:
184 print 'No peaks to fit.'
187 fitplot=copy.deepcopy(displayed_plot)
189 #Pick up force baseline
191 self.basepoints=self.baseline_points(peak_location, displayed_plot)
193 boundaries=[self.basepoints[0].index, self.basepoints[1].index]
195 to_average=displayed_plot.vectors[1][1][boundaries[0]:boundaries[1]] #y points to average
196 avg=np.mean(to_average)
198 clicks=self.config['baseline_clicks']
201 avg=displayed_plot.vectors[1][1][contact_point_index]
203 avg=displayed_plot.vectors[1][1][cindex]
205 for peak in peak_location:
209 fit_points=self.fit_interval_nm(peak, displayed_plot, self.config['auto_fit_nm'], True)
210 peak_point=self._clickize(displayed_plot.vectors[1][0],displayed_plot.vectors[1][1],peak)
211 other_fit_point=self._clickize(displayed_plot.vectors[1][0],displayed_plot.vectors[1][1],peak-fit_points)
214 points=[contact_point, peak_point, other_fit_point]
216 if abs(peak_point.index-other_fit_point.index) < 2:
219 if self.config['fit_function']=='wlc':
221 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)
222 elif self.config['fit_function']=='fjc':
223 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)
225 print 'Unknown fit function'
226 print 'Please set fit_function as wlc or fjc'
231 delta_to_measure=displayed_plot.vectors[1][1][peak-delta_force:peak+delta_force]
232 y=min(delta_to_measure)
233 #save force values (pN)
235 slope=self.linefit_between(peak-slope_span,peak)[0]
238 #check fitted data and, if right, add peak to the measurement
239 #FIXME: code duplication
240 if len(params)==1: #if we did choose 1-value fit
241 p_lengths.append(pl_value)
242 c_lengths.append(params[0]*(1.0e+9))
243 sigma_p_lengths.append(0)
244 sigma_c_lengths.append(fit_errors[0]*(1.0e+9))
245 forces.append(abs(y-avg)*(1.0e+12))
247 #Add WLC fit lines to plot
248 fitplot.add_set(xfit,yfit)
249 if len(fitplot.styles)==0:
253 fitplot.styles.append(None)
254 fitplot.colors.append(None)
256 p_leng=params[1]*(1.0e+9)
257 #check if persistent length makes sense. otherwise, discard peak.
258 if p_leng>self.config['auto_min_p'] and p_leng<self.config['auto_max_p']:
259 p_lengths.append(p_leng)
260 c_lengths.append(params[0]*(1.0e+9))
261 sigma_c_lengths.append(fit_errors[0]*(1.0e+9))
262 sigma_p_lengths.append(fit_errors[1]*(1.0e+9))
263 forces.append(abs(y-avg)*(1.0e+12))
266 #Add WLC fit lines to plot
267 fitplot.add_set(xfit,yfit)
268 if len(fitplot.styles)==0:
272 fitplot.styles.append(None)
273 fitplot.colors.append(None)
278 #add basepoints to fitplot
279 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]])
280 fitplot.styles.append('scatter')
281 fitplot.colors.append(None)
283 #Show wlc fits and peak locations
284 self._send_plot([fitplot])
286 print 'Using fit function: ',self.config['fit_function']
287 print 'Measurements for all peaks detected:'
288 print 'contour (nm)', c_lengths
289 print 'sigma contour (nm)',sigma_c_lengths
290 print 'p (nm)',p_lengths
291 print 'sigma p (nm)',sigma_p_lengths
292 print 'forces (pN)',forces
293 print 'slopes (N/m)',slopes
296 while controller==False:
297 #Ask the user what peaks to ignore from analysis.
298 print 'Peaks to ignore (0,1...n from contact point,return to take all)'
299 print 'N to discard measurement'
300 exclude_raw=raw_input('Input:')
306 if not exclude_raw=='':
307 exclude=exclude_raw.split(',')
309 exclude=[int(item) for item in exclude]
315 sigma_c_lengths[i]=None
316 sigma_p_lengths[i]=None
322 #Clean data vectors from ignored peaks
323 #FIXME:code duplication
324 c_lengths=[item for item in c_lengths if item != None]
325 p_lengths=[item for item in p_lengths if item != None]
326 forces=[item for item in forces if item != None]
327 slopes=[item for item in slopes if item != None]
328 sigma_c_lengths=[item for item in sigma_c_lengths if item != None]
329 sigma_p_lengths=[item for item in sigma_p_lengths if item != None]
331 print 'Measurements for chosen peaks:'
332 print 'contour (nm)',c_lengths
333 print 'sigma contour (nm)',sigma_c_lengths
334 print 'p (nm)',p_lengths
335 print 'sigma p (nm)',sigma_p_lengths
336 print 'forces (pN)',forces
337 print 'slopes (N/m)',slopes
340 if self.autofile=='':
341 self.autofile=raw_input('Autopeak filename? (return to ignore) ')
342 if self.autofile=='':
346 if not os.path.exists(self.autofile):
347 f=open(self.autofile,'w+')
348 f.write('Analysis started '+time.asctime()+'\n')
349 f.write('----------------------------------------\n')
350 f.write('; Contour length (nm) ; Persistence length (nm) ; Max.Force (pN) ; Slope (N/m) ; Sigma contour (nm) ; Sigma persistence (nm)\n')
354 f=open(self.autofile,'a+')
356 f.write(self.current.path+'\n')
357 for i in range(len(c_lengths)):
358 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')
361 self.do_note('autopeak')