1 # Copyright (C) 2010-2012 W. Trevor King <wking@drexel.edu>
3 # This file is part of Hooke.
5 # Hooke is free software: you can redistribute it and/or modify it under the
6 # terms of the GNU Lesser General Public License as published by the Free
7 # Software Foundation, either version 3 of the License, or (at your option) any
10 # Hooke is distributed in the hope that it will be useful, but WITHOUT ANY
11 # WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12 # A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
15 # You should have received a copy of the GNU Lesser General Public License
16 # along with Hooke. If not, see <http://www.gnu.org/licenses/>.
18 """The ``curve`` module provides :class:`CurvePlugin` and several
19 associated :class:`hooke.command.Command`\s for handling
20 :mod:`hooke.curve` classes.
30 from ..command import Command, Argument, Failure
31 from ..command_stack import CommandStack
32 from ..curve import Data
33 from ..engine import CommandMessage
34 from ..util.calculus import derivative
35 from ..util.fft import unitary_avg_power_spectrum
36 from ..util.si import ppSI, join_data_label, split_data_label
38 from .playlist import current_playlist_callback
41 # Define common or complicated arguments
43 def current_curve_callback(hooke, command, argument, value, load=True):
46 playlist = current_playlist_callback(hooke, command, argument, value)
47 curve = playlist.current(load=load)
49 raise Failure('No curves in %s' % playlist)
52 def unloaded_current_curve_callback(hooke, command, argument, value):
53 return current_curve_callback(
54 hooke=hooke, command=command, argument=argument, value=value,
57 CurveArgument = Argument(
58 name='curve', type='curve', callback=current_curve_callback,
60 :class:`hooke.curve.Curve` to act on. Defaults to the current curve
61 of the current playlist.
64 def _name_argument(name, default, help):
67 return Argument(name=name, type='string', default=default, help=help)
69 def block_argument(*args, **kwargs):
72 return _name_argument(*args, **kwargs)
74 def column_argument(*args, **kwargs):
77 return _name_argument(*args, **kwargs)
80 # Define useful command subclasses
82 class CurveCommand (Command):
83 """A :class:`~hooke.command.Command` operating on a
84 :class:`~hooke.curve.Curve`.
86 def __init__(self, **kwargs):
87 if 'arguments' in kwargs:
88 kwargs['arguments'].insert(0, CurveArgument)
90 kwargs['arguments'] = [CurveArgument]
91 super(CurveCommand, self).__init__(**kwargs)
93 def _curve(self, hooke, params):
94 """Get the selected curve.
98 `hooke` is intended to attach the selected curve to the local
99 playlist; the returned curve should not be effected by the
100 state of `hooke`. This is important for reliable
101 :class:`~hooke.command_stack.CommandStack`\s.
103 # HACK? rely on params['curve'] being bound to the local hooke
104 # playlist (i.e. not a copy, as you would get by passing a
105 # curve through the queue). Ugh. Stupid queues. As an
106 # alternative, we could pass lookup information through the
108 return params['curve']
110 def _add_to_command_stack(self, params):
111 """Store the command name and current `params` values in the
112 curve's `.command_stack`.
114 If this would duplicate the command currently on top of the
115 stack, no action is taken. Call early on, or watch out for
116 repeated param processing.
118 Recommended practice is to *not* lock in argument values that
119 are loaded from the plugin's :attr:`.config`.
123 Perhaps we should subclass :meth:`_run` and use :func:`super`,
124 or embed this in :meth:`run` to avoid subclasses calling this
125 method explicitly, with all the tedium and brittality that
126 implies. On the other hand, the current implemtnation allows
127 CurveCommands that don't effect the curve itself
128 (e.g. :class:`GetCommand`) to avoid adding themselves to the
131 if params['stack'] == True:
132 curve = self._curve(hooke=None, params=params)
133 if (len(curve.command_stack) > 0
134 and curve.command_stack[-1].command == self.name
135 and curve.command_stack[-1].arguments == params):
136 pass # no need to place duplicate calls on the stack.
138 curve.command_stack.append(CommandMessage(
139 self.name, dict(params)))
142 class BlockCommand (CurveCommand):
143 """A :class:`CurveCommand` operating on a :class:`~hooke.curve.Data` block.
145 def __init__(self, blocks=None, **kwargs):
147 blocks = [('block', None, 'Name of the data block to act on.')]
149 for name,default,help in blocks:
150 block_args.append(block_argument(name, default, help))
151 self._block_arguments = block_args
152 if 'arguments' not in kwargs:
153 kwargs['arguments'] = []
154 kwargs['arguments'] = block_args + kwargs['arguments']
155 super(BlockCommand, self).__init__(**kwargs)
157 def _block_names(self, hooke, params):
158 curve = self._curve(hooke, params)
159 return [b.info['name'] for b in curve.data]
161 def _block_index(self, hooke, params, name=None):
163 name = self._block_arguments[0].name
164 block_name = params[name]
165 if block_name == None:
166 curve = self._curve(hooke=hooke, params=params)
167 if len(curve.data) == 0:
168 raise Failure('no blocks in %s' % curve)
169 block_name = curve.data[0].info['name']
170 names = self._block_names(hooke=hooke, params=params)
172 return names.index(block_name)
173 except ValueError, e:
174 curve = self._curve(hooke, params)
175 raise Failure('no block named %s in %s (%s): %s'
176 % (block_name, curve, names, e))
178 def _block(self, hooke, params, name=None):
179 # HACK? rely on params['block'] being bound to the local hooke
180 # playlist (i.e. not a copy, as you would get by passing a
181 # block through the queue). Ugh. Stupid queues. As an
182 # alternative, we could pass lookup information through the
184 curve = self._curve(hooke, params)
185 index = self._block_index(hooke, params, name)
186 return curve.data[index]
189 class ColumnAccessCommand (BlockCommand):
190 """A :class:`BlockCommand` accessing a :class:`~hooke.curve.Data`
193 def __init__(self, columns=None, **kwargs):
195 columns = [('column', None, 'Name of the data column to act on.')]
197 for name,default,help in columns:
198 column_args.append(column_argument(name, default, help))
199 self._column_arguments = column_args
200 if 'arguments' not in kwargs:
201 kwargs['arguments'] = []
202 kwargs['arguments'] = column_args + kwargs['arguments']
203 super(ColumnAccessCommand, self).__init__(**kwargs)
205 def _get_column(self, hooke, params, block_name=None, column_name=None):
206 if column_name == None:
207 column_name = self._column_arguments[0].name
208 column_name = params[column_name]
209 if column_name is None:
211 block = self._block(hooke, params, block_name)
212 columns = block.info['columns']
214 column_index = columns.index(column_name)
215 except ValueError, e:
216 raise Failure('%s not in %s (%s): %s'
217 % (column_name, block.info['name'], columns, e))
218 return block[:,column_index]
221 class ColumnAddingCommand (ColumnAccessCommand):
222 """A :class:`ColumnAccessCommand` that also adds columns.
224 def __init__(self, new_columns=None, **kwargs):
225 if new_columns == None:
228 for name,default,help in new_columns:
229 column_args.append(column_argument(name, default, help))
230 self._new_column_arguments = column_args
231 if 'arguments' not in kwargs:
232 kwargs['arguments'] = []
233 kwargs['arguments'] = column_args + kwargs['arguments']
234 super(ColumnAddingCommand, self).__init__(**kwargs)
236 def _get_column(self, hooke, params, block_name=None, column_name=None):
237 if column_name == None and len(self._column_arguments) == 0:
238 column_name = self._new_column_arguments[0].name
239 return super(ColumnAddingCommand, self)._get_column(
240 hooke=hooke, params=params, block_name=block_name,
241 column_name=column_name)
243 def _set_column(self, hooke, params, block_name=None, column_name=None,
245 if column_name == None:
246 column_name = self._column_arguments[0].name
247 column_name = params[column_name]
248 block = self._block(hooke=hooke, params=params, name=block_name)
249 if column_name not in block.info['columns']:
250 new = Data((block.shape[0], block.shape[1]+1), dtype=block.dtype)
251 new.info = copy.deepcopy(block.info)
253 new.info['columns'].append(column_name)
255 block_index = self._block_index(hooke, params, name=block_name)
256 self._curve(hooke, params).data[block_index] = block
257 column_index = block.info['columns'].index(column_name)
258 block[:,column_index] = values
263 class CurvePlugin (Builtin):
265 super(CurvePlugin, self).__init__(name='curve')
267 GetCommand(self), InfoCommand(self), BlockInfoCommand(self),
268 DeltaCommand(self), ExportCommand(self), DifferenceCommand(self),
269 DerivativeCommand(self), PowerSpectrumCommand(self),
270 ScaledColumnAdditionCommand(self), ClearStackCommand(self)]
275 class GetCommand (CurveCommand):
276 """Return a :class:`hooke.curve.Curve`.
278 def __init__(self, plugin):
279 super(GetCommand, self).__init__(
280 name='get curve', help=self.__doc__, plugin=plugin)
282 def _run(self, hooke, inqueue, outqueue, params):
283 outqueue.put(self._curve(hooke, params))
286 class InfoCommand (CurveCommand):
287 """Get selected information about a :class:`hooke.curve.Curve`.
289 def __init__(self, plugin):
291 Argument(name='all', type='bool', default=False, count=1,
292 help='Get all curve information.'),
294 self.fields = ['name', 'path', 'driver', 'note', 'command stack',
295 'blocks', 'block names', 'block sizes']
296 for field in self.fields:
297 args.append(Argument(
298 name=field, type='bool', default=False, count=1,
299 help='Get curve %s' % field))
300 super(InfoCommand, self).__init__(
301 name='curve info', arguments=args,
302 help=self.__doc__, plugin=plugin)
304 def _run(self, hooke, inqueue, outqueue, params):
305 curve = self._curve(hooke, params)
307 for key in self.fields:
308 fields[key] = params[key]
309 if reduce(lambda x,y: x or y, fields.values()) == False:
310 params['all'] = True # No specific fields set, default to 'all'
311 if params['all'] == True:
312 for key in self.fields:
315 for key in self.fields:
316 if fields[key] == True:
317 get = getattr(self, '_get_%s' % key.replace(' ', '_'))
318 lines.append('%s: %s' % (key, get(curve)))
319 outqueue.put('\n'.join(lines))
321 def _get_name(self, curve):
324 def _get_path(self, curve):
327 def _get_driver(self, curve):
330 def _get_note(self, curve):
331 return curve.info.get('note', None)
333 def _get_command_stack(self, curve):
334 return curve.command_stack
336 def _get_blocks(self, curve):
337 return len(curve.data)
339 def _get_block_names(self, curve):
340 return [block.info['name'] for block in curve.data]
342 def _get_block_sizes(self, curve):
343 return [block.shape for block in curve.data]
346 class BlockInfoCommand (BlockCommand):
347 """Get selected information about a :class:`hooke.curve.Curve` data block.
349 def __init__(self, plugin):
350 super(BlockInfoCommand, self).__init__(
351 name='block info', arguments=[
353 name='key', count=-1, optional=False,
354 help='Dot-separted (.) key selection regexp.'),
358 File name for the output (appended).
361 help=self.__doc__, plugin=plugin)
363 def _run(self, hooke, inqueue, outqueue, params):
364 block = self._block(hooke, params)
365 values = {'index': self._block_index(hooke, params)}
366 for key in params['key']:
367 keys = [(0, key.split('.'), block.info)]
369 index,key_stack,info = keys.pop(0)
370 regexp = re.compile(key_stack[index])
372 for k,v in info.items():
375 new_stack = copy.copy(key_stack)
377 if index+1 == len(key_stack):
379 for k in new_stack[:-1]:
383 vals[new_stack[-1]] = v
385 keys.append((index+1, new_stack, v))
388 'no match found for %s (%s) in %s'
389 % (key_stack[index], key, sorted(info.keys())))
390 if params['output'] != None:
391 curve = self._curve(hooke, params)
392 with open(params['output'], 'a') as f:
393 yaml.dump({curve.name:{
395 block.info['name']: values
400 class DeltaCommand (BlockCommand):
401 """Get distance information between two points.
403 With two points A and B, the returned distances are A-B.
405 def __init__(self, plugin):
406 super(DeltaCommand, self).__init__(
409 Argument(name='point', type='point', optional=False, count=2,
411 Indicies of points bounding the selected data.
413 Argument(name='SI', type='bool', default=False,
415 Return distances in SI notation.
418 help=self.__doc__, plugin=plugin)
420 def _run(self, hooke, inqueue, outqueue, params):
421 data = self._block(hooke, params)
422 As = data[params['point'][0],:]
423 Bs = data[params['point'][1],:]
424 ds = [A-B for A,B in zip(As, Bs)]
425 if params['SI'] == False:
426 out = [(name, d) for name,d in zip(data.info['columns'], ds)]
429 for name,d in zip(data.info['columns'], ds):
430 n,units = split_data_label(name)
432 (n, ppSI(value=d, unit=units, decimals=2)))
436 class ExportCommand (BlockCommand):
437 """Export a :class:`hooke.curve.Curve` data block as TAB-delimeted
440 A "#" prefixed header will optionally appear at the beginning of
441 the file naming the columns.
443 def __init__(self, plugin):
444 super(ExportCommand, self).__init__(
447 Argument(name='output', type='file', default='curve.dat',
449 File name for the output data. Defaults to 'curve.dat'
451 Argument(name='header', type='bool', default=True,
453 True if you want the column-naming header line.
456 help=self.__doc__, plugin=plugin)
458 def _run(self, hooke, inqueue, outqueue, params):
459 data = self._block(hooke, params)
461 with open(os.path.expanduser(params['output']), 'w') as f:
462 if params['header'] == True:
463 f.write('# %s \n' % ('\t'.join(data.info['columns'])))
464 numpy.savetxt(f, data, delimiter='\t')
467 class DifferenceCommand (ColumnAddingCommand):
468 """Calculate the difference between two columns of data.
470 The difference is added to block A as a new column.
472 Note that the command will fail if the columns have different
473 lengths, so be careful when differencing columns from different
476 def __init__(self, plugin):
477 super(DifferenceCommand, self).__init__(
481 'Name of block A in A-B. Defaults to the first block'),
483 'Name of block B in A-B. Defaults to matching `block A`.'),
488 Column of data from block A to difference. Defaults to the first column.
492 Column of data from block B to difference. Defaults to matching `column A`.
496 ('output column', None,
498 Name of the new column for storing the difference (without units, defaults to
499 `difference of <block A> <column A> and <block B> <column B>`).
502 help=self.__doc__, plugin=plugin)
504 def _run(self, hooke, inqueue, outqueue, params):
505 self._add_to_command_stack(params)
506 params = self._setup_params(hooke=hooke, params=params)
507 data_A = self._get_column(hooke=hooke, params=params,
508 block_name='block A',
509 column_name='column A')
510 data_B = self._get_column(hooke=hooke, params=params,
511 block_name='block B',
512 column_name='column B')
513 out = data_A - data_B
514 self._set_column(hooke=hooke, params=params,
515 block_name='block A',
516 column_name='output column',
519 def _setup_params(self, hooke, params):
520 curve = self._curve(hooke, params)
521 if params['block A'] == None:
522 params['block A'] = curve.data[0].info['name']
523 if params['block B'] == None:
524 params['block B'] = params['block A']
525 block_A = self._block(hooke, params=params, name='block A')
526 block_B = self._block(hooke, params=params, name='block B')
527 if params['column A'] == None:
528 params['column A'] = block.info['columns'][0]
529 if params['column B'] == None:
530 params['column B'] = params['column A']
531 a_name,a_unit = split_data_label(params['column A'])
532 b_name,b_unit = split_data_label(params['column B'])
534 raise Failure('Unit missmatch: %s != %s' % (a_unit, b_unit))
535 if params['output column'] == None:
536 params['output column'] = join_data_label(
537 'difference of %s %s and %s %s' % (
538 block_A.info['name'], a_name,
539 block_B.info['name'], b_name),
542 params['output column'] = join_data_label(
543 params['output column'], a_unit)
547 class DerivativeCommand (ColumnAddingCommand):
548 """Calculate the derivative (actually, the discrete differentiation)
551 See :func:`hooke.util.calculus.derivative` for implementation
554 def __init__(self, plugin):
555 super(DerivativeCommand, self).__init__(
559 'Column of data block to differentiate with respect to.'),
561 'Column of data block to differentiate.'),
564 ('output column', None,
566 Name of the new column for storing the derivative (without units, defaults to
567 `derivative of <f column> with respect to <x column>`).
571 Argument(name='weights', type='dict', default={-1:-0.5, 1:0.5},
573 Weighting scheme dictionary for finite differencing. Defaults to
574 central differencing.
577 help=self.__doc__, plugin=plugin)
579 def _run(self, hooke, inqueue, outqueue, params):
580 self._add_to_command_stack(params)
581 params = self._setup_params(hooke=hooke, params=params)
582 x_data = self._get_column(hooke=hooke, params=params,
583 column_name='x column')
584 f_data = self._get_column(hooke=hooke, params=params,
585 column_name='f column')
587 x_data=x_data, f_data=f_data, weights=params['weights'])
588 self._set_column(hooke=hooke, params=params,
589 column_name='output column',
592 def _setup_params(self, hooke, params):
593 curve = self._curve(hooke, params)
594 x_name,x_unit = split_data_label(params['x column'])
595 f_name,f_unit = split_data_label(params['f column'])
596 d_unit = '%s/%s' % (f_unit, x_unit)
597 if params['output column'] == None:
598 params['output column'] = join_data_label(
599 'derivative of %s with respect to %s' % (
603 params['output column'] = join_data_label(
604 params['output column'], d_unit)
608 class PowerSpectrumCommand (ColumnAddingCommand):
609 """Calculate the power spectrum of a data column.
611 def __init__(self, plugin):
612 super(PowerSpectrumCommand, self).__init__(
613 name='power spectrum',
615 Argument(name='output block', type='string',
617 Name of the new data block for storing the power spectrum (defaults to
618 `power spectrum of <source block name> <source column name>`).
620 Argument(name='bounds', type='point', optional=True, count=2,
622 Indicies of points bounding the selected data.
624 Argument(name='freq', type='float', default=1.0,
628 Argument(name='freq units', type='string', default='Hz',
630 Units for the sampling frequency.
632 Argument(name='chunk size', type='int', default=2048,
634 Number of samples per chunk. Use a power of two.
636 Argument(name='overlap', type='bool', default=False,
638 If `True`, each chunk overlaps the previous chunk by half its length.
639 Otherwise, the chunks are end-to-end, and not overlapping.
642 help=self.__doc__, plugin=plugin)
644 def _run(self, hooke, inqueue, outqueue, params):
645 self._add_to_command_stack(params)
646 params = self._setup_params(hooke=hooke, params=params)
647 data = self._get_column(hooke=hooke, params=params)
648 bounds = params['bounds']
650 data = data[bounds[0]:bounds[1]]
651 freq_axis,power = unitary_avg_power_spectrum(
652 data, freq=params['freq'],
653 chunk_size=params['chunk size'],
654 overlap=params['overlap'])
655 b = Data((len(freq_axis),2), dtype=data.dtype)
656 b.info['name'] = params['output block']
657 b.info['columns'] = [
658 params['output freq column'],
659 params['output power column'],
661 self._curve(hooke, params).data.append(b)
662 self._set_column(hooke, params, block_name='output block',
663 column_name='output freq column',
665 self._set_column(hooke, params, block_name='output block',
666 column_name='output power column',
670 def _setup_params(self, hooke, params):
671 if params['output block'] in self._block_names(hooke, params):
672 raise Failure('output block %s already exists in %s.'
673 % (params['output block'],
674 self._curve(hooke, params)))
675 data = self._get_column(hooke=hooke, params=params)
676 d_name,d_unit = split_data_label(data.info['name'])
677 if params['output block'] == None:
678 params['output block'] = 'power spectrum of %s %s' % (
679 data.info['name'], params['column'])
680 self.params['output freq column'] = join_data_label(
681 'frequency axis', params['freq units'])
682 self.params['output power column'] = join_data_label(
683 'power density', '%s^2/%s' % (data_units, params['freq units']))
687 class ScaledColumnAdditionCommand (ColumnAddingCommand):
688 """Add one affine transformed column to another: `o=A*i1+B*i2+C`.
690 def __init__(self, plugin):
691 super(ScaledColumnAdditionCommand, self).__init__(
692 name='scaled column addition',
694 ('input column 1', 'input column (m)', """
695 Name of the first column to use as the transform input.
697 ('input column 2', None, """
698 Name of the second column to use as the transform input.
702 ('output column', 'output column (m)', """
703 Name of the column to use as the transform output.
707 Argument(name='scale 1', type='float', default=None,
709 A float value for the first scale constant.
711 Argument(name='scale 1 name', type='string', default=None,
713 The name of the first scale constant in the `.info` dictionary.
715 Argument(name='scale 2', type='float', default=None,
717 A float value for the second scale constant.
719 Argument(name='scale 2 name', type='string', default=None,
721 The name of the second scale constant in the `.info` dictionary.
723 Argument(name='constant', type='float', default=None,
725 A float value for the offset constant.
727 Argument(name='constant name', type='string', default=None,
729 The name of the offset constant in the `.info` dictionary.
732 help=self.__doc__, plugin=plugin)
734 def _run(self, hooke, inqueue, outqueue, params):
735 self._add_to_command_stack(params)
736 i1 = self._get_column(hooke=hooke, params=params,
737 column_name='input column 1')
738 i2 = self._get_column(hooke=hooke, params=params,
739 column_name='input column 2')
744 # what if both i1 and i2 are None?
745 a = self._get_constant(params, i1.info, 'scale 1')
746 b = self._get_constant(params, i2.info, 'scale 2')
747 c = self._get_constant(params, i1.info, 'constant')
748 out = a*i1 + b*i2 + c
749 self._set_column(hooke=hooke, params=params,
750 column_name='output column', values=out)
752 def _get_constant(self, params, info, name):
754 pname = params[name + ' name']
756 if pname is not None:
757 pname_entries = pname.split('|')
759 for entry in pname_entries:
761 if a is None and b is None:
770 class ClearStackCommand (CurveCommand):
771 """Empty a curve's command stack.
773 def __init__(self, plugin):
774 super(ClearStackCommand, self).__init__(
775 name='clear curve command stack',
776 help=self.__doc__, plugin=plugin)
777 i,arg = [(i,arg) for i,arg in enumerate(self.arguments)
778 if arg.name == 'curve'][0]
780 arg.callback = unloaded_current_curve_callback
781 self.arguments[i] = arg
783 def _run(self, hooke, inqueue, outqueue, params):
784 curve = self._curve(hooke, params)
785 curve.command_stack = CommandStack()
788 class OldCruft (object):
790 def do_forcebase(self,args):
794 Measures the difference in force (in pN) between a point and a baseline
795 took as the average between two points.
797 The baseline is fixed once for a given curve and different force measurements,
798 unless the user wants it to be recalculated
800 Syntax: forcebase [rebase]
801 rebase: Forces forcebase to ask again the baseline
802 max: Instead of asking for a point to measure, asks for two points and use
803 the maximum peak in between
805 rebase=False #if true=we select rebase
806 maxpoint=False #if true=we measure the maximum peak
808 plot=self._get_displayed_plot()
809 whatset=1 #fixme: for all sets
810 if 'rebase' in args or (self.basecurrent != self.current.path):
816 print 'Select baseline'
817 self.basepoints=self._measure_N_points(N=2, whatset=whatset)
818 self.basecurrent=self.current.path
821 print 'Select two points'
822 points=self._measure_N_points(N=2, whatset=whatset)
823 boundpoints=[points[0].index, points[1].index]
826 y=min(plot.vectors[whatset][1][boundpoints[0]:boundpoints[1]])
828 print 'Chosen interval not valid. Try picking it again. Did you pick the same point as begin and end of interval?'
830 print 'Select point to measure'
831 points=self._measure_N_points(N=1, whatset=whatset)
832 #whatplot=points[0].dest
833 y=points[0].graph_coords[1]
835 #fixme: code duplication
836 boundaries=[self.basepoints[0].index, self.basepoints[1].index]
838 to_average=plot.vectors[whatset][1][boundaries[0]:boundaries[1]] #y points to average
840 avg=np.mean(to_average)
842 print str(forcebase*(10**12))+' pN'
843 to_dump='forcebase '+self.current.path+' '+str(forcebase*(10**12))+' pN'
844 self.outlet.push(to_dump)
847 def do_slope(self,args):
851 Measures the slope of a delimited chunk on the return trace.
852 The chunk can be delimited either by two manual clicks, or have
853 a fixed width, given as an argument.
855 Syntax: slope [width]
856 The facultative [width] parameter specifies how many
857 points will be considered for the fit. If [width] is
858 specified, only one click will be required.
859 (c) Marco Brucale, Massimo Sandal 2008
862 # Reads the facultative width argument
868 # Decides between the two forms of user input, as per (args)
870 # Gets the Xs of two clicked points as indexes on the current curve vector
871 print 'Click twice to delimit chunk'
872 points=self._measure_N_points(N=2,whatset=1)
874 print 'Click once on the leftmost point of the chunk (i.e.usually the peak)'
875 points=self._measure_N_points(N=1,whatset=1)
877 slope=self._slope(points,fitspan)
879 # Outputs the relevant slope parameter
882 to_dump='slope '+self.current.path+' '+str(slope)
883 self.outlet.push(to_dump)
885 def _slope(self,points,fitspan):
886 # Calls the function linefit_between
887 parameters=[0,0,[],[]]
889 clickedpoints=[points[0].index,points[1].index]
892 clickedpoints=[points[0].index-fitspan,points[0].index]
895 parameters=self.linefit_between(clickedpoints[0],clickedpoints[1])
897 print 'Cannot fit. Did you click twice the same point?'
900 # Outputs the relevant slope parameter
902 print str(parameters[0])
903 to_dump='slope '+self.curve.path+' '+str(parameters[0])
904 self.outlet.push(to_dump)
906 # Makes a vector with the fitted parameters and sends it to the GUI
907 xtoplot=parameters[2]
911 ytoplot.append((x*parameters[0])+parameters[1])
913 clickvector_x, clickvector_y=[], []
915 clickvector_x.append(item.graph_coords[0])
916 clickvector_y.append(item.graph_coords[1])
918 lineplot=self._get_displayed_plot(0) #get topmost displayed plot
920 lineplot.add_set(xtoplot,ytoplot)
921 lineplot.add_set(clickvector_x, clickvector_y)
924 if lineplot.styles==[]:
925 lineplot.styles=[None,None,None,'scatter']
927 lineplot.styles+=[None,'scatter']
928 if lineplot.colors==[]:
929 lineplot.colors=[None,None,'black',None]
931 lineplot.colors+=['black',None]
934 self._send_plot([lineplot])
939 def linefit_between(self,index1,index2,whatset=1):
941 Creates two vectors (xtofit,ytofit) slicing out from the
942 current return trace a portion delimited by the two indexes
944 Then does a least squares linear fit on that slice.
945 Finally returns [0]=the slope, [1]=the intercept of the
946 fitted 1st grade polynomial, and [2,3]=the actual (x,y) vectors
948 (c) Marco Brucale, Massimo Sandal 2008
950 # Translates the indexes into two vectors containing the x,y data to fit
951 xtofit=self.plots[0].vectors[whatset][0][index1:index2]
952 ytofit=self.plots[0].vectors[whatset][1][index1:index2]
954 # Does the actual linear fitting (simple least squares with numpy.polyfit)
956 linefit=np.polyfit(xtofit,ytofit,1)
958 return (linefit[0],linefit[1],xtofit,ytofit)