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