hooke/plugin/curve.py

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