hooke/plugin/curve.py

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