hooke/plugin/curve.py

   1 # Copyright (C) 2010-2012 W. Trevor King <wking@drexel.edu>
   2 #
   3 # This file is part of Hooke.
   4 #
   5 # Hooke is free software: you can redistribute it and/or modify it
   6 # under the terms of the GNU Lesser General Public License as
   7 # published by the Free Software Foundation, either version 3 of the
   8 # License, or (at your option) any later version.
   9 #
  10 # Hooke is distributed in the hope that it will be useful, but WITHOUT
  11 # ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  12 # or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General
  13 # Public License for more details.
  14 #
  15 # You should have received a copy of the GNU Lesser General Public
  16 # License along with Hooke.  If not, see
  17 # <http://www.gnu.org/licenses/>.
  18
  19 """The ``curve`` module provides :class:`CurvePlugin` and several
  20 associated :class:`hooke.command.Command`\s for handling
  21 :mod:`hooke.curve` classes.
  22 """
  23
  24 import copy
  25 import os.path
  26 import re
  27
  28 import numpy
  29 import yaml
  30
  31 from ..command import Command, Argument, Failure
  32 from ..command_stack import CommandStack
  33 from ..curve import Data
  34 from ..engine import CommandMessage
  35 from ..util.calculus import derivative
  36 from ..util.fft import unitary_avg_power_spectrum
  37 from ..util.si import ppSI, join_data_label, split_data_label
  38 from . import Builtin
  39 from .playlist import current_playlist_callback
  40
  41
  42 # Define common or complicated arguments
  43
  44 def current_curve_callback(hooke, command, argument, value, load=True):
  45     if value != None:
  46         return value
  47     playlist = current_playlist_callback(hooke, command, argument, value)
  48     curve = playlist.current(load=load)
  49     if curve == None:
  50         raise Failure('No curves in %s' % playlist)
  51     return curve
  52
  53 def unloaded_current_curve_callback(hooke, command, argument, value):
  54     return current_curve_callback(
  55         hooke=hooke, command=command, argument=argument, value=value,
  56         load=False)
  57
  58 CurveArgument = Argument(
  59     name='curve', type='curve', callback=current_curve_callback,
  60     help="""
  61 :class:`hooke.curve.Curve` to act on.  Defaults to the current curve
  62 of the current playlist.
  63 """.strip())
  64
  65 def _name_argument(name, default, help):
  66     """TODO
  67     """
  68     return Argument(name=name, type='string', default=default, help=help)
  69
  70 def block_argument(*args, **kwargs):
  71     """TODO
  72     """
  73     return _name_argument(*args, **kwargs)
  74
  75 def column_argument(*args, **kwargs):
  76     """TODO
  77     """
  78     return _name_argument(*args, **kwargs)
  79
  80
  81 # Define useful command subclasses
  82
  83 class CurveCommand (Command):
  84     """A :class:`~hooke.command.Command` operating on a
  85     :class:`~hooke.curve.Curve`.
  86     """
  87     def __init__(self, **kwargs):
  88         if 'arguments' in kwargs:
  89             kwargs['arguments'].insert(0, CurveArgument)
  90         else:
  91             kwargs['arguments'] = [CurveArgument]
  92         super(CurveCommand, self).__init__(**kwargs)
  93
  94     def _curve(self, hooke, params):
  95         """Get the selected curve.
  96
  97         Notes
  98         -----
  99         `hooke` is intended to attach the selected curve to the local
 100         playlist; the returned curve should not be effected by the
 101         state of `hooke`.  This is important for reliable
 102         :class:`~hooke.command_stack.CommandStack`\s.
 103         """
 104         # HACK? rely on params['curve'] being bound to the local hooke
 105         # playlist (i.e. not a copy, as you would get by passing a
 106         # curve through the queue).  Ugh.  Stupid queues.  As an
 107         # alternative, we could pass lookup information through the
 108         # queue...
 109         return params['curve']
 110
 111     def _add_to_command_stack(self, params):
 112         """Store the command name and current `params` values in the
 113         curve's `.command_stack`.
 114
 115         If this would duplicate the command currently on top of the
 116         stack, no action is taken.  Call early on, or watch out for
 117         repeated param processing.
 118
 119         Recommended practice is to *not* lock in argument values that
 120         are loaded from the plugin's :attr:`.config`.
 121
 122         Notes
 123         -----
 124         Perhaps we should subclass :meth:`_run` and use :func:`super`,
 125         or embed this in :meth:`run` to avoid subclasses calling this
 126         method explicitly, with all the tedium and brittality that
 127         implies.  On the other hand, the current implemtnation allows
 128         CurveCommands that don't effect the curve itself
 129         (e.g. :class:`GetCommand`) to avoid adding themselves to the
 130         stack entirely.
 131         """
 132         if params['stack'] == True:
 133             curve = self._curve(hooke=None, params=params)
 134             if (len(curve.command_stack) > 0
 135                 and curve.command_stack[-1].command == self.name
 136                 and curve.command_stack[-1].arguments == params):
 137                 pass  # no need to place duplicate calls on the stack.
 138             else:
 139                 curve.command_stack.append(CommandMessage(
 140                         self.name, dict(params)))
 141
 142
 143 class BlockCommand (CurveCommand):
 144     """A :class:`CurveCommand` operating on a :class:`~hooke.curve.Data` block.
 145     """
 146     def __init__(self, blocks=None, **kwargs):
 147         if blocks == None:
 148             blocks = [('block', None, 'Name of the data block to act on.')]
 149         block_args = []
 150         for name,default,help in blocks:
 151             block_args.append(block_argument(name, default, help))
 152         self._block_arguments = block_args
 153         if 'arguments' not in kwargs:
 154             kwargs['arguments'] = []
 155         kwargs['arguments'] = block_args + kwargs['arguments']
 156         super(BlockCommand, self).__init__(**kwargs)
 157
 158     def _block_names(self, hooke, params):
 159         curve = self._curve(hooke, params)
 160         return [b.info['name'] for b in curve.data]
 161
 162     def _block_index(self, hooke, params, name=None):
 163         if name == None:
 164             name = self._block_arguments[0].name
 165         block_name = params[name]
 166         if block_name == None:
 167             curve = self._curve(hooke=hooke, params=params)
 168             if len(curve.data) == 0:
 169                 raise Failure('no blocks in %s' % curve)
 170             block_name = curve.data[0].info['name']
 171         names = self._block_names(hooke=hooke, params=params)
 172         try:
 173             return names.index(block_name)
 174         except ValueError, e:
 175             curve = self._curve(hooke, params)
 176             raise Failure('no block named %s in %s (%s): %s'
 177                           % (block_name, curve, names, e))
 178
 179     def _block(self, hooke, params, name=None):
 180         # HACK? rely on params['block'] being bound to the local hooke
 181         # playlist (i.e. not a copy, as you would get by passing a
 182         # block through the queue).  Ugh.  Stupid queues.  As an
 183         # alternative, we could pass lookup information through the
 184         # queue...
 185         curve = self._curve(hooke, params)
 186         index = self._block_index(hooke, params, name)
 187         return curve.data[index]
 188
 189
 190 class ColumnAccessCommand (BlockCommand):
 191     """A :class:`BlockCommand` accessing a :class:`~hooke.curve.Data`
 192     block column.
 193     """
 194     def __init__(self, columns=None, **kwargs):
 195         if columns == None:
 196             columns = [('column', None, 'Name of the data column to act on.')]
 197         column_args = []
 198         for name,default,help in columns:
 199             column_args.append(column_argument(name, default, help))
 200         self._column_arguments = column_args
 201         if 'arguments' not in kwargs:
 202             kwargs['arguments'] = []
 203         kwargs['arguments'] = column_args + kwargs['arguments']
 204         super(ColumnAccessCommand, self).__init__(**kwargs)
 205
 206     def _get_column(self, hooke, params, block_name=None, column_name=None):
 207         if column_name == None:
 208             column_name = self._column_arguments[0].name
 209         column_name = params[column_name]
 210         if column_name is None:
 211             return None
 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             ScaledColumnAdditionCommand(self), 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(os.path.expanduser(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 ScaledColumnAdditionCommand (ColumnAddingCommand):
 689     """Add one affine transformed column to another: `o=A*i1+B*i2+C`.
 690     """
 691     def __init__(self, plugin):
 692         super(ScaledColumnAdditionCommand, self).__init__(
 693             name='scaled column addition',
 694             columns=[
 695                 ('input column 1', 'input column (m)', """
 696 Name of the first column to use as the transform input.
 697 """.strip()),
 698                 ('input column 2', None, """
 699 Name of the second column to use as the transform input.
 700 """.strip()),
 701                 ],
 702             new_columns=[
 703                 ('output column', 'output column (m)', """
 704 Name of the column to use as the transform output.
 705 """.strip()),
 706                 ],
 707             arguments=[
 708                 Argument(name='scale 1', type='float', default=None,
 709                          help="""
 710 A float value for the first scale constant.
 711 """.strip()),
 712                 Argument(name='scale 1 name', type='string', default=None,
 713                          help="""
 714 The name of the first scale constant in the `.info` dictionary.
 715 """.strip()),
 716                 Argument(name='scale 2', type='float', default=None,
 717                          help="""
 718 A float value for the second scale constant.
 719 """.strip()),
 720                 Argument(name='scale 2 name', type='string', default=None,
 721                          help="""
 722 The name of the second scale constant in the `.info` dictionary.
 723 """.strip()),
 724                 Argument(name='constant', type='float', default=None,
 725                          help="""
 726 A float value for the offset constant.
 727 """.strip()),
 728                 Argument(name='constant name', type='string', default=None,
 729                          help="""
 730 The name of the offset constant in the `.info` dictionary.
 731 """.strip()),
 732                 ],
 733             help=self.__doc__, plugin=plugin)
 734
 735     def _run(self, hooke, inqueue, outqueue, params):
 736         self._add_to_command_stack(params)
 737         i1 = self._get_column(hooke=hooke, params=params,
 738                                     column_name='input column 1')
 739         i2 = self._get_column(hooke=hooke, params=params,
 740                                     column_name='input column 2')
 741         if i1 is None:
 742             i1 = 0
 743         if i2 is None:
 744             i2 = 0
 745         # what if both i1 and i2 are None?
 746         a = self._get_constant(params, i1.info, 'scale 1')
 747         b = self._get_constant(params, i2.info, 'scale 2')
 748         c = self._get_constant(params, i1.info, 'constant')
 749         out = a*i1 + b*i2 + c
 750         self._set_column(hooke=hooke, params=params,
 751                          column_name='output column', values=out)
 752
 753     def _get_constant(self, params, info, name):
 754         a = params[name]
 755         pname = params[name + ' name']
 756         b = None
 757         if pname is not None:
 758             pname_entries = pname.split('|')
 759             b = info
 760             for entry in pname_entries:
 761                 b = b[entry]
 762         if a is None and b is None:
 763             return 0
 764         if a is None:
 765             a = 1
 766         if b is None:
 767             b = 1
 768         return a*b
 769
 770
 771 class ClearStackCommand (CurveCommand):
 772     """Empty a curve's command stack.
 773     """
 774     def __init__(self, plugin):
 775         super(ClearStackCommand, self).__init__(
 776             name='clear curve command stack',
 777             help=self.__doc__, plugin=plugin)
 778         i,arg = [(i,arg) for i,arg in enumerate(self.arguments)
 779                  if arg.name == 'curve'][0]
 780         arg = copy.copy(arg)
 781         arg.callback = unloaded_current_curve_callback
 782         self.arguments[i] = arg
 783
 784     def _run(self, hooke, inqueue, outqueue, params):
 785         curve = self._curve(hooke, params)
 786         curve.command_stack = CommandStack()
 787
 788
 789 class OldCruft (object):
 790
 791     def do_forcebase(self,args):
 792         '''
 793         FORCEBASE
 794         (generalvclamp.py)
 795         Measures the difference in force (in pN) between a point and a baseline
 796         took as the average between two points.
 797
 798         The baseline is fixed once for a given curve and different force measurements,
 799         unless the user wants it to be recalculated
 800         ------------
 801         Syntax: forcebase [rebase]
 802                 rebase: Forces forcebase to ask again the baseline
 803                 max: Instead of asking for a point to measure, asks for two points and use
 804                      the maximum peak in between
 805         '''
 806         rebase=False #if true=we select rebase
 807         maxpoint=False #if true=we measure the maximum peak
 808
 809         plot=self._get_displayed_plot()
 810         whatset=1 #fixme: for all sets
 811         if 'rebase' in args or (self.basecurrent != self.current.path):
 812             rebase=True
 813         if 'max' in args:
 814             maxpoint=True
 815
 816         if rebase:
 817             print 'Select baseline'
 818             self.basepoints=self._measure_N_points(N=2, whatset=whatset)
 819             self.basecurrent=self.current.path
 820
 821         if maxpoint:
 822             print 'Select two points'
 823             points=self._measure_N_points(N=2, whatset=whatset)
 824             boundpoints=[points[0].index, points[1].index]
 825             boundpoints.sort()
 826             try:
 827                 y=min(plot.vectors[whatset][1][boundpoints[0]:boundpoints[1]])
 828             except ValueError:
 829                 print 'Chosen interval not valid. Try picking it again. Did you pick the same point as begin and end of interval?'
 830         else:
 831             print 'Select point to measure'
 832             points=self._measure_N_points(N=1, whatset=whatset)
 833             #whatplot=points[0].dest
 834             y=points[0].graph_coords[1]
 835
 836         #fixme: code duplication
 837         boundaries=[self.basepoints[0].index, self.basepoints[1].index]
 838         boundaries.sort()
 839         to_average=plot.vectors[whatset][1][boundaries[0]:boundaries[1]] #y points to average
 840
 841         avg=np.mean(to_average)
 842         forcebase=abs(y-avg)
 843         print str(forcebase*(10**12))+' pN'
 844         to_dump='forcebase '+self.current.path+' '+str(forcebase*(10**12))+' pN'
 845         self.outlet.push(to_dump)
 846
 847     #---SLOPE---
 848     def do_slope(self,args):
 849         '''
 850         SLOPE
 851         (generalvclamp.py)
 852         Measures the slope of a delimited chunk on the return trace.
 853         The chunk can be delimited either by two manual clicks, or have
 854         a fixed width, given as an argument.
 855         ---------------
 856         Syntax: slope [width]
 857                 The facultative [width] parameter specifies how many
 858                 points will be considered for the fit. If [width] is
 859                 specified, only one click will be required.
 860         (c) Marco Brucale, Massimo Sandal 2008
 861         '''
 862
 863         # Reads the facultative width argument
 864         try:
 865             fitspan=int(args)
 866         except:
 867             fitspan=0
 868
 869         # Decides between the two forms of user input, as per (args)
 870         if fitspan == 0:
 871             # Gets the Xs of two clicked points as indexes on the current curve vector
 872             print 'Click twice to delimit chunk'
 873             points=self._measure_N_points(N=2,whatset=1)
 874         else:
 875             print 'Click once on the leftmost point of the chunk (i.e.usually the peak)'
 876             points=self._measure_N_points(N=1,whatset=1)
 877
 878         slope=self._slope(points,fitspan)
 879
 880         # Outputs the relevant slope parameter
 881         print 'Slope:'
 882         print str(slope)
 883         to_dump='slope '+self.current.path+' '+str(slope)
 884         self.outlet.push(to_dump)
 885
 886     def _slope(self,points,fitspan):
 887         # Calls the function linefit_between
 888         parameters=[0,0,[],[]]
 889         try:
 890             clickedpoints=[points[0].index,points[1].index]
 891             clickedpoints.sort()
 892         except:
 893             clickedpoints=[points[0].index-fitspan,points[0].index]
 894
 895         try:
 896             parameters=self.linefit_between(clickedpoints[0],clickedpoints[1])
 897         except:
 898             print 'Cannot fit. Did you click twice the same point?'
 899             return
 900
 901         # Outputs the relevant slope parameter
 902         print 'Slope:'
 903         print str(parameters[0])
 904         to_dump='slope '+self.curve.path+' '+str(parameters[0])
 905         self.outlet.push(to_dump)
 906
 907         # Makes a vector with the fitted parameters and sends it to the GUI
 908         xtoplot=parameters[2]
 909         ytoplot=[]
 910         x=0
 911         for x in xtoplot:
 912             ytoplot.append((x*parameters[0])+parameters[1])
 913
 914         clickvector_x, clickvector_y=[], []
 915         for item in points:
 916             clickvector_x.append(item.graph_coords[0])
 917             clickvector_y.append(item.graph_coords[1])
 918
 919         lineplot=self._get_displayed_plot(0) #get topmost displayed plot
 920
 921         lineplot.add_set(xtoplot,ytoplot)
 922         lineplot.add_set(clickvector_x, clickvector_y)
 923
 924
 925         if lineplot.styles==[]:
 926             lineplot.styles=[None,None,None,'scatter']
 927         else:
 928             lineplot.styles+=[None,'scatter']
 929         if lineplot.colors==[]:
 930             lineplot.colors=[None,None,'black',None]
 931         else:
 932             lineplot.colors+=['black',None]
 933
 934
 935         self._send_plot([lineplot])
 936
 937         return parameters[0]
 938
 939
 940     def linefit_between(self,index1,index2,whatset=1):
 941         '''
 942         Creates two vectors (xtofit,ytofit) slicing out from the
 943         current return trace a portion delimited by the two indexes
 944         given as arguments.
 945         Then does a least squares linear fit on that slice.
 946         Finally returns [0]=the slope, [1]=the intercept of the
 947         fitted 1st grade polynomial, and [2,3]=the actual (x,y) vectors
 948         used for the fit.
 949         (c) Marco Brucale, Massimo Sandal 2008
 950         '''
 951         # Translates the indexes into two vectors containing the x,y data to fit
 952         xtofit=self.plots[0].vectors[whatset][0][index1:index2]
 953         ytofit=self.plots[0].vectors[whatset][1][index1:index2]
 954
 955         # Does the actual linear fitting (simple least squares with numpy.polyfit)
 956         linefit=[]
 957         linefit=np.polyfit(xtofit,ytofit,1)
 958
 959         return (linefit[0],linefit[1],xtofit,ytofit)