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 When you first issue the command, it will ask for the filename. If you are giving the filename
322 of an existing file, autopeak will resume it and append measurements to it. If you are giving
323 a new filename, it will create the file and append to it until you close Hooke.
325 Useful variables (to set with SET command):
327 temperature= temperature of the system for wlc fit (in K)
328 auto_fit_points = number of points to fit before the peak maximum, for wlc
329 auto_slope_span = number of points on which measure the slope, for slope
333 T=self.config['temperature']
335 fit_points=int(self.config['auto_fit_points'])
336 slope_span=int(self.config['auto_slope_span'])
339 rebase=False #if true=we select rebase
342 displayed_plot=self._get_displayed_plot(0)
344 if self.current.curve.experiment != 'smfs' or len(displayed_plot.vectors) < 2:
345 print 'Cannot work on this curve.'
348 #Look for custom persistent length / custom temperature
349 for arg in args.split():
350 #look for a persistent length argument.
352 pl_expression=arg.split('=')
353 pl_value=float(pl_expression[1]) #actual value
354 #look for a T argument. FIXME: spaces are not allowed between 'pl' and value
355 if ('t=' in arg[0:2]) or ('T=' in arg[0:2]):
356 t_expression=arg.split('=')
357 T=float(t_expression[1])
359 #Handle contact point arguments
360 #FIXME: code duplication
361 if 'reclick' in args.split():
362 print 'Click contact point'
363 contact_point=self._measure_N_points(N=1, whatset=1)[0]
364 contact_point_index=contact_point.index
365 self.wlccontact_point=contact_point
366 self.wlccontact_index=contact_point.index
367 self.wlccurrent=self.current.path
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=self._measure_N_points(N=1, whatset=1)[0]
372 contact_point_index=contact_point.index
373 self.wlccontact_point=contact_point
374 self.wlccontact_index=contact_point.index
375 self.wlccurrent=self.current.path
377 contact_point=self.wlccontact_point
378 contact_point_index=self.wlccontact_index
380 #Automatically find contact point
381 cindex=self.find_contact_point()
382 contact_point=ClickedPoint()
383 contact_point.absolute_coords=displayed_plot.vectors[1][0][cindex], displayed_plot.vectors[1][1][cindex]
384 contact_point.find_graph_coords(displayed_plot.vectors[1][0], displayed_plot.vectors[1][1])
385 contact_point.is_marker=True
388 #Pick up force baseline
389 whatset=1 #fixme: for all sets
390 if 'rebase' in args or (self.basecurrent != self.current.path):
393 print 'Select baseline'
394 self.basepoints=self._measure_N_points(N=2, whatset=whatset)
395 self.basecurrent=self.current.path
396 boundaries=[self.basepoints[0].index, self.basepoints[1].index]
398 to_average=displayed_plot.vectors[1][1][boundaries[0]:boundaries[1]] #y points to average
399 avg=np.mean(to_average)
402 defplot=self.current.curve.default_plots()[0]
403 flatten=self._find_plotmanip('flatten') #Extract flatten plotmanip
404 defplot=flatten(defplot, self.current, customvalue=1) #Flatten curve before feeding it to has_peaks
405 peak_location,peak_size=self.has_peaks(defplot, self.convfilt_config['mindeviation'])
407 #Create a new plot to send
408 fitplot=copy.deepcopy(displayed_plot)
410 #Initialize data vectors
416 #Cycle between peaks and do analysis.
417 for peak in peak_location:
419 #FIXME: clean wlc fitting, to avoid this clickedpoint hell
420 #-create a clicked point for the peak point
421 peak_point=ClickedPoint()
422 peak_point.absolute_coords=displayed_plot.vectors[1][0][peak], displayed_plot.vectors[1][1][peak]
423 peak_point.find_graph_coords(displayed_plot.vectors[1][0], displayed_plot.vectors[1][1])
424 #-create a clicked point for the other fit point
425 other_fit_point=ClickedPoint()
426 other_fit_point.absolute_coords=displayed_plot.vectors[1][0][peak-fit_points], displayed_plot.vectors[1][1][peak-fit_points]
427 other_fit_point.find_graph_coords(displayed_plot.vectors[1][0], displayed_plot.vectors[1][1])
429 points=[contact_point, peak_point, other_fit_point]
430 params, yfit, xfit = self.wlc_fit(points, displayed_plot.vectors[1][0], displayed_plot.vectors[1][1],pl_value,T)
431 #save wlc values (nm)
432 c_lengths.append(params[0]*(1.0e+9))
433 if len(params)==2: #if we did choose 2-value fit
434 p_lengths.append(params[1]*(1.0e+9))
436 p_lengths.append(pl_value)
437 #Add WLC fit lines to plot
438 fitplot.add_set(xfit,yfit)
439 if len(fitplot.styles)==0:
442 fitplot.styles.append(None)
445 delta_to_measure=displayed_plot.vectors[1][1][peak-delta_force:peak+delta_force]
446 y=min(delta_to_measure)
447 #save force values (pN)
448 forces.append(abs(y-avg)*(1.0e+12))
451 slope=self.linefit_between(peak-slope_span,peak)[0]
454 #Show wlc fits and peak locations
455 self._send_plot([fitplot])
458 #Ask the user what peaks to ignore from analysis.
459 print 'Peaks to ignore (0,1...n from contact point,return to take all)'
460 print 'N to discard measurement'
461 exclude_raw=raw_input('Input:')
465 if not exclude_raw=='':
466 exclude=exclude_raw.split(',')
468 exclude=[int(item) for item in exclude]
475 print 'Bad input, taking all...'
476 #Clean data vectors from ignored peaks
477 c_lengths=[item for item in c_lengths if item != None]
478 p_lengths=[item for item in p_lengths if item != None]
479 forces=[item for item in forces if item != None]
480 slopes=[item for item in slopes if item != None]
481 print 'contour (nm)',c_lengths
482 print 'p (nm)',p_lengths
483 print 'forces (pN)',forces
484 print 'slopes (N/m)',slopes
487 if self.autofile=='':
488 self.autofile=raw_input('Filename? (return to ignore) ')
489 if self.autofile=='':
493 if not os.path.exists(self.autofile):
494 f=open(self.autofile,'w+')
495 f.write('Analysis started '+time.asctime()+'\n')
496 f.write('----------------------------------------\n')
497 f.write('; Contour length (nm) ; Persistence length (nm) ; Max.Force (pN) ; Slope (N/m) \n')
501 f=open(self.autofile,'a+')
503 f.write(self.current.path+'\n')
504 for i in range(len(c_lengths)):
505 f.write(' ; '+str(c_lengths[i])+' ; '+str(p_lengths[i])+' ; '+str(forces[i])+' ; '+str(slopes[i])+'\n')
508 self.do_note('autopeak')