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