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', 'experiment', 'driver', 'filetype',
 296                        'note', 'command stack', 'blocks', '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_experiment(self, curve):
 329         return curve.info.get('experiment', None)
 330
 331     def _get_driver(self, curve):
 332         return curve.driver
 333
 334     def _get_filetype(self, curve):
 335         return curve.info.get('filetype', None)
 336
 337     def _get_note(self, curve):
 338         return curve.info.get('note', None)
 339
 340     def _get_command_stack(self, curve):
 341         return curve.command_stack
 342
 343     def _get_blocks(self, curve):
 344         return len(curve.data)
 345
 346     def _get_block_sizes(self, curve):
 347         return [block.shape for block in curve.data]
 348
 349
 350 class BlockInfoCommand (BlockCommand):
 351     """Get selected information about a :class:`hooke.curve.Curve` data block.
 352     """
 353     def __init__(self, plugin):
 354         super(BlockInfoCommand, self).__init__(
 355             name='block info', arguments=[
 356                 Argument(
 357                     name='key', count=-1, optional=False,
 358                     help='Dot-separted (.) key selection regexp.'),
 359                 Argument(
 360                     name='output',
 361                     help="""
 362 File name for the output (appended).
 363 """.strip()),
 364                 ],
 365             help=self.__doc__, plugin=plugin)
 366
 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)]
 372             while len(keys) > 0:
 373                 index,key_stack,info = keys.pop(0)
 374                 regexp = re.compile(key_stack[index])
 375                 matched = False
 376                 for k,v in info.items():
 377                     if regexp.match(k):
 378                         matched = True
 379                         new_stack = copy.copy(key_stack)
 380                         new_stack[index] = k
 381                         if index+1 == len(key_stack):
 382                             vals = values
 383                             for k in new_stack[:-1]:
 384                                 if k not in vals:
 385                                     vals[k] = {}
 386                                 vals = vals[k]
 387                             vals[new_stack[-1]] = v
 388                         else:
 389                             keys.append((index+1, new_stack, v))
 390                 if matched == False:
 391                     raise ValueError('no match found for %s in %s'
 392                                      % (key_stack[index], key))
 393         if params['output'] != None:
 394             curve = self._curve(hooke, params)
 395             with open(params['output'], 'a') as f:
 396                 yaml.dump({curve.name:{
 397                             'path': curve.path,
 398                             block.info['name']: values
 399                             }}, f)
 400         outqueue.put(values)
 401
 402
 403 class DeltaCommand (BlockCommand):
 404     """Get distance information between two points.
 405
 406     With two points A and B, the returned distances are A-B.
 407     """
 408     def __init__(self, plugin):
 409         super(DeltaCommand, self).__init__(
 410             name='delta',
 411             arguments=[
 412                 Argument(name='point', type='point', optional=False, count=2,
 413                          help="""
 414 Indicies of points bounding the selected data.
 415 """.strip()),
 416                 Argument(name='SI', type='bool', default=False,
 417                          help="""
 418 Return distances in SI notation.
 419 """.strip())
 420                 ],
 421             help=self.__doc__, plugin=plugin)
 422
 423     def _run(self, hooke, inqueue, outqueue, params):
 424         data = self._block(hooke, params)
 425         As = data[params['point'][0],:]
 426         Bs = data[params['point'][1],:]
 427         ds = [A-B for A,B in zip(As, Bs)]
 428         if params['SI'] == False:
 429             out = [(name, d) for name,d in zip(data.info['columns'], ds)]
 430         else:
 431             out = []
 432             for name,d in zip(data.info['columns'], ds):
 433                 n,units = split_data_label(name)
 434                 out.append(
 435                   (n, ppSI(value=d, unit=units, decimals=2)))
 436         outqueue.put(out)
 437
 438
 439 class ExportCommand (BlockCommand):
 440     """Export a :class:`hooke.curve.Curve` data block as TAB-delimeted
 441     ASCII text.
 442
 443     A "#" prefixed header will optionally appear at the beginning of
 444     the file naming the columns.
 445     """
 446     def __init__(self, plugin):
 447         super(ExportCommand, self).__init__(
 448             name='export block',
 449             arguments=[
 450                 Argument(name='output', type='file', default='curve.dat',
 451                          help="""
 452 File name for the output data.  Defaults to 'curve.dat'
 453 """.strip()),
 454                 Argument(name='header', type='bool', default=True,
 455                          help="""
 456 True if you want the column-naming header line.
 457 """.strip()),
 458                 ],
 459             help=self.__doc__, plugin=plugin)
 460
 461     def _run(self, hooke, inqueue, outqueue, params):
 462         data = self._block(hooke, params)
 463
 464         with open(params['output'], 'w') as f:
 465             if params['header'] == True:
 466                 f.write('# %s \n' % ('\t'.join(data.info['columns'])))
 467             numpy.savetxt(f, data, delimiter='\t')
 468
 469
 470 class DifferenceCommand (ColumnAddingCommand):
 471     """Calculate the difference between two columns of data.
 472
 473     The difference is added to block A as a new column.
 474
 475     Note that the command will fail if the columns have different
 476     lengths, so be careful when differencing columns from different
 477     blocks.
 478     """
 479     def __init__(self, plugin):
 480         super(DifferenceCommand, self).__init__(
 481             name='difference',
 482             blocks=[
 483                 ('block A', None,
 484                  'Name of block A in A-B.  Defaults to the first block'),
 485                 ('block B', None,
 486                  'Name of block B in A-B.  Defaults to matching `block A`.'),
 487                 ],
 488             columns=[
 489                 ('column A', None,
 490                  """
 491 Column of data from block A to difference.  Defaults to the first column.
 492 """.strip()),
 493                 ('column B', None,
 494                  """
 495 Column of data from block B to difference.  Defaults to matching `column A`.
 496 """.strip()),
 497                 ],
 498             new_columns=[
 499                 ('output column', None,
 500                  """
 501 Name of the new column for storing the difference (without units, defaults to
 502 `difference of <block A> <column A> and <block B> <column B>`).
 503 """.strip()),
 504                 ],
 505             help=self.__doc__, plugin=plugin)
 506
 507     def _run(self, hooke, inqueue, outqueue, params):
 508         self._add_to_command_stack(params)
 509         params = self._setup_params(hooke=hooke, params=params)
 510         data_A = self._get_column(hooke=hooke, params=params,
 511                                   block_name='block A',
 512                                   column_name='column A')
 513         data_B = self._get_column(hooke=hooke, params=params,
 514                                   block_name='block B',
 515                                   column_name='column B')
 516         out = data_A - data_B
 517         self._set_column(hooke=hooke, params=params,
 518                          block_name='block A',
 519                          column_name='output column',
 520                          values=out)
 521
 522     def _setup_params(self, hooke, params):
 523         curve = self._curve(hooke, params)
 524         if params['block A'] == None:
 525             params['block A'] = curve.data[0].info['name']
 526         if params['block B'] == None:
 527             params['block B'] = params['block A']
 528         block_A = self._block(hooke, params=params, name='block A')
 529         block_B = self._block(hooke, params=params, name='block B')
 530         if params['column A'] == None:
 531             params['column A'] = block.info['columns'][0]
 532         if params['column B'] == None:
 533             params['column B'] = params['column A']
 534         a_name,a_unit = split_data_label(params['column A'])
 535         b_name,b_unit = split_data_label(params['column B'])
 536         if a_unit != b_unit:
 537             raise Failure('Unit missmatch: %s != %s' % (a_unit, b_unit))
 538         if params['output column'] == None:
 539             params['output column'] = join_data_label(
 540                 'difference of %s %s and %s %s' % (
 541                     block_A.info['name'], a_name,
 542                     block_B.info['name'], b_name),
 543                 a_unit)
 544         else:
 545             params['output column'] = join_data_label(
 546                 params['output column'], a_unit)
 547         return params
 548
 549
 550 class DerivativeCommand (ColumnAddingCommand):
 551     """Calculate the derivative (actually, the discrete differentiation)
 552     of a data column.
 553
 554     See :func:`hooke.util.calculus.derivative` for implementation
 555     details.
 556     """
 557     def __init__(self, plugin):
 558         super(DerivativeCommand, self).__init__(
 559             name='derivative',
 560             columns=[
 561                 ('x column', None,
 562                  'Column of data block to differentiate with respect to.'),
 563                 ('f column', None,
 564                  'Column of data block to differentiate.'),
 565                 ],
 566             new_columns=[
 567                 ('output column', None,
 568                  """
 569 Name of the new column for storing the derivative (without units, defaults to
 570 `derivative of <f column> with respect to <x column>`).
 571 """.strip()),
 572                 ],
 573             arguments=[
 574                 Argument(name='weights', type='dict', default={-1:-0.5, 1:0.5},
 575                          help="""
 576 Weighting scheme dictionary for finite differencing.  Defaults to
 577 central differencing.
 578 """.strip()),
 579                 ],
 580             help=self.__doc__, plugin=plugin)
 581
 582     def _run(self, hooke, inqueue, outqueue, params):
 583         self._add_to_command_stack(params)
 584         params = self._setup_params(hooke=hooke, params=params)
 585         x_data = self._get_column(hooke=hooke, params=params,
 586                                   column_name='x column')
 587         f_data = self._get_column(hooke=hooke, params=params,
 588                                   column_name='f column')
 589         d = derivative(
 590             x_data=x_data, f_data=f_data, weights=params['weights'])
 591         self._set_column(hooke=hooke, params=params,
 592                          column_name='output column',
 593                          values=d)
 594
 595     def _setup_params(self, hooke, params):
 596         curve = self._curve(hooke, params)
 597         x_name,x_unit = split_data_label(params['x column'])
 598         f_name,f_unit = split_data_label(params['f column'])
 599         d_unit = '%s/%s' % (f_unit, x_unit)
 600         if params['output column'] == None:
 601             params['output column'] = join_data_label(
 602                 'derivative of %s with respect to %s' % (
 603                     f_name, x_name),
 604                 d_unit)
 605         else:
 606             params['output column'] = join_data_label(
 607                 params['output column'], d_unit)
 608         return params
 609
 610
 611 class PowerSpectrumCommand (ColumnAddingCommand):
 612     """Calculate the power spectrum of a data column.
 613     """
 614     def __init__(self, plugin):
 615         super(PowerSpectrumCommand, self).__init__(
 616             name='power spectrum',
 617             arguments=[
 618                 Argument(name='output block', type='string',
 619                          help="""
 620 Name of the new data block for storing the power spectrum (defaults to
 621 `power spectrum of <source block name> <source column name>`).
 622 """.strip()),
 623                 Argument(name='bounds', type='point', optional=True, count=2,
 624                          help="""
 625 Indicies of points bounding the selected data.
 626 """.strip()),
 627                 Argument(name='freq', type='float', default=1.0,
 628                          help="""
 629 Sampling frequency.
 630 """.strip()),
 631                 Argument(name='freq units', type='string', default='Hz',
 632                          help="""
 633 Units for the sampling frequency.
 634 """.strip()),
 635                 Argument(name='chunk size', type='int', default=2048,
 636                          help="""
 637 Number of samples per chunk.  Use a power of two.
 638 """.strip()),
 639                 Argument(name='overlap', type='bool', default=False,
 640                          help="""
 641 If `True`, each chunk overlaps the previous chunk by half its length.
 642 Otherwise, the chunks are end-to-end, and not overlapping.
 643 """.strip()),
 644                 ],
 645             help=self.__doc__, plugin=plugin)
 646
 647     def _run(self, hooke, inqueue, outqueue, params):
 648         self._add_to_command_stack(params)
 649         params = self._setup_params(hooke=hooke, params=params)
 650         data = self._get_column(hooke=hooke, params=params)
 651         bounds = params['bounds']
 652         if bounds != None:
 653             data = data[bounds[0]:bounds[1]]
 654         freq_axis,power = unitary_avg_power_spectrum(
 655             data, freq=params['freq'],
 656             chunk_size=params['chunk size'],
 657             overlap=params['overlap'])
 658         b = Data((len(freq_axis),2), dtype=data.dtype)
 659         b.info['name'] = params['output block']
 660         b.info['columns'] = [
 661             params['output freq column'],
 662             params['output power column'],
 663             ]
 664         self._curve(hooke, params).data.append(b)
 665         self._set_column(hooke, params, block_name='output block',
 666                          column_name='output freq column',
 667                          values=freq_axis)
 668         self._set_column(hooke, params, block_name='output block',
 669                          column_name='output power column',
 670                          values=power)
 671         outqueue.put(b)
 672
 673     def _setup_params(self, hooke, params):
 674         if params['output block'] in self._block_names(hooke, params):
 675             raise Failure('output block %s already exists in %s.'
 676                           % (params['output block'],
 677                              self._curve(hooke, params)))
 678         data = self._get_column(hooke=hooke, params=params)
 679         d_name,d_unit = split_data_label(data.info['name'])
 680         if params['output block'] == None:
 681             params['output block'] = 'power spectrum of %s %s' % (
 682                 data.info['name'], params['column'])
 683         self.params['output freq column'] = join_data_label(
 684             'frequency axis', params['freq units'])
 685         self.params['output power column'] = join_data_label(
 686             'power density', '%s^2/%s' % (data_units, params['freq units']))
 687         return params
 688
 689
 690 class ClearStackCommand (CurveCommand):
 691     """Empty a curve's command stack.
 692     """
 693     def __init__(self, plugin):
 694         super(ClearStackCommand, self).__init__(
 695             name='clear curve command stack',
 696             help=self.__doc__, plugin=plugin)
 697         i,arg = [(i,arg) for i,arg in enumerate(self.arguments)
 698                  if arg.name == 'curve'][0]
 699         arg = copy.copy(arg)
 700         arg.callback = unloaded_current_curve_callback
 701         self.arguments[i] = arg
 702
 703     def _run(self, hooke, inqueue, outqueue, params):
 704         curve = self._curve(hooke, params)
 705         curve.command_stack = CommandStack()
 706
 707
 708 class OldCruft (object):
 709
 710     def do_forcebase(self,args):
 711         '''
 712         FORCEBASE
 713         (generalvclamp.py)
 714         Measures the difference in force (in pN) between a point and a baseline
 715         took as the average between two points.
 716
 717         The baseline is fixed once for a given curve and different force measurements,
 718         unless the user wants it to be recalculated
 719         ------------
 720         Syntax: forcebase [rebase]
 721                 rebase: Forces forcebase to ask again the baseline
 722                 max: Instead of asking for a point to measure, asks for two points and use
 723                      the maximum peak in between
 724         '''
 725         rebase=False #if true=we select rebase
 726         maxpoint=False #if true=we measure the maximum peak
 727
 728         plot=self._get_displayed_plot()
 729         whatset=1 #fixme: for all sets
 730         if 'rebase' in args or (self.basecurrent != self.current.path):
 731             rebase=True
 732         if 'max' in args:
 733             maxpoint=True
 734
 735         if rebase:
 736             print 'Select baseline'
 737             self.basepoints=self._measure_N_points(N=2, whatset=whatset)
 738             self.basecurrent=self.current.path
 739
 740         if maxpoint:
 741             print 'Select two points'
 742             points=self._measure_N_points(N=2, whatset=whatset)
 743             boundpoints=[points[0].index, points[1].index]
 744             boundpoints.sort()
 745             try:
 746                 y=min(plot.vectors[whatset][1][boundpoints[0]:boundpoints[1]])
 747             except ValueError:
 748                 print 'Chosen interval not valid. Try picking it again. Did you pick the same point as begin and end of interval?'
 749         else:
 750             print 'Select point to measure'
 751             points=self._measure_N_points(N=1, whatset=whatset)
 752             #whatplot=points[0].dest
 753             y=points[0].graph_coords[1]
 754
 755         #fixme: code duplication
 756         boundaries=[self.basepoints[0].index, self.basepoints[1].index]
 757         boundaries.sort()
 758         to_average=plot.vectors[whatset][1][boundaries[0]:boundaries[1]] #y points to average
 759
 760         avg=np.mean(to_average)
 761         forcebase=abs(y-avg)
 762         print str(forcebase*(10**12))+' pN'
 763         to_dump='forcebase '+self.current.path+' '+str(forcebase*(10**12))+' pN'
 764         self.outlet.push(to_dump)
 765
 766     #---SLOPE---
 767     def do_slope(self,args):
 768         '''
 769         SLOPE
 770         (generalvclamp.py)
 771         Measures the slope of a delimited chunk on the return trace.
 772         The chunk can be delimited either by two manual clicks, or have
 773         a fixed width, given as an argument.
 774         ---------------
 775         Syntax: slope [width]
 776                 The facultative [width] parameter specifies how many
 777                 points will be considered for the fit. If [width] is
 778                 specified, only one click will be required.
 779         (c) Marco Brucale, Massimo Sandal 2008
 780         '''
 781
 782         # Reads the facultative width argument
 783         try:
 784             fitspan=int(args)
 785         except:
 786             fitspan=0
 787
 788         # Decides between the two forms of user input, as per (args)
 789         if fitspan == 0:
 790             # Gets the Xs of two clicked points as indexes on the current curve vector
 791             print 'Click twice to delimit chunk'
 792             points=self._measure_N_points(N=2,whatset=1)
 793         else:
 794             print 'Click once on the leftmost point of the chunk (i.e.usually the peak)'
 795             points=self._measure_N_points(N=1,whatset=1)
 796
 797         slope=self._slope(points,fitspan)
 798
 799         # Outputs the relevant slope parameter
 800         print 'Slope:'
 801         print str(slope)
 802         to_dump='slope '+self.current.path+' '+str(slope)
 803         self.outlet.push(to_dump)
 804
 805     def _slope(self,points,fitspan):
 806         # Calls the function linefit_between
 807         parameters=[0,0,[],[]]
 808         try:
 809             clickedpoints=[points[0].index,points[1].index]
 810             clickedpoints.sort()
 811         except:
 812             clickedpoints=[points[0].index-fitspan,points[0].index]
 813
 814         try:
 815             parameters=self.linefit_between(clickedpoints[0],clickedpoints[1])
 816         except:
 817             print 'Cannot fit. Did you click twice the same point?'
 818             return
 819
 820         # Outputs the relevant slope parameter
 821         print 'Slope:'
 822         print str(parameters[0])
 823         to_dump='slope '+self.curve.path+' '+str(parameters[0])
 824         self.outlet.push(to_dump)
 825
 826         # Makes a vector with the fitted parameters and sends it to the GUI
 827         xtoplot=parameters[2]
 828         ytoplot=[]
 829         x=0
 830         for x in xtoplot:
 831             ytoplot.append((x*parameters[0])+parameters[1])
 832
 833         clickvector_x, clickvector_y=[], []
 834         for item in points:
 835             clickvector_x.append(item.graph_coords[0])
 836             clickvector_y.append(item.graph_coords[1])
 837
 838         lineplot=self._get_displayed_plot(0) #get topmost displayed plot
 839
 840         lineplot.add_set(xtoplot,ytoplot)
 841         lineplot.add_set(clickvector_x, clickvector_y)
 842
 843
 844         if lineplot.styles==[]:
 845             lineplot.styles=[None,None,None,'scatter']
 846         else:
 847             lineplot.styles+=[None,'scatter']
 848         if lineplot.colors==[]:
 849             lineplot.colors=[None,None,'black',None]
 850         else:
 851             lineplot.colors+=['black',None]
 852
 853
 854         self._send_plot([lineplot])
 855
 856         return parameters[0]
 857
 858
 859     def linefit_between(self,index1,index2,whatset=1):
 860         '''
 861         Creates two vectors (xtofit,ytofit) slicing out from the
 862         current return trace a portion delimited by the two indexes
 863         given as arguments.
 864         Then does a least squares linear fit on that slice.
 865         Finally returns [0]=the slope, [1]=the intercept of the
 866         fitted 1st grade polynomial, and [2,3]=the actual (x,y) vectors
 867         used for the fit.
 868         (c) Marco Brucale, Massimo Sandal 2008
 869         '''
 870         # Translates the indexes into two vectors containing the x,y data to fit
 871         xtofit=self.plots[0].vectors[whatset][0][index1:index2]
 872         ytofit=self.plots[0].vectors[whatset][1][index1:index2]
 873
 874         # Does the actual linear fitting (simple least squares with numpy.polyfit)
 875         linefit=[]
 876         linefit=np.polyfit(xtofit,ytofit,1)
 877
 878         return (linefit[0],linefit[1],xtofit,ytofit)