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 modify it
10 # under the terms of the GNU Lesser General Public License as
11 # published by the Free Software Foundation, either version 3 of the
12 # License, or (at your option) any later version.
14 # Hooke is distributed in the hope that it will be useful, but WITHOUT
15 # ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
16 # or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General
17 # 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.
33 from ..command import Command, Argument, Failure, NullQueue
34 from ..config import Setting
35 from ..curve import Data
36 from ..plugin import Plugin
37 from ..util.fit import PoorFit, ModelFitter
38 from .curve import CurveArgument
41 def scale(hooke, curve, block=None):
42 """Run 'add block force array' on `block`.
44 Runs 'zero block surface contact point' first, if necessary. Does
45 not run either command if the columns they add to the block are
48 If `block` is `None`, scale all blocks in `curve`.
50 commands = hooke.commands
51 contact = [c for c in hooke.commands
52 if c.name == 'zero block surface contact point'][0]
53 force = [c for c in hooke.commands if c.name == 'add block force array'][0]
55 outqueue = NullQueue()
57 for i in range(len(curve.data)):
58 scale(hooke, curve, block=i)
60 params = {'curve':curve, 'block':block}
62 if ('surface distance (m)' not in b.info.columns
63 or 'surface adjusted deflection (m)' not in b.info.columns):
65 contact._run(hooke, inqueue, outqueue, params)
67 raise PoorFit('Could not fit %s %s: %s'
68 % (curve.path, block, str(e)))
69 if ('deflection (N)' not in b.info.columns):
70 force._run(hooke, inqueue, outqueue, params)
73 class SurfacePositionModel (ModelFitter):
76 The bilinear model is symmetric, but the parameter guessing and
77 sanity checks assume the contact region occurs for lower indicies
78 ("left of") the non-contact region. We also assume that
79 tip-surface attractions produce positive deflections.
83 Algorithm borrowed from WTK's `piezo package`_, specifically
84 from :func:`piezo.z_piezo_utils.analyzeSurfPosData`.
87 http://www.physics.drexel.edu/~wking/code/git/git.php?p=piezo.git
89 Fits the data to the bilinear :method:`model`.
91 In order for this model to produce a satisfactory fit, there
92 should be enough data in the off-surface region that interactions
93 due to proteins, etc. will not seriously skew the fit in the
99 def model(self, params):
100 """A continuous, bilinear model.
107 p_0 + p_1 x & \text{if $x <= p_2$}, \\
108 p_0 + p_1 p_2 + p_3 (x-p_2) & \text{if $x >= p_2$}.
111 Where :math:`p_0` is a vertical offset, :math:`p_1` is the slope
112 of the first region, :math:`p_2` is the transition location, and
113 :math:`p_3` is the slope of the second region.
115 p = params # convenient alias
116 if self.info.get('force zero non-contact slope', None) == True:
118 p.append(0.) # restore the non-contact slope parameter
119 r2 = numpy.round(abs(p[2]))
121 self._model_data[:r2] = p[0] + p[1] * numpy.arange(r2)
122 if r2 < len(self._data)-1:
123 self._model_data[r2:] = \
124 p[0] + p[1]*p[2] + p[3] * numpy.arange(len(self._data)-r2)
125 return self._model_data
127 def set_data(self, data, info=None):
128 super(SurfacePositionModel, self).set_data(data, info)
132 self.info['min position'] = 0
133 self.info['max position'] = len(data)
134 self.info['max deflection'] = data.max()
135 self.info['min deflection'] = data.min()
136 self.info['position range'] = self.info['max position'] - self.info['min position']
137 self.info['deflection range'] = self.info['max deflection'] - self.info['min deflection']
139 def guess_initial_params(self, outqueue=None):
140 """Guess the initial parameters.
144 We guess initial parameters such that the offset (:math:`p_1`)
145 matches the minimum deflection, the kink (:math:`p_2`) occurs in
146 the middle of the data, the initial (contact) slope (:math:`p_0`)
147 produces the maximum deflection at the left-most point, and the
148 final (non-contact) slope (:math:`p_3`) is zero.
150 left_offset = self.info['min deflection']
151 left_slope = 2*(self.info['deflection range']
152 /self.info['position range'])
153 kink_position = (self.info['max position']
154 +self.info['min position'])/2.0
156 self.info['guessed contact slope'] = left_slope
157 params = [left_offset, left_slope, kink_position, right_slope]
158 if self.info.get('force zero non-contact slope', None) == True:
162 def guess_scale(self, params, outqueue=None):
163 """Guess the parameter scales.
168 We guess offset scale (:math:`p_0`) as one tenth of the total
169 deflection range, the kink scale (:math:`p_2`) as one tenth of
170 the total index range, the initial (contact) slope scale
171 (:math:`p_1`) as one tenth of the contact slope estimation,
172 and the final (non-contact) slope scale (:math:`p_3`) is as
173 one tenth of the initial slope scale.
175 offset_scale = self.info['deflection range']/10.
176 left_slope_scale = abs(params[1])/10.
177 kink_scale = self.info['position range']/10.
178 right_slope_scale = left_slope_scale/10.
179 scale = [offset_scale, left_slope_scale, kink_scale, right_slope_scale]
180 if self.info.get('force zero non-contact slope', None) == True:
184 def fit(self, *args, **kwargs):
185 self.info['guessed contact slope'] = None
186 params = super(SurfacePositionModel, self).fit(*args, **kwargs)
187 params[2] = abs(params[2])
188 if self.info.get('force zero non-contact slope', None) == True:
189 params = list(params)
190 params.append(0.) # restore the non-contact slope parameter
192 # check that the fit is reasonable, see the :meth:`model` docstring
193 # for parameter descriptions. HACK: hardcoded cutoffs.
194 if abs(params[3]*10) > abs(params[1]) :
195 raise PoorFit('Slope in non-contact region, or no slope in contact')
196 if params[2] < self.info['min position']+0.02*self.info['position range']:
198 'No kink (kink %g less than %g, need more space to left)'
200 self.info['min position']+0.02*self.info['position range']))
201 if params[2] > self.info['max position']-0.02*self.info['position range']:
203 'No kink (kink %g more than %g, need more space to right)'
205 self.info['max position']-0.02*self.info['position range']))
206 if (self.info['guessed contact slope'] != None
207 and abs(params[1]) < 0.5 * abs(self.info['guessed contact slope'])):
208 raise PoorFit('Too far (contact slope %g, but expected ~%g'
209 % (params[3], self.info['guessed contact slope']))
212 class VelocityClampPlugin (Plugin):
214 super(VelocityClampPlugin, self).__init__(name='vclamp')
216 SurfaceContactCommand(self), ForceCommand(self),
219 def default_settings(self):
221 Setting(section=self.setting_section, help=self.__doc__),
222 Setting(section=self.setting_section,
223 option='surface contact point algorithm',
225 help='Select the surface contact point algorithm. See the documentation for descriptions of available algorithms.')
229 class SurfaceContactCommand (Command):
230 """Automatically determine a block's surface contact point.
232 Uses the block's `z piezo (m)` and `deflection (m)` arrays.
233 Stores the contact parameters in `block.info`'s `surface distance
234 offset (m)` and `surface deflection offset (m)`. Model-specific
235 fitting information is stored in `surface detection parameters`.
237 The adjusted data columns `surface distance (m)` and `surface
238 adjusted deflection (m)` are also added to the block.
240 You can select the contact point algorithm with the creatively
241 named `surface contact point algorithm` configuration setting.
242 Currently available options are:
244 * fmms (:meth:`find_contact_point_fmms`)
245 * ms (:meth:`find_contact_point_ms`)
246 * wtk (:meth:`find_contact_point_wtk`)
248 def __init__(self, plugin):
249 super(SurfaceContactCommand, self).__init__(
250 name='zero block surface contact point',
253 Argument(name='block', aliases=['set'], type='int', default=0,
255 Data block for which the force should be calculated. For an
256 approach/retract force curve, `0` selects the approaching curve and `1`
257 selects the retracting curve.
260 help=self.__doc__, plugin=plugin)
262 def _run(self, hooke, inqueue, outqueue, params):
263 data = params['curve'].data[params['block']]
264 # HACK? rely on params['curve'] being bound to the local hooke
265 # playlist (i.e. not a copy, as you would get by passing a
266 # curve through the queue). Ugh. Stupid queues. As an
267 # alternative, we could pass lookup information through the
269 new = Data((data.shape[0], data.shape[1]+2), dtype=data.dtype)
270 new.info = copy.deepcopy(data.info)
272 new.info['columns'].extend(
273 ['surface distance (m)', 'surface adjusted deflection (m)'])
274 z_data = data[:,data.info['columns'].index('z piezo (m)')]
275 d_data = data[:,data.info['columns'].index('deflection (m)')]
276 i,deflection_offset,ps = self.find_contact_point(
277 params['curve'], z_data, d_data, outqueue)
278 surface_offset = z_data[i]
279 new.info['surface distance offset (m)'] = surface_offset
280 new.info['surface deflection offset (m)'] = deflection_offset
281 new.info['surface detection parameters'] = ps
282 new[:,-2] = z_data - surface_offset
283 new[:,-1] = d_data - deflection_offset
284 params['curve'].data[params['block']] = new
286 def find_contact_point(self, curve, z_data, d_data, outqueue=None):
287 """Railyard for the `find_contact_point_*` family.
289 Uses the `surface contact point algorithm` configuration
290 setting to call the appropriate backend algorithm.
292 fn = getattr(self, 'find_contact_point_%s'
293 % self.plugin.config['surface contact point algorithm'])
294 return fn(curve, z_data, d_data, outqueue)
296 def find_contact_point_fmms(self, curve, z_data, d_data, outqueue=None):
297 """Algorithm by Francesco Musiani and Massimo Sandal.
303 0) Driver-specific workarounds, e.g. deal with the PicoForce
304 trigger bug by excluding retraction portions with excessive
306 1) Select the second half (non-contact side) of the retraction
308 2) Fit the selection to a line.
309 3) If the fit is not almost horizontal, halve the selection
311 4) Average the selection and use it as a baseline.
312 5) Slide in from the start (contact side) of the retraction
313 curve, until you find a point with greater than baseline
314 deflection. That point is the contact point.
316 if curve.info['filetype'] == 'picoforce':
317 # Take care of the picoforce trigger bug (TODO: example
318 # data file demonstrating the bug). We exclude portions
319 # of the curve that have too much standard deviation.
320 # Yes, a lot of magic is here.
321 check_start = len(d_data)-len(d_data)/20
322 monster_start = len(d_data)
324 # look at the non-contact tail
325 non_monster = d_data[check_start:monster_start]
326 if non_monster.std() < 2e-10: # HACK: hardcoded cutoff
328 else: # move further away from the monster
329 check_start -= len(d_data)/50
330 monster_start -= len(d_data)/50
331 z_data = z_data[:monster_start]
332 d_data = d_data[:monster_start]
334 # take half of the thing to start
335 selection_start = len(d_data)/2
337 z_chunk = z_data[selection_start:]
338 d_chunk = d_data[selection_start:]
339 slope,intercept,r,two_tailed_prob,stderr_of_the_estimate = \
340 scipy.stats.linregress(z_chunk, d_chunk)
341 # We stop if we found an almost-horizontal fit or if we're
342 # getting to small a selection. FIXME: 0.1 and 5./6 here
343 # are "magic numbers" (although reasonable)
344 if (abs(slope) < 0.1 # deflection (m) / surface (m)
345 or selection_start > 5./6*len(d_data)):
347 selection_start += 10
349 d_baseline = d_chunk.mean()
351 # find the first point above the calculated baseline
353 while i < len(d_data) and d_data[i] < ymean:
355 return (i, d_baseline, {})
357 def find_contact_point_ms(self, curve, z_data, d_data, outqueue=None):
358 """Algorithm by Massimo Sandal.
362 WTK: At least the commits are by Massimo, and I see no notes
363 attributing the algorithm to anyone else.
369 xext=raw_plot.vectors[0][0]
370 yext=raw_plot.vectors[0][1]
371 xret2=raw_plot.vectors[1][0]
372 yret=raw_plot.vectors[1][1]
374 first_point=[xext[0], yext[0]]
375 last_point=[xext[-1], yext[-1]]
377 #regr=scipy.polyfit(first_point, last_point,1)[0:2]
378 diffx=abs(first_point[0]-last_point[0])
379 diffy=abs(first_point[1]-last_point[1])
381 #using polyfit results in numerical errors. good old algebra.
383 b=first_point[1]-(a*first_point[0])
384 baseline=scipy.polyval((a,b), xext)
386 ysub=[item-basitem for item,basitem in zip(yext,baseline)]
388 contact=ysub.index(min(ysub))
390 return xext,ysub,contact
392 #now, exploit a ClickedPoint instance to calculate index...
394 dummy.absolute_coords=(x_intercept,y_intercept)
395 dummy.find_graph_coords(xret2,yret)
398 return dummy.index, regr, regr_contact
402 def find_contact_point_wtk(self, curve, z_data, d_data, outqueue=None):
403 """Algorithm by W. Trevor King.
407 Uses :func:`analyze_surf_pos_data_wtk` internally.
409 reverse = z_data[0] > z_data[-1]
410 if reverse == True: # approaching, contact region on the right
411 d_data = d_data[::-1]
412 s = SurfacePositionModel(d_data)
413 s.info['force zero non-contact slope'] = True
414 offset,contact_slope,surface_index,non_contact_slope = s.fit(
418 'contact slope': contact_slope,
419 'surface index': surface_index,
420 'non-contact slope': non_contact_slope,
423 deflection_offset = offset + contact_slope*surface_index,
425 surface_index = len(d_data)-1-surface_index
426 return (numpy.round(surface_index), deflection_offset, info)
428 class ForceCommand (Command):
429 """Calculate a block's `deflection (N)` array.
431 Uses the block's `deflection (m)` array and `spring constant
434 def __init__(self, plugin):
435 super(ForceCommand, self).__init__(
436 name='add block force array',
439 Argument(name='block', aliases=['set'], type='int', default=0,
441 Data block for which the force should be calculated. For an
442 approach/retract force curve, `0` selects the approaching curve and `1`
443 selects the retracting curve.
446 help=self.__doc__, plugin=plugin)
448 def _run(self, hooke, inqueue, outqueue, params):
449 data = params['curve'].data[params['block']]
450 # HACK? rely on params['curve'] being bound to the local hooke
451 # playlist (i.e. not a copy, as you would get by passing a
452 # curve through the queue). Ugh. Stupid queues. As an
453 # alternative, we could pass lookup information through the
455 new = Data((data.shape[0], data.shape[1]+1), dtype=data.dtype)
456 new.info = copy.deepcopy(data.info)
458 new.info['columns'].append('deflection (N)')
459 d_data = data[:,data.info['columns'].index('surface adjusted deflection (m)')]
460 new[:,-1] = d_data * data.info['spring constant (N/m)']
461 params['curve'].data[params['block']] = new
464 class generalvclampCommands(object):
466 def do_subtplot(self, args):
469 (procplots.py plugin)
470 Plots the difference between ret and ext current curve
474 #FIXME: sub_filter and sub_order must be args
476 if len(self.plots[0].vectors) != 2:
477 print 'This command only works on a curve with two different plots.'
480 outplot=self.subtract_curves(sub_order=1)
482 plot_graph=self.list_of_events['plot_graph']
483 wx.PostEvent(self.frame,plot_graph(plots=[outplot]))
485 def _plug_init(self):
486 self.basecurrent=None
490 def do_distance(self,args):
494 Measure the distance (in nm) between two points.
495 For a standard experiment this is the delta X distance.
496 For a force clamp experiment this is the delta Y distance (actually becomes
501 if self.current.curve.experiment == 'clamp':
502 print 'You wanted to use zpiezo perhaps?'
505 dx,unitx,dy,unity=self._delta(set=1)
506 print str(dx*(10**9))+' nm'
507 to_dump='distance '+self.current.path+' '+str(dx*(10**9))+' nm'
508 self.outlet.push(to_dump)
511 def do_force(self,args):
515 Measure the force difference (in pN) between two points
519 if self.current.curve.experiment == 'clamp':
520 print 'This command makes no sense for a force clamp experiment.'
522 dx,unitx,dy,unity=self._delta(set=1)
523 print str(dy*(10**12))+' pN'
524 to_dump='force '+self.current.path+' '+str(dy*(10**12))+' pN'
525 self.outlet.push(to_dump)
528 def do_forcebase(self,args):
532 Measures the difference in force (in pN) between a point and a baseline
533 took as the average between two points.
535 The baseline is fixed once for a given curve and different force measurements,
536 unless the user wants it to be recalculated
538 Syntax: forcebase [rebase]
539 rebase: Forces forcebase to ask again the baseline
540 max: Instead of asking for a point to measure, asks for two points and use
541 the maximum peak in between
543 rebase=False #if true=we select rebase
544 maxpoint=False #if true=we measure the maximum peak
546 plot=self._get_displayed_plot()
547 whatset=1 #fixme: for all sets
548 if 'rebase' in args or (self.basecurrent != self.current.path):
554 print 'Select baseline'
555 self.basepoints=self._measure_N_points(N=2, whatset=whatset)
556 self.basecurrent=self.current.path
559 print 'Select two points'
560 points=self._measure_N_points(N=2, whatset=whatset)
561 boundpoints=[points[0].index, points[1].index]
564 y=min(plot.vectors[whatset][1][boundpoints[0]:boundpoints[1]])
566 print 'Chosen interval not valid. Try picking it again. Did you pick the same point as begin and end of interval?'
568 print 'Select point to measure'
569 points=self._measure_N_points(N=1, whatset=whatset)
570 #whatplot=points[0].dest
571 y=points[0].graph_coords[1]
573 #fixme: code duplication
574 boundaries=[self.basepoints[0].index, self.basepoints[1].index]
576 to_average=plot.vectors[whatset][1][boundaries[0]:boundaries[1]] #y points to average
578 avg=np.mean(to_average)
580 print str(forcebase*(10**12))+' pN'
581 to_dump='forcebase '+self.current.path+' '+str(forcebase*(10**12))+' pN'
582 self.outlet.push(to_dump)
584 def plotmanip_multiplier(self, plot, current):
586 Multiplies all the Y values of an SMFS curve by a value stored in the 'force_multiplier'
587 configuration variable. Useful for calibrations and other stuff.
591 if current.curve.experiment != 'smfs':
594 #only one set is present...
595 if len(self.plots[0].vectors) != 2:
599 if (self.config['force_multiplier']==1):
602 for i in range(len(plot.vectors[0][1])):
603 plot.vectors[0][1][i]=plot.vectors[0][1][i]*self.config['force_multiplier']
605 for i in range(len(plot.vectors[1][1])):
606 plot.vectors[1][1][i]=plot.vectors[1][1][i]*self.config['force_multiplier']
611 def plotmanip_flatten(self, plot, current, customvalue=False):
613 Subtracts a polynomial fit to the non-contact part of the curve, as to flatten it.
614 the best polynomial fit is chosen among polynomials of degree 1 to n, where n is
615 given by the configuration file or by the customvalue.
617 customvalue= int (>0) --> starts the function even if config says no (default=False)
621 if current.curve.experiment != 'smfs':
624 #only one set is present...
625 if len(self.plots[0].vectors) != 2:
628 #config is not flatten, and customvalue flag is false too
629 if (not self.config['flatten']) and (not customvalue):
636 max_cycles=customvalue
638 max_cycles=self.config['flatten'] #Using > 1 usually doesn't help and can give artefacts. However, it could be useful too.
640 contact_index=self.find_contact_point()
642 valn=[[] for item in range(max_exponent)]
643 yrn=[0.0 for item in range(max_exponent)]
644 errn=[0.0 for item in range(max_exponent)]
646 #Check if we have a proper numerical value
650 #Loudly and annoyingly complain if it's not a number, then fallback to zero
651 print '''Warning: flatten value is not a number!
652 Use "set flatten" or edit hooke.conf to set it properly
656 for i in range(int(max_cycles)):
658 x_ext=plot.vectors[0][0][contact_index+delta_contact:]
659 y_ext=plot.vectors[0][1][contact_index+delta_contact:]
660 x_ret=plot.vectors[1][0][contact_index+delta_contact:]
661 y_ret=plot.vectors[1][1][contact_index+delta_contact:]
662 for exponent in range(max_exponent):
664 valn[exponent]=sp.polyfit(x_ext,y_ext,exponent)
665 yrn[exponent]=sp.polyval(valn[exponent],x_ret)
666 errn[exponent]=sp.sqrt(sum((yrn[exponent]-y_ext)**2)/float(len(y_ext)))
668 print 'Cannot flatten!'
672 best_exponent=errn.index(min(errn))
675 ycorr_ext=y_ext-yrn[best_exponent]+y_ext[0] #noncontact part
676 yjoin_ext=np.array(plot.vectors[0][1][0:contact_index+delta_contact]) #contact part
678 ycorr_ret=y_ret-yrn[best_exponent]+y_ext[0] #noncontact part
679 yjoin_ret=np.array(plot.vectors[1][1][0:contact_index+delta_contact]) #contact part
681 ycorr_ext=np.concatenate((yjoin_ext, ycorr_ext))
682 ycorr_ret=np.concatenate((yjoin_ret, ycorr_ret))
684 plot.vectors[0][1]=list(ycorr_ext)
685 plot.vectors[1][1]=list(ycorr_ret)
690 def do_slope(self,args):
694 Measures the slope of a delimited chunk on the return trace.
695 The chunk can be delimited either by two manual clicks, or have
696 a fixed width, given as an argument.
698 Syntax: slope [width]
699 The facultative [width] parameter specifies how many
700 points will be considered for the fit. If [width] is
701 specified, only one click will be required.
702 (c) Marco Brucale, Massimo Sandal 2008
705 # Reads the facultative width argument
711 # Decides between the two forms of user input, as per (args)
713 # Gets the Xs of two clicked points as indexes on the current curve vector
714 print 'Click twice to delimit chunk'
715 points=self._measure_N_points(N=2,whatset=1)
717 print 'Click once on the leftmost point of the chunk (i.e.usually the peak)'
718 points=self._measure_N_points(N=1,whatset=1)
720 slope=self._slope(points,fitspan)
722 # Outputs the relevant slope parameter
725 to_dump='slope '+self.current.path+' '+str(slope)
726 self.outlet.push(to_dump)
728 def _slope(self,points,fitspan):
729 # Calls the function linefit_between
730 parameters=[0,0,[],[]]
732 clickedpoints=[points[0].index,points[1].index]
735 clickedpoints=[points[0].index-fitspan,points[0].index]
738 parameters=self.linefit_between(clickedpoints[0],clickedpoints[1])
740 print 'Cannot fit. Did you click twice the same point?'
743 # Outputs the relevant slope parameter
745 print str(parameters[0])
746 to_dump='slope '+self.curve.path+' '+str(parameters[0])
747 self.outlet.push(to_dump)
749 # Makes a vector with the fitted parameters and sends it to the GUI
750 xtoplot=parameters[2]
754 ytoplot.append((x*parameters[0])+parameters[1])
756 clickvector_x, clickvector_y=[], []
758 clickvector_x.append(item.graph_coords[0])
759 clickvector_y.append(item.graph_coords[1])
761 lineplot=self._get_displayed_plot(0) #get topmost displayed plot
763 lineplot.add_set(xtoplot,ytoplot)
764 lineplot.add_set(clickvector_x, clickvector_y)
767 if lineplot.styles==[]:
768 lineplot.styles=[None,None,None,'scatter']
770 lineplot.styles+=[None,'scatter']
771 if lineplot.colors==[]:
772 lineplot.colors=[None,None,'black',None]
774 lineplot.colors+=['black',None]
777 self._send_plot([lineplot])
782 def linefit_between(self,index1,index2,whatset=1):
784 Creates two vectors (xtofit,ytofit) slicing out from the
785 current return trace a portion delimited by the two indexes
787 Then does a least squares linear fit on that slice.
788 Finally returns [0]=the slope, [1]=the intercept of the
789 fitted 1st grade polynomial, and [2,3]=the actual (x,y) vectors
791 (c) Marco Brucale, Massimo Sandal 2008
793 # Translates the indexes into two vectors containing the x,y data to fit
794 xtofit=self.plots[0].vectors[whatset][0][index1:index2]
795 ytofit=self.plots[0].vectors[whatset][1][index1:index2]
797 # Does the actual linear fitting (simple least squares with numpy.polyfit)
799 linefit=np.polyfit(xtofit,ytofit,1)
801 return (linefit[0],linefit[1],xtofit,ytofit)
804 def fit_interval_nm(self,start_index,plot,nm,backwards):
806 Calculates the number of points to fit, given a fit interval in nm
807 start_index: index of point
809 backwards: if true, finds a point backwards.
811 whatset=1 #FIXME: should be decidable
812 x_vect=plot.vectors[1][0]
816 start=x_vect[start_index]
818 while abs(x_vect[i]-x_vect[start_index])*(10**9) < nm:
819 if i==0 or i==maxlen-1: #we reached boundaries of vector!
831 def find_current_peaks(self,noflatten, a=True, maxpeak=True):
834 a=self.convfilt_config['mindeviation']
838 print "Bad input, using default."
839 abs_devs=self.convfilt_config['mindeviation']
841 defplot=self.current.curve.default_plots()[0]
843 flatten=self._find_plotmanip('flatten') #Extract flatten plotmanip
844 defplot=flatten(defplot, self.current, customvalue=1) #Flatten curve before feeding it to has_peaks
845 pk_location,peak_size=self.has_peaks(defplot, abs_devs, maxpeak)
846 return pk_location, peak_size
849 def pickup_contact_point(self,N=1,whatset=1):
850 '''macro to pick up the contact point by clicking'''
851 contact_point=self._measure_N_points(N=1, whatset=1)[0]
852 contact_point_index=contact_point.index
853 self.wlccontact_point=contact_point
854 self.wlccontact_index=contact_point.index
855 self.wlccurrent=self.current.path
856 return contact_point, contact_point_index
859 def baseline_points(self,peak_location, displayed_plot):
860 clicks=self.config['baseline_clicks']
863 base_index_0=peak_location[-1]+self.fit_interval_nm(peak_location[-1], displayed_plot, self.config['auto_right_baseline'],False)
864 self.basepoints.append(self._clickize(displayed_plot.vectors[1][0],displayed_plot.vectors[1][1],base_index_0))
865 base_index_1=self.basepoints[0].index+self.fit_interval_nm(self.basepoints[0].index, displayed_plot, self.config['auto_left_baseline'],False)
866 self.basepoints.append(self._clickize(displayed_plot.vectors[1][0],displayed_plot.vectors[1][1],base_index_1))
868 print 'Select baseline'
870 self.basepoints=self._measure_N_points(N=1, whatset=1)
871 base_index_1=self.basepoints[0].index+self.fit_interval_nm(self.basepoints[0].index, displayed_plot, self.config['auto_left_baseline'], False)
872 self.basepoints.append(self._clickize(displayed_plot.vectors[1][0],displayed_plot.vectors[1][1],base_index_1))
874 self.basepoints=self._measure_N_points(N=2, whatset=1)
876 self.basecurrent=self.current.path
877 return self.basepoints