1 # Copyright (C) 2008-2010 Alberto Gomez-Casado
4 # Massimo Sandal <devicerandom@gmail.com>
5 # W. Trevor King <wking@drexel.edu>
7 # This file is part of Hooke.
9 # Hooke is free software: you can redistribute it and/or
10 # modify it under the terms of the GNU Lesser General Public
11 # License as published by the Free Software Foundation, either
12 # version 3 of the License, or (at your option) any later version.
14 # Hooke is distributed in the hope that it will be useful,
15 # but WITHOUT ANY WARRANTY; without even the implied warranty of
16 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 # GNU Lesser General Public License for more details.
19 # You should have received a copy of the GNU Lesser General Public
20 # License along with Hooke. If not, see
21 # <http://www.gnu.org/licenses/>.
23 """The ``vclamp`` module provides :class:`VelocityClampPlugin` and
24 several associated :class:`hooke.command.Command`\s for handling
25 common velocity clamp analysis tasks.
30 from ..command import Command, Argument, Failure
31 from ..curve import Data
32 from ..plugin import Builtin
38 for i,block in enumerate(curve.data):
39 data = Data((block.shape[0], block.shape[1]+2), dtype=block.dtype)
40 data.info = copy.deepcopy(block.info)
42 data.info['columns'].extend(['surface z piezo (m)', 'deflection (N)'])
43 z_data = data[:,data.info['columns'].index('z piezo (m)')]
44 d_data = data[:,data.info['columns'].index('deflection (m)')]
45 surface_offset = min(z_data) # TODO
46 data.info['surface offset (m)'] = surface_offset
47 data[:,-2] = z_data - surface_offset
48 data[:,-1] = d_data * data.info['spring constant (N/m)']
53 class VelocityClampPlugin (Builtin):
55 super(VelocityClampPlugin, self).__init__(name='vclamp')
57 # NextCommand(self), PreviousCommand(self), JumpCommand(self),
60 class generalvclampCommands(object):
62 def do_subtplot(self, args):
66 Plots the difference between ret and ext current curve
70 #FIXME: sub_filter and sub_order must be args
72 if len(self.plots[0].vectors) != 2:
73 print 'This command only works on a curve with two different plots.'
76 outplot=self.subtract_curves(sub_order=1)
78 plot_graph=self.list_of_events['plot_graph']
79 wx.PostEvent(self.frame,plot_graph(plots=[outplot]))
86 def do_distance(self,args):
90 Measure the distance (in nm) between two points.
91 For a standard experiment this is the delta X distance.
92 For a force clamp experiment this is the delta Y distance (actually becomes
97 if self.current.curve.experiment == 'clamp':
98 print 'You wanted to use zpiezo perhaps?'
101 dx,unitx,dy,unity=self._delta(set=1)
102 print str(dx*(10**9))+' nm'
103 to_dump='distance '+self.current.path+' '+str(dx*(10**9))+' nm'
104 self.outlet.push(to_dump)
107 def do_force(self,args):
111 Measure the force difference (in pN) between two points
115 if self.current.curve.experiment == 'clamp':
116 print 'This command makes no sense for a force clamp experiment.'
118 dx,unitx,dy,unity=self._delta(set=1)
119 print str(dy*(10**12))+' pN'
120 to_dump='force '+self.current.path+' '+str(dy*(10**12))+' pN'
121 self.outlet.push(to_dump)
124 def do_forcebase(self,args):
128 Measures the difference in force (in pN) between a point and a baseline
129 took as the average between two points.
131 The baseline is fixed once for a given curve and different force measurements,
132 unless the user wants it to be recalculated
134 Syntax: forcebase [rebase]
135 rebase: Forces forcebase to ask again the baseline
136 max: Instead of asking for a point to measure, asks for two points and use
137 the maximum peak in between
139 rebase=False #if true=we select rebase
140 maxpoint=False #if true=we measure the maximum peak
142 plot=self._get_displayed_plot()
143 whatset=1 #fixme: for all sets
144 if 'rebase' in args or (self.basecurrent != self.current.path):
150 print 'Select baseline'
151 self.basepoints=self._measure_N_points(N=2, whatset=whatset)
152 self.basecurrent=self.current.path
155 print 'Select two points'
156 points=self._measure_N_points(N=2, whatset=whatset)
157 boundpoints=[points[0].index, points[1].index]
160 y=min(plot.vectors[whatset][1][boundpoints[0]:boundpoints[1]])
162 print 'Chosen interval not valid. Try picking it again. Did you pick the same point as begin and end of interval?'
164 print 'Select point to measure'
165 points=self._measure_N_points(N=1, whatset=whatset)
166 #whatplot=points[0].dest
167 y=points[0].graph_coords[1]
169 #fixme: code duplication
170 boundaries=[self.basepoints[0].index, self.basepoints[1].index]
172 to_average=plot.vectors[whatset][1][boundaries[0]:boundaries[1]] #y points to average
174 avg=np.mean(to_average)
176 print str(forcebase*(10**12))+' pN'
177 to_dump='forcebase '+self.current.path+' '+str(forcebase*(10**12))+' pN'
178 self.outlet.push(to_dump)
180 def plotmanip_multiplier(self, plot, current):
182 Multiplies all the Y values of an SMFS curve by a value stored in the 'force_multiplier'
183 configuration variable. Useful for calibrations and other stuff.
187 if current.curve.experiment != 'smfs':
190 #only one set is present...
191 if len(self.plots[0].vectors) != 2:
195 if (self.config['force_multiplier']==1):
198 for i in range(len(plot.vectors[0][1])):
199 plot.vectors[0][1][i]=plot.vectors[0][1][i]*self.config['force_multiplier']
201 for i in range(len(plot.vectors[1][1])):
202 plot.vectors[1][1][i]=plot.vectors[1][1][i]*self.config['force_multiplier']
207 def plotmanip_flatten(self, plot, current, customvalue=False):
209 Subtracts a polynomial fit to the non-contact part of the curve, as to flatten it.
210 the best polynomial fit is chosen among polynomials of degree 1 to n, where n is
211 given by the configuration file or by the customvalue.
213 customvalue= int (>0) --> starts the function even if config says no (default=False)
217 if current.curve.experiment != 'smfs':
220 #only one set is present...
221 if len(self.plots[0].vectors) != 2:
224 #config is not flatten, and customvalue flag is false too
225 if (not self.config['flatten']) and (not customvalue):
232 max_cycles=customvalue
234 max_cycles=self.config['flatten'] #Using > 1 usually doesn't help and can give artefacts. However, it could be useful too.
236 contact_index=self.find_contact_point()
238 valn=[[] for item in range(max_exponent)]
239 yrn=[0.0 for item in range(max_exponent)]
240 errn=[0.0 for item in range(max_exponent)]
242 #Check if we have a proper numerical value
246 #Loudly and annoyingly complain if it's not a number, then fallback to zero
247 print '''Warning: flatten value is not a number!
248 Use "set flatten" or edit hooke.conf to set it properly
252 for i in range(int(max_cycles)):
254 x_ext=plot.vectors[0][0][contact_index+delta_contact:]
255 y_ext=plot.vectors[0][1][contact_index+delta_contact:]
256 x_ret=plot.vectors[1][0][contact_index+delta_contact:]
257 y_ret=plot.vectors[1][1][contact_index+delta_contact:]
258 for exponent in range(max_exponent):
260 valn[exponent]=sp.polyfit(x_ext,y_ext,exponent)
261 yrn[exponent]=sp.polyval(valn[exponent],x_ret)
262 errn[exponent]=sp.sqrt(sum((yrn[exponent]-y_ext)**2)/float(len(y_ext)))
264 print 'Cannot flatten!'
268 best_exponent=errn.index(min(errn))
271 ycorr_ext=y_ext-yrn[best_exponent]+y_ext[0] #noncontact part
272 yjoin_ext=np.array(plot.vectors[0][1][0:contact_index+delta_contact]) #contact part
274 ycorr_ret=y_ret-yrn[best_exponent]+y_ext[0] #noncontact part
275 yjoin_ret=np.array(plot.vectors[1][1][0:contact_index+delta_contact]) #contact part
277 ycorr_ext=np.concatenate((yjoin_ext, ycorr_ext))
278 ycorr_ret=np.concatenate((yjoin_ret, ycorr_ret))
280 plot.vectors[0][1]=list(ycorr_ext)
281 plot.vectors[1][1]=list(ycorr_ret)
286 def do_slope(self,args):
290 Measures the slope of a delimited chunk on the return trace.
291 The chunk can be delimited either by two manual clicks, or have
292 a fixed width, given as an argument.
294 Syntax: slope [width]
295 The facultative [width] parameter specifies how many
296 points will be considered for the fit. If [width] is
297 specified, only one click will be required.
298 (c) Marco Brucale, Massimo Sandal 2008
301 # Reads the facultative width argument
307 # Decides between the two forms of user input, as per (args)
309 # Gets the Xs of two clicked points as indexes on the current curve vector
310 print 'Click twice to delimit chunk'
311 points=self._measure_N_points(N=2,whatset=1)
313 print 'Click once on the leftmost point of the chunk (i.e.usually the peak)'
314 points=self._measure_N_points(N=1,whatset=1)
316 slope=self._slope(points,fitspan)
318 # Outputs the relevant slope parameter
321 to_dump='slope '+self.current.path+' '+str(slope)
322 self.outlet.push(to_dump)
324 def _slope(self,points,fitspan):
325 # Calls the function linefit_between
326 parameters=[0,0,[],[]]
328 clickedpoints=[points[0].index,points[1].index]
331 clickedpoints=[points[0].index-fitspan,points[0].index]
334 parameters=self.linefit_between(clickedpoints[0],clickedpoints[1])
336 print 'Cannot fit. Did you click twice the same point?'
339 # Outputs the relevant slope parameter
341 print str(parameters[0])
342 to_dump='slope '+self.curve.path+' '+str(parameters[0])
343 self.outlet.push(to_dump)
345 # Makes a vector with the fitted parameters and sends it to the GUI
346 xtoplot=parameters[2]
350 ytoplot.append((x*parameters[0])+parameters[1])
352 clickvector_x, clickvector_y=[], []
354 clickvector_x.append(item.graph_coords[0])
355 clickvector_y.append(item.graph_coords[1])
357 lineplot=self._get_displayed_plot(0) #get topmost displayed plot
359 lineplot.add_set(xtoplot,ytoplot)
360 lineplot.add_set(clickvector_x, clickvector_y)
363 if lineplot.styles==[]:
364 lineplot.styles=[None,None,None,'scatter']
366 lineplot.styles+=[None,'scatter']
367 if lineplot.colors==[]:
368 lineplot.colors=[None,None,'black',None]
370 lineplot.colors+=['black',None]
373 self._send_plot([lineplot])
378 def linefit_between(self,index1,index2,whatset=1):
380 Creates two vectors (xtofit,ytofit) slicing out from the
381 current return trace a portion delimited by the two indexes
383 Then does a least squares linear fit on that slice.
384 Finally returns [0]=the slope, [1]=the intercept of the
385 fitted 1st grade polynomial, and [2,3]=the actual (x,y) vectors
387 (c) Marco Brucale, Massimo Sandal 2008
389 # Translates the indexes into two vectors containing the x,y data to fit
390 xtofit=self.plots[0].vectors[whatset][0][index1:index2]
391 ytofit=self.plots[0].vectors[whatset][1][index1:index2]
393 # Does the actual linear fitting (simple least squares with numpy.polyfit)
395 linefit=np.polyfit(xtofit,ytofit,1)
397 return (linefit[0],linefit[1],xtofit,ytofit)
400 def fit_interval_nm(self,start_index,plot,nm,backwards):
402 Calculates the number of points to fit, given a fit interval in nm
403 start_index: index of point
405 backwards: if true, finds a point backwards.
407 whatset=1 #FIXME: should be decidable
408 x_vect=plot.vectors[1][0]
412 start=x_vect[start_index]
414 while abs(x_vect[i]-x_vect[start_index])*(10**9) < nm:
415 if i==0 or i==maxlen-1: #we reached boundaries of vector!
427 def find_current_peaks(self,noflatten, a=True, maxpeak=True):
430 a=self.convfilt_config['mindeviation']
434 print "Bad input, using default."
435 abs_devs=self.convfilt_config['mindeviation']
437 defplot=self.current.curve.default_plots()[0]
439 flatten=self._find_plotmanip('flatten') #Extract flatten plotmanip
440 defplot=flatten(defplot, self.current, customvalue=1) #Flatten curve before feeding it to has_peaks
441 pk_location,peak_size=self.has_peaks(defplot, abs_devs, maxpeak)
442 return pk_location, peak_size
445 def pickup_contact_point(self,N=1,whatset=1):
446 '''macro to pick up the contact point by clicking'''
447 contact_point=self._measure_N_points(N=1, whatset=1)[0]
448 contact_point_index=contact_point.index
449 self.wlccontact_point=contact_point
450 self.wlccontact_index=contact_point.index
451 self.wlccurrent=self.current.path
452 return contact_point, contact_point_index
456 def baseline_points(self,peak_location, displayed_plot):
457 clicks=self.config['baseline_clicks']
460 base_index_0=peak_location[-1]+self.fit_interval_nm(peak_location[-1], displayed_plot, self.config['auto_right_baseline'],False)
461 self.basepoints.append(self._clickize(displayed_plot.vectors[1][0],displayed_plot.vectors[1][1],base_index_0))
462 base_index_1=self.basepoints[0].index+self.fit_interval_nm(self.basepoints[0].index, displayed_plot, self.config['auto_left_baseline'],False)
463 self.basepoints.append(self._clickize(displayed_plot.vectors[1][0],displayed_plot.vectors[1][1],base_index_1))
465 print 'Select baseline'
467 self.basepoints=self._measure_N_points(N=1, whatset=1)
468 base_index_1=self.basepoints[0].index+self.fit_interval_nm(self.basepoints[0].index, displayed_plot, self.config['auto_left_baseline'], False)
469 self.basepoints.append(self._clickize(displayed_plot.vectors[1][0],displayed_plot.vectors[1][1],base_index_1))
471 self.basepoints=self._measure_N_points(N=2, whatset=1)
473 self.basecurrent=self.current.path
474 return self.basepoints