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