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 #====================
288 def do_autopeak(self,args):
289 #FIXME: this function is too long. split it and make it rational.
290 #FIXME: also, *generalize fits* to allow FJC and any other model in the future!
292 def fit_interval_nm(start_index,plot,nm,backwards):
296 Calculates the number of points to fit, given a fit interval in nm
297 start_index: index of point
299 backwards: if true, finds a point backwards.
303 x_vect=plot.vectors[1][0]
307 start=x_vect[start_index]
309 while abs(x_vect[i]-x_vect[start_index])*(10**9) < nm:
310 if i==0 or i==maxlen-1: #we reached boundaries of vector!
322 T=self.config['temperature']
324 if 'usepoints' in args.split():
325 fit_points=int(self.config['auto_fit_points'])
331 slope_span=int(self.config['auto_slope_span'])
333 rebase=False #if true=we select rebase
336 displayed_plot=self._get_displayed_plot(0)
338 if self.current.curve.experiment != 'smfs' or len(displayed_plot.vectors) < 2:
339 print 'Cannot work on this curve.'
342 #Look for custom persistent length / custom temperature
343 for arg in args.split():
344 #look for a persistent length argument.
346 pl_expression=arg.split('=')
347 pl_value=float(pl_expression[1]) #actual value
348 #look for a T argument. FIXME: spaces are not allowed between 'pl' and value
349 if ('t=' in arg[0:2]) or ('T=' in arg[0:2]):
350 t_expression=arg.split('=')
351 T=float(t_expression[1])
353 #Handle contact point arguments
354 def pickup_contact_point():
355 #macro to pick up the contact point by clicking
356 contact_point=self._measure_N_points(N=1, whatset=1)[0]
357 contact_point_index=contact_point.index
358 self.wlccontact_point=contact_point
359 self.wlccontact_index=contact_point.index
360 self.wlccurrent=self.current.path
361 return contact_point, contact_point_index
363 if 'reclick' in args.split():
364 print 'Click contact point'
365 contact_point, contact_point_index = pickup_contact_point()
366 elif 'noauto' in args.split():
367 if self.wlccontact_index==None or self.wlccurrent != self.current.path:
368 print 'Click contact point'
369 contact_point , contact_point_index = pickup_contact_point()
371 contact_point=self.wlccontact_point
372 contact_point_index=self.wlccontact_index
374 #Automatically find contact point
375 cindex=self.find_contact_point()
376 contact_point=ClickedPoint()
377 contact_point.absolute_coords=displayed_plot.vectors[1][0][cindex], displayed_plot.vectors[1][1][cindex]
378 contact_point.find_graph_coords(displayed_plot.vectors[1][0], displayed_plot.vectors[1][1])
379 contact_point.is_marker=True
383 defplot=self.current.curve.default_plots()[0]
384 flatten=self._find_plotmanip('flatten') #Extract flatten plotmanip
385 defplot=flatten(defplot, self.current, customvalue=1) #Flatten curve before feeding it to has_peaks
386 peak_location,peak_size=self.has_peaks(defplot, self.convfilt_config['mindeviation'])
388 #Create a new plot to send
389 fitplot=copy.deepcopy(displayed_plot)
391 #Pick up force baseline
392 whatset=1 #fixme: for all sets
393 if 'rebase' in args or (self.basecurrent != self.current.path):
396 clicks=self.config['baseline_clicks']
399 self.basepoints.append(ClickedPoint())
400 self.basepoints.append(ClickedPoint())
401 self.basepoints[0].index=peak_location[-1]+fit_interval_nm(peak_location[-1], displayed_plot, self.config['auto_right_baseline'],False)
402 self.basepoints[1].index=self.basepoints[0].index+fit_interval_nm(self.basepoints[0].index, displayed_plot, self.config['auto_left_baseline'],False)
403 for point in self.basepoints:
404 #for graphing etc. purposes, fill-in with coordinates
405 point.absolute_coords=displayed_plot.vectors[1][0][point.index], displayed_plot.vectors[1][1][point.index]
406 point.find_graph_coords(displayed_plot.vectors[1][0], displayed_plot.vectors[1][1])
408 print 'Select baseline'
410 self.basepoints=self._measure_N_points(N=1, whatset=whatset)
411 self.basepoints.append(ClickedPoint())
412 self.basepoints[1].index=self.basepoints[0].index+fit_interval_nm(self.basepoints[0].index, displayed_plot, self.config['auto_left_baseline'], False)
413 #for graphing etc. purposes, fill-in with coordinates
414 self.basepoints[1].absolute_coords=displayed_plot.vectors[1][0][self.basepoints[1].index], displayed_plot.vectors[1][1][self.basepoints[1].index]
415 self.basepoints[1].find_graph_coords(displayed_plot.vectors[1][0], displayed_plot.vectors[1][1])
417 self.basepoints=self._measure_N_points(N=2, whatset=whatset)
419 self.basecurrent=self.current.path
421 boundaries=[self.basepoints[0].index, self.basepoints[1].index]
423 to_average=displayed_plot.vectors[1][1][boundaries[0]:boundaries[1]] #y points to average
424 avg=np.mean(to_average)
427 #Initialize data vectors
435 #Cycle between peaks and do analysis.
436 for peak in peak_location:
438 #FIXME: clean wlc fitting, to avoid this clickedpoint hell
439 #-create a clicked point for the peak point
440 peak_point=ClickedPoint()
441 peak_point.absolute_coords=displayed_plot.vectors[1][0][peak], displayed_plot.vectors[1][1][peak]
442 peak_point.find_graph_coords(displayed_plot.vectors[1][0], displayed_plot.vectors[1][1])
445 fit_points=fit_interval_nm(peak, displayed_plot, self.config['auto_fit_nm'], True)
447 #-create a clicked point for the other fit point
448 other_fit_point=ClickedPoint()
449 other_fit_point.absolute_coords=displayed_plot.vectors[1][0][peak-fit_points], displayed_plot.vectors[1][1][peak-fit_points]
450 other_fit_point.find_graph_coords(displayed_plot.vectors[1][0], displayed_plot.vectors[1][1])
452 points=[contact_point, peak_point, other_fit_point]
454 #Check if we have enough points for a fit. If not, wlc_fit could crash
455 if abs(peak_point.index-other_fit_point.index) < 2:
458 params, yfit, xfit = self.wlc_fit(points, displayed_plot.vectors[1][0], displayed_plot.vectors[1][1],pl_value,T)
459 #save wlc values (nm)
460 c_lengths.append(params[0]*(1.0e+9))
461 if len(params)==2: #if we did choose 2-value fit
462 p_lengths.append(params[1]*(1.0e+9))
464 p_lengths.append(pl_value)
465 #Add WLC fit lines to plot
466 fitplot.add_set(xfit,yfit)
468 if len(fitplot.styles)==0:
471 fitplot.styles.append(None)
474 delta_to_measure=displayed_plot.vectors[1][1][peak-delta_force:peak+delta_force]
475 y=min(delta_to_measure)
476 #save force values (pN)
477 forces.append(abs(y-avg)*(1.0e+12))
480 slope=self.linefit_between(peak-slope_span,peak)[0]
484 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]])
485 fitplot.styles.append('scatter')
487 #Show wlc fits and peak locations
488 self._send_plot([fitplot])
491 #Ask the user what peaks to ignore from analysis.
492 print 'Peaks to ignore (0,1...n from contact point,return to take all)'
493 print 'N to discard measurement'
494 exclude_raw=raw_input('Input:')
498 if not exclude_raw=='':
499 exclude=exclude_raw.split(',')
501 exclude=[int(item) for item in exclude]
508 print 'Bad input, taking all...'
509 #Clean data vectors from ignored peaks
510 c_lengths=[item for item in c_lengths if item != None]
511 p_lengths=[item for item in p_lengths if item != None]
512 forces=[item for item in forces if item != None]
513 slopes=[item for item in slopes if item != None]
514 print 'contour (nm)',c_lengths
515 print 'p (nm)',p_lengths
516 print 'forces (pN)',forces
517 print 'slopes (N/m)',slopes
520 if self.autofile=='':
521 self.autofile=raw_input('Autopeak filename? (return to ignore) ')
522 if self.autofile=='':
526 if not os.path.exists(self.autofile):
527 f=open(self.autofile,'w+')
528 f.write('Analysis started '+time.asctime()+'\n')
529 f.write('----------------------------------------\n')
530 f.write('; Contour length (nm) ; Persistence length (nm) ; Max.Force (pN) ; Slope (N/m) \n')
534 f=open(self.autofile,'a+')
536 f.write(self.current.path+'\n')
537 for i in range(len(c_lengths)):
538 f.write(' ; '+str(c_lengths[i])+' ; '+str(p_lengths[i])+' ; '+str(forces[i])+' ; '+str(slopes[i])+'\n')
541 self.do_note('autopeak')