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]
54 cant_adjust = [c for c in hooke.commands
55 if c.name == 'add block cantilever adjusted extension array'][0]
57 outqueue = NullQueue()
59 for i in range(len(curve.data)):
60 scale(hooke, curve, block=i)
62 params = {'curve':curve, 'block':block}
64 if ('surface distance (m)' not in b.info['columns']
65 or 'surface adjusted deflection (m)' not in b.info['columns']):
67 contact._run(hooke, inqueue, outqueue, params)
69 raise PoorFit('Could not fit %s %s: %s'
70 % (curve.path, block, str(e)))
71 if ('deflection (N)' not in b.info['columns']):
72 force._run(hooke, inqueue, outqueue, params)
73 if ('cantilever adjusted extension (m)' not in b.info['columns']):
74 cant_adjust._run(hooke, inqueue, outqueue, params)
77 class SurfacePositionModel (ModelFitter):
80 The bilinear model is symmetric, but the parameter guessing and
81 sanity checks assume the contact region occurs for lower indicies
82 ("left of") the non-contact region. We also assume that
83 tip-surface attractions produce positive deflections.
87 Algorithm borrowed from WTK's `piezo package`_, specifically
88 from :func:`piezo.z_piezo_utils.analyzeSurfPosData`.
91 http://www.physics.drexel.edu/~wking/code/git/git.php?p=piezo.git
93 Fits the data to the bilinear :method:`model`.
95 In order for this model to produce a satisfactory fit, there
96 should be enough data in the off-surface region that interactions
97 due to proteins, etc. will not seriously skew the fit in the
103 def model(self, params):
104 """A continuous, bilinear model.
111 p_0 + p_1 x & \text{if $x <= p_2$}, \\
112 p_0 + p_1 p_2 + p_3 (x-p_2) & \text{if $x >= p_2$}.
115 Where :math:`p_0` is a vertical offset, :math:`p_1` is the slope
116 of the first region, :math:`p_2` is the transition location, and
117 :math:`p_3` is the slope of the second region.
119 p = params # convenient alias
120 if self.info.get('force zero non-contact slope', None) == True:
122 p.append(0.) # restore the non-contact slope parameter
123 r2 = numpy.round(abs(p[2]))
125 self._model_data[:r2] = p[0] + p[1] * numpy.arange(r2)
126 if r2 < len(self._data)-1:
127 self._model_data[r2:] = \
128 p[0] + p[1]*p[2] + p[3] * numpy.arange(len(self._data)-r2)
129 return self._model_data
131 def set_data(self, data, info=None):
132 super(SurfacePositionModel, self).set_data(data, info)
136 self.info['min position'] = 0
137 self.info['max position'] = len(data)
138 self.info['max deflection'] = data.max()
139 self.info['min deflection'] = data.min()
140 self.info['position range'] = self.info['max position'] - self.info['min position']
141 self.info['deflection range'] = self.info['max deflection'] - self.info['min deflection']
143 def guess_initial_params(self, outqueue=None):
144 """Guess the initial parameters.
148 We guess initial parameters such that the offset (:math:`p_1`)
149 matches the minimum deflection, the kink (:math:`p_2`) occurs in
150 the middle of the data, the initial (contact) slope (:math:`p_0`)
151 produces the maximum deflection at the left-most point, and the
152 final (non-contact) slope (:math:`p_3`) is zero.
154 left_offset = self.info['min deflection']
155 left_slope = 2*(self.info['deflection range']
156 /self.info['position range'])
157 kink_position = (self.info['max position']
158 +self.info['min position'])/2.0
160 self.info['guessed contact slope'] = left_slope
161 params = [left_offset, left_slope, kink_position, right_slope]
162 if self.info.get('force zero non-contact slope', None) == True:
166 def guess_scale(self, params, outqueue=None):
167 """Guess the parameter scales.
172 We guess offset scale (:math:`p_0`) as one tenth of the total
173 deflection range, the kink scale (:math:`p_2`) as one tenth of
174 the total index range, the initial (contact) slope scale
175 (:math:`p_1`) as one tenth of the contact slope estimation,
176 and the final (non-contact) slope scale (:math:`p_3`) is as
177 one tenth of the initial slope scale.
179 offset_scale = self.info['deflection range']/10.
180 left_slope_scale = abs(params[1])/10.
181 kink_scale = self.info['position range']/10.
182 right_slope_scale = left_slope_scale/10.
183 scale = [offset_scale, left_slope_scale, kink_scale, right_slope_scale]
184 if self.info.get('force zero non-contact slope', None) == True:
188 def fit(self, *args, **kwargs):
189 self.info['guessed contact slope'] = None
190 params = super(SurfacePositionModel, self).fit(*args, **kwargs)
191 params[2] = abs(params[2])
192 if self.info.get('force zero non-contact slope', None) == True:
193 params = list(params)
194 params.append(0.) # restore the non-contact slope parameter
196 # check that the fit is reasonable, see the :meth:`model` docstring
197 # for parameter descriptions. HACK: hardcoded cutoffs.
198 if abs(params[3]*10) > abs(params[1]) :
199 raise PoorFit('Slope in non-contact region, or no slope in contact')
200 if params[2] < self.info['min position']+0.02*self.info['position range']:
202 'No kink (kink %g less than %g, need more space to left)'
204 self.info['min position']+0.02*self.info['position range']))
205 if params[2] > self.info['max position']-0.02*self.info['position range']:
207 'No kink (kink %g more than %g, need more space to right)'
209 self.info['max position']-0.02*self.info['position range']))
210 if (self.info['guessed contact slope'] != None
211 and abs(params[1]) < 0.5 * abs(self.info['guessed contact slope'])):
212 raise PoorFit('Too far (contact slope %g, but expected ~%g'
213 % (params[3], self.info['guessed contact slope']))
216 class VelocityClampPlugin (Plugin):
218 super(VelocityClampPlugin, self).__init__(name='vclamp')
220 SurfaceContactCommand(self), ForceCommand(self),
221 CantileverAdjustedExtensionCommand(self),
224 def default_settings(self):
226 Setting(section=self.setting_section, help=self.__doc__),
227 Setting(section=self.setting_section,
228 option='surface contact point algorithm',
230 help='Select the surface contact point algorithm. See the documentation for descriptions of available algorithms.')
234 class SurfaceContactCommand (Command):
235 """Automatically determine a block's surface contact point.
237 Uses the block's `z piezo (m)` and `deflection (m)` arrays.
238 Stores the contact parameters in `block.info`'s `surface distance
239 offset (m)` and `surface deflection offset (m)`. Model-specific
240 fitting information is stored in `surface detection parameters`.
242 The adjusted data columns `surface distance (m)` and `surface
243 adjusted deflection (m)` are also added to the block.
245 You can select the contact point algorithm with the creatively
246 named `surface contact point algorithm` configuration setting.
247 Currently available options are:
249 * fmms (:meth:`find_contact_point_fmms`)
250 * ms (:meth:`find_contact_point_ms`)
251 * wtk (:meth:`find_contact_point_wtk`)
253 def __init__(self, plugin):
254 super(SurfaceContactCommand, self).__init__(
255 name='zero block surface contact point',
258 Argument(name='block', aliases=['set'], type='int', default=0,
260 Data block for which the force should be calculated. For an
261 approach/retract force curve, `0` selects the approaching curve and `1`
262 selects the retracting curve.
265 help=self.__doc__, plugin=plugin)
267 def _run(self, hooke, inqueue, outqueue, params):
268 data = params['curve'].data[params['block']]
269 # HACK? rely on params['curve'] being bound to the local hooke
270 # playlist (i.e. not a copy, as you would get by passing a
271 # curve through the queue). Ugh. Stupid queues. As an
272 # alternative, we could pass lookup information through the
274 new = Data((data.shape[0], data.shape[1]+2), dtype=data.dtype)
275 new.info = copy.deepcopy(data.info)
277 new.info['columns'].extend(
278 ['surface distance (m)', 'surface adjusted deflection (m)'])
279 z_data = data[:,data.info['columns'].index('z piezo (m)')]
280 d_data = data[:,data.info['columns'].index('deflection (m)')]
281 i,deflection_offset,ps = self.find_contact_point(
282 params['curve'], z_data, d_data, outqueue)
283 surface_offset = z_data[i]
284 new.info['surface distance offset (m)'] = surface_offset
285 new.info['surface deflection offset (m)'] = deflection_offset
286 new.info['surface detection parameters'] = ps
287 new[:,-2] = z_data - surface_offset
288 new[:,-1] = d_data - deflection_offset
289 params['curve'].data[params['block']] = new
291 def find_contact_point(self, curve, z_data, d_data, outqueue=None):
292 """Railyard for the `find_contact_point_*` family.
294 Uses the `surface contact point algorithm` configuration
295 setting to call the appropriate backend algorithm.
297 fn = getattr(self, 'find_contact_point_%s'
298 % self.plugin.config['surface contact point algorithm'])
299 return fn(curve, z_data, d_data, outqueue)
301 def find_contact_point_fmms(self, curve, z_data, d_data, outqueue=None):
302 """Algorithm by Francesco Musiani and Massimo Sandal.
308 0) Driver-specific workarounds, e.g. deal with the PicoForce
309 trigger bug by excluding retraction portions with excessive
311 1) Select the second half (non-contact side) of the retraction
313 2) Fit the selection to a line.
314 3) If the fit is not almost horizontal, halve the selection
316 4) Average the selection and use it as a baseline.
317 5) Slide in from the start (contact side) of the retraction
318 curve, until you find a point with greater than baseline
319 deflection. That point is the contact point.
321 if curve.info['filetype'] == 'picoforce':
322 # Take care of the picoforce trigger bug (TODO: example
323 # data file demonstrating the bug). We exclude portions
324 # of the curve that have too much standard deviation.
325 # Yes, a lot of magic is here.
326 check_start = len(d_data)-len(d_data)/20
327 monster_start = len(d_data)
329 # look at the non-contact tail
330 non_monster = d_data[check_start:monster_start]
331 if non_monster.std() < 2e-10: # HACK: hardcoded cutoff
333 else: # move further away from the monster
334 check_start -= len(d_data)/50
335 monster_start -= len(d_data)/50
336 z_data = z_data[:monster_start]
337 d_data = d_data[:monster_start]
339 # take half of the thing to start
340 selection_start = len(d_data)/2
342 z_chunk = z_data[selection_start:]
343 d_chunk = d_data[selection_start:]
344 slope,intercept,r,two_tailed_prob,stderr_of_the_estimate = \
345 scipy.stats.linregress(z_chunk, d_chunk)
346 # We stop if we found an almost-horizontal fit or if we're
347 # getting to small a selection. FIXME: 0.1 and 5./6 here
348 # are "magic numbers" (although reasonable)
349 if (abs(slope) < 0.1 # deflection (m) / surface (m)
350 or selection_start > 5./6*len(d_data)):
352 selection_start += 10
354 d_baseline = d_chunk.mean()
356 # find the first point above the calculated baseline
358 while i < len(d_data) and d_data[i] < ymean:
360 return (i, d_baseline, {})
362 def find_contact_point_ms(self, curve, z_data, d_data, outqueue=None):
363 """Algorithm by Massimo Sandal.
367 WTK: At least the commits are by Massimo, and I see no notes
368 attributing the algorithm to anyone else.
374 xext=raw_plot.vectors[0][0]
375 yext=raw_plot.vectors[0][1]
376 xret2=raw_plot.vectors[1][0]
377 yret=raw_plot.vectors[1][1]
379 first_point=[xext[0], yext[0]]
380 last_point=[xext[-1], yext[-1]]
382 #regr=scipy.polyfit(first_point, last_point,1)[0:2]
383 diffx=abs(first_point[0]-last_point[0])
384 diffy=abs(first_point[1]-last_point[1])
386 #using polyfit results in numerical errors. good old algebra.
388 b=first_point[1]-(a*first_point[0])
389 baseline=scipy.polyval((a,b), xext)
391 ysub=[item-basitem for item,basitem in zip(yext,baseline)]
393 contact=ysub.index(min(ysub))
395 return xext,ysub,contact
397 #now, exploit a ClickedPoint instance to calculate index...
399 dummy.absolute_coords=(x_intercept,y_intercept)
400 dummy.find_graph_coords(xret2,yret)
403 return dummy.index, regr, regr_contact
407 def find_contact_point_wtk(self, curve, z_data, d_data, outqueue=None):
408 """Algorithm by W. Trevor King.
412 Uses :func:`analyze_surf_pos_data_wtk` internally.
414 reverse = z_data[0] > z_data[-1]
415 if reverse == True: # approaching, contact region on the right
416 d_data = d_data[::-1]
417 s = SurfacePositionModel(d_data)
418 s.info['force zero non-contact slope'] = True
419 offset,contact_slope,surface_index,non_contact_slope = s.fit(
423 'contact slope': contact_slope,
424 'surface index': surface_index,
425 'non-contact slope': non_contact_slope,
428 deflection_offset = offset + contact_slope*surface_index,
430 surface_index = len(d_data)-1-surface_index
431 return (numpy.round(surface_index), deflection_offset, info)
433 class ForceCommand (Command):
434 """Calculate a block's `deflection (N)` array.
436 Uses the block's `deflection (m)` array and
437 `spring constant (N/m)`.
439 def __init__(self, plugin):
440 super(ForceCommand, self).__init__(
441 name='add block force array',
444 Argument(name='block', aliases=['set'], type='int', default=0,
446 Data block for which the force should be calculated. For an
447 approach/retract force curve, `0` selects the approaching curve and `1`
448 selects the retracting curve.
451 help=self.__doc__, plugin=plugin)
453 def _run(self, hooke, inqueue, outqueue, params):
454 data = params['curve'].data[params['block']]
455 # HACK? rely on params['curve'] being bound to the local hooke
456 # playlist (i.e. not a copy, as you would get by passing a
457 # curve through the queue). Ugh. Stupid queues. As an
458 # alternative, we could pass lookup information through the
460 new = Data((data.shape[0], data.shape[1]+1), dtype=data.dtype)
461 new.info = copy.deepcopy(data.info)
463 new.info['columns'].append('deflection (N)')
464 d_data = data[:,data.info['columns'].index('surface adjusted deflection (m)')]
465 new[:,-1] = d_data * data.info['spring constant (N/m)']
466 params['curve'].data[params['block']] = new
469 class CantileverAdjustedExtensionCommand (Command):
470 """Calculate a block's `cantilever adjusted extension (m)` array.
472 Uses the block's `deflection (m)` and `surface distance offset (m)`
473 arrays and `spring constant (N/m)`.
475 def __init__(self, plugin):
476 super(CantileverAdjustedExtensionCommand, self).__init__(
477 name='add block cantilever adjusted extension array',
480 Argument(name='block', aliases=['set'], type='int', default=0,
482 Data block for which the adjusted extension should be calculated. For
483 an approach/retract force curve, `0` selects the approaching curve and
484 `1` selects the retracting curve.
487 help=self.__doc__, plugin=plugin)
489 def _run(self, hooke, inqueue, outqueue, params):
490 data = params['curve'].data[params['block']]
491 # HACK? rely on params['curve'] being bound to the local hooke
492 # playlist (i.e. not a copy, as you would get by passing a
493 # curve through the queue). Ugh. Stupid queues. As an
494 # alternative, we could pass lookup information through the
496 new = Data((data.shape[0], data.shape[1]+1), dtype=data.dtype)
497 new.info = copy.deepcopy(data.info)
499 new.info['columns'].append('cantilever adjusted extension (m)')
500 z_data = data[:,data.info['columns'].index('surface distance (m)')]
501 d_data = data[:,data.info['columns'].index('deflection (N)')]
502 new[:,-1] = z_data - d_data / data.info['spring constant (N/m)']
503 params['curve'].data[params['block']] = new
506 class generalvclampCommands(object):
508 def do_subtplot(self, args):
511 (procplots.py plugin)
512 Plots the difference between ret and ext current curve
516 #FIXME: sub_filter and sub_order must be args
518 if len(self.plots[0].vectors) != 2:
519 print 'This command only works on a curve with two different plots.'
522 outplot=self.subtract_curves(sub_order=1)
524 plot_graph=self.list_of_events['plot_graph']
525 wx.PostEvent(self.frame,plot_graph(plots=[outplot]))
527 def _plug_init(self):
528 self.basecurrent=None
532 def do_distance(self,args):
536 Measure the distance (in nm) between two points.
537 For a standard experiment this is the delta X distance.
538 For a force clamp experiment this is the delta Y distance (actually becomes
543 if self.current.curve.experiment == 'clamp':
544 print 'You wanted to use zpiezo perhaps?'
547 dx,unitx,dy,unity=self._delta(set=1)
548 print str(dx*(10**9))+' nm'
549 to_dump='distance '+self.current.path+' '+str(dx*(10**9))+' nm'
550 self.outlet.push(to_dump)
553 def do_force(self,args):
557 Measure the force difference (in pN) between two points
561 if self.current.curve.experiment == 'clamp':
562 print 'This command makes no sense for a force clamp experiment.'
564 dx,unitx,dy,unity=self._delta(set=1)
565 print str(dy*(10**12))+' pN'
566 to_dump='force '+self.current.path+' '+str(dy*(10**12))+' pN'
567 self.outlet.push(to_dump)
570 def do_forcebase(self,args):
574 Measures the difference in force (in pN) between a point and a baseline
575 took as the average between two points.
577 The baseline is fixed once for a given curve and different force measurements,
578 unless the user wants it to be recalculated
580 Syntax: forcebase [rebase]
581 rebase: Forces forcebase to ask again the baseline
582 max: Instead of asking for a point to measure, asks for two points and use
583 the maximum peak in between
585 rebase=False #if true=we select rebase
586 maxpoint=False #if true=we measure the maximum peak
588 plot=self._get_displayed_plot()
589 whatset=1 #fixme: for all sets
590 if 'rebase' in args or (self.basecurrent != self.current.path):
596 print 'Select baseline'
597 self.basepoints=self._measure_N_points(N=2, whatset=whatset)
598 self.basecurrent=self.current.path
601 print 'Select two points'
602 points=self._measure_N_points(N=2, whatset=whatset)
603 boundpoints=[points[0].index, points[1].index]
606 y=min(plot.vectors[whatset][1][boundpoints[0]:boundpoints[1]])
608 print 'Chosen interval not valid. Try picking it again. Did you pick the same point as begin and end of interval?'
610 print 'Select point to measure'
611 points=self._measure_N_points(N=1, whatset=whatset)
612 #whatplot=points[0].dest
613 y=points[0].graph_coords[1]
615 #fixme: code duplication
616 boundaries=[self.basepoints[0].index, self.basepoints[1].index]
618 to_average=plot.vectors[whatset][1][boundaries[0]:boundaries[1]] #y points to average
620 avg=np.mean(to_average)
622 print str(forcebase*(10**12))+' pN'
623 to_dump='forcebase '+self.current.path+' '+str(forcebase*(10**12))+' pN'
624 self.outlet.push(to_dump)
626 def plotmanip_multiplier(self, plot, current):
628 Multiplies all the Y values of an SMFS curve by a value stored in the 'force_multiplier'
629 configuration variable. Useful for calibrations and other stuff.
633 if current.curve.experiment != 'smfs':
636 #only one set is present...
637 if len(self.plots[0].vectors) != 2:
641 if (self.config['force_multiplier']==1):
644 for i in range(len(plot.vectors[0][1])):
645 plot.vectors[0][1][i]=plot.vectors[0][1][i]*self.config['force_multiplier']
647 for i in range(len(plot.vectors[1][1])):
648 plot.vectors[1][1][i]=plot.vectors[1][1][i]*self.config['force_multiplier']
653 def plotmanip_flatten(self, plot, current, customvalue=False):
655 Subtracts a polynomial fit to the non-contact part of the curve, as to flatten it.
656 the best polynomial fit is chosen among polynomials of degree 1 to n, where n is
657 given by the configuration file or by the customvalue.
659 customvalue= int (>0) --> starts the function even if config says no (default=False)
663 if current.curve.experiment != 'smfs':
666 #only one set is present...
667 if len(self.plots[0].vectors) != 2:
670 #config is not flatten, and customvalue flag is false too
671 if (not self.config['flatten']) and (not customvalue):
678 max_cycles=customvalue
680 max_cycles=self.config['flatten'] #Using > 1 usually doesn't help and can give artefacts. However, it could be useful too.
682 contact_index=self.find_contact_point()
684 valn=[[] for item in range(max_exponent)]
685 yrn=[0.0 for item in range(max_exponent)]
686 errn=[0.0 for item in range(max_exponent)]
688 #Check if we have a proper numerical value
692 #Loudly and annoyingly complain if it's not a number, then fallback to zero
693 print '''Warning: flatten value is not a number!
694 Use "set flatten" or edit hooke.conf to set it properly
698 for i in range(int(max_cycles)):
700 x_ext=plot.vectors[0][0][contact_index+delta_contact:]
701 y_ext=plot.vectors[0][1][contact_index+delta_contact:]
702 x_ret=plot.vectors[1][0][contact_index+delta_contact:]
703 y_ret=plot.vectors[1][1][contact_index+delta_contact:]
704 for exponent in range(max_exponent):
706 valn[exponent]=sp.polyfit(x_ext,y_ext,exponent)
707 yrn[exponent]=sp.polyval(valn[exponent],x_ret)
708 errn[exponent]=sp.sqrt(sum((yrn[exponent]-y_ext)**2)/float(len(y_ext)))
710 print 'Cannot flatten!'
714 best_exponent=errn.index(min(errn))
717 ycorr_ext=y_ext-yrn[best_exponent]+y_ext[0] #noncontact part
718 yjoin_ext=np.array(plot.vectors[0][1][0:contact_index+delta_contact]) #contact part
720 ycorr_ret=y_ret-yrn[best_exponent]+y_ext[0] #noncontact part
721 yjoin_ret=np.array(plot.vectors[1][1][0:contact_index+delta_contact]) #contact part
723 ycorr_ext=np.concatenate((yjoin_ext, ycorr_ext))
724 ycorr_ret=np.concatenate((yjoin_ret, ycorr_ret))
726 plot.vectors[0][1]=list(ycorr_ext)
727 plot.vectors[1][1]=list(ycorr_ret)
732 def do_slope(self,args):
736 Measures the slope of a delimited chunk on the return trace.
737 The chunk can be delimited either by two manual clicks, or have
738 a fixed width, given as an argument.
740 Syntax: slope [width]
741 The facultative [width] parameter specifies how many
742 points will be considered for the fit. If [width] is
743 specified, only one click will be required.
744 (c) Marco Brucale, Massimo Sandal 2008
747 # Reads the facultative width argument
753 # Decides between the two forms of user input, as per (args)
755 # Gets the Xs of two clicked points as indexes on the current curve vector
756 print 'Click twice to delimit chunk'
757 points=self._measure_N_points(N=2,whatset=1)
759 print 'Click once on the leftmost point of the chunk (i.e.usually the peak)'
760 points=self._measure_N_points(N=1,whatset=1)
762 slope=self._slope(points,fitspan)
764 # Outputs the relevant slope parameter
767 to_dump='slope '+self.current.path+' '+str(slope)
768 self.outlet.push(to_dump)
770 def _slope(self,points,fitspan):
771 # Calls the function linefit_between
772 parameters=[0,0,[],[]]
774 clickedpoints=[points[0].index,points[1].index]
777 clickedpoints=[points[0].index-fitspan,points[0].index]
780 parameters=self.linefit_between(clickedpoints[0],clickedpoints[1])
782 print 'Cannot fit. Did you click twice the same point?'
785 # Outputs the relevant slope parameter
787 print str(parameters[0])
788 to_dump='slope '+self.curve.path+' '+str(parameters[0])
789 self.outlet.push(to_dump)
791 # Makes a vector with the fitted parameters and sends it to the GUI
792 xtoplot=parameters[2]
796 ytoplot.append((x*parameters[0])+parameters[1])
798 clickvector_x, clickvector_y=[], []
800 clickvector_x.append(item.graph_coords[0])
801 clickvector_y.append(item.graph_coords[1])
803 lineplot=self._get_displayed_plot(0) #get topmost displayed plot
805 lineplot.add_set(xtoplot,ytoplot)
806 lineplot.add_set(clickvector_x, clickvector_y)
809 if lineplot.styles==[]:
810 lineplot.styles=[None,None,None,'scatter']
812 lineplot.styles+=[None,'scatter']
813 if lineplot.colors==[]:
814 lineplot.colors=[None,None,'black',None]
816 lineplot.colors+=['black',None]
819 self._send_plot([lineplot])
824 def linefit_between(self,index1,index2,whatset=1):
826 Creates two vectors (xtofit,ytofit) slicing out from the
827 current return trace a portion delimited by the two indexes
829 Then does a least squares linear fit on that slice.
830 Finally returns [0]=the slope, [1]=the intercept of the
831 fitted 1st grade polynomial, and [2,3]=the actual (x,y) vectors
833 (c) Marco Brucale, Massimo Sandal 2008
835 # Translates the indexes into two vectors containing the x,y data to fit
836 xtofit=self.plots[0].vectors[whatset][0][index1:index2]
837 ytofit=self.plots[0].vectors[whatset][1][index1:index2]
839 # Does the actual linear fitting (simple least squares with numpy.polyfit)
841 linefit=np.polyfit(xtofit,ytofit,1)
843 return (linefit[0],linefit[1],xtofit,ytofit)
846 def fit_interval_nm(self,start_index,plot,nm,backwards):
848 Calculates the number of points to fit, given a fit interval in nm
849 start_index: index of point
851 backwards: if true, finds a point backwards.
853 whatset=1 #FIXME: should be decidable
854 x_vect=plot.vectors[1][0]
858 start=x_vect[start_index]
860 while abs(x_vect[i]-x_vect[start_index])*(10**9) < nm:
861 if i==0 or i==maxlen-1: #we reached boundaries of vector!
873 def find_current_peaks(self,noflatten, a=True, maxpeak=True):
876 a=self.convfilt_config['mindeviation']
880 print "Bad input, using default."
881 abs_devs=self.convfilt_config['mindeviation']
883 defplot=self.current.curve.default_plots()[0]
885 flatten=self._find_plotmanip('flatten') #Extract flatten plotmanip
886 defplot=flatten(defplot, self.current, customvalue=1) #Flatten curve before feeding it to has_peaks
887 pk_location,peak_size=self.has_peaks(defplot, abs_devs, maxpeak)
888 return pk_location, peak_size
891 def pickup_contact_point(self,N=1,whatset=1):
892 '''macro to pick up the contact point by clicking'''
893 contact_point=self._measure_N_points(N=1, whatset=1)[0]
894 contact_point_index=contact_point.index
895 self.wlccontact_point=contact_point
896 self.wlccontact_index=contact_point.index
897 self.wlccurrent=self.current.path
898 return contact_point, contact_point_index
901 def baseline_points(self,peak_location, displayed_plot):
902 clicks=self.config['baseline_clicks']
905 base_index_0=peak_location[-1]+self.fit_interval_nm(peak_location[-1], displayed_plot, self.config['auto_right_baseline'],False)
906 self.basepoints.append(self._clickize(displayed_plot.vectors[1][0],displayed_plot.vectors[1][1],base_index_0))
907 base_index_1=self.basepoints[0].index+self.fit_interval_nm(self.basepoints[0].index, displayed_plot, self.config['auto_left_baseline'],False)
908 self.basepoints.append(self._clickize(displayed_plot.vectors[1][0],displayed_plot.vectors[1][1],base_index_1))
910 print 'Select baseline'
912 self.basepoints=self._measure_N_points(N=1, whatset=1)
913 base_index_1=self.basepoints[0].index+self.fit_interval_nm(self.basepoints[0].index, displayed_plot, self.config['auto_left_baseline'], False)
914 self.basepoints.append(self._clickize(displayed_plot.vectors[1][0],displayed_plot.vectors[1][1],base_index_1))
916 self.basepoints=self._measure_N_points(N=2, whatset=1)
918 self.basecurrent=self.current.path
919 return self.basepoints