1 # -*- coding: utf-8 -*-
3 from libhooke import WX_GOOD, ClickedPoint
5 wxversion.select(WX_GOOD)
6 from wx import PostEvent
14 warnings.simplefilter('ignore',np.RankWarning)
17 class autopeakCommands:
19 def do_autopeak(self,args):
23 Automatically performs a number of analyses on the peaks of the given curve.
24 Currently it automatically:
25 - fits peaks with WLC or FJC function (depending on how the fit_function variable is set)
26 - measures peak maximum forces with a baseline
27 - measures slope in proximity of peak maximum
28 Requires flatten plotmanipulator , fit.py plugin , flatfilts.py plugin with convfilt
31 autopeak [rebase] [pl=value] [t=value] [noauto] [reclick]
33 rebase : Re-asks baseline interval
35 pl=[value] : Use a fixed persistent length for the fit. If pl is not given,
36 the fit will be a 2-variable
37 fit. DO NOT put spaces between 'pl', '=' and the value.
38 The value must be in meters.
39 Scientific notation like 0.35e-9 is fine.
41 t=[value] : Use a user-defined temperature. The value must be in
42 kelvins; by default it is 293 K.
43 DO NOT put spaces between 't', '=' and the value.
45 noauto : allows for clicking the contact point by
46 hand (otherwise it is automatically estimated) the first time.
47 If subsequent measurements are made, the same contact point
48 clicked the first time is used
50 reclick : redefines by hand the contact point, if noauto has been used before
51 but the user is unsatisfied of the previously choosen contact point.
53 usepoints : fit interval by number of points instead than by nanometers
55 noflatten : does not use the "flatten" plot manipulator
57 When you first issue the command, it will ask for the filename. If you are giving the filename
58 of an existing file, autopeak will resume it and append measurements to it. If you are giving
59 a new filename, it will create the file and append to it until you close Hooke.
62 Useful variables (to set with SET command):
64 fit_function = type of function to use for elasticity. If "wlc" worm-like chain is used, if "fjc" freely jointed
67 temperature= temperature of the system for wlc/fjc fit (in K)
69 auto_slope_span = number of points on which measure the slope, for slope
71 auto_fit_nm = number of nm to fit before the peak maximum, for WLC/FJC (if usepoints false)
72 auto_fit_points = number of points to fit before the peak maximum, for WLC/FJC (if usepoints true)
74 baseline_clicks = -1: no baseline, f=0 at the contact point (whether hand-picked or automatically found)
76 1: decide baseline with a single click and length defined in auto_left_baseline
77 2: let user click points of baseline
78 auto_left_baseline = length in nm to use as baseline from the right point (if baseline_clicks=0 , 1)
79 auto_right_baseline = distance in nm of peak-most baseline point from last peak (if baseline_clicks = 0)
81 auto_min_p ; auto_max_p = Minimum and maximum persistence length (if using WLC) or Kuhn length (if using FJC)
82 outside of which the peak is automatically discarded (in nm)
86 #FIXME: to move outside function
87 def fit_interval_nm(start_index,plot,nm,backwards):
89 Calculates the number of points to fit, given a fit interval in nm
90 start_index: index of point
92 backwards: if true, finds a point backwards.
94 whatset=1 #FIXME: should be decidable
95 x_vect=plot.vectors[1][0]
99 start=x_vect[start_index]
101 while abs(x_vect[i]-x_vect[start_index])*(10**9) < nm:
102 if i==0 or i==maxlen-1: #we reached boundaries of vector!
112 def pickup_contact_point():
113 '''macro to pick up the contact point by clicking'''
114 contact_point=self._measure_N_points(N=1, whatset=1)[0]
115 contact_point_index=contact_point.index
116 self.wlccontact_point=contact_point
117 self.wlccontact_index=contact_point.index
118 self.wlccurrent=self.current.path
119 return contact_point, contact_point_index
121 def find_current_peaks(noflatten):
123 defplot=self.current.curve.default_plots()[0]
125 flatten=self._find_plotmanip('flatten') #Extract flatten plotmanip
126 defplot=flatten(defplot, self.current, customvalue=1) #Flatten curve before feeding it to has_peaks
127 peak_location,peak_size=self.has_peaks(defplot, self.convfilt_config['mindeviation'])
128 return peak_location, peak_size
130 #default fit etc. variables
132 T=self.config['temperature']
134 slope_span=int(self.config['auto_slope_span'])
136 rebase=False #if true=we select rebase
137 noflatten=False #if true=we avoid flattening
139 #initialize output data vectors
148 displayed_plot=self._get_displayed_plot(0)
150 #COMMAND LINE PARSING
151 #--Using points instead of nm interval
152 if 'usepoints' in args.split():
153 fit_points=int(self.config['auto_fit_points'])
158 #--Recalculate baseline
159 if 'rebase' in args or (self.basecurrent != self.current.path):
162 if 'noflatten' in args:
165 #--Custom persistent length / custom temperature
166 for arg in args.split():
167 #look for a persistent length argument.
169 pl_expression=arg.split('=')
170 pl_value=float(pl_expression[1]) #actual value
171 #look for a T argument. FIXME: spaces are not allowed between 'pl' and value
172 if ('t=' in arg[0:2]) or ('T=' in arg[0:2]):
173 t_expression=arg.split('=')
174 T=float(t_expression[1])
175 #--Contact point arguments
176 if 'reclick' in args.split():
177 print 'Click contact point'
178 contact_point, contact_point_index = pickup_contact_point()
179 elif 'noauto' in args.split():
180 if self.wlccontact_index==None or self.wlccurrent != self.current.path:
181 print 'Click contact point'
182 contact_point , contact_point_index = pickup_contact_point()
184 contact_point=self.wlccontact_point
185 contact_point_index=self.wlccontact_index
187 #Automatically find contact point
188 cindex=self.find_contact_point()
189 contact_point=self._clickize(displayed_plot.vectors[1][0], displayed_plot.vectors[1][1], cindex)
190 #--END COMMAND LINE PARSING--
193 peak_location, peak_size = find_current_peaks(noflatten)
195 if len(peak_location) == 0:
196 print 'No peaks to fit.'
199 fitplot=copy.deepcopy(displayed_plot)
201 #Pick up force baseline
203 clicks=self.config['baseline_clicks']
206 base_index_0=peak_location[-1]+fit_interval_nm(peak_location[-1], displayed_plot, self.config['auto_right_baseline'],False)
207 self.basepoints.append(self._clickize(displayed_plot.vectors[1][0],displayed_plot.vectors[1][1],base_index_0))
208 base_index_1=self.basepoints[0].index+fit_interval_nm(self.basepoints[0].index, displayed_plot, self.config['auto_left_baseline'],False)
209 self.basepoints.append(self._clickize(displayed_plot.vectors[1][0],displayed_plot.vectors[1][1],base_index_1))
211 print 'Select baseline'
213 self.basepoints=self._measure_N_points(N=1, whatset=whatset)
214 base_index_1=self.basepoints[0].index+fit_interval_nm(self.basepoints[0].index, displayed_plot, self.config['auto_left_baseline'], False)
215 self.basepoints.append(self._clickize(displayed_plot.vectors[1][0],displayed_plot.vectors[1][1],base_index_1))
217 self.basepoints=self._measure_N_points(N=2, whatset=whatset)
219 self.basecurrent=self.current.path
221 boundaries=[self.basepoints[0].index, self.basepoints[1].index]
223 to_average=displayed_plot.vectors[1][1][boundaries[0]:boundaries[1]] #y points to average
224 avg=np.mean(to_average)
226 clicks=self.config['baseline_clicks']
229 avg=displayed_plot.vectors[1][1][contact_point_index]
231 avg=displayed_plot.vectors[1][1][cindex]
233 for peak in peak_location:
237 fit_points=fit_interval_nm(peak, displayed_plot, self.config['auto_fit_nm'], True)
238 peak_point=self._clickize(displayed_plot.vectors[1][0],displayed_plot.vectors[1][1],peak)
239 other_fit_point=self._clickize(displayed_plot.vectors[1][0],displayed_plot.vectors[1][1],peak-fit_points)
242 points=[contact_point, peak_point, other_fit_point]
244 if abs(peak_point.index-other_fit_point.index) < 2:
247 if self.config['fit_function']=='wlc':
249 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)
250 elif self.config['fit_function']=='fjc':
251 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)
253 print 'Unknown fit function'
254 print 'Please set fit_function as wlc or fjc'
259 delta_to_measure=displayed_plot.vectors[1][1][peak-delta_force:peak+delta_force]
260 y=min(delta_to_measure)
261 #save force values (pN)
263 slope=self.linefit_between(peak-slope_span,peak)[0]
266 #check fitted data and, if right, add peak to the measurement
267 #FIXME: code duplication
268 if len(params)==1: #if we did choose 1-value fit
269 p_lengths.append(pl_value)
270 c_lengths.append(params[0]*(1.0e+9))
271 sigma_p_lengths.append(0)
272 sigma_c_lengths.append(fit_errors[0]*(1.0e+9))
273 forces.append(abs(y-avg)*(1.0e+12))
275 #Add WLC fit lines to plot
276 fitplot.add_set(xfit,yfit)
277 if len(fitplot.styles)==0:
281 fitplot.styles.append(None)
282 fitplot.colors.append(None)
284 p_leng=params[1]*(1.0e+9)
285 #check if persistent length makes sense. otherwise, discard peak.
286 if p_leng>self.config['auto_min_p'] and p_leng<self.config['auto_max_p']:
287 p_lengths.append(p_leng)
288 c_lengths.append(params[0]*(1.0e+9))
289 sigma_c_lengths.append(fit_errors[0]*(1.0e+9))
290 sigma_p_lengths.append(fit_errors[1]*(1.0e+9))
291 forces.append(abs(y-avg)*(1.0e+12))
294 #Add WLC fit lines to plot
295 fitplot.add_set(xfit,yfit)
296 if len(fitplot.styles)==0:
300 fitplot.styles.append(None)
301 fitplot.colors.append(None)
306 #add basepoints to fitplot
307 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]])
308 fitplot.styles.append('scatter')
309 fitplot.colors.append(None)
311 #Show wlc fits and peak locations
312 self._send_plot([fitplot])
315 print 'Using fit function: ',self.config['fit_function']
316 print 'Measurements for all peaks detected:'
317 print 'contour (nm)', c_lengths
318 print 'sigma contour (nm)',sigma_c_lengths
319 print 'p (nm)',p_lengths
320 print 'sigma p (nm)',sigma_p_lengths
321 print 'forces (pN)',forces
322 print 'slopes (N/m)',slopes
325 while controller==False:
326 #Ask the user what peaks to ignore from analysis.
327 print 'Peaks to ignore (0,1...n from contact point,return to take all)'
328 print 'N to discard measurement'
329 exclude_raw=raw_input('Input:')
335 if not exclude_raw=='':
336 exclude=exclude_raw.split(',')
338 exclude=[int(item) for item in exclude]
344 sigma_c_lengths[i]=None
345 sigma_p_lengths[i]=None
351 #Clean data vectors from ignored peaks
352 #FIXME:code duplication
353 c_lengths=[item for item in c_lengths if item != None]
354 p_lengths=[item for item in p_lengths if item != None]
355 forces=[item for item in forces if item != None]
356 slopes=[item for item in slopes if item != None]
357 sigma_c_lengths=[item for item in sigma_c_lengths if item != None]
358 sigma_p_lengths=[item for item in sigma_p_lengths if item != None]
360 print 'Measurements for chosen peaks:'
361 print 'contour (nm)',c_lengths
362 print 'sigma contour (nm)',sigma_c_lengths
363 print 'p (nm)',p_lengths
364 print 'sigma p (nm)',sigma_p_lengths
365 print 'forces (pN)',forces
366 print 'slopes (N/m)',slopes
369 if self.autofile=='':
370 self.autofile=raw_input('Autopeak filename? (return to ignore) ')
371 if self.autofile=='':
375 if not os.path.exists(self.autofile):
376 f=open(self.autofile,'w+')
377 f.write('Analysis started '+time.asctime()+'\n')
378 f.write('----------------------------------------\n')
379 f.write('; Contour length (nm) ; Persistence length (nm) ; Max.Force (pN) ; Slope (N/m) ; Sigma contour (nm) ; Sigma persistence (nm)\n')
383 f=open(self.autofile,'a+')
385 f.write(self.current.path+'\n')
386 for i in range(len(c_lengths)):
387 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')
390 self.do_note('autopeak')