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)
286 #====================
288 #====================
290 def do_autopeak(self,args):
291 #FIXME: this function is too long. split it and make it rational.
292 #FIXME: also, *generalize fits* to allow FJC and any other model in the future!
294 def fit_interval_nm(start_index,plot,nm,backwards):
298 Calculates the number of points to fit, given a fit interval in nm
299 start_index: index of point
301 backwards: if true, finds a point backwards.
305 x_vect=plot.vectors[1][0]
309 start=x_vect[start_index]
311 while abs(x_vect[i]-x_vect[start_index])*(10**9) < nm:
312 if i==0 or i==maxlen-1: #we reached boundaries of vector!
324 T=self.config['temperature']
326 if 'usepoints' in args.split():
327 fit_points=int(self.config['auto_fit_points'])
333 slope_span=int(self.config['auto_slope_span'])
335 rebase=False #if true=we select rebase
338 displayed_plot=self._get_displayed_plot(0)
340 if self.current.curve.experiment != 'smfs' or len(displayed_plot.vectors) < 2:
341 print 'Cannot work on this curve.'
344 #Look for custom persistent length / custom temperature
345 for arg in args.split():
346 #look for a persistent length argument.
348 pl_expression=arg.split('=')
349 pl_value=float(pl_expression[1]) #actual value
350 #look for a T argument. FIXME: spaces are not allowed between 'pl' and value
351 if ('t=' in arg[0:2]) or ('T=' in arg[0:2]):
352 t_expression=arg.split('=')
353 T=float(t_expression[1])
355 #Handle contact point arguments
356 def pickup_contact_point():
357 #macro to pick up the contact point by clicking
358 contact_point=self._measure_N_points(N=1, whatset=1)[0]
359 contact_point_index=contact_point.index
360 self.wlccontact_point=contact_point
361 self.wlccontact_index=contact_point.index
362 self.wlccurrent=self.current.path
363 return contact_point, contact_point_index
365 if 'reclick' in args.split():
366 print 'Click contact point'
367 contact_point, contact_point_index = pickup_contact_point()
368 elif 'noauto' in args.split():
369 if self.wlccontact_index==None or self.wlccurrent != self.current.path:
370 print 'Click contact point'
371 contact_point , contact_point_index = pickup_contact_point()
373 contact_point=self.wlccontact_point
374 contact_point_index=self.wlccontact_index
376 #Automatically find contact point
377 cindex=self.find_contact_point()
378 contact_point=ClickedPoint()
379 contact_point.absolute_coords=displayed_plot.vectors[1][0][cindex], displayed_plot.vectors[1][1][cindex]
380 contact_point.find_graph_coords(displayed_plot.vectors[1][0], displayed_plot.vectors[1][1])
381 contact_point.is_marker=True
385 defplot=self.current.curve.default_plots()[0]
386 flatten=self._find_plotmanip('flatten') #Extract flatten plotmanip
387 defplot=flatten(defplot, self.current, customvalue=1) #Flatten curve before feeding it to has_peaks
388 peak_location,peak_size=self.has_peaks(defplot, self.convfilt_config['mindeviation'])
390 #Create a new plot to send
391 fitplot=copy.deepcopy(displayed_plot)
393 #Pick up force baseline
394 whatset=1 #fixme: for all sets
395 if 'rebase' in args or (self.basecurrent != self.current.path):
398 clicks=self.config['baseline_clicks']
401 self.basepoints.append(ClickedPoint())
402 self.basepoints.append(ClickedPoint())
403 self.basepoints[0].index=peak_location[-1]+fit_interval_nm(peak_location[-1], displayed_plot, self.config['auto_right_baseline'],False)
404 self.basepoints[1].index=self.basepoints[0].index+fit_interval_nm(self.basepoints[0].index, displayed_plot, self.config['auto_left_baseline'],False)
405 for point in self.basepoints:
406 #for graphing etc. purposes, fill-in with coordinates
407 point.absolute_coords=displayed_plot.vectors[1][0][point.index], displayed_plot.vectors[1][1][point.index]
408 point.find_graph_coords(displayed_plot.vectors[1][0], displayed_plot.vectors[1][1])
410 print 'Select baseline'
412 self.basepoints=self._measure_N_points(N=1, whatset=whatset)
413 self.basepoints.append(ClickedPoint())
414 self.basepoints[1].index=self.basepoints[0].index+fit_interval_nm(self.basepoints[0].index, displayed_plot, self.config['auto_left_baseline'], False)
415 #for graphing etc. purposes, fill-in with coordinates
416 self.basepoints[1].absolute_coords=displayed_plot.vectors[1][0][self.basepoints[1].index], displayed_plot.vectors[1][1][self.basepoints[1].index]
417 self.basepoints[1].find_graph_coords(displayed_plot.vectors[1][0], displayed_plot.vectors[1][1])
419 self.basepoints=self._measure_N_points(N=2, whatset=whatset)
421 self.basecurrent=self.current.path
423 boundaries=[self.basepoints[0].index, self.basepoints[1].index]
425 to_average=displayed_plot.vectors[1][1][boundaries[0]:boundaries[1]] #y points to average
426 avg=np.mean(to_average)
429 #Initialize data vectors
437 #Cycle between peaks and do analysis.
438 for peak in peak_location:
440 #FIXME: clean wlc fitting, to avoid this clickedpoint hell
441 #-create a clicked point for the peak point
442 peak_point=ClickedPoint()
443 peak_point.absolute_coords=displayed_plot.vectors[1][0][peak], displayed_plot.vectors[1][1][peak]
444 peak_point.find_graph_coords(displayed_plot.vectors[1][0], displayed_plot.vectors[1][1])
447 fit_points=fit_interval_nm(peak, displayed_plot, self.config['auto_fit_nm'], True)
449 #-create a clicked point for the other fit point
450 other_fit_point=ClickedPoint()
451 other_fit_point.absolute_coords=displayed_plot.vectors[1][0][peak-fit_points], displayed_plot.vectors[1][1][peak-fit_points]
452 other_fit_point.find_graph_coords(displayed_plot.vectors[1][0], displayed_plot.vectors[1][1])
454 points=[contact_point, peak_point, other_fit_point]
456 #Check if we have enough points for a fit. If not, wlc_fit could crash
457 if abs(peak_point.index-other_fit_point.index) < 2:
460 params, yfit, xfit = self.wlc_fit(points, displayed_plot.vectors[1][0], displayed_plot.vectors[1][1],pl_value,T)
461 #save wlc values (nm)
462 c_lengths.append(params[0]*(1.0e+9))
463 if len(params)==2: #if we did choose 2-value fit
464 p_lengths.append(params[1]*(1.0e+9))
466 p_lengths.append(pl_value)
467 #Add WLC fit lines to plot
468 fitplot.add_set(xfit,yfit)
470 if len(fitplot.styles)==0:
473 fitplot.styles.append(None)
476 delta_to_measure=displayed_plot.vectors[1][1][peak-delta_force:peak+delta_force]
477 y=min(delta_to_measure)
478 #save force values (pN)
479 forces.append(abs(y-avg)*(1.0e+12))
482 slope=self.linefit_between(peak-slope_span,peak)[0]
486 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]])
487 fitplot.styles.append('scatter')
489 #Show wlc fits and peak locations
490 self._send_plot([fitplot])
493 #Ask the user what peaks to ignore from analysis.
494 print 'Peaks to ignore (0,1...n from contact point,return to take all)'
495 print 'N to discard measurement'
496 exclude_raw=raw_input('Input:')
500 if not exclude_raw=='':
501 exclude=exclude_raw.split(',')
503 exclude=[int(item) for item in exclude]
510 print 'Bad input, taking all...'
511 #Clean data vectors from ignored peaks
512 c_lengths=[item for item in c_lengths if item != None]
513 p_lengths=[item for item in p_lengths if item != None]
514 forces=[item for item in forces if item != None]
515 slopes=[item for item in slopes if item != None]
516 print 'contour (nm)',c_lengths
517 print 'p (nm)',p_lengths
518 print 'forces (pN)',forces
519 print 'slopes (N/m)',slopes
522 if self.autofile=='':
523 self.autofile=raw_input('Autopeak filename? (return to ignore) ')
524 if self.autofile=='':
528 if not os.path.exists(self.autofile):
529 f=open(self.autofile,'w+')
530 f.write('Analysis started '+time.asctime()+'\n')
531 f.write('----------------------------------------\n')
532 f.write('; Contour length (nm) ; Persistence length (nm) ; Max.Force (pN) ; Slope (N/m) \n')
536 f=open(self.autofile,'a+')
538 f.write(self.current.path+'\n')
539 for i in range(len(c_lengths)):
540 f.write(' ; '+str(c_lengths[i])+' ; '+str(p_lengths[i])+' ; '+str(forces[i])+' ; '+str(slopes[i])+'\n')
543 self.do_note('autopeak')