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 ..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
39 from .playlist import current_playlist_callback
42 # Define common or complicated arguments
44 def current_curve_callback(hooke, command, argument, value):
47 playlist = current_playlist_callback(hooke, command, argument, value)
48 curve = playlist.current()
50 raise Failure('No curves in %s' % playlist)
53 CurveArgument = Argument(
54 name='curve', type='curve', callback=current_curve_callback,
56 :class:`hooke.curve.Curve` to act on. Defaults to the current curve
57 of the current playlist.
60 def _name_argument(name, default, help):
63 return Argument(name=name, type='string', default=default, help=help)
65 def block_argument(*args, **kwargs):
68 return _name_argument(*args, **kwargs)
70 def column_argument(*args, **kwargs):
73 return _name_argument(*args, **kwargs)
76 # Define useful command subclasses
78 class CurveCommand (Command):
79 """A :class:`~hooke.command.Command` operating on a
80 :class:`~hooke.curve.Curve`.
82 def __init__(self, **kwargs):
83 if 'arguments' in kwargs:
84 kwargs['arguments'].insert(0, CurveArgument)
86 kwargs['arguments'] = [CurveArgument]
87 super(CurveCommand, self).__init__(**kwargs)
89 def _curve(self, hooke, params):
90 """Get the selected curve.
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.
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
104 return params['curve']
106 def _add_to_command_stack(self, params):
107 """Store the command name and current `params` values in the
108 curve's `.command_stack`.
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.
114 Recommended practice is to *not* lock in argument values that
115 are loaded from the plugin's :attr:`.config`.
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
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.
134 curve.command_stack.append(CommandMessage(
138 class BlockCommand (CurveCommand):
139 """A :class:`CurveCommand` operating on a :class:`~hooke.curve.Data` block.
141 def __init__(self, blocks=None, **kwargs):
143 blocks = [('block', None, 'Name of the data block to act on.')]
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)
153 def _block_names(self, hooke, params):
154 curve = self._curve(hooke, params)
155 return [b.info['name'] for b in curve.data]
157 def _block_index(self, hooke, params, 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)
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))
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
180 curve = self._curve(hooke, params)
181 index = self._block_index(hooke, params, name)
182 return curve.data[index]
185 class ColumnAccessCommand (BlockCommand):
186 """A :class:`BlockCommand` accessing a :class:`~hooke.curve.Data`
189 def __init__(self, columns=None, **kwargs):
191 columns = [('column', None, 'Name of the data column to act on.')]
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)
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']
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]
215 class ColumnAddingCommand (ColumnAccessCommand):
216 """A :class:`ColumnAccessCommand` that also adds columns.
218 def __init__(self, new_columns=None, **kwargs):
219 if new_columns == None:
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)
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)
237 def _set_column(self, hooke, params, block_name=None, column_name=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)
247 new.info['columns'].append(column_name)
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
257 class CurvePlugin (Builtin):
259 super(CurvePlugin, self).__init__(name='curve')
261 GetCommand(self), InfoCommand(self), DeltaCommand(self),
262 ExportCommand(self), DifferenceCommand(self),
263 DerivativeCommand(self), PowerSpectrumCommand(self),
264 ClearStackCommand(self)]
269 class GetCommand (CurveCommand):
270 """Return a :class:`hooke.curve.Curve`.
272 def __init__(self, plugin):
273 super(GetCommand, self).__init__(
274 name='get curve', help=self.__doc__, plugin=plugin)
276 def _run(self, hooke, inqueue, outqueue, params):
277 outqueue.put(self._curve(hooke, params))
280 class InfoCommand (CurveCommand):
281 """Get selected information about a :class:`hooke.curve.Curve`.
283 def __init__(self, plugin):
285 Argument(name='all', type='bool', default=False, count=1,
286 help='Get all curve information.'),
288 self.fields = ['name', 'path', 'experiment', 'driver', 'filetype', 'note',
289 '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)
298 def _run(self, hooke, inqueue, outqueue, params):
299 curve = self._curve(hooke, params)
301 for key in self.fields:
302 fields[key] = params[key]
303 if reduce(lambda x,y: x and 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:
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))
315 def _get_name(self, curve):
318 def _get_path(self, curve):
321 def _get_experiment(self, curve):
322 return curve.info.get('experiment', None)
324 def _get_driver(self, curve):
327 def _get_filetype(self, curve):
328 return curve.info.get('filetype', None)
330 def _get_note(self, curve):
331 return curve.info.get('note', None)
333 def _get_blocks(self, curve):
334 return len(curve.data)
336 def _get_block_sizes(self, curve):
337 return [block.shape for block in curve.data]
340 class DeltaCommand (BlockCommand):
341 """Get distance information between two points.
343 With two points A and B, the returned distances are A-B.
345 def __init__(self, plugin):
346 super(DeltaCommand, self).__init__(
349 Argument(name='point', type='point', optional=False, count=2,
351 Indicies of points bounding the selected data.
353 Argument(name='SI', type='bool', default=False,
355 Return distances in SI notation.
358 help=self.__doc__, plugin=plugin)
360 def _run(self, hooke, inqueue, outqueue, params):
361 data = self._block(hooke, params)
362 As = data[params['point'][0],:]
363 Bs = data[params['point'][1],:]
364 ds = [A-B for A,B in zip(As, Bs)]
365 if params['SI'] == False:
366 out = [(name, d) for name,d in zip(data.info['columns'], ds)]
369 for name,d in zip(data.info['columns'], ds):
370 n,units = split_data_label(name)
372 (n, ppSI(value=d, unit=units, decimals=2)))
376 class ExportCommand (BlockCommand):
377 """Export a :class:`hooke.curve.Curve` data block as TAB-delimeted
380 A "#" prefixed header will optionally appear at the beginning of
381 the file naming the columns.
383 def __init__(self, plugin):
384 super(ExportCommand, self).__init__(
387 Argument(name='output', type='file', default='curve.dat',
389 File name for the output data. Defaults to 'curve.dat'
391 Argument(name='header', type='bool', default=True,
393 True if you want the column-naming header line.
396 help=self.__doc__, plugin=plugin)
398 def _run(self, hooke, inqueue, outqueue, params):
399 data = self._block(hooke, params)
401 with open(params['output'], 'w') as f:
402 if params['header'] == True:
403 f.write('# %s \n' % ('\t'.join(data.info['columns'])))
404 numpy.savetxt(f, data, delimiter='\t')
407 class DifferenceCommand (ColumnAddingCommand):
408 """Calculate the difference between two columns of data.
410 The difference is added to block A as a new column.
412 Note that the command will fail if the columns have different
413 lengths, so be careful when differencing columns from different
416 def __init__(self, plugin):
417 super(DifferenceCommand, self).__init__(
421 'Name of block A in A-B. Defaults to the first block'),
423 'Name of block B in A-B. Defaults to matching `block A`.'),
428 Column of data from block A to difference. Defaults to the first column.
432 Column of data from block B to difference. Defaults to matching `column A`.
436 ('output column', None,
438 Name of the new column for storing the difference (without units, defaults to
439 `difference of <block A> <column A> and <block B> <column B>`).
442 help=self.__doc__, plugin=plugin)
444 def _run(self, hooke, inqueue, outqueue, params):
445 self._add_to_command_stack(params)
446 params = self.__setup_params(hooke=hooke, params=params)
447 data_A = self._get_column(hooke=hooke, params=params,
448 block_name='block A',
449 column_name='column A')
450 data_B = self._get_column(hooke=hooke, params=params,
451 block_name='block B',
452 column_name='column B')
453 out = data_A - data_B
454 self._set_column(hooke=hooke, params=params,
455 block_name='block A',
456 column_name='output column',
459 def __setup_params(self, hooke, params):
460 curve = self._curve(hooke, params)
461 if params['block A'] == None:
462 params['block A'] = curve.data[0].info['name']
463 if params['block B'] == None:
464 params['block B'] = params['block A']
465 block_A = self._block(hooke, params=params, name='block A')
466 block_B = self._block(hooke, params=params, name='block B')
467 if params['column A'] == None:
468 params['column A'] = block.info['columns'][0]
469 if params['column B'] == None:
470 params['column B'] = params['column A']
471 a_name,a_unit = split_data_label(params['column A'])
472 b_name,b_unit = split_data_label(params['column B'])
474 raise Failure('Unit missmatch: %s != %s' % (a_unit, b_unit))
475 if params['output column'] == None:
476 params['output column'] = join_data_label(
477 'difference of %s %s and %s %s' % (
478 block_A.info['name'], a_name,
479 block_B.info['name'], b_name),
482 params['output column'] = join_data_label(
483 params['output column'], a_unit)
487 class DerivativeCommand (ColumnAddingCommand):
488 """Calculate the derivative (actually, the discrete differentiation)
491 See :func:`hooke.util.calculus.derivative` for implementation
494 def __init__(self, plugin):
495 super(DerivativeCommand, self).__init__(
499 'Column of data block to differentiate with respect to.'),
501 'Column of data block to differentiate.'),
504 ('output column', None,
506 Name of the new column for storing the derivative (without units, defaults to
507 `derivative of <f column> with respect to <x column>`).
511 Argument(name='weights', type='dict', default={-1:-0.5, 1:0.5},
513 Weighting scheme dictionary for finite differencing. Defaults to
514 central differencing.
517 help=self.__doc__, plugin=plugin)
519 def _run(self, hooke, inqueue, outqueue, params):
520 self._add_to_command_stack(params)
521 params = self.__setup_params(hooke=hooke, params=params)
522 x_data = self._get_column(hooke=hooke, params=params,
523 column_name='x column')
524 f_data = self._get_column(hooke=hooke, params=params,
525 column_name='f column')
527 x_data=x_data, f_data=f_data, weights=params['weights'])
528 self._set_column(hooke=hooke, params=params,
529 column_name='output column',
532 def __setup_params(self, hooke, params):
533 curve = self._curve(hooke, params)
534 x_name,x_unit = split_data_label(params['x column'])
535 f_name,f_unit = split_data_label(params['f column'])
536 d_unit = '%s/%s' % (f_unit, x_unit)
537 if params['output column'] == None:
538 params['output column'] = join_data_label(
539 'derivative of %s with respect to %s' % (
543 params['output column'] = join_data_label(
544 params['output column'], d_unit)
548 class PowerSpectrumCommand (ColumnAddingCommand):
549 """Calculate the power spectrum of a data column.
551 def __init__(self, plugin):
552 super(PowerSpectrumCommand, self).__init__(
553 name='power spectrum',
555 Argument(name='output block', type='string',
557 Name of the new data block for storing the power spectrum (defaults to
558 `power spectrum of <source block name> <source column name>`).
560 Argument(name='bounds', type='point', optional=True, count=2,
562 Indicies of points bounding the selected data.
564 Argument(name='freq', type='float', default=1.0,
568 Argument(name='freq units', type='string', default='Hz',
570 Units for the sampling frequency.
572 Argument(name='chunk size', type='int', default=2048,
574 Number of samples per chunk. Use a power of two.
576 Argument(name='overlap', type='bool', default=False,
578 If `True`, each chunk overlaps the previous chunk by half its length.
579 Otherwise, the chunks are end-to-end, and not overlapping.
582 help=self.__doc__, plugin=plugin)
584 def _run(self, hooke, inqueue, outqueue, params):
585 self._add_to_command_stack(params)
586 params = self.__setup_params(hooke=hooke, params=params)
587 data = self._get_column(hooke=hooke, params=params)
588 bounds = params['bounds']
590 data = data[bounds[0]:bounds[1]]
591 freq_axis,power = unitary_avg_power_spectrum(
592 data, freq=params['freq'],
593 chunk_size=params['chunk size'],
594 overlap=params['overlap'])
595 b = Data((len(freq_axis),2), dtype=data.dtype)
596 b.info['name'] = params['output block']
597 b.info['columns'] = [
598 params['output freq column'],
599 params['output power column'],
601 self._curve(hooke, params).data.append(b)
602 self._set_column(hooke, params, block_name='output block',
603 column_name='output freq column',
605 self._set_column(hooke, params, block_name='output block',
606 column_name='output power column',
610 def __setup_params(self, hooke, params):
611 if params['output block'] in self._block_names(hooke, params):
612 raise Failure('output block %s already exists in %s.'
613 % (params['output block'],
614 self._curve(hooke, params)))
615 data = self._get_column(hooke=hooke, params=params)
616 d_name,d_unit = split_data_label(data.info['name'])
617 if params['output block'] == None:
618 params['output block'] = 'power spectrum of %s %s' % (
619 data.info['name'], params['column'])
620 self.params['output freq column'] = join_data_label(
621 'frequency axis', params['freq units'])
622 self.params['output power column'] = join_data_label(
623 'power density', '%s^2/%s' % (data_units, params['freq units']))
627 class ClearStackCommand (CurveCommand):
628 """Empty a curve's command stack.
630 def __init__(self, plugin):
631 super(ClearStackCommand, self).__init__(
632 name='clear curve command stack',
633 help=self.__doc__, plugin=plugin)
635 def _run(self, hooke, inqueue, outqueue, params):
636 curve = self._curve(hooke, params)
637 curve.command_stack = CommandStack()
640 class OldCruft (object):
642 def do_forcebase(self,args):
646 Measures the difference in force (in pN) between a point and a baseline
647 took as the average between two points.
649 The baseline is fixed once for a given curve and different force measurements,
650 unless the user wants it to be recalculated
652 Syntax: forcebase [rebase]
653 rebase: Forces forcebase to ask again the baseline
654 max: Instead of asking for a point to measure, asks for two points and use
655 the maximum peak in between
657 rebase=False #if true=we select rebase
658 maxpoint=False #if true=we measure the maximum peak
660 plot=self._get_displayed_plot()
661 whatset=1 #fixme: for all sets
662 if 'rebase' in args or (self.basecurrent != self.current.path):
668 print 'Select baseline'
669 self.basepoints=self._measure_N_points(N=2, whatset=whatset)
670 self.basecurrent=self.current.path
673 print 'Select two points'
674 points=self._measure_N_points(N=2, whatset=whatset)
675 boundpoints=[points[0].index, points[1].index]
678 y=min(plot.vectors[whatset][1][boundpoints[0]:boundpoints[1]])
680 print 'Chosen interval not valid. Try picking it again. Did you pick the same point as begin and end of interval?'
682 print 'Select point to measure'
683 points=self._measure_N_points(N=1, whatset=whatset)
684 #whatplot=points[0].dest
685 y=points[0].graph_coords[1]
687 #fixme: code duplication
688 boundaries=[self.basepoints[0].index, self.basepoints[1].index]
690 to_average=plot.vectors[whatset][1][boundaries[0]:boundaries[1]] #y points to average
692 avg=np.mean(to_average)
694 print str(forcebase*(10**12))+' pN'
695 to_dump='forcebase '+self.current.path+' '+str(forcebase*(10**12))+' pN'
696 self.outlet.push(to_dump)
699 def do_slope(self,args):
703 Measures the slope of a delimited chunk on the return trace.
704 The chunk can be delimited either by two manual clicks, or have
705 a fixed width, given as an argument.
707 Syntax: slope [width]
708 The facultative [width] parameter specifies how many
709 points will be considered for the fit. If [width] is
710 specified, only one click will be required.
711 (c) Marco Brucale, Massimo Sandal 2008
714 # Reads the facultative width argument
720 # Decides between the two forms of user input, as per (args)
722 # Gets the Xs of two clicked points as indexes on the current curve vector
723 print 'Click twice to delimit chunk'
724 points=self._measure_N_points(N=2,whatset=1)
726 print 'Click once on the leftmost point of the chunk (i.e.usually the peak)'
727 points=self._measure_N_points(N=1,whatset=1)
729 slope=self._slope(points,fitspan)
731 # Outputs the relevant slope parameter
734 to_dump='slope '+self.current.path+' '+str(slope)
735 self.outlet.push(to_dump)
737 def _slope(self,points,fitspan):
738 # Calls the function linefit_between
739 parameters=[0,0,[],[]]
741 clickedpoints=[points[0].index,points[1].index]
744 clickedpoints=[points[0].index-fitspan,points[0].index]
747 parameters=self.linefit_between(clickedpoints[0],clickedpoints[1])
749 print 'Cannot fit. Did you click twice the same point?'
752 # Outputs the relevant slope parameter
754 print str(parameters[0])
755 to_dump='slope '+self.curve.path+' '+str(parameters[0])
756 self.outlet.push(to_dump)
758 # Makes a vector with the fitted parameters and sends it to the GUI
759 xtoplot=parameters[2]
763 ytoplot.append((x*parameters[0])+parameters[1])
765 clickvector_x, clickvector_y=[], []
767 clickvector_x.append(item.graph_coords[0])
768 clickvector_y.append(item.graph_coords[1])
770 lineplot=self._get_displayed_plot(0) #get topmost displayed plot
772 lineplot.add_set(xtoplot,ytoplot)
773 lineplot.add_set(clickvector_x, clickvector_y)
776 if lineplot.styles==[]:
777 lineplot.styles=[None,None,None,'scatter']
779 lineplot.styles+=[None,'scatter']
780 if lineplot.colors==[]:
781 lineplot.colors=[None,None,'black',None]
783 lineplot.colors+=['black',None]
786 self._send_plot([lineplot])
791 def linefit_between(self,index1,index2,whatset=1):
793 Creates two vectors (xtofit,ytofit) slicing out from the
794 current return trace a portion delimited by the two indexes
796 Then does a least squares linear fit on that slice.
797 Finally returns [0]=the slope, [1]=the intercept of the
798 fitted 1st grade polynomial, and [2,3]=the actual (x,y) vectors
800 (c) Marco Brucale, Massimo Sandal 2008
802 # Translates the indexes into two vectors containing the x,y data to fit
803 xtofit=self.plots[0].vectors[whatset][0][index1:index2]
804 ytofit=self.plots[0].vectors[whatset][1][index1:index2]
806 # Does the actual linear fitting (simple least squares with numpy.polyfit)
808 linefit=np.polyfit(xtofit,ytofit,1)
810 return (linefit[0],linefit[1],xtofit,ytofit)