6 Plugin regarding general velocity clamp measurements
9 from libhooke import WX_GOOD, ClickedPoint
11 wxversion.select(WX_GOOD)
12 from wx import PostEvent
20 warnings.simplefilter('ignore',np.RankWarning)
23 class generalvclampCommands:
30 def do_distance(self,args):
34 Measure the distance (in nm) between two points.
35 For a standard experiment this is the delta X distance.
36 For a force clamp experiment this is the delta Y distance (actually becomes
41 if self.current.curve.experiment == 'clamp':
42 print 'You wanted to use zpiezo perhaps?'
45 dx,unitx,dy,unity=self._delta(set=1)
46 print str(dx*(10**9))+' nm'
47 to_dump='distance '+self.current.path+' '+str(dx*(10**9))+' nm'
48 self.outlet.push(to_dump)
51 def do_force(self,args):
55 Measure the force difference (in pN) between two points
59 if self.current.curve.experiment == 'clamp':
60 print 'This command makes no sense for a force clamp experiment.'
62 dx,unitx,dy,unity=self._delta(set=1)
63 print str(dy*(10**12))+' pN'
64 to_dump='force '+self.current.path+' '+str(dy*(10**12))+' pN'
65 self.outlet.push(to_dump)
68 def do_forcebase(self,args):
72 Measures the difference in force (in pN) between a point and a baseline
73 took as the average between two points.
75 The baseline is fixed once for a given curve and different force measurements,
76 unless the user wants it to be recalculated
78 Syntax: forcebase [rebase]
79 rebase: Forces forcebase to ask again the baseline
80 max: Instead of asking for a point to measure, asks for two points and use
81 the maximum peak in between
83 rebase=False #if true=we select rebase
84 maxpoint=False #if true=we measure the maximum peak
86 plot=self._get_displayed_plot()
87 whatset=1 #fixme: for all sets
88 if 'rebase' in args or (self.basecurrent != self.current.path):
94 print 'Select baseline'
95 self.basepoints=self._measure_N_points(N=2, whatset=whatset)
96 self.basecurrent=self.current.path
99 print 'Select two points'
100 points=self._measure_N_points(N=2, whatset=whatset)
101 boundpoints=[points[0].index, points[1].index]
104 y=min(plot.vectors[whatset][1][boundpoints[0]:boundpoints[1]])
106 print 'Chosen interval not valid. Try picking it again. Did you pick the same point as begin and end of interval?'
108 print 'Select point to measure'
109 points=self._measure_N_points(N=1, whatset=whatset)
110 #whatplot=points[0].dest
111 y=points[0].graph_coords[1]
113 #fixme: code duplication
114 boundaries=[self.basepoints[0].index, self.basepoints[1].index]
116 to_average=plot.vectors[whatset][1][boundaries[0]:boundaries[1]] #y points to average
118 avg=np.mean(to_average)
120 print str(forcebase*(10**12))+' pN'
121 to_dump='forcebase '+self.current.path+' '+str(forcebase*(10**12))+' pN'
122 self.outlet.push(to_dump)
125 def plotmanip_flatten(self, plot, current, customvalue=False):
127 Subtracts a polynomial fit to the non-contact part of the curve, as to flatten it.
128 the best polynomial fit is chosen among polynomials of degree 1 to n, where n is
129 given by the configuration file or by the customvalue.
131 customvalue= int (>0) --> starts the function even if config says no (default=False)
135 if current.curve.experiment != 'smfs':
138 #only one set is present...
139 if len(self.plots[0].vectors) != 2:
142 #config is not flatten, and customvalue flag is false too
143 if (not self.config['flatten']) and (not customvalue):
150 max_cycles=customvalue
152 max_cycles=self.config['flatten'] #Using > 1 usually doesn't help and can give artefacts. However, it could be useful too.
154 contact_index=self.find_contact_point()
156 valn=[[] for item in range(max_exponent)]
157 yrn=[0.0 for item in range(max_exponent)]
158 errn=[0.0 for item in range(max_exponent)]
160 for i in range(int(max_cycles)):
162 x_ext=plot.vectors[0][0][contact_index+delta_contact:]
163 y_ext=plot.vectors[0][1][contact_index+delta_contact:]
164 x_ret=plot.vectors[1][0][contact_index+delta_contact:]
165 y_ret=plot.vectors[1][1][contact_index+delta_contact:]
166 for exponent in range(max_exponent):
168 valn[exponent]=sp.polyfit(x_ext,y_ext,exponent)
169 yrn[exponent]=sp.polyval(valn[exponent],x_ret)
170 errn[exponent]=sp.sqrt(sum((yrn[exponent]-y_ext)**2)/float(len(y_ext)))
172 print 'Cannot flatten!'
176 best_exponent=errn.index(min(errn))
179 ycorr_ext=y_ext-yrn[best_exponent]+y_ext[0] #noncontact part
180 yjoin_ext=np.array(plot.vectors[0][1][0:contact_index+delta_contact]) #contact part
182 ycorr_ret=y_ret-yrn[best_exponent]+y_ext[0] #noncontact part
183 yjoin_ret=np.array(plot.vectors[1][1][0:contact_index+delta_contact]) #contact part
185 ycorr_ext=np.concatenate((yjoin_ext, ycorr_ext))
186 ycorr_ret=np.concatenate((yjoin_ret, ycorr_ret))
188 plot.vectors[0][1]=list(ycorr_ext)
189 plot.vectors[1][1]=list(ycorr_ret)
194 def do_slope(self,args):
198 Measures the slope of a delimited chunk on the return trace.
199 The chunk can be delimited either by two manual clicks, or have
200 a fixed width, given as an argument.
202 Syntax: slope [width]
203 The facultative [width] parameter specifies how many
204 points will be considered for the fit. If [width] is
205 specified, only one click will be required.
206 (c) Marco Brucale, Massimo Sandal 2008
209 # Reads the facultative width argument
215 # Decides between the two forms of user input, as per (args)
217 # Gets the Xs of two clicked points as indexes on the current curve vector
218 print 'Click twice to delimit chunk'
220 points=self._measure_N_points(N=2,whatset=1)
221 clickedpoints=[points[0].index,points[1].index]
224 print 'Click once on the leftmost point of the chunk (i.e.usually the peak)'
226 points=self._measure_N_points(N=1,whatset=1)
227 clickedpoints=[points[0].index-fitspan,points[0].index]
229 # Calls the function linefit_between
230 parameters=[0,0,[],[]]
231 parameters=self.linefit_between(clickedpoints[0],clickedpoints[1])
233 # Outputs the relevant slope parameter
235 print str(parameters[0])
236 to_dump='slope '+self.current.path+' '+str(parameters[0])
237 self.outlet.push(to_dump)
239 # Makes a vector with the fitted parameters and sends it to the GUI
240 xtoplot=parameters[2]
244 ytoplot.append((x*parameters[0])+parameters[1])
246 clickvector_x, clickvector_y=[], []
248 clickvector_x.append(item.graph_coords[0])
249 clickvector_y.append(item.graph_coords[1])
251 lineplot=self._get_displayed_plot(0) #get topmost displayed plot
253 lineplot.add_set(xtoplot,ytoplot)
254 lineplot.add_set(clickvector_x, clickvector_y)
256 if lineplot.styles==[]:
257 lineplot.styles=[None,None,None,'scatter']
259 lineplot.styles+=[None,'scatter']
261 self._send_plot([lineplot])
263 def linefit_between(self,index1,index2):
265 Creates two vectors (xtofit,ytofit) slicing out from the
266 current return trace a portion delimited by the two indexes
268 Then does a least squares linear fit on that slice.
269 Finally returns [0]=the slope, [1]=the intercept of the
270 fitted 1st grade polynomial, and [2,3]=the actual (x,y) vectors
272 (c) Marco Brucale, Massimo Sandal 2008
274 # Translates the indexes into two vectors containing the x,y data to fit
275 xtofit=self.plots[0].vectors[1][0][index1:index2]
276 ytofit=self.plots[0].vectors[1][1][index1:index2]
278 # Does the actual linear fitting (simple least squares with numpy.polyfit)
280 linefit=np.polyfit(xtofit,ytofit,1)
282 return (linefit[0],linefit[1],xtofit,ytofit)
284 #====================
286 #====================
287 def do_autopeak(self,args):
291 Automatically performs a number of analyses on the peaks of the given curve.
292 Currently it automatically:
293 - fits peaks with WLC function
294 - measures peak maximum forces with a baseline
295 - measures slope in proximity of peak maximum
296 Requires flatten plotmanipulator , fit.py plugin , flatfilts.py plugin with convfilt
299 autopeak [rebase] [pl=value] [t=value] [noauto] [reclick]
301 rebase : Re-asks baseline interval
303 pl=[value] : Use a fixed persistent length for the fit. If pl is not given,
304 the fit will be a 2-variable
305 fit. DO NOT put spaces between 'pl', '=' and the value.
306 The value must be in meters.
307 Scientific notation like 0.35e-9 is fine.
309 t=[value] : Use a user-defined temperature. The value must be in
310 kelvins; by default it is 293 K.
311 DO NOT put spaces between 't', '=' and the value.
313 noauto : allows for clicking the contact point by
314 hand (otherwise it is automatically estimated) the first time.
315 If subsequent measurements are made, the same contact point
316 clicked the first time is used
318 reclick : redefines by hand the contact point, if noauto has been used before
319 but the user is unsatisfied of the previously choosen contact point.
321 usepoints : fit interval by number of points instead than by nanometers
323 When you first issue the command, it will ask for the filename. If you are giving the filename
324 of an existing file, autopeak will resume it and append measurements to it. If you are giving
325 a new filename, it will create the file and append to it until you close Hooke.
328 Useful variables (to set with SET command):
330 temperature= temperature of the system for wlc fit (in K)
332 auto_slope_span = number of points on which measure the slope, for slope
334 auto_fit_nm = number of nm to fit before the peak maximum, for WLC (if usepoints false)
335 auto_fit_points = number of points to fit before the peak maximum, for WLC (if usepoints true)
337 baseline_clicks = 0: automatic baseline
338 1: decide baseline with a single click and length defined in auto_left_baseline
339 2: let user click points of baseline
340 auto_left_baseline = length in nm to use as baseline from the right point (if baseline_clicks=0 , 1)
341 auto_right_baseline = distance in nm of peak-most baseline point from last peak (if baseline_clicks = 0)
344 def fit_interval_nm(start_index,plot,nm,backwards):
346 Calculates the number of points to fit, given a fit interval in nm
347 start_index: index of point
349 backwards: if true, finds a point backwards.
351 x_vect=plot.vectors[1][0]
355 start=x_vect[start_index]
357 while abs(x_vect[i]-x_vect[start_index])*(10**9) < nm:
359 if i==0 or i==maxlen: #we reached boundaries of vector!
371 T=self.config['temperature']
373 if 'usepoints' in args.split():
374 fit_points=int(self.config['auto_fit_points'])
380 slope_span=int(self.config['auto_slope_span'])
382 rebase=False #if true=we select rebase
385 displayed_plot=self._get_displayed_plot(0)
387 if self.current.curve.experiment != 'smfs' or len(displayed_plot.vectors) < 2:
388 print 'Cannot work on this curve.'
391 #Look for custom persistent length / custom temperature
392 for arg in args.split():
393 #look for a persistent length argument.
395 pl_expression=arg.split('=')
396 pl_value=float(pl_expression[1]) #actual value
397 #look for a T argument. FIXME: spaces are not allowed between 'pl' and value
398 if ('t=' in arg[0:2]) or ('T=' in arg[0:2]):
399 t_expression=arg.split('=')
400 T=float(t_expression[1])
403 #Handle contact point arguments
404 def pickup_contact_point():
405 '''macro to pick up the contact point by clicking'''
406 contact_point=self._measure_N_points(N=1, whatset=1)[0]
407 contact_point_index=contact_point.index
408 self.wlccontact_point=contact_point
409 self.wlccontact_index=contact_point.index
410 self.wlccurrent=self.current.path
411 return contact_point, contact_point_index
413 if 'reclick' in args.split():
414 print 'Click contact point'
415 contact_point, contact_point_index = pickup_contact_point()
416 elif 'noauto' in args.split():
417 if self.wlccontact_index==None or self.wlccurrent != self.current.path:
418 print 'Click contact point'
419 contact_point , contact_point_index = pickup_contact_point()
421 contact_point=self.wlccontact_point
422 contact_point_index=self.wlccontact_index
424 #Automatically find contact point
425 cindex=self.find_contact_point()
426 contact_point=ClickedPoint()
427 contact_point.absolute_coords=displayed_plot.vectors[1][0][cindex], displayed_plot.vectors[1][1][cindex]
428 contact_point.find_graph_coords(displayed_plot.vectors[1][0], displayed_plot.vectors[1][1])
429 contact_point.is_marker=True
433 defplot=self.current.curve.default_plots()[0]
434 flatten=self._find_plotmanip('flatten') #Extract flatten plotmanip
435 defplot=flatten(defplot, self.current, customvalue=1) #Flatten curve before feeding it to has_peaks
436 peak_location,peak_size=self.has_peaks(defplot, self.convfilt_config['mindeviation'])
438 #Create a new plot to send
439 fitplot=copy.deepcopy(displayed_plot)
441 #Pick up force baseline
442 whatset=1 #fixme: for all sets
443 if 'rebase' in args or (self.basecurrent != self.current.path):
446 clicks=self.config['baseline_clicks']
449 self.basepoints.append(ClickedPoint())
450 self.basepoints.append(ClickedPoint())
451 self.basepoints[0].index=peak_location[-1]+fit_interval_nm(peak_location[-1], displayed_plot, self.config['auto_right_baseline'],False)
452 self.basepoints[1].index=self.basepoints[0].index+fit_interval_nm(self.basepoints[0].index, displayed_plot, self.config['auto_left_baseline'],False)
453 for point in self.basepoints:
454 #for graphing etc. purposes, fill-in with coordinates
455 point.absolute_coords=displayed_plot.vectors[1][0][point.index], displayed_plot.vectors[1][1][point.index]
456 point.find_graph_coords(displayed_plot.vectors[1][0], displayed_plot.vectors[1][1])
458 print 'Select baseline'
460 self.basepoints=self._measure_N_points(N=1, whatset=whatset)
461 self.basepoints.append(ClickedPoint())
462 self.basepoints[1].index=self.basepoints[0].index+fit_interval_nm(self.basepoints[0].index, displayed_plot, self.config['auto_left_baseline'], False)
463 #for graphing etc. purposes, fill-in with coordinates
464 self.basepoints[1].absolute_coords=displayed_plot.vectors[1][0][self.basepoints[1].index], displayed_plot.vectors[1][1][self.basepoints[1].index]
465 self.basepoints[1].find_graph_coords(displayed_plot.vectors[1][0], displayed_plot.vectors[1][1])
467 self.basepoints=self._measure_N_points(N=2, whatset=whatset)
469 self.basecurrent=self.current.path
471 boundaries=[self.basepoints[0].index, self.basepoints[1].index]
473 to_average=displayed_plot.vectors[1][1][boundaries[0]:boundaries[1]] #y points to average
474 avg=np.mean(to_average)
477 #Initialize data vectors
485 #Cycle between peaks and do analysis.
486 for peak in peak_location:
488 #FIXME: clean wlc fitting, to avoid this clickedpoint hell
489 #-create a clicked point for the peak point
490 peak_point=ClickedPoint()
491 peak_point.absolute_coords=displayed_plot.vectors[1][0][peak], displayed_plot.vectors[1][1][peak]
492 peak_point.find_graph_coords(displayed_plot.vectors[1][0], displayed_plot.vectors[1][1])
495 fit_points=fit_interval_nm(peak, displayed_plot, self.config['auto_fit_nm'], True)
497 #-create a clicked point for the other fit point
498 other_fit_point=ClickedPoint()
499 other_fit_point.absolute_coords=displayed_plot.vectors[1][0][peak-fit_points], displayed_plot.vectors[1][1][peak-fit_points]
500 other_fit_point.find_graph_coords(displayed_plot.vectors[1][0], displayed_plot.vectors[1][1])
502 points=[contact_point, peak_point, other_fit_point]
504 #Check if we have enough points for a fit. If not, wlc_fit could crash
505 if abs(peak_point.index-other_fit_point.index) < 2:
508 params, yfit, xfit = self.wlc_fit(points, displayed_plot.vectors[1][0], displayed_plot.vectors[1][1],pl_value,T)
509 #save wlc values (nm)
510 c_lengths.append(params[0]*(1.0e+9))
511 if len(params)==2: #if we did choose 2-value fit
512 p_lengths.append(params[1]*(1.0e+9))
514 p_lengths.append(pl_value)
515 #Add WLC fit lines to plot
516 fitplot.add_set(xfit,yfit)
518 if len(fitplot.styles)==0:
521 fitplot.styles.append(None)
524 delta_to_measure=displayed_plot.vectors[1][1][peak-delta_force:peak+delta_force]
525 y=min(delta_to_measure)
526 #save force values (pN)
527 forces.append(abs(y-avg)*(1.0e+12))
530 slope=self.linefit_between(peak-slope_span,peak)[0]
534 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]])
535 fitplot.styles.append('scatter')
537 #Show wlc fits and peak locations
538 self._send_plot([fitplot])
541 #Ask the user what peaks to ignore from analysis.
542 print 'Peaks to ignore (0,1...n from contact point,return to take all)'
543 print 'N to discard measurement'
544 exclude_raw=raw_input('Input:')
548 if not exclude_raw=='':
549 exclude=exclude_raw.split(',')
551 exclude=[int(item) for item in exclude]
558 print 'Bad input, taking all...'
559 #Clean data vectors from ignored peaks
560 c_lengths=[item for item in c_lengths if item != None]
561 p_lengths=[item for item in p_lengths if item != None]
562 forces=[item for item in forces if item != None]
563 slopes=[item for item in slopes if item != None]
564 print 'contour (nm)',c_lengths
565 print 'p (nm)',p_lengths
566 print 'forces (pN)',forces
567 print 'slopes (N/m)',slopes
570 if self.autofile=='':
571 self.autofile=raw_input('Autopeak filename? (return to ignore) ')
572 if self.autofile=='':
576 if not os.path.exists(self.autofile):
577 f=open(self.autofile,'w+')
578 f.write('Analysis started '+time.asctime()+'\n')
579 f.write('----------------------------------------\n')
580 f.write('; Contour length (nm) ; Persistence length (nm) ; Max.Force (pN) ; Slope (N/m) \n')
584 f=open(self.autofile,'a+')
586 f.write(self.current.path+'\n')
587 for i in range(len(c_lengths)):
588 f.write(' ; '+str(c_lengths[i])+' ; '+str(p_lengths[i])+' ; '+str(forces[i])+' ; '+str(slopes[i])+'\n')
591 self.do_note('autopeak')