3 """The autopeak module provides :class:`Autopeak`, a
4 :class:`hooke.plugin.Plugin` for automatically extracting force peaks
5 (unfolding events) from force curves.
8 from hooke.libhooke import WX_GOOD
11 wxversion.select(WX_GOOD)
12 from wx import PostEvent
20 warnings.simplefilter('ignore',np.RankWarning)
22 #from .. import ui.gui.results as results
25 class autopeakCommands(object):
27 TODO: autopeak docstring.
29 def do_autopeak(self, args):
33 Automatically performs a number of analyses on the peaks of the given curve.
34 Currently it automatically:
35 - fits peaks with WLC or FJC function (depending on how the fit_function variable is set)
36 - measures peak maximum forces with a baseline
37 - measures slope in proximity of peak maximum
38 Requires flatten plotmanipulator , fit.py plugin , flatfilts.py plugin with convfilt
41 autopeak [rebase] [pl=value] [t=value] [noauto] [reclick]
43 rebase : Re-asks baseline interval
45 pl=[value] : Use a fixed persistent length for the fit. If pl is not given,
46 the fit will be a 2-variable
47 fit. DO NOT put spaces between 'pl', '=' and the value.
48 The value must be in meters.
49 Scientific notation like 0.35e-9 is fine.
51 t=[value] : Use a user-defined temperature. The value must be in
52 kelvins; by default it is 293 K.
53 DO NOT put spaces between 't', '=' and the value.
55 noauto : allows for clicking the contact point by
56 hand (otherwise it is automatically estimated) the first time.
57 If subsequent measurements are made, the same contact point
58 clicked the first time is used
60 reclick : redefines by hand the contact point, if noauto has been used before
61 but the user is unsatisfied of the previously choosen contact point.
63 usepoints : fit interval by number of points instead than by nanometers
65 noflatten : does not use the "flatten" plot manipulator
67 When you first issue the command, it will ask for the filename. If you are giving the filename
68 of an existing file, autopeak will resume it and append measurements to it. If you are giving
69 a new filename, it will create the file and append to it until you close Hooke.
72 Useful variables (to set with SET command):
74 fit_function = type of function to use for elasticity. If "wlc" worm-like chain is used, if "fjc" freely jointed
77 temperature= temperature of the system for wlc/fjc fit (in K)
79 auto_slope_span = number of points on which measure the slope, for slope
81 auto_fit_nm = number of nm to fit before the peak maximum, for WLC/FJC (if usepoints false)
82 auto_fit_points = number of points to fit before the peak maximum, for WLC/FJC (if usepoints true)
84 baseline_clicks = -1: no baseline, f=0 at the contact point (whether hand-picked or automatically found)
86 1: decide baseline with a single click and length defined in auto_left_baseline
87 2: let user click points of baseline
88 auto_left_baseline = length in nm to use as baseline from the right point (if baseline_clicks=0 , 1)
89 auto_right_baseline = distance in nm of peak-most baseline point from last peak (if baseline_clicks = 0)
91 auto_min_p ; auto_max_p = Minimum and maximum persistence length (if using WLC) or Kuhn length (if using FJC)
92 outside of which the peak is automatically discarded (in nm)
95 #default fit etc. variables
97 T=self.config['temperature']
99 slope_span=int(self.config['auto_slope_span'])
101 rebase=False #if true=we select rebase
102 noflatten=False #if true=we avoid flattening
104 #initialize output data vectors
113 displayed_plot=self._get_displayed_plot(0)
115 #COMMAND LINE PARSING
116 #--Using points instead of nm interval
117 if 'usepoints' in args.split():
118 fit_points=int(self.config['auto_fit_points'])
123 #--Recalculate baseline
124 if 'rebase' in args or (self.basecurrent != self.current.path):
127 if 'noflatten' in args:
130 #--Custom persistent length / custom temperature
131 for arg in args.split():
132 #look for a persistent length argument.
134 pl_expression=arg.split('=')
135 pl_value=float(pl_expression[1]) #actual value
136 #look for a T argument. FIXME: spaces are not allowed between 'pl' and value
137 if ('t=' in arg[0:2]) or ('T=' in arg[0:2]):
138 t_expression=arg.split('=')
139 T=float(t_expression[1])
140 #--Contact point arguments
141 if 'reclick' in args.split():
142 print 'Click contact point'
143 contact_point, contact_point_index = self.pickup_contact_point()
144 elif 'noauto' in args.split():
145 if self.wlccontact_index==None or self.wlccurrent != self.current.path:
146 print 'Click contact point'
147 contact_point , contact_point_index = self.pickup_contact_point()
149 contact_point=self.wlccontact_point
150 contact_point_index=self.wlccontact_index
152 #Automatically find contact point
153 cindex=self.find_contact_point()
154 contact_point=self._clickize(displayed_plot.vectors[1][0], displayed_plot.vectors[1][1], cindex)
155 #--END COMMAND LINE PARSING--
158 peak_location, peak_size = self.find_current_peaks(noflatten)
160 if len(peak_location) == 0:
161 print 'No peaks to fit.'
164 fitplot=copy.deepcopy(displayed_plot)
166 #Pick up force baseline
168 self.basepoints=self.baseline_points(peak_location, displayed_plot)
170 boundaries=[self.basepoints[0].index, self.basepoints[1].index]
172 to_average=displayed_plot.vectors[1][1][boundaries[0]:boundaries[1]] #y points to average
173 avg=np.mean(to_average)
175 clicks=self.config['baseline_clicks']
178 avg=displayed_plot.vectors[1][1][contact_point_index]
180 avg=displayed_plot.vectors[1][1][cindex]
182 for peak in peak_location:
186 fit_points=self.fit_interval_nm(peak, displayed_plot, self.config['auto_fit_nm'], True)
187 peak_point=self._clickize(displayed_plot.vectors[1][0],displayed_plot.vectors[1][1],peak)
188 other_fit_point=self._clickize(displayed_plot.vectors[1][0],displayed_plot.vectors[1][1],peak-fit_points)
191 points=[contact_point, peak_point, other_fit_point]
193 if abs(peak_point.index-other_fit_point.index) < 2:
196 if self.config['fit_function']=='wlc':
198 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)
199 elif self.config['fit_function']=='fjc':
200 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)
202 print 'Unknown fit function'
203 print 'Please set fit_function as wlc or fjc'
208 delta_to_measure=displayed_plot.vectors[1][1][peak-delta_force:peak+delta_force]
209 y=min(delta_to_measure)
210 #save force values (pN)
212 slope=self.linefit_between(peak-slope_span,peak)[0]
215 #check fitted data and, if right, add peak to the measurement
216 #FIXME: code duplication
217 if len(params)==1: #if we did choose 1-value fit
218 p_lengths.append(pl_value)
219 c_lengths.append(params[0]*(1.0e+9))
220 sigma_p_lengths.append(0)
221 sigma_c_lengths.append(fit_errors[0]*(1.0e+9))
222 forces.append(abs(y-avg)*(1.0e+12))
224 #Add WLC fit lines to plot
225 fitplot.add_set(xfit,yfit)
226 if len(fitplot.styles)==0:
230 fitplot.styles.append(None)
231 fitplot.colors.append(None)
233 p_leng=params[1]*(1.0e+9)
234 #check if persistent length makes sense. otherwise, discard peak.
235 if p_leng>self.config['auto_min_p'] and p_leng<self.config['auto_max_p']:
236 p_lengths.append(p_leng)
237 c_lengths.append(params[0]*(1.0e+9))
238 sigma_c_lengths.append(fit_errors[0]*(1.0e+9))
239 sigma_p_lengths.append(fit_errors[1]*(1.0e+9))
240 forces.append(abs(y-avg)*(1.0e+12))
243 #Add WLC fit lines to plot
244 fitplot.add_set(xfit,yfit)
245 if len(fitplot.styles)==0:
249 fitplot.styles.append(None)
250 fitplot.colors.append(None)
255 #add basepoints to fitplot
256 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]])
257 fitplot.styles.append('scatter')
258 fitplot.colors.append(None)
260 #Show wlc fits and peak locations
261 self._send_plot([fitplot])
263 print 'Using fit function: ',self.config['fit_function']
264 print 'Measurements for all peaks detected:'
265 print 'contour (nm)', c_lengths
266 print 'sigma contour (nm)',sigma_c_lengths
267 print 'p (nm)',p_lengths
268 print 'sigma p (nm)',sigma_p_lengths
269 print 'forces (pN)',forces
270 print 'slopes (N/m)',slopes
273 while controller==False:
274 #Ask the user what peaks to ignore from analysis.
275 print 'Peaks to ignore (0,1...n from contact point,return to take all)'
276 print 'N to discard measurement'
277 exclude_raw=raw_input('Input:')
283 if not exclude_raw=='':
284 exclude=exclude_raw.split(',')
286 exclude=[int(item) for item in exclude]
292 sigma_c_lengths[i]=None
293 sigma_p_lengths[i]=None
299 #Clean data vectors from ignored peaks
300 #FIXME:code duplication
301 c_lengths=[item for item in c_lengths if item != None]
302 p_lengths=[item for item in p_lengths if item != None]
303 forces=[item for item in forces if item != None]
304 slopes=[item for item in slopes if item != None]
305 sigma_c_lengths=[item for item in sigma_c_lengths if item != None]
306 sigma_p_lengths=[item for item in sigma_p_lengths if item != None]
308 print 'Measurements for chosen peaks:'
309 print 'contour (nm)',c_lengths
310 print 'sigma contour (nm)',sigma_c_lengths
311 print 'p (nm)',p_lengths
312 print 'sigma p (nm)',sigma_p_lengths
313 print 'forces (pN)',forces
314 print 'slopes (N/m)',slopes
317 if self.autofile=='':
318 self.autofile=raw_input('Autopeak filename? (return to ignore) ')
319 if self.autofile=='':
323 if not os.path.exists(self.autofile):
324 f=open(self.autofile,'w+')
325 f.write('Analysis started '+time.asctime()+'\n')
326 f.write('----------------------------------------\n')
327 f.write('; Contour length (nm) ; Persistence length (nm) ; Max.Force (pN) ; Slope (N/m) ; Sigma contour (nm) ; Sigma persistence (nm)\n')
331 f=open(self.autofile,'a+')
333 f.write(self.current.path+'\n')
334 for i in range(len(c_lengths)):
335 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')
338 self.do_note('autopeak')