hooke/plugin/curve.py

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