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