1 # Copyright (C) 2008-2011 Alberto Gomez-Casado
3 # Massimo Sandal <devicerandom@gmail.com>
4 # W. Trevor King <wking@drexel.edu>
6 # This file is part of Hooke.
8 # Hooke is free software: you can redistribute it and/or modify it
9 # under the terms of the GNU Lesser General Public License as
10 # published by the Free Software Foundation, either version 3 of the
11 # License, or (at your option) any later version.
13 # Hooke is distributed in the hope that it will be useful, but WITHOUT
14 # ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
15 # or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General
16 # Public License for more details.
18 # You should have received a copy of the GNU Lesser General Public
19 # License along with Hooke. If not, see
20 # <http://www.gnu.org/licenses/>.
22 """The ``curve`` module provides :class:`CurvePlugin` and several
23 associated :class:`hooke.command.Command`\s for handling
24 :mod:`hooke.curve` classes.
34 from ..command import Command, Argument, Failure
35 from ..command_stack import CommandStack
36 from ..curve import Data
37 from ..engine import CommandMessage
38 from ..util.calculus import derivative
39 from ..util.fft import unitary_avg_power_spectrum
40 from ..util.si import ppSI, join_data_label, split_data_label
42 from .playlist import current_playlist_callback
45 # Define common or complicated arguments
47 def current_curve_callback(hooke, command, argument, value, load=True):
50 playlist = current_playlist_callback(hooke, command, argument, value)
51 curve = playlist.current(load=load)
53 raise Failure('No curves in %s' % playlist)
56 def unloaded_current_curve_callback(hooke, command, argument, value):
57 return current_curve_callback(
58 hooke=hooke, command=command, argument=argument, value=value,
61 CurveArgument = Argument(
62 name='curve', type='curve', callback=current_curve_callback,
64 :class:`hooke.curve.Curve` to act on. Defaults to the current curve
65 of the current playlist.
68 def _name_argument(name, default, help):
71 return Argument(name=name, type='string', default=default, help=help)
73 def block_argument(*args, **kwargs):
76 return _name_argument(*args, **kwargs)
78 def column_argument(*args, **kwargs):
81 return _name_argument(*args, **kwargs)
84 # Define useful command subclasses
86 class CurveCommand (Command):
87 """A :class:`~hooke.command.Command` operating on a
88 :class:`~hooke.curve.Curve`.
90 def __init__(self, **kwargs):
91 if 'arguments' in kwargs:
92 kwargs['arguments'].insert(0, CurveArgument)
94 kwargs['arguments'] = [CurveArgument]
95 super(CurveCommand, self).__init__(**kwargs)
97 def _curve(self, hooke, params):
98 """Get the selected curve.
102 `hooke` is intended to attach the selected curve to the local
103 playlist; the returned curve should not be effected by the
104 state of `hooke`. This is important for reliable
105 :class:`~hooke.command_stack.CommandStack`\s.
107 # HACK? rely on params['curve'] being bound to the local hooke
108 # playlist (i.e. not a copy, as you would get by passing a
109 # curve through the queue). Ugh. Stupid queues. As an
110 # alternative, we could pass lookup information through the
112 return params['curve']
114 def _add_to_command_stack(self, params):
115 """Store the command name and current `params` values in the
116 curve's `.command_stack`.
118 If this would duplicate the command currently on top of the
119 stack, no action is taken. Call early on, or watch out for
120 repeated param processing.
122 Recommended practice is to *not* lock in argument values that
123 are loaded from the plugin's :attr:`.config`.
127 Perhaps we should subclass :meth:`_run` and use :func:`super`,
128 or embed this in :meth:`run` to avoid subclasses calling this
129 method explicitly, with all the tedium and brittality that
130 implies. On the other hand, the current implemtnation allows
131 CurveCommands that don't effect the curve itself
132 (e.g. :class:`GetCommand`) to avoid adding themselves to the
135 if params['stack'] == True:
136 curve = self._curve(hooke=None, params=params)
137 if (len(curve.command_stack) > 0
138 and curve.command_stack[-1].command == self.name
139 and curve.command_stack[-1].arguments == params):
140 pass # no need to place duplicate calls on the stack.
142 curve.command_stack.append(CommandMessage(
143 self.name, dict(params)))
146 class BlockCommand (CurveCommand):
147 """A :class:`CurveCommand` operating on a :class:`~hooke.curve.Data` block.
149 def __init__(self, blocks=None, **kwargs):
151 blocks = [('block', None, 'Name of the data block to act on.')]
153 for name,default,help in blocks:
154 block_args.append(block_argument(name, default, help))
155 self._block_arguments = block_args
156 if 'arguments' not in kwargs:
157 kwargs['arguments'] = []
158 kwargs['arguments'] = block_args + kwargs['arguments']
159 super(BlockCommand, self).__init__(**kwargs)
161 def _block_names(self, hooke, params):
162 curve = self._curve(hooke, params)
163 return [b.info['name'] for b in curve.data]
165 def _block_index(self, hooke, params, name=None):
167 name = self._block_arguments[0].name
168 block_name = params[name]
169 if block_name == None:
170 curve = self._curve(hooke=hooke, params=params)
171 if len(curve.data) == 0:
172 raise Failure('no blocks in %s' % curve)
173 block_name = curve.data[0].info['name']
174 names = self._block_names(hooke=hooke, params=params)
176 return names.index(block_name)
177 except ValueError, e:
178 curve = self._curve(hooke, params)
179 raise Failure('no block named %s in %s (%s): %s'
180 % (block_name, curve, names, e))
182 def _block(self, hooke, params, name=None):
183 # HACK? rely on params['block'] being bound to the local hooke
184 # playlist (i.e. not a copy, as you would get by passing a
185 # block through the queue). Ugh. Stupid queues. As an
186 # alternative, we could pass lookup information through the
188 curve = self._curve(hooke, params)
189 index = self._block_index(hooke, params, name)
190 return curve.data[index]
193 class ColumnAccessCommand (BlockCommand):
194 """A :class:`BlockCommand` accessing a :class:`~hooke.curve.Data`
197 def __init__(self, columns=None, **kwargs):
199 columns = [('column', None, 'Name of the data column to act on.')]
201 for name,default,help in columns:
202 column_args.append(column_argument(name, default, help))
203 self._column_arguments = column_args
204 if 'arguments' not in kwargs:
205 kwargs['arguments'] = []
206 kwargs['arguments'] = column_args + kwargs['arguments']
207 super(ColumnAccessCommand, self).__init__(**kwargs)
209 def _get_column(self, hooke, params, block_name=None, column_name=None):
210 if column_name == None:
211 column_name = self._column_arguments[0].name
212 column_name = params[column_name]
213 if column_name is None:
215 block = self._block(hooke, params, block_name)
216 columns = block.info['columns']
218 column_index = columns.index(column_name)
219 except ValueError, e:
220 raise Failure('%s not in %s (%s): %s'
221 % (column_name, block.info['name'], columns, e))
222 return block[:,column_index]
225 class ColumnAddingCommand (ColumnAccessCommand):
226 """A :class:`ColumnAccessCommand` that also adds columns.
228 def __init__(self, new_columns=None, **kwargs):
229 if new_columns == None:
232 for name,default,help in new_columns:
233 column_args.append(column_argument(name, default, help))
234 self._new_column_arguments = column_args
235 if 'arguments' not in kwargs:
236 kwargs['arguments'] = []
237 kwargs['arguments'] = column_args + kwargs['arguments']
238 super(ColumnAddingCommand, self).__init__(**kwargs)
240 def _get_column(self, hooke, params, block_name=None, column_name=None):
241 if column_name == None and len(self._column_arguments) == 0:
242 column_name = self._new_column_arguments[0].name
243 return super(ColumnAddingCommand, self)._get_column(
244 hooke=hooke, params=params, block_name=block_name,
245 column_name=column_name)
247 def _set_column(self, hooke, params, block_name=None, column_name=None,
249 if column_name == None:
250 column_name = self._column_arguments[0].name
251 column_name = params[column_name]
252 block = self._block(hooke=hooke, params=params, name=block_name)
253 if column_name not in block.info['columns']:
254 new = Data((block.shape[0], block.shape[1]+1), dtype=block.dtype)
255 new.info = copy.deepcopy(block.info)
257 new.info['columns'].append(column_name)
259 block_index = self._block_index(hooke, params, name=block_name)
260 self._curve(hooke, params).data[block_index] = block
261 column_index = block.info['columns'].index(column_name)
262 block[:,column_index] = values
267 class CurvePlugin (Builtin):
269 super(CurvePlugin, self).__init__(name='curve')
271 GetCommand(self), InfoCommand(self), BlockInfoCommand(self),
272 DeltaCommand(self), ExportCommand(self), DifferenceCommand(self),
273 DerivativeCommand(self), PowerSpectrumCommand(self),
274 ScaledColumnAdditionCommand(self), ClearStackCommand(self)]
279 class GetCommand (CurveCommand):
280 """Return a :class:`hooke.curve.Curve`.
282 def __init__(self, plugin):
283 super(GetCommand, self).__init__(
284 name='get curve', help=self.__doc__, plugin=plugin)
286 def _run(self, hooke, inqueue, outqueue, params):
287 outqueue.put(self._curve(hooke, params))
290 class InfoCommand (CurveCommand):
291 """Get selected information about a :class:`hooke.curve.Curve`.
293 def __init__(self, plugin):
295 Argument(name='all', type='bool', default=False, count=1,
296 help='Get all curve information.'),
298 self.fields = ['name', 'path', 'driver', 'note', 'command stack',
299 'blocks', 'block names', 'block sizes']
300 for field in self.fields:
301 args.append(Argument(
302 name=field, type='bool', default=False, count=1,
303 help='Get curve %s' % field))
304 super(InfoCommand, self).__init__(
305 name='curve info', arguments=args,
306 help=self.__doc__, plugin=plugin)
308 def _run(self, hooke, inqueue, outqueue, params):
309 curve = self._curve(hooke, params)
311 for key in self.fields:
312 fields[key] = params[key]
313 if reduce(lambda x,y: x or y, fields.values()) == False:
314 params['all'] = True # No specific fields set, default to 'all'
315 if params['all'] == True:
316 for key in self.fields:
319 for key in self.fields:
320 if fields[key] == True:
321 get = getattr(self, '_get_%s' % key.replace(' ', '_'))
322 lines.append('%s: %s' % (key, get(curve)))
323 outqueue.put('\n'.join(lines))
325 def _get_name(self, curve):
328 def _get_path(self, curve):
331 def _get_driver(self, curve):
334 def _get_note(self, curve):
335 return curve.info.get('note', None)
337 def _get_command_stack(self, curve):
338 return curve.command_stack
340 def _get_blocks(self, curve):
341 return len(curve.data)
343 def _get_block_names(self, curve):
344 return [block.info['name'] for block in curve.data]
346 def _get_block_sizes(self, curve):
347 return [block.shape for block in curve.data]
350 class BlockInfoCommand (BlockCommand):
351 """Get selected information about a :class:`hooke.curve.Curve` data block.
353 def __init__(self, plugin):
354 super(BlockInfoCommand, self).__init__(
355 name='block info', arguments=[
357 name='key', count=-1, optional=False,
358 help='Dot-separted (.) key selection regexp.'),
362 File name for the output (appended).
365 help=self.__doc__, plugin=plugin)
367 def _run(self, hooke, inqueue, outqueue, params):
368 block = self._block(hooke, params)
369 values = {'index': self._block_index(hooke, params)}
370 for key in params['key']:
371 keys = [(0, key.split('.'), block.info)]
373 index,key_stack,info = keys.pop(0)
374 regexp = re.compile(key_stack[index])
376 for k,v in info.items():
379 new_stack = copy.copy(key_stack)
381 if index+1 == len(key_stack):
383 for k in new_stack[:-1]:
387 vals[new_stack[-1]] = v
389 keys.append((index+1, new_stack, v))
392 'no match found for %s (%s) in %s'
393 % (key_stack[index], key, sorted(info.keys())))
394 if params['output'] != None:
395 curve = self._curve(hooke, params)
396 with open(params['output'], 'a') as f:
397 yaml.dump({curve.name:{
399 block.info['name']: values
404 class DeltaCommand (BlockCommand):
405 """Get distance information between two points.
407 With two points A and B, the returned distances are A-B.
409 def __init__(self, plugin):
410 super(DeltaCommand, self).__init__(
413 Argument(name='point', type='point', optional=False, count=2,
415 Indicies of points bounding the selected data.
417 Argument(name='SI', type='bool', default=False,
419 Return distances in SI notation.
422 help=self.__doc__, plugin=plugin)
424 def _run(self, hooke, inqueue, outqueue, params):
425 data = self._block(hooke, params)
426 As = data[params['point'][0],:]
427 Bs = data[params['point'][1],:]
428 ds = [A-B for A,B in zip(As, Bs)]
429 if params['SI'] == False:
430 out = [(name, d) for name,d in zip(data.info['columns'], ds)]
433 for name,d in zip(data.info['columns'], ds):
434 n,units = split_data_label(name)
436 (n, ppSI(value=d, unit=units, decimals=2)))
440 class ExportCommand (BlockCommand):
441 """Export a :class:`hooke.curve.Curve` data block as TAB-delimeted
444 A "#" prefixed header will optionally appear at the beginning of
445 the file naming the columns.
447 def __init__(self, plugin):
448 super(ExportCommand, self).__init__(
451 Argument(name='output', type='file', default='curve.dat',
453 File name for the output data. Defaults to 'curve.dat'
455 Argument(name='header', type='bool', default=True,
457 True if you want the column-naming header line.
460 help=self.__doc__, plugin=plugin)
462 def _run(self, hooke, inqueue, outqueue, params):
463 data = self._block(hooke, params)
465 with open(os.path.expanduser(params['output']), 'w') as f:
466 if params['header'] == True:
467 f.write('# %s \n' % ('\t'.join(data.info['columns'])))
468 numpy.savetxt(f, data, delimiter='\t')
471 class DifferenceCommand (ColumnAddingCommand):
472 """Calculate the difference between two columns of data.
474 The difference is added to block A as a new column.
476 Note that the command will fail if the columns have different
477 lengths, so be careful when differencing columns from different
480 def __init__(self, plugin):
481 super(DifferenceCommand, self).__init__(
485 'Name of block A in A-B. Defaults to the first block'),
487 'Name of block B in A-B. Defaults to matching `block A`.'),
492 Column of data from block A to difference. Defaults to the first column.
496 Column of data from block B to difference. Defaults to matching `column A`.
500 ('output column', None,
502 Name of the new column for storing the difference (without units, defaults to
503 `difference of <block A> <column A> and <block B> <column B>`).
506 help=self.__doc__, plugin=plugin)
508 def _run(self, hooke, inqueue, outqueue, params):
509 self._add_to_command_stack(params)
510 params = self._setup_params(hooke=hooke, params=params)
511 data_A = self._get_column(hooke=hooke, params=params,
512 block_name='block A',
513 column_name='column A')
514 data_B = self._get_column(hooke=hooke, params=params,
515 block_name='block B',
516 column_name='column B')
517 out = data_A - data_B
518 self._set_column(hooke=hooke, params=params,
519 block_name='block A',
520 column_name='output column',
523 def _setup_params(self, hooke, params):
524 curve = self._curve(hooke, params)
525 if params['block A'] == None:
526 params['block A'] = curve.data[0].info['name']
527 if params['block B'] == None:
528 params['block B'] = params['block A']
529 block_A = self._block(hooke, params=params, name='block A')
530 block_B = self._block(hooke, params=params, name='block B')
531 if params['column A'] == None:
532 params['column A'] = block.info['columns'][0]
533 if params['column B'] == None:
534 params['column B'] = params['column A']
535 a_name,a_unit = split_data_label(params['column A'])
536 b_name,b_unit = split_data_label(params['column B'])
538 raise Failure('Unit missmatch: %s != %s' % (a_unit, b_unit))
539 if params['output column'] == None:
540 params['output column'] = join_data_label(
541 'difference of %s %s and %s %s' % (
542 block_A.info['name'], a_name,
543 block_B.info['name'], b_name),
546 params['output column'] = join_data_label(
547 params['output column'], a_unit)
551 class DerivativeCommand (ColumnAddingCommand):
552 """Calculate the derivative (actually, the discrete differentiation)
555 See :func:`hooke.util.calculus.derivative` for implementation
558 def __init__(self, plugin):
559 super(DerivativeCommand, self).__init__(
563 'Column of data block to differentiate with respect to.'),
565 'Column of data block to differentiate.'),
568 ('output column', None,
570 Name of the new column for storing the derivative (without units, defaults to
571 `derivative of <f column> with respect to <x column>`).
575 Argument(name='weights', type='dict', default={-1:-0.5, 1:0.5},
577 Weighting scheme dictionary for finite differencing. Defaults to
578 central differencing.
581 help=self.__doc__, plugin=plugin)
583 def _run(self, hooke, inqueue, outqueue, params):
584 self._add_to_command_stack(params)
585 params = self._setup_params(hooke=hooke, params=params)
586 x_data = self._get_column(hooke=hooke, params=params,
587 column_name='x column')
588 f_data = self._get_column(hooke=hooke, params=params,
589 column_name='f column')
591 x_data=x_data, f_data=f_data, weights=params['weights'])
592 self._set_column(hooke=hooke, params=params,
593 column_name='output column',
596 def _setup_params(self, hooke, params):
597 curve = self._curve(hooke, params)
598 x_name,x_unit = split_data_label(params['x column'])
599 f_name,f_unit = split_data_label(params['f column'])
600 d_unit = '%s/%s' % (f_unit, x_unit)
601 if params['output column'] == None:
602 params['output column'] = join_data_label(
603 'derivative of %s with respect to %s' % (
607 params['output column'] = join_data_label(
608 params['output column'], d_unit)
612 class PowerSpectrumCommand (ColumnAddingCommand):
613 """Calculate the power spectrum of a data column.
615 def __init__(self, plugin):
616 super(PowerSpectrumCommand, self).__init__(
617 name='power spectrum',
619 Argument(name='output block', type='string',
621 Name of the new data block for storing the power spectrum (defaults to
622 `power spectrum of <source block name> <source column name>`).
624 Argument(name='bounds', type='point', optional=True, count=2,
626 Indicies of points bounding the selected data.
628 Argument(name='freq', type='float', default=1.0,
632 Argument(name='freq units', type='string', default='Hz',
634 Units for the sampling frequency.
636 Argument(name='chunk size', type='int', default=2048,
638 Number of samples per chunk. Use a power of two.
640 Argument(name='overlap', type='bool', default=False,
642 If `True`, each chunk overlaps the previous chunk by half its length.
643 Otherwise, the chunks are end-to-end, and not overlapping.
646 help=self.__doc__, plugin=plugin)
648 def _run(self, hooke, inqueue, outqueue, params):
649 self._add_to_command_stack(params)
650 params = self._setup_params(hooke=hooke, params=params)
651 data = self._get_column(hooke=hooke, params=params)
652 bounds = params['bounds']
654 data = data[bounds[0]:bounds[1]]
655 freq_axis,power = unitary_avg_power_spectrum(
656 data, freq=params['freq'],
657 chunk_size=params['chunk size'],
658 overlap=params['overlap'])
659 b = Data((len(freq_axis),2), dtype=data.dtype)
660 b.info['name'] = params['output block']
661 b.info['columns'] = [
662 params['output freq column'],
663 params['output power column'],
665 self._curve(hooke, params).data.append(b)
666 self._set_column(hooke, params, block_name='output block',
667 column_name='output freq column',
669 self._set_column(hooke, params, block_name='output block',
670 column_name='output power column',
674 def _setup_params(self, hooke, params):
675 if params['output block'] in self._block_names(hooke, params):
676 raise Failure('output block %s already exists in %s.'
677 % (params['output block'],
678 self._curve(hooke, params)))
679 data = self._get_column(hooke=hooke, params=params)
680 d_name,d_unit = split_data_label(data.info['name'])
681 if params['output block'] == None:
682 params['output block'] = 'power spectrum of %s %s' % (
683 data.info['name'], params['column'])
684 self.params['output freq column'] = join_data_label(
685 'frequency axis', params['freq units'])
686 self.params['output power column'] = join_data_label(
687 'power density', '%s^2/%s' % (data_units, params['freq units']))
691 class ScaledColumnAdditionCommand (ColumnAddingCommand):
692 """Add one affine transformed column to another: `o=A*i1+B*i2+C`.
694 def __init__(self, plugin):
695 super(ScaledColumnAdditionCommand, self).__init__(
696 name='scaled column addition',
698 ('input column 1', 'input column (m)', """
699 Name of the first column to use as the transform input.
701 ('input column 2', None, """
702 Name of the second column to use as the transform input.
706 ('output column', 'output column (m)', """
707 Name of the column to use as the transform output.
711 Argument(name='scale 1', type='float', default=None,
713 A float value for the first scale constant.
715 Argument(name='scale 1 name', type='string', default=None,
717 The name of the first scale constant in the `.info` dictionary.
719 Argument(name='scale 2', type='float', default=None,
721 A float value for the second scale constant.
723 Argument(name='scale 2 name', type='string', default=None,
725 The name of the second scale constant in the `.info` dictionary.
727 Argument(name='constant', type='float', default=None,
729 A float value for the offset constant.
731 Argument(name='constant name', type='string', default=None,
733 The name of the offset constant in the `.info` dictionary.
736 help=self.__doc__, plugin=plugin)
738 def _run(self, hooke, inqueue, outqueue, params):
739 self._add_to_command_stack(params)
740 i1 = self._get_column(hooke=hooke, params=params,
741 column_name='input column 1')
742 i2 = self._get_column(hooke=hooke, params=params,
743 column_name='input column 2')
748 # what if both i1 and i2 are None?
749 a = self._get_constant(params, i1.info, 'scale 1')
750 b = self._get_constant(params, i2.info, 'scale 2')
751 c = self._get_constant(params, i1.info, 'constant')
752 out = a*i1 + b*i2 + c
753 self._set_column(hooke=hooke, params=params,
754 column_name='output column', values=out)
756 def _get_constant(self, params, info, name):
758 pname = params[name + ' name']
760 if pname is not None:
761 pname_entries = pname.split('|')
763 for entry in pname_entries:
765 if a is None and b is None:
774 class ClearStackCommand (CurveCommand):
775 """Empty a curve's command stack.
777 def __init__(self, plugin):
778 super(ClearStackCommand, self).__init__(
779 name='clear curve command stack',
780 help=self.__doc__, plugin=plugin)
781 i,arg = [(i,arg) for i,arg in enumerate(self.arguments)
782 if arg.name == 'curve'][0]
784 arg.callback = unloaded_current_curve_callback
785 self.arguments[i] = arg
787 def _run(self, hooke, inqueue, outqueue, params):
788 curve = self._curve(hooke, params)
789 curve.command_stack = CommandStack()
792 class OldCruft (object):
794 def do_forcebase(self,args):
798 Measures the difference in force (in pN) between a point and a baseline
799 took as the average between two points.
801 The baseline is fixed once for a given curve and different force measurements,
802 unless the user wants it to be recalculated
804 Syntax: forcebase [rebase]
805 rebase: Forces forcebase to ask again the baseline
806 max: Instead of asking for a point to measure, asks for two points and use
807 the maximum peak in between
809 rebase=False #if true=we select rebase
810 maxpoint=False #if true=we measure the maximum peak
812 plot=self._get_displayed_plot()
813 whatset=1 #fixme: for all sets
814 if 'rebase' in args or (self.basecurrent != self.current.path):
820 print 'Select baseline'
821 self.basepoints=self._measure_N_points(N=2, whatset=whatset)
822 self.basecurrent=self.current.path
825 print 'Select two points'
826 points=self._measure_N_points(N=2, whatset=whatset)
827 boundpoints=[points[0].index, points[1].index]
830 y=min(plot.vectors[whatset][1][boundpoints[0]:boundpoints[1]])
832 print 'Chosen interval not valid. Try picking it again. Did you pick the same point as begin and end of interval?'
834 print 'Select point to measure'
835 points=self._measure_N_points(N=1, whatset=whatset)
836 #whatplot=points[0].dest
837 y=points[0].graph_coords[1]
839 #fixme: code duplication
840 boundaries=[self.basepoints[0].index, self.basepoints[1].index]
842 to_average=plot.vectors[whatset][1][boundaries[0]:boundaries[1]] #y points to average
844 avg=np.mean(to_average)
846 print str(forcebase*(10**12))+' pN'
847 to_dump='forcebase '+self.current.path+' '+str(forcebase*(10**12))+' pN'
848 self.outlet.push(to_dump)
851 def do_slope(self,args):
855 Measures the slope of a delimited chunk on the return trace.
856 The chunk can be delimited either by two manual clicks, or have
857 a fixed width, given as an argument.
859 Syntax: slope [width]
860 The facultative [width] parameter specifies how many
861 points will be considered for the fit. If [width] is
862 specified, only one click will be required.
863 (c) Marco Brucale, Massimo Sandal 2008
866 # Reads the facultative width argument
872 # Decides between the two forms of user input, as per (args)
874 # Gets the Xs of two clicked points as indexes on the current curve vector
875 print 'Click twice to delimit chunk'
876 points=self._measure_N_points(N=2,whatset=1)
878 print 'Click once on the leftmost point of the chunk (i.e.usually the peak)'
879 points=self._measure_N_points(N=1,whatset=1)
881 slope=self._slope(points,fitspan)
883 # Outputs the relevant slope parameter
886 to_dump='slope '+self.current.path+' '+str(slope)
887 self.outlet.push(to_dump)
889 def _slope(self,points,fitspan):
890 # Calls the function linefit_between
891 parameters=[0,0,[],[]]
893 clickedpoints=[points[0].index,points[1].index]
896 clickedpoints=[points[0].index-fitspan,points[0].index]
899 parameters=self.linefit_between(clickedpoints[0],clickedpoints[1])
901 print 'Cannot fit. Did you click twice the same point?'
904 # Outputs the relevant slope parameter
906 print str(parameters[0])
907 to_dump='slope '+self.curve.path+' '+str(parameters[0])
908 self.outlet.push(to_dump)
910 # Makes a vector with the fitted parameters and sends it to the GUI
911 xtoplot=parameters[2]
915 ytoplot.append((x*parameters[0])+parameters[1])
917 clickvector_x, clickvector_y=[], []
919 clickvector_x.append(item.graph_coords[0])
920 clickvector_y.append(item.graph_coords[1])
922 lineplot=self._get_displayed_plot(0) #get topmost displayed plot
924 lineplot.add_set(xtoplot,ytoplot)
925 lineplot.add_set(clickvector_x, clickvector_y)
928 if lineplot.styles==[]:
929 lineplot.styles=[None,None,None,'scatter']
931 lineplot.styles+=[None,'scatter']
932 if lineplot.colors==[]:
933 lineplot.colors=[None,None,'black',None]
935 lineplot.colors+=['black',None]
938 self._send_plot([lineplot])
943 def linefit_between(self,index1,index2,whatset=1):
945 Creates two vectors (xtofit,ytofit) slicing out from the
946 current return trace a portion delimited by the two indexes
948 Then does a least squares linear fit on that slice.
949 Finally returns [0]=the slope, [1]=the intercept of the
950 fitted 1st grade polynomial, and [2,3]=the actual (x,y) vectors
952 (c) Marco Brucale, Massimo Sandal 2008
954 # Translates the indexes into two vectors containing the x,y data to fit
955 xtofit=self.plots[0].vectors[whatset][0][index1:index2]
956 ytofit=self.plots[0].vectors[whatset][1][index1:index2]
958 # Does the actual linear fitting (simple least squares with numpy.polyfit)
960 linefit=np.polyfit(xtofit,ytofit,1)
962 return (linefit[0],linefit[1],xtofit,ytofit)