hooke/plugin/curve.py

   1 # Copyright (C) 2008-2010 Alberto Gomez-Casado
   2 #                         Fabrizio Benedetti
   3 #                         Massimo Sandal <devicerandom@gmail.com>
   4 #                         W. Trevor King <wking@drexel.edu>
   5 #
   6 # This file is part of Hooke.
   7 #
   8 # Hooke is free software: you can redistribute it and/or modify it
   9 # under the terms of the GNU Lesser General Public License as
  10 # published by the Free Software Foundation, either version 3 of the
  11 # License, or (at your option) any later version.
  12 #
  13 # Hooke is distributed in the hope that it will be useful, but WITHOUT
  14 # ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  15 # or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General
  16 # Public License for more details.
  17 #
  18 # You should have received a copy of the GNU Lesser General Public
  19 # License along with Hooke.  If not, see
  20 # <http://www.gnu.org/licenses/>.
  21
  22 """The ``curve`` module provides :class:`CurvePlugin` and several
  23 associated :class:`hooke.command.Command`\s for handling
  24 :mod:`hooke.curve` classes.
  25 """
  26
  27 import copy
  28
  29 import numpy
  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, 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         block = self._block(hooke, params, block_name)
 211         columns = block.info['columns']
 212         try:
 213             column_index = columns.index(column_name)
 214         except ValueError, e:
 215             raise Failure('%s not in %s (%s): %s'
 216                           % (column_name, block.info['name'], columns, e))
 217         return block[:,column_index]
 218
 219
 220 class ColumnAddingCommand (ColumnAccessCommand):
 221     """A :class:`ColumnAccessCommand` that also adds columns.
 222     """
 223     def __init__(self, new_columns=None, **kwargs):
 224         if new_columns == None:
 225             new_columns = []
 226         column_args = []
 227         for name,default,help in new_columns:
 228             column_args.append(column_argument(name, default, help))
 229         self._new_column_arguments = column_args
 230         if 'arguments' not in kwargs:
 231             kwargs['arguments'] = []
 232         kwargs['arguments'] = column_args + kwargs['arguments']
 233         super(ColumnAddingCommand, self).__init__(**kwargs)
 234
 235     def _get_column(self, hooke, params, block_name=None, column_name=None):
 236         if column_name == None and len(self._column_arguments) == 0:
 237             column_name = self._new_column_arguments[0].name
 238         return super(ColumnAddingCommand, self)._get_column(
 239             hooke=hooke, params=params, block_name=block_name,
 240             column_name=column_name)
 241
 242     def _set_column(self, hooke, params, block_name=None, column_name=None,
 243                     values=None):
 244         if column_name == None:
 245             column_name = self._column_arguments[0].name
 246         column_name = params[column_name]
 247         block = self._block(hooke=hooke, params=params, name=block_name)
 248         if column_name not in block.info['columns']:
 249             new = Data((block.shape[0], block.shape[1]+1), dtype=block.dtype)
 250             new.info = copy.deepcopy(block.info)
 251             new[:,:-1] = block
 252             new.info['columns'].append(column_name)
 253             block = new
 254             block_index = self._block_index(hooke, params, name=block_name)
 255             self._curve(hooke, params).data[block_index] = block
 256         column_index = block.info['columns'].index(column_name)
 257         block[:,column_index] = values
 258
 259
 260 # The plugin itself
 261
 262 class CurvePlugin (Builtin):
 263     def __init__(self):
 264         super(CurvePlugin, self).__init__(name='curve')
 265         self._commands = [
 266             GetCommand(self), InfoCommand(self), DeltaCommand(self),
 267             ExportCommand(self), DifferenceCommand(self),
 268             DerivativeCommand(self), PowerSpectrumCommand(self),
 269             ClearStackCommand(self)]
 270
 271
 272 # Define commands
 273
 274 class GetCommand (CurveCommand):
 275     """Return a :class:`hooke.curve.Curve`.
 276     """
 277     def __init__(self, plugin):
 278         super(GetCommand, self).__init__(
 279             name='get curve', help=self.__doc__, plugin=plugin)
 280
 281     def _run(self, hooke, inqueue, outqueue, params):
 282         outqueue.put(self._curve(hooke, params))
 283
 284
 285 class InfoCommand (CurveCommand):
 286     """Get selected information about a :class:`hooke.curve.Curve`.
 287     """
 288     def __init__(self, plugin):
 289         args = [
 290             Argument(name='all', type='bool', default=False, count=1,
 291                      help='Get all curve information.'),
 292             ]
 293         self.fields = ['name', 'path', 'experiment', 'driver', 'filetype',
 294                        'note', 'command stack', 'blocks', 'block sizes']
 295         for field in self.fields:
 296             args.append(Argument(
 297                     name=field, type='bool', default=False, count=1,
 298                     help='Get curve %s' % field))
 299         super(InfoCommand, self).__init__(
 300             name='curve info', arguments=args,
 301             help=self.__doc__, plugin=plugin)
 302
 303     def _run(self, hooke, inqueue, outqueue, params):
 304         curve = self._curve(hooke, params)
 305         fields = {}
 306         for key in self.fields:
 307             fields[key] = params[key]
 308         if reduce(lambda x,y: x or y, fields.values()) == False:
 309             params['all'] = True # No specific fields set, default to 'all'
 310         if params['all'] == True:
 311             for key in self.fields:
 312                 fields[key] = True
 313         lines = []
 314         for key in self.fields:
 315             if fields[key] == True:
 316                 get = getattr(self, '_get_%s' % key.replace(' ', '_'))
 317                 lines.append('%s: %s' % (key, get(curve)))
 318         outqueue.put('\n'.join(lines))
 319
 320     def _get_name(self, curve):
 321         return curve.name
 322
 323     def _get_path(self, curve):
 324         return curve.path
 325
 326     def _get_experiment(self, curve):
 327         return curve.info.get('experiment', None)
 328
 329     def _get_driver(self, curve):
 330         return curve.driver
 331
 332     def _get_filetype(self, curve):
 333         return curve.info.get('filetype', None)
 334
 335     def _get_note(self, curve):
 336         return curve.info.get('note', None)
 337
 338     def _get_command_stack(self, curve):
 339         return curve.command_stack
 340
 341     def _get_blocks(self, curve):
 342         return len(curve.data)
 343
 344     def _get_block_sizes(self, curve):
 345         return [block.shape for block in curve.data]
 346
 347
 348 class DeltaCommand (BlockCommand):
 349     """Get distance information between two points.
 350
 351     With two points A and B, the returned distances are A-B.
 352     """
 353     def __init__(self, plugin):
 354         super(DeltaCommand, self).__init__(
 355             name='delta',
 356             arguments=[
 357                 Argument(name='point', type='point', optional=False, count=2,
 358                          help="""
 359 Indicies of points bounding the selected data.
 360 """.strip()),
 361                 Argument(name='SI', type='bool', default=False,
 362                          help="""
 363 Return distances in SI notation.
 364 """.strip())
 365                 ],
 366             help=self.__doc__, plugin=plugin)
 367
 368     def _run(self, hooke, inqueue, outqueue, params):
 369         data = self._block(hooke, params)
 370         As = data[params['point'][0],:]
 371         Bs = data[params['point'][1],:]
 372         ds = [A-B for A,B in zip(As, Bs)]
 373         if params['SI'] == False:
 374             out = [(name, d) for name,d in zip(data.info['columns'], ds)]
 375         else:
 376             out = []
 377             for name,d in zip(data.info['columns'], ds):
 378                 n,units = split_data_label(name)
 379                 out.append(
 380                   (n, ppSI(value=d, unit=units, decimals=2)))
 381         outqueue.put(out)
 382
 383
 384 class ExportCommand (BlockCommand):
 385     """Export a :class:`hooke.curve.Curve` data block as TAB-delimeted
 386     ASCII text.
 387
 388     A "#" prefixed header will optionally appear at the beginning of
 389     the file naming the columns.
 390     """
 391     def __init__(self, plugin):
 392         super(ExportCommand, self).__init__(
 393             name='export block',
 394             arguments=[
 395                 Argument(name='output', type='file', default='curve.dat',
 396                          help="""
 397 File name for the output data.  Defaults to 'curve.dat'
 398 """.strip()),
 399                 Argument(name='header', type='bool', default=True,
 400                          help="""
 401 True if you want the column-naming header line.
 402 """.strip()),
 403                 ],
 404             help=self.__doc__, plugin=plugin)
 405
 406     def _run(self, hooke, inqueue, outqueue, params):
 407         data = self._block(hooke, params)
 408
 409         with open(params['output'], 'w') as f:
 410             if params['header'] == True:
 411                 f.write('# %s \n' % ('\t'.join(data.info['columns'])))
 412             numpy.savetxt(f, data, delimiter='\t')
 413
 414
 415 class DifferenceCommand (ColumnAddingCommand):
 416     """Calculate the difference between two columns of data.
 417
 418     The difference is added to block A as a new column.
 419
 420     Note that the command will fail if the columns have different
 421     lengths, so be careful when differencing columns from different
 422     blocks.
 423     """
 424     def __init__(self, plugin):
 425         super(DifferenceCommand, self).__init__(
 426             name='difference',
 427             blocks=[
 428                 ('block A', None,
 429                  'Name of block A in A-B.  Defaults to the first block'),
 430                 ('block B', None,
 431                  'Name of block B in A-B.  Defaults to matching `block A`.'),
 432                 ],
 433             columns=[
 434                 ('column A', None,
 435                  """
 436 Column of data from block A to difference.  Defaults to the first column.
 437 """.strip()),
 438                 ('column B', None,
 439                  """
 440 Column of data from block B to difference.  Defaults to matching `column A`.
 441 """.strip()),
 442                 ],
 443             new_columns=[
 444                 ('output column', None,
 445                  """
 446 Name of the new column for storing the difference (without units, defaults to
 447 `difference of <block A> <column A> and <block B> <column B>`).
 448 """.strip()),
 449                 ],
 450             help=self.__doc__, plugin=plugin)
 451
 452     def _run(self, hooke, inqueue, outqueue, params):
 453         self._add_to_command_stack(params)
 454         params = self._setup_params(hooke=hooke, params=params)
 455         data_A = self._get_column(hooke=hooke, params=params,
 456                                   block_name='block A',
 457                                   column_name='column A')
 458         data_B = self._get_column(hooke=hooke, params=params,
 459                                   block_name='block B',
 460                                   column_name='column B')
 461         out = data_A - data_B
 462         self._set_column(hooke=hooke, params=params,
 463                          block_name='block A',
 464                          column_name='output column',
 465                          values=out)
 466
 467     def _setup_params(self, hooke, params):
 468         curve = self._curve(hooke, params)
 469         if params['block A'] == None:
 470             params['block A'] = curve.data[0].info['name']
 471         if params['block B'] == None:
 472             params['block B'] = params['block A']
 473         block_A = self._block(hooke, params=params, name='block A')
 474         block_B = self._block(hooke, params=params, name='block B')
 475         if params['column A'] == None:
 476             params['column A'] = block.info['columns'][0]
 477         if params['column B'] == None:
 478             params['column B'] = params['column A']
 479         a_name,a_unit = split_data_label(params['column A'])
 480         b_name,b_unit = split_data_label(params['column B'])
 481         if a_unit != b_unit:
 482             raise Failure('Unit missmatch: %s != %s' % (a_unit, b_unit))
 483         if params['output column'] == None:
 484             params['output column'] = join_data_label(
 485                 'difference of %s %s and %s %s' % (
 486                     block_A.info['name'], a_name,
 487                     block_B.info['name'], b_name),
 488                 a_unit)
 489         else:
 490             params['output column'] = join_data_label(
 491                 params['output column'], a_unit)
 492         return params
 493
 494
 495 class DerivativeCommand (ColumnAddingCommand):
 496     """Calculate the derivative (actually, the discrete differentiation)
 497     of a data column.
 498
 499     See :func:`hooke.util.calculus.derivative` for implementation
 500     details.
 501     """
 502     def __init__(self, plugin):
 503         super(DerivativeCommand, self).__init__(
 504             name='derivative',
 505             columns=[
 506                 ('x column', None,
 507                  'Column of data block to differentiate with respect to.'),
 508                 ('f column', None,
 509                  'Column of data block to differentiate.'),
 510                 ],
 511             new_columns=[
 512                 ('output column', None,
 513                  """
 514 Name of the new column for storing the derivative (without units, defaults to
 515 `derivative of <f column> with respect to <x column>`).
 516 """.strip()),
 517                 ],
 518             arguments=[
 519                 Argument(name='weights', type='dict', default={-1:-0.5, 1:0.5},
 520                          help="""
 521 Weighting scheme dictionary for finite differencing.  Defaults to
 522 central differencing.
 523 """.strip()),
 524                 ],
 525             help=self.__doc__, plugin=plugin)
 526
 527     def _run(self, hooke, inqueue, outqueue, params):
 528         self._add_to_command_stack(params)
 529         params = self._setup_params(hooke=hooke, params=params)
 530         x_data = self._get_column(hooke=hooke, params=params,
 531                                   column_name='x column')
 532         f_data = self._get_column(hooke=hooke, params=params,
 533                                   column_name='f column')
 534         d = derivative(
 535             x_data=x_data, f_data=f_data, weights=params['weights'])
 536         self._set_column(hooke=hooke, params=params,
 537                          column_name='output column',
 538                          values=d)
 539
 540     def _setup_params(self, hooke, params):
 541         curve = self._curve(hooke, params)
 542         x_name,x_unit = split_data_label(params['x column'])
 543         f_name,f_unit = split_data_label(params['f column'])
 544         d_unit = '%s/%s' % (f_unit, x_unit)
 545         if params['output column'] == None:
 546             params['output column'] = join_data_label(
 547                 'derivative of %s with respect to %s' % (
 548                     f_name, x_name),
 549                 d_unit)
 550         else:
 551             params['output column'] = join_data_label(
 552                 params['output column'], d_unit)
 553         return params
 554
 555
 556 class PowerSpectrumCommand (ColumnAddingCommand):
 557     """Calculate the power spectrum of a data column.
 558     """
 559     def __init__(self, plugin):
 560         super(PowerSpectrumCommand, self).__init__(
 561             name='power spectrum',
 562             arguments=[
 563                 Argument(name='output block', type='string',
 564                          help="""
 565 Name of the new data block for storing the power spectrum (defaults to
 566 `power spectrum of <source block name> <source column name>`).
 567 """.strip()),
 568                 Argument(name='bounds', type='point', optional=True, count=2,
 569                          help="""
 570 Indicies of points bounding the selected data.
 571 """.strip()),
 572                 Argument(name='freq', type='float', default=1.0,
 573                          help="""
 574 Sampling frequency.
 575 """.strip()),
 576                 Argument(name='freq units', type='string', default='Hz',
 577                          help="""
 578 Units for the sampling frequency.
 579 """.strip()),
 580                 Argument(name='chunk size', type='int', default=2048,
 581                          help="""
 582 Number of samples per chunk.  Use a power of two.
 583 """.strip()),
 584                 Argument(name='overlap', type='bool', default=False,
 585                          help="""
 586 If `True`, each chunk overlaps the previous chunk by half its length.
 587 Otherwise, the chunks are end-to-end, and not overlapping.
 588 """.strip()),
 589                 ],
 590             help=self.__doc__, plugin=plugin)
 591
 592     def _run(self, hooke, inqueue, outqueue, params):
 593         self._add_to_command_stack(params)
 594         params = self._setup_params(hooke=hooke, params=params)
 595         data = self._get_column(hooke=hooke, params=params)
 596         bounds = params['bounds']
 597         if bounds != None:
 598             data = data[bounds[0]:bounds[1]]
 599         freq_axis,power = unitary_avg_power_spectrum(
 600             data, freq=params['freq'],
 601             chunk_size=params['chunk size'],
 602             overlap=params['overlap'])
 603         b = Data((len(freq_axis),2), dtype=data.dtype)
 604         b.info['name'] = params['output block']
 605         b.info['columns'] = [
 606             params['output freq column'],
 607             params['output power column'],
 608             ]
 609         self._curve(hooke, params).data.append(b)
 610         self._set_column(hooke, params, block_name='output block',
 611                          column_name='output freq column',
 612                          values=freq_axis)
 613         self._set_column(hooke, params, block_name='output block',
 614                          column_name='output power column',
 615                          values=power)
 616         outqueue.put(b)
 617
 618     def _setup_params(self, hooke, params):
 619         if params['output block'] in self._block_names(hooke, params):
 620             raise Failure('output block %s already exists in %s.'
 621                           % (params['output block'],
 622                              self._curve(hooke, params)))
 623         data = self._get_column(hooke=hooke, params=params)
 624         d_name,d_unit = split_data_label(data.info['name'])
 625         if params['output block'] == None:
 626             params['output block'] = 'power spectrum of %s %s' % (
 627                 data.info['name'], params['column'])
 628         self.params['output freq column'] = join_data_label(
 629             'frequency axis', params['freq units'])
 630         self.params['output power column'] = join_data_label(
 631             'power density', '%s^2/%s' % (data_units, params['freq units']))
 632         return params
 633
 634
 635 class ClearStackCommand (CurveCommand):
 636     """Empty a curve's command stack.
 637     """
 638     def __init__(self, plugin):
 639         super(ClearStackCommand, self).__init__(
 640             name='clear curve command stack',
 641             help=self.__doc__, plugin=plugin)
 642         i,arg = [(i,arg) for i,arg in enumerate(self.arguments)
 643                  if arg.name == 'curve'][0]
 644         arg = copy.copy(arg)
 645         arg.callback = unloaded_current_curve_callback
 646         self.arguments[i] = arg
 647
 648     def _run(self, hooke, inqueue, outqueue, params):
 649         curve = self._curve(hooke, params)
 650         curve.command_stack = CommandStack()
 651
 652
 653 class OldCruft (object):
 654
 655     def do_forcebase(self,args):
 656         '''
 657         FORCEBASE
 658         (generalvclamp.py)
 659         Measures the difference in force (in pN) between a point and a baseline
 660         took as the average between two points.
 661
 662         The baseline is fixed once for a given curve and different force measurements,
 663         unless the user wants it to be recalculated
 664         ------------
 665         Syntax: forcebase [rebase]
 666                 rebase: Forces forcebase to ask again the baseline
 667                 max: Instead of asking for a point to measure, asks for two points and use
 668                      the maximum peak in between
 669         '''
 670         rebase=False #if true=we select rebase
 671         maxpoint=False #if true=we measure the maximum peak
 672
 673         plot=self._get_displayed_plot()
 674         whatset=1 #fixme: for all sets
 675         if 'rebase' in args or (self.basecurrent != self.current.path):
 676             rebase=True
 677         if 'max' in args:
 678             maxpoint=True
 679
 680         if rebase:
 681             print 'Select baseline'
 682             self.basepoints=self._measure_N_points(N=2, whatset=whatset)
 683             self.basecurrent=self.current.path
 684
 685         if maxpoint:
 686             print 'Select two points'
 687             points=self._measure_N_points(N=2, whatset=whatset)
 688             boundpoints=[points[0].index, points[1].index]
 689             boundpoints.sort()
 690             try:
 691                 y=min(plot.vectors[whatset][1][boundpoints[0]:boundpoints[1]])
 692             except ValueError:
 693                 print 'Chosen interval not valid. Try picking it again. Did you pick the same point as begin and end of interval?'
 694         else:
 695             print 'Select point to measure'
 696             points=self._measure_N_points(N=1, whatset=whatset)
 697             #whatplot=points[0].dest
 698             y=points[0].graph_coords[1]
 699
 700         #fixme: code duplication
 701         boundaries=[self.basepoints[0].index, self.basepoints[1].index]
 702         boundaries.sort()
 703         to_average=plot.vectors[whatset][1][boundaries[0]:boundaries[1]] #y points to average
 704
 705         avg=np.mean(to_average)
 706         forcebase=abs(y-avg)
 707         print str(forcebase*(10**12))+' pN'
 708         to_dump='forcebase '+self.current.path+' '+str(forcebase*(10**12))+' pN'
 709         self.outlet.push(to_dump)
 710
 711     #---SLOPE---
 712     def do_slope(self,args):
 713         '''
 714         SLOPE
 715         (generalvclamp.py)
 716         Measures the slope of a delimited chunk on the return trace.
 717         The chunk can be delimited either by two manual clicks, or have
 718         a fixed width, given as an argument.
 719         ---------------
 720         Syntax: slope [width]
 721                 The facultative [width] parameter specifies how many
 722                 points will be considered for the fit. If [width] is
 723                 specified, only one click will be required.
 724         (c) Marco Brucale, Massimo Sandal 2008
 725         '''
 726
 727         # Reads the facultative width argument
 728         try:
 729             fitspan=int(args)
 730         except:
 731             fitspan=0
 732
 733         # Decides between the two forms of user input, as per (args)
 734         if fitspan == 0:
 735             # Gets the Xs of two clicked points as indexes on the current curve vector
 736             print 'Click twice to delimit chunk'
 737             points=self._measure_N_points(N=2,whatset=1)
 738         else:
 739             print 'Click once on the leftmost point of the chunk (i.e.usually the peak)'
 740             points=self._measure_N_points(N=1,whatset=1)
 741
 742         slope=self._slope(points,fitspan)
 743
 744         # Outputs the relevant slope parameter
 745         print 'Slope:'
 746         print str(slope)
 747         to_dump='slope '+self.current.path+' '+str(slope)
 748         self.outlet.push(to_dump)
 749
 750     def _slope(self,points,fitspan):
 751         # Calls the function linefit_between
 752         parameters=[0,0,[],[]]
 753         try:
 754             clickedpoints=[points[0].index,points[1].index]
 755             clickedpoints.sort()
 756         except:
 757             clickedpoints=[points[0].index-fitspan,points[0].index]
 758
 759         try:
 760             parameters=self.linefit_between(clickedpoints[0],clickedpoints[1])
 761         except:
 762             print 'Cannot fit. Did you click twice the same point?'
 763             return
 764
 765         # Outputs the relevant slope parameter
 766         print 'Slope:'
 767         print str(parameters[0])
 768         to_dump='slope '+self.curve.path+' '+str(parameters[0])
 769         self.outlet.push(to_dump)
 770
 771         # Makes a vector with the fitted parameters and sends it to the GUI
 772         xtoplot=parameters[2]
 773         ytoplot=[]
 774         x=0
 775         for x in xtoplot:
 776             ytoplot.append((x*parameters[0])+parameters[1])
 777
 778         clickvector_x, clickvector_y=[], []
 779         for item in points:
 780             clickvector_x.append(item.graph_coords[0])
 781             clickvector_y.append(item.graph_coords[1])
 782
 783         lineplot=self._get_displayed_plot(0) #get topmost displayed plot
 784
 785         lineplot.add_set(xtoplot,ytoplot)
 786         lineplot.add_set(clickvector_x, clickvector_y)
 787
 788
 789         if lineplot.styles==[]:
 790             lineplot.styles=[None,None,None,'scatter']
 791         else:
 792             lineplot.styles+=[None,'scatter']
 793         if lineplot.colors==[]:
 794             lineplot.colors=[None,None,'black',None]
 795         else:
 796             lineplot.colors+=['black',None]
 797
 798
 799         self._send_plot([lineplot])
 800
 801         return parameters[0]
 802
 803
 804     def linefit_between(self,index1,index2,whatset=1):
 805         '''
 806         Creates two vectors (xtofit,ytofit) slicing out from the
 807         current return trace a portion delimited by the two indexes
 808         given as arguments.
 809         Then does a least squares linear fit on that slice.
 810         Finally returns [0]=the slope, [1]=the intercept of the
 811         fitted 1st grade polynomial, and [2,3]=the actual (x,y) vectors
 812         used for the fit.
 813         (c) Marco Brucale, Massimo Sandal 2008
 814         '''
 815         # Translates the indexes into two vectors containing the x,y data to fit
 816         xtofit=self.plots[0].vectors[whatset][0][index1:index2]
 817         ytofit=self.plots[0].vectors[whatset][1][index1:index2]
 818
 819         # Does the actual linear fitting (simple least squares with numpy.polyfit)
 820         linefit=[]
 821         linefit=np.polyfit(xtofit,ytofit,1)
 822
 823         return (linefit[0],linefit[1],xtofit,ytofit)