1 # Copyright (C) 2008-2010 Alberto Gomez-Kasai
3 # Massimo Sandal <devicerandom@gmail.com>
4 # W. Trevor King <wking@drexel.edu>
6 # This file is part of Hooke.
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.
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.
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/>.
22 """The ``curve`` module provides :class:`CurvePlugin` and several
23 associated :class:`hooke.command.Command`\s for handling
24 :mod:`hooke.curve` classes.
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
38 from .playlist import current_playlist_callback
41 # Define common or complicated arguments
43 def current_curve_callback(hooke, command, argument, value):
46 playlist = current_playlist_callback(hooke, command, argument, value)
47 curve = playlist.current()
49 raise Failure('No curves in %s' % playlist)
52 CurveArgument = Argument(
53 name='curve', type='curve', callback=current_curve_callback,
55 :class:`hooke.curve.Curve` to act on. Defaults to the current curve
56 of the current playlist.
59 def _name_argument(name, default, help):
62 return Argument(name=name, type='string', default=default, help=help)
64 def block_argument(*args, **kwargs):
67 return _name_argument(*args, **kwargs)
69 def column_argument(*args, **kwargs):
72 return _name_argument(*args, **kwargs)
75 # Define useful command subclasses
77 class CurveCommand (Command):
78 """A :class:`~hooke.command.Command` operating on a
79 :class:`~hooke.curve.Curve`.
81 def __init__(self, **kwargs):
82 if 'arguments' in kwargs:
83 kwargs['arguments'].insert(0, CurveArgument)
85 kwargs['arguments'] = [CurveArgument]
86 super(CurveCommand, self).__init__(**kwargs)
88 def _curve(self, hooke, params):
89 """Get the selected curve.
93 `hooke` is intended to attach the selected curve to the local
94 playlist; 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.
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
103 return params['curve']
105 def _add_to_command_stack(self, params):
106 """Store the command name and current `params` values in the
107 curve's `.command_stack`.
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.
113 Recommended practice is to *not* lock in argument values that
114 are loaded from the plugin's :attr:`.config`.
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
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.
132 curve.command_stack.append(CommandMessage(
136 class BlockCommand (CurveCommand):
137 """A :class:`CurveCommand` operating on a :class:`~hooke.curve.Data` block.
139 def __init__(self, blocks=None, **kwargs):
141 blocks = [('block', None, 'Name of the data block to act on.')]
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)
151 def _block_names(self, hooke, params):
152 curve = self._curve(hooke, params)
153 return [b.info['name'] for b in curve.data]
155 def _block_index(self, hooke, params, 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)
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))
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
178 curve = self._curve(hooke, params)
179 index = self._block_index(hooke, params, name)
180 return curve.data[index]
183 class ColumnAccessCommand (BlockCommand):
184 """A :class:`BlockCommand` accessing a :class:`~hooke.curve.Data`
187 def __init__(self, columns=None, **kwargs):
189 columns = [('column', None, 'Name of the data column to act on.')]
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)
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']
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]
213 class ColumnAddingCommand (ColumnAccessCommand):
214 """A :class:`ColumnAccessCommand` that also adds columns.
216 def __init__(self, new_columns=None, **kwargs):
217 if new_columns == None:
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)
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)
235 def _set_column(self, hooke, params, block_name=None, column_name=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)
245 new.info['columns'].append(column_name)
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
255 class CurvePlugin (Builtin):
257 super(CurvePlugin, self).__init__(name='curve')
259 GetCommand(self), InfoCommand(self), DeltaCommand(self),
260 ExportCommand(self), DifferenceCommand(self),
261 DerivativeCommand(self), PowerSpectrumCommand(self)]
266 class GetCommand (CurveCommand):
267 """Return a :class:`hooke.curve.Curve`.
269 def __init__(self, plugin):
270 super(GetCommand, self).__init__(
271 name='get curve', help=self.__doc__, plugin=plugin)
273 def _run(self, hooke, inqueue, outqueue, params):
274 outqueue.put(self._curve(hooke, params))
277 class InfoCommand (CurveCommand):
278 """Get selected information about a :class:`hooke.curve.Curve`.
280 def __init__(self, plugin):
282 Argument(name='all', type='bool', default=False, count=1,
283 help='Get all curve information.'),
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)
295 def _run(self, hooke, inqueue, outqueue, params):
296 curve = self._curve(hooke, params)
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:
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))
312 def _get_name(self, curve):
315 def _get_path(self, curve):
318 def _get_experiment(self, curve):
319 return curve.info.get('experiment', None)
321 def _get_driver(self, curve):
324 def _get_filetype(self, curve):
325 return curve.info.get('filetype', None)
327 def _get_note(self, curve):
328 return curve.info.get('note', None)
330 def _get_blocks(self, curve):
331 return len(curve.data)
333 def _get_block_sizes(self, curve):
334 return [block.shape for block in curve.data]
337 class DeltaCommand (BlockCommand):
338 """Get distance information between two points.
340 With two points A and B, the returned distances are A-B.
342 def __init__(self, plugin):
343 super(DeltaCommand, self).__init__(
346 Argument(name='point', type='point', optional=False, count=2,
348 Indicies of points bounding the selected data.
350 Argument(name='SI', type='bool', default=False,
352 Return distances in SI notation.
355 help=self.__doc__, plugin=plugin)
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)]
366 for name,d in zip(data.info['columns'], ds):
367 n,units = split_data_label(name)
369 (n, ppSI(value=d, unit=units, decimals=2)))
373 class ExportCommand (BlockCommand):
374 """Export a :class:`hooke.curve.Curve` data block as TAB-delimeted
377 A "#" prefixed header will optionally appear at the beginning of
378 the file naming the columns.
380 def __init__(self, plugin):
381 super(ExportCommand, self).__init__(
384 Argument(name='output', type='file', default='curve.dat',
386 File name for the output data. Defaults to 'curve.dat'
388 Argument(name='header', type='bool', default=True,
390 True if you want the column-naming header line.
393 help=self.__doc__, plugin=plugin)
395 def _run(self, hooke, inqueue, outqueue, params):
396 data = self._block(hooke, params)
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')
404 class DifferenceCommand (ColumnAddingCommand):
405 """Calculate the difference between two columns of data.
407 The difference is added to block A as a new column.
409 Note that the command will fail if the columns have different
410 lengths, so be careful when differencing columns from different
413 def __init__(self, plugin):
414 super(DifferenceCommand, self).__init__(
418 'Name of block A in A-B. Defaults to the first block'),
420 'Name of block B in A-B. Defaults to matching `block A`.'),
425 Column of data from block A to difference. Defaults to the first column.
429 Column of data from block B to difference. Defaults to matching `column A`.
433 ('output column', None,
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>`).
439 help=self.__doc__, plugin=plugin)
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',
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'])
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),
479 params['output column'] = join_data_label(
480 params['output column'], a_unit)
484 class DerivativeCommand (ColumnAddingCommand):
485 """Calculate the derivative (actually, the discrete differentiation)
488 See :func:`hooke.util.calculus.derivative` for implementation
491 def __init__(self, plugin):
492 super(DerivativeCommand, self).__init__(
496 'Column of data block to differentiate with respect to.'),
498 'Column of data block to differentiate.'),
501 ('output column', None,
503 Name of the new column for storing the derivative (without units, defaults to
504 `derivative of <f column> with respect to <x column>`).
508 Argument(name='weights', type='dict', default={-1:-0.5, 1:0.5},
510 Weighting scheme dictionary for finite differencing. Defaults to
511 central differencing.
514 help=self.__doc__, plugin=plugin)
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')
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',
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' % (
540 params['output column'] = join_data_label(
541 params['output column'], d_unit)
545 class PowerSpectrumCommand (ColumnAddingCommand):
546 """Calculate the power spectrum of a data column.
548 def __init__(self, plugin):
549 super(PowerSpectrumCommand, self).__init__(
550 name='power spectrum',
552 Argument(name='output block', type='string',
554 Name of the new data block for storing the power spectrum (defaults to
555 `power spectrum of <source block name> <source column name>`).
557 Argument(name='bounds', type='point', optional=True, count=2,
559 Indicies of points bounding the selected data.
561 Argument(name='freq', type='float', default=1.0,
565 Argument(name='freq units', type='string', default='Hz',
567 Units for the sampling frequency.
569 Argument(name='chunk size', type='int', default=2048,
571 Number of samples per chunk. Use a power of two.
573 Argument(name='overlap', type='bool', default=False,
575 If `True`, each chunk overlaps the previous chunk by half its length.
576 Otherwise, the chunks are end-to-end, and not overlapping.
579 help=self.__doc__, plugin=plugin)
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']
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'],
598 self._curve(hooke, params).data.append(b)
599 self._set_column(hooke, params, block_name='output block',
600 column_name='output freq column',
602 self._set_column(hooke, params, block_name='output block',
603 column_name='output power column',
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']))
624 class OldCruft (object):
626 def do_forcebase(self,args):
630 Measures the difference in force (in pN) between a point and a baseline
631 took as the average between two points.
633 The baseline is fixed once for a given curve and different force measurements,
634 unless the user wants it to be recalculated
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
641 rebase=False #if true=we select rebase
642 maxpoint=False #if true=we measure the maximum peak
644 plot=self._get_displayed_plot()
645 whatset=1 #fixme: for all sets
646 if 'rebase' in args or (self.basecurrent != self.current.path):
652 print 'Select baseline'
653 self.basepoints=self._measure_N_points(N=2, whatset=whatset)
654 self.basecurrent=self.current.path
657 print 'Select two points'
658 points=self._measure_N_points(N=2, whatset=whatset)
659 boundpoints=[points[0].index, points[1].index]
662 y=min(plot.vectors[whatset][1][boundpoints[0]:boundpoints[1]])
664 print 'Chosen interval not valid. Try picking it again. Did you pick the same point as begin and end of interval?'
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]
671 #fixme: code duplication
672 boundaries=[self.basepoints[0].index, self.basepoints[1].index]
674 to_average=plot.vectors[whatset][1][boundaries[0]:boundaries[1]] #y points to average
676 avg=np.mean(to_average)
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)
683 def do_slope(self,args):
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.
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
698 # Reads the facultative width argument
704 # Decides between the two forms of user input, as per (args)
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)
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)
713 slope=self._slope(points,fitspan)
715 # Outputs the relevant slope parameter
718 to_dump='slope '+self.current.path+' '+str(slope)
719 self.outlet.push(to_dump)
721 def _slope(self,points,fitspan):
722 # Calls the function linefit_between
723 parameters=[0,0,[],[]]
725 clickedpoints=[points[0].index,points[1].index]
728 clickedpoints=[points[0].index-fitspan,points[0].index]
731 parameters=self.linefit_between(clickedpoints[0],clickedpoints[1])
733 print 'Cannot fit. Did you click twice the same point?'
736 # Outputs the relevant slope parameter
738 print str(parameters[0])
739 to_dump='slope '+self.curve.path+' '+str(parameters[0])
740 self.outlet.push(to_dump)
742 # Makes a vector with the fitted parameters and sends it to the GUI
743 xtoplot=parameters[2]
747 ytoplot.append((x*parameters[0])+parameters[1])
749 clickvector_x, clickvector_y=[], []
751 clickvector_x.append(item.graph_coords[0])
752 clickvector_y.append(item.graph_coords[1])
754 lineplot=self._get_displayed_plot(0) #get topmost displayed plot
756 lineplot.add_set(xtoplot,ytoplot)
757 lineplot.add_set(clickvector_x, clickvector_y)
760 if lineplot.styles==[]:
761 lineplot.styles=[None,None,None,'scatter']
763 lineplot.styles+=[None,'scatter']
764 if lineplot.colors==[]:
765 lineplot.colors=[None,None,'black',None]
767 lineplot.colors+=['black',None]
770 self._send_plot([lineplot])
775 def linefit_between(self,index1,index2,whatset=1):
777 Creates two vectors (xtofit,ytofit) slicing out from the
778 current return trace a portion delimited by the two indexes
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
784 (c) Marco Brucale, Massimo Sandal 2008
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]
790 # Does the actual linear fitting (simple least squares with numpy.polyfit)
792 linefit=np.polyfit(xtofit,ytofit,1)
794 return (linefit[0],linefit[1],xtofit,ytofit)