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