3ef141ce0b14aa2d5d707f9108dccfd132c9e5cf
[hooke.git] / hooke / plugin / curve.py
1 # Copyright (C) 2008-2010 Alberto Gomez-Casado
2 #                         Fabrizio Benedetti
3 #                         Massimo Sandal <devicerandom@gmail.com>
4 #                         W. Trevor King <wking@drexel.edu>
5 #
6 # This file is part of Hooke.
7 #
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.
12 #
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.
17 #
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/>.
21
22 """The ``curve`` module provides :class:`CurvePlugin` and several
23 associated :class:`hooke.command.Command`\s for handling
24 :mod:`hooke.curve` classes.
25 """
26
27 import copy
28 import re
29
30 import numpy
31 import yaml
32
33 from ..command import Command, Argument, Failure
34 from ..command_stack import CommandStack
35 from ..curve import Data
36 from ..engine import CommandMessage
37 from ..util.calculus import derivative
38 from ..util.fft import unitary_avg_power_spectrum
39 from ..util.si import ppSI, join_data_label, split_data_label
40 from . import Builtin
41 from .playlist import current_playlist_callback
42
43
44 # Define common or complicated arguments
45
46 def current_curve_callback(hooke, command, argument, value, load=True):
47     if value != None:
48         return value
49     playlist = current_playlist_callback(hooke, command, argument, value)
50     curve = playlist.current(load=load)
51     if curve == None:
52         raise Failure('No curves in %s' % playlist)
53     return curve
54
55 def unloaded_current_curve_callback(hooke, command, argument, value):
56     return current_curve_callback(
57         hooke=hooke, command=command, argument=argument, value=value,
58         load=False)
59
60 CurveArgument = Argument(
61     name='curve', type='curve', callback=current_curve_callback,
62     help="""
63 :class:`hooke.curve.Curve` to act on.  Defaults to the current curve
64 of the current playlist.
65 """.strip())
66
67 def _name_argument(name, default, help):
68     """TODO
69     """
70     return Argument(name=name, type='string', default=default, help=help)
71
72 def block_argument(*args, **kwargs):
73     """TODO
74     """
75     return _name_argument(*args, **kwargs)
76
77 def column_argument(*args, **kwargs):
78     """TODO
79     """
80     return _name_argument(*args, **kwargs)
81
82
83 # Define useful command subclasses
84
85 class CurveCommand (Command):
86     """A :class:`~hooke.command.Command` operating on a
87     :class:`~hooke.curve.Curve`.
88     """
89     def __init__(self, **kwargs):
90         if 'arguments' in kwargs:
91             kwargs['arguments'].insert(0, CurveArgument)
92         else:
93             kwargs['arguments'] = [CurveArgument]
94         super(CurveCommand, self).__init__(**kwargs)
95
96     def _curve(self, hooke, params):
97         """Get the selected curve.
98
99         Notes
100         -----
101         `hooke` is intended to attach the selected curve to the local
102         playlist; the returned curve should not be effected by the
103         state of `hooke`.  This is important for reliable
104         :class:`~hooke.command_stack.CommandStack`\s.
105         """
106         # HACK? rely on params['curve'] being bound to the local hooke
107         # playlist (i.e. not a copy, as you would get by passing a
108         # curve through the queue).  Ugh.  Stupid queues.  As an
109         # alternative, we could pass lookup information through the
110         # queue...
111         return params['curve']
112
113     def _add_to_command_stack(self, params):
114         """Store the command name and current `params` values in the
115         curve's `.command_stack`.
116
117         If this would duplicate the command currently on top of the
118         stack, no action is taken.  Call early on, or watch out for
119         repeated param processing.
120
121         Recommended practice is to *not* lock in argument values that
122         are loaded from the plugin's :attr:`.config`.
123
124         Notes
125         -----
126         Perhaps we should subclass :meth:`_run` and use :func:`super`,
127         or embed this in :meth:`run` to avoid subclasses calling this
128         method explicitly, with all the tedium and brittality that
129         implies.  On the other hand, the current implemtnation allows
130         CurveCommands that don't effect the curve itself
131         (e.g. :class:`GetCommand`) to avoid adding themselves to the
132         stack entirely.
133         """
134         if params['stack'] == True:
135             curve = self._curve(hooke=None, params=params)
136             if (len(curve.command_stack) > 0
137                 and curve.command_stack[-1].command == self.name
138                 and curve.command_stack[-1].arguments == params):
139                 pass  # no need to place duplicate calls on the stack.
140             else:
141                 curve.command_stack.append(CommandMessage(
142                         self.name, dict(params)))
143
144
145 class BlockCommand (CurveCommand):
146     """A :class:`CurveCommand` operating on a :class:`~hooke.curve.Data` block.
147     """
148     def __init__(self, blocks=None, **kwargs):
149         if blocks == None:
150             blocks = [('block', None, 'Name of the data block to act on.')]
151         block_args = []
152         for name,default,help in blocks:
153             block_args.append(block_argument(name, default, help))
154         self._block_arguments = block_args
155         if 'arguments' not in kwargs:
156             kwargs['arguments'] = []
157         kwargs['arguments'] = block_args + kwargs['arguments']
158         super(BlockCommand, self).__init__(**kwargs)
159
160     def _block_names(self, hooke, params):
161         curve = self._curve(hooke, params)
162         return [b.info['name'] for b in curve.data]
163
164     def _block_index(self, hooke, params, name=None):
165         if name == None:
166             name = self._block_arguments[0].name
167         block_name = params[name]
168         if block_name == None:
169             curve = self._curve(hooke=hooke, params=params)
170             if len(curve.data) == 0:
171                 raise Failure('no blocks in %s' % curve)
172             block_name = curve.data[0].info['name']
173         names = self._block_names(hooke=hooke, params=params)
174         try:
175             return names.index(block_name)
176         except ValueError, e:
177             curve = self._curve(hooke, params)
178             raise Failure('no block named %s in %s (%s): %s'
179                           % (block_name, curve, names, e))
180
181     def _block(self, hooke, params, name=None):
182         # HACK? rely on params['block'] being bound to the local hooke
183         # playlist (i.e. not a copy, as you would get by passing a
184         # block through the queue).  Ugh.  Stupid queues.  As an
185         # alternative, we could pass lookup information through the
186         # queue...
187         curve = self._curve(hooke, params)
188         index = self._block_index(hooke, params, name)
189         return curve.data[index]
190
191
192 class ColumnAccessCommand (BlockCommand):
193     """A :class:`BlockCommand` accessing a :class:`~hooke.curve.Data`
194     block column.
195     """
196     def __init__(self, columns=None, **kwargs):
197         if columns == None:
198             columns = [('column', None, 'Name of the data column to act on.')]
199         column_args = []
200         for name,default,help in columns:
201             column_args.append(column_argument(name, default, help))
202         self._column_arguments = column_args
203         if 'arguments' not in kwargs:
204             kwargs['arguments'] = []
205         kwargs['arguments'] = column_args + kwargs['arguments']
206         super(ColumnAccessCommand, self).__init__(**kwargs)
207
208     def _get_column(self, hooke, params, block_name=None, column_name=None):
209         if column_name == None:
210             column_name = self._column_arguments[0].name
211         column_name = params[column_name]
212         block = self._block(hooke, params, block_name)
213         columns = block.info['columns']
214         try:
215             column_index = columns.index(column_name)
216         except ValueError, e:
217             raise Failure('%s not in %s (%s): %s'
218                           % (column_name, block.info['name'], columns, e))
219         return block[:,column_index]
220
221
222 class ColumnAddingCommand (ColumnAccessCommand):
223     """A :class:`ColumnAccessCommand` that also adds columns.
224     """
225     def __init__(self, new_columns=None, **kwargs):
226         if new_columns == None:
227             new_columns = []
228         column_args = []
229         for name,default,help in new_columns:
230             column_args.append(column_argument(name, default, help))
231         self._new_column_arguments = column_args
232         if 'arguments' not in kwargs:
233             kwargs['arguments'] = []
234         kwargs['arguments'] = column_args + kwargs['arguments']
235         super(ColumnAddingCommand, self).__init__(**kwargs)
236
237     def _get_column(self, hooke, params, block_name=None, column_name=None):
238         if column_name == None and len(self._column_arguments) == 0:
239             column_name = self._new_column_arguments[0].name
240         return super(ColumnAddingCommand, self)._get_column(
241             hooke=hooke, params=params, block_name=block_name,
242             column_name=column_name)
243
244     def _set_column(self, hooke, params, block_name=None, column_name=None,
245                     values=None):
246         if column_name == None:
247             column_name = self._column_arguments[0].name
248         column_name = params[column_name]
249         block = self._block(hooke=hooke, params=params, name=block_name)
250         if column_name not in block.info['columns']:
251             new = Data((block.shape[0], block.shape[1]+1), dtype=block.dtype)
252             new.info = copy.deepcopy(block.info)
253             new[:,:-1] = block
254             new.info['columns'].append(column_name)
255             block = new
256             block_index = self._block_index(hooke, params, name=block_name)
257             self._curve(hooke, params).data[block_index] = block
258         column_index = block.info['columns'].index(column_name)
259         block[:,column_index] = values
260
261
262 # The plugin itself
263
264 class CurvePlugin (Builtin):
265     def __init__(self):
266         super(CurvePlugin, self).__init__(name='curve')
267         self._commands = [
268             GetCommand(self), InfoCommand(self), BlockInfoCommand(self),
269             DeltaCommand(self), ExportCommand(self), DifferenceCommand(self),
270             DerivativeCommand(self), PowerSpectrumCommand(self),
271             ClearStackCommand(self)]
272
273
274 # Define commands
275
276 class GetCommand (CurveCommand):
277     """Return a :class:`hooke.curve.Curve`.
278     """
279     def __init__(self, plugin):
280         super(GetCommand, self).__init__(
281             name='get curve', help=self.__doc__, plugin=plugin)
282
283     def _run(self, hooke, inqueue, outqueue, params):
284         outqueue.put(self._curve(hooke, params))
285
286
287 class InfoCommand (CurveCommand):
288     """Get selected information about a :class:`hooke.curve.Curve`.
289     """
290     def __init__(self, plugin):
291         args = [
292             Argument(name='all', type='bool', default=False, count=1,
293                      help='Get all curve information.'),
294             ]
295         self.fields = ['name', 'path', 'driver', 'note', 'command stack',
296                        'blocks', 'block names', 'block sizes']
297         for field in self.fields:
298             args.append(Argument(
299                     name=field, type='bool', default=False, count=1,
300                     help='Get curve %s' % field))
301         super(InfoCommand, self).__init__(
302             name='curve info', arguments=args,
303             help=self.__doc__, plugin=plugin)
304
305     def _run(self, hooke, inqueue, outqueue, params):
306         curve = self._curve(hooke, params)
307         fields = {}
308         for key in self.fields:
309             fields[key] = params[key]
310         if reduce(lambda x,y: x or y, fields.values()) == False:
311             params['all'] = True # No specific fields set, default to 'all'
312         if params['all'] == True:
313             for key in self.fields:
314                 fields[key] = True
315         lines = []
316         for key in self.fields:
317             if fields[key] == True:
318                 get = getattr(self, '_get_%s' % key.replace(' ', '_'))
319                 lines.append('%s: %s' % (key, get(curve)))
320         outqueue.put('\n'.join(lines))
321
322     def _get_name(self, curve):
323         return curve.name
324
325     def _get_path(self, curve):
326         return curve.path
327
328     def _get_driver(self, curve):
329         return curve.driver
330
331     def _get_note(self, curve):
332         return curve.info.get('note', None)
333                               
334     def _get_command_stack(self, curve):
335         return curve.command_stack
336
337     def _get_blocks(self, curve):
338         return len(curve.data)
339
340     def _get_block_names(self, curve):
341         return [block.info['name'] for block in curve.data]
342
343     def _get_block_sizes(self, curve):
344         return [block.shape for block in curve.data]
345
346
347 class BlockInfoCommand (BlockCommand):
348     """Get selected information about a :class:`hooke.curve.Curve` data block.
349     """
350     def __init__(self, plugin):
351         super(BlockInfoCommand, self).__init__(
352             name='block info', arguments=[
353                 Argument(
354                     name='key', count=-1, optional=False,
355                     help='Dot-separted (.) key selection regexp.'),
356                 Argument(
357                     name='output',
358                     help="""
359 File name for the output (appended).
360 """.strip()),
361                 ],
362             help=self.__doc__, plugin=plugin)
363
364     def _run(self, hooke, inqueue, outqueue, params):
365         block = self._block(hooke, params)
366         values = {'index': self._block_index(hooke, params)}
367         for key in params['key']:
368             keys = [(0, key.split('.'), block.info)]
369             while len(keys) > 0:
370                 index,key_stack,info = keys.pop(0)
371                 regexp = re.compile(key_stack[index])
372                 matched = False
373                 for k,v in info.items():
374                     if regexp.match(k):
375                         matched = True
376                         new_stack = copy.copy(key_stack)
377                         new_stack[index] = k
378                         if index+1 == len(key_stack):
379                             vals = values
380                             for k in new_stack[:-1]:
381                                 if k not in vals:
382                                     vals[k] = {}
383                                 vals = vals[k]
384                             vals[new_stack[-1]] = v
385                         else:
386                             keys.append((index+1, new_stack, v))
387                 if matched == False:
388                     raise ValueError(
389                         'no match found for %s (%s) in %s'
390                         % (key_stack[index], key, sorted(info.keys())))
391         if params['output'] != None:
392             curve = self._curve(hooke, params)
393             with open(params['output'], 'a') as f:
394                 yaml.dump({curve.name:{
395                             'path': curve.path,
396                             block.info['name']: values
397                             }}, f)
398         outqueue.put(values)
399
400
401 class DeltaCommand (BlockCommand):
402     """Get distance information between two points.
403
404     With two points A and B, the returned distances are A-B.
405     """
406     def __init__(self, plugin):
407         super(DeltaCommand, self).__init__(
408             name='delta',
409             arguments=[
410                 Argument(name='point', type='point', optional=False, count=2,
411                          help="""
412 Indicies of points bounding the selected data.
413 """.strip()),
414                 Argument(name='SI', type='bool', default=False,
415                          help="""
416 Return distances in SI notation.
417 """.strip())
418                 ],
419             help=self.__doc__, plugin=plugin)
420
421     def _run(self, hooke, inqueue, outqueue, params):
422         data = self._block(hooke, params)
423         As = data[params['point'][0],:]
424         Bs = data[params['point'][1],:]
425         ds = [A-B for A,B in zip(As, Bs)]
426         if params['SI'] == False:
427             out = [(name, d) for name,d in zip(data.info['columns'], ds)]
428         else:
429             out = []
430             for name,d in zip(data.info['columns'], ds):
431                 n,units = split_data_label(name)
432                 out.append(
433                   (n, ppSI(value=d, unit=units, decimals=2)))
434         outqueue.put(out)
435
436
437 class ExportCommand (BlockCommand):
438     """Export a :class:`hooke.curve.Curve` data block as TAB-delimeted
439     ASCII text.
440
441     A "#" prefixed header will optionally appear at the beginning of
442     the file naming the columns.
443     """
444     def __init__(self, plugin):
445         super(ExportCommand, self).__init__(
446             name='export block',
447             arguments=[
448                 Argument(name='output', type='file', default='curve.dat',
449                          help="""
450 File name for the output data.  Defaults to 'curve.dat'
451 """.strip()),
452                 Argument(name='header', type='bool', default=True,
453                          help="""
454 True if you want the column-naming header line.
455 """.strip()),
456                 ],
457             help=self.__doc__, plugin=plugin)
458
459     def _run(self, hooke, inqueue, outqueue, params):
460         data = self._block(hooke, params)
461
462         with open(params['output'], 'w') as f:
463             if params['header'] == True:
464                 f.write('# %s \n' % ('\t'.join(data.info['columns'])))
465             numpy.savetxt(f, data, delimiter='\t')
466
467
468 class DifferenceCommand (ColumnAddingCommand):
469     """Calculate the difference between two columns of data.
470
471     The difference is added to block A as a new column.
472
473     Note that the command will fail if the columns have different
474     lengths, so be careful when differencing columns from different
475     blocks.
476     """
477     def __init__(self, plugin):
478         super(DifferenceCommand, self).__init__(
479             name='difference',
480             blocks=[
481                 ('block A', None,
482                  'Name of block A in A-B.  Defaults to the first block'),
483                 ('block B', None,
484                  'Name of block B in A-B.  Defaults to matching `block A`.'),
485                 ],
486             columns=[
487                 ('column A', None,
488                  """
489 Column of data from block A to difference.  Defaults to the first column.
490 """.strip()),
491                 ('column B', None,
492                  """
493 Column of data from block B to difference.  Defaults to matching `column A`.
494 """.strip()),
495                 ],
496             new_columns=[
497                 ('output column', None,
498                  """
499 Name of the new column for storing the difference (without units, defaults to
500 `difference of <block A> <column A> and <block B> <column B>`).
501 """.strip()),
502                 ],
503             help=self.__doc__, plugin=plugin)
504
505     def _run(self, hooke, inqueue, outqueue, params):
506         self._add_to_command_stack(params)
507         params = self._setup_params(hooke=hooke, params=params)
508         data_A = self._get_column(hooke=hooke, params=params,
509                                   block_name='block A',
510                                   column_name='column A')
511         data_B = self._get_column(hooke=hooke, params=params,
512                                   block_name='block B',
513                                   column_name='column B')
514         out = data_A - data_B
515         self._set_column(hooke=hooke, params=params,
516                          block_name='block A',
517                          column_name='output column',
518                          values=out)
519
520     def _setup_params(self, hooke, params):
521         curve = self._curve(hooke, params)
522         if params['block A'] == None:
523             params['block A'] = curve.data[0].info['name']
524         if params['block B'] == None:
525             params['block B'] = params['block A']
526         block_A = self._block(hooke, params=params, name='block A')
527         block_B = self._block(hooke, params=params, name='block B')
528         if params['column A'] == None:
529             params['column A'] = block.info['columns'][0]
530         if params['column B'] == None:
531             params['column B'] = params['column A']
532         a_name,a_unit = split_data_label(params['column A'])
533         b_name,b_unit = split_data_label(params['column B'])
534         if a_unit != b_unit:
535             raise Failure('Unit missmatch: %s != %s' % (a_unit, b_unit))
536         if params['output column'] == None:
537             params['output column'] = join_data_label(
538                 'difference of %s %s and %s %s' % (
539                     block_A.info['name'], a_name,
540                     block_B.info['name'], b_name),
541                 a_unit)
542         else:
543             params['output column'] = join_data_label(
544                 params['output column'], a_unit)
545         return params
546
547
548 class DerivativeCommand (ColumnAddingCommand):
549     """Calculate the derivative (actually, the discrete differentiation)
550     of a data column.
551
552     See :func:`hooke.util.calculus.derivative` for implementation
553     details.
554     """
555     def __init__(self, plugin):
556         super(DerivativeCommand, self).__init__(
557             name='derivative',
558             columns=[
559                 ('x column', None,
560                  'Column of data block to differentiate with respect to.'),
561                 ('f column', None,
562                  'Column of data block to differentiate.'),
563                 ],
564             new_columns=[
565                 ('output column', None,
566                  """
567 Name of the new column for storing the derivative (without units, defaults to
568 `derivative of <f column> with respect to <x column>`).
569 """.strip()),
570                 ],
571             arguments=[
572                 Argument(name='weights', type='dict', default={-1:-0.5, 1:0.5},
573                          help="""
574 Weighting scheme dictionary for finite differencing.  Defaults to
575 central differencing.
576 """.strip()),
577                 ],
578             help=self.__doc__, plugin=plugin)
579
580     def _run(self, hooke, inqueue, outqueue, params):
581         self._add_to_command_stack(params)
582         params = self._setup_params(hooke=hooke, params=params)
583         x_data = self._get_column(hooke=hooke, params=params,
584                                   column_name='x column')
585         f_data = self._get_column(hooke=hooke, params=params,
586                                   column_name='f column')
587         d = derivative(
588             x_data=x_data, f_data=f_data, weights=params['weights'])
589         self._set_column(hooke=hooke, params=params,
590                          column_name='output column',
591                          values=d)
592
593     def _setup_params(self, hooke, params):
594         curve = self._curve(hooke, params)
595         x_name,x_unit = split_data_label(params['x column'])
596         f_name,f_unit = split_data_label(params['f column'])
597         d_unit = '%s/%s' % (f_unit, x_unit)
598         if params['output column'] == None:
599             params['output column'] = join_data_label(
600                 'derivative of %s with respect to %s' % (
601                     f_name, x_name),
602                 d_unit)
603         else:
604             params['output column'] = join_data_label(
605                 params['output column'], d_unit)
606         return params
607
608
609 class PowerSpectrumCommand (ColumnAddingCommand):
610     """Calculate the power spectrum of a data column.
611     """
612     def __init__(self, plugin):
613         super(PowerSpectrumCommand, self).__init__(
614             name='power spectrum',
615             arguments=[
616                 Argument(name='output block', type='string',
617                          help="""
618 Name of the new data block for storing the power spectrum (defaults to
619 `power spectrum of <source block name> <source column name>`).
620 """.strip()),
621                 Argument(name='bounds', type='point', optional=True, count=2,
622                          help="""
623 Indicies of points bounding the selected data.
624 """.strip()),
625                 Argument(name='freq', type='float', default=1.0,
626                          help="""
627 Sampling frequency.
628 """.strip()),
629                 Argument(name='freq units', type='string', default='Hz',
630                          help="""
631 Units for the sampling frequency.
632 """.strip()),
633                 Argument(name='chunk size', type='int', default=2048,
634                          help="""
635 Number of samples per chunk.  Use a power of two.
636 """.strip()),
637                 Argument(name='overlap', type='bool', default=False,
638                          help="""
639 If `True`, each chunk overlaps the previous chunk by half its length.
640 Otherwise, the chunks are end-to-end, and not overlapping.
641 """.strip()),
642                 ],
643             help=self.__doc__, plugin=plugin)
644
645     def _run(self, hooke, inqueue, outqueue, params):
646         self._add_to_command_stack(params)
647         params = self._setup_params(hooke=hooke, params=params)
648         data = self._get_column(hooke=hooke, params=params)
649         bounds = params['bounds']
650         if bounds != None:
651             data = data[bounds[0]:bounds[1]]
652         freq_axis,power = unitary_avg_power_spectrum(
653             data, freq=params['freq'],
654             chunk_size=params['chunk size'],
655             overlap=params['overlap'])
656         b = Data((len(freq_axis),2), dtype=data.dtype)
657         b.info['name'] = params['output block']
658         b.info['columns'] = [
659             params['output freq column'],
660             params['output power column'],
661             ]
662         self._curve(hooke, params).data.append(b)
663         self._set_column(hooke, params, block_name='output block',
664                          column_name='output freq column',
665                          values=freq_axis)
666         self._set_column(hooke, params, block_name='output block',
667                          column_name='output power column',
668                          values=power)
669         outqueue.put(b)
670
671     def _setup_params(self, hooke, params):
672         if params['output block'] in self._block_names(hooke, params):
673             raise Failure('output block %s already exists in %s.'
674                           % (params['output block'],
675                              self._curve(hooke, params)))
676         data = self._get_column(hooke=hooke, params=params)
677         d_name,d_unit = split_data_label(data.info['name'])
678         if params['output block'] == None:
679             params['output block'] = 'power spectrum of %s %s' % (
680                 data.info['name'], params['column'])
681         self.params['output freq column'] = join_data_label(
682             'frequency axis', params['freq units'])
683         self.params['output power column'] = join_data_label(
684             'power density', '%s^2/%s' % (data_units, params['freq units']))
685         return params
686
687
688 class ClearStackCommand (CurveCommand):
689     """Empty a curve's command stack.
690     """
691     def __init__(self, plugin):
692         super(ClearStackCommand, self).__init__(
693             name='clear curve command stack',
694             help=self.__doc__, plugin=plugin)
695         i,arg = [(i,arg) for i,arg in enumerate(self.arguments)
696                  if arg.name == 'curve'][0]
697         arg = copy.copy(arg)
698         arg.callback = unloaded_current_curve_callback
699         self.arguments[i] = arg
700
701     def _run(self, hooke, inqueue, outqueue, params):
702         curve = self._curve(hooke, params)
703         curve.command_stack = CommandStack()
704
705
706 class OldCruft (object):
707
708     def do_forcebase(self,args):
709         '''
710         FORCEBASE
711         (generalvclamp.py)
712         Measures the difference in force (in pN) between a point and a baseline
713         took as the average between two points.
714
715         The baseline is fixed once for a given curve and different force measurements,
716         unless the user wants it to be recalculated
717         ------------
718         Syntax: forcebase [rebase]
719                 rebase: Forces forcebase to ask again the baseline
720                 max: Instead of asking for a point to measure, asks for two points and use
721                      the maximum peak in between
722         '''
723         rebase=False #if true=we select rebase
724         maxpoint=False #if true=we measure the maximum peak
725
726         plot=self._get_displayed_plot()
727         whatset=1 #fixme: for all sets
728         if 'rebase' in args or (self.basecurrent != self.current.path):
729             rebase=True
730         if 'max' in args:
731             maxpoint=True
732
733         if rebase:
734             print 'Select baseline'
735             self.basepoints=self._measure_N_points(N=2, whatset=whatset)
736             self.basecurrent=self.current.path
737
738         if maxpoint:
739             print 'Select two points'
740             points=self._measure_N_points(N=2, whatset=whatset)
741             boundpoints=[points[0].index, points[1].index]
742             boundpoints.sort()
743             try:
744                 y=min(plot.vectors[whatset][1][boundpoints[0]:boundpoints[1]])
745             except ValueError:
746                 print 'Chosen interval not valid. Try picking it again. Did you pick the same point as begin and end of interval?'
747         else:
748             print 'Select point to measure'
749             points=self._measure_N_points(N=1, whatset=whatset)
750             #whatplot=points[0].dest
751             y=points[0].graph_coords[1]
752
753         #fixme: code duplication
754         boundaries=[self.basepoints[0].index, self.basepoints[1].index]
755         boundaries.sort()
756         to_average=plot.vectors[whatset][1][boundaries[0]:boundaries[1]] #y points to average
757
758         avg=np.mean(to_average)
759         forcebase=abs(y-avg)
760         print str(forcebase*(10**12))+' pN'
761         to_dump='forcebase '+self.current.path+' '+str(forcebase*(10**12))+' pN'
762         self.outlet.push(to_dump)
763
764     #---SLOPE---
765     def do_slope(self,args):
766         '''
767         SLOPE
768         (generalvclamp.py)
769         Measures the slope of a delimited chunk on the return trace.
770         The chunk can be delimited either by two manual clicks, or have
771         a fixed width, given as an argument.
772         ---------------
773         Syntax: slope [width]
774                 The facultative [width] parameter specifies how many
775                 points will be considered for the fit. If [width] is
776                 specified, only one click will be required.
777         (c) Marco Brucale, Massimo Sandal 2008
778         '''
779
780         # Reads the facultative width argument
781         try:
782             fitspan=int(args)
783         except:
784             fitspan=0
785
786         # Decides between the two forms of user input, as per (args)
787         if fitspan == 0:
788             # Gets the Xs of two clicked points as indexes on the current curve vector
789             print 'Click twice to delimit chunk'
790             points=self._measure_N_points(N=2,whatset=1)
791         else:
792             print 'Click once on the leftmost point of the chunk (i.e.usually the peak)'
793             points=self._measure_N_points(N=1,whatset=1)
794             
795         slope=self._slope(points,fitspan)
796
797         # Outputs the relevant slope parameter
798         print 'Slope:'
799         print str(slope)
800         to_dump='slope '+self.current.path+' '+str(slope)
801         self.outlet.push(to_dump)
802
803     def _slope(self,points,fitspan):
804         # Calls the function linefit_between
805         parameters=[0,0,[],[]]
806         try:
807             clickedpoints=[points[0].index,points[1].index]
808             clickedpoints.sort()
809         except:
810             clickedpoints=[points[0].index-fitspan,points[0].index]        
811
812         try:
813             parameters=self.linefit_between(clickedpoints[0],clickedpoints[1])
814         except:
815             print 'Cannot fit. Did you click twice the same point?'
816             return
817              
818         # Outputs the relevant slope parameter
819         print 'Slope:'
820         print str(parameters[0])
821         to_dump='slope '+self.curve.path+' '+str(parameters[0])
822         self.outlet.push(to_dump)
823
824         # Makes a vector with the fitted parameters and sends it to the GUI
825         xtoplot=parameters[2]
826         ytoplot=[]
827         x=0
828         for x in xtoplot:
829             ytoplot.append((x*parameters[0])+parameters[1])
830
831         clickvector_x, clickvector_y=[], []
832         for item in points:
833             clickvector_x.append(item.graph_coords[0])
834             clickvector_y.append(item.graph_coords[1])
835
836         lineplot=self._get_displayed_plot(0) #get topmost displayed plot
837
838         lineplot.add_set(xtoplot,ytoplot)
839         lineplot.add_set(clickvector_x, clickvector_y)
840
841
842         if lineplot.styles==[]:
843             lineplot.styles=[None,None,None,'scatter']
844         else:
845             lineplot.styles+=[None,'scatter']
846         if lineplot.colors==[]:
847             lineplot.colors=[None,None,'black',None]
848         else:
849             lineplot.colors+=['black',None]
850         
851         
852         self._send_plot([lineplot])
853
854         return parameters[0]
855
856
857     def linefit_between(self,index1,index2,whatset=1):
858         '''
859         Creates two vectors (xtofit,ytofit) slicing out from the
860         current return trace a portion delimited by the two indexes
861         given as arguments.
862         Then does a least squares linear fit on that slice.
863         Finally returns [0]=the slope, [1]=the intercept of the
864         fitted 1st grade polynomial, and [2,3]=the actual (x,y) vectors
865         used for the fit.
866         (c) Marco Brucale, Massimo Sandal 2008
867         '''
868         # Translates the indexes into two vectors containing the x,y data to fit
869         xtofit=self.plots[0].vectors[whatset][0][index1:index2]
870         ytofit=self.plots[0].vectors[whatset][1][index1:index2]
871
872         # Does the actual linear fitting (simple least squares with numpy.polyfit)
873         linefit=[]
874         linefit=np.polyfit(xtofit,ytofit,1)
875
876         return (linefit[0],linefit[1],xtofit,ytofit)