1 # Copyright (C) 2008-2010 Alberto Gomez-Casado
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, load=True):
47 playlist = current_playlist_callback(hooke, command, argument, value)
48 curve = playlist.current(load=load)
50 raise Failure('No curves in %s' % playlist)
53 def unloaded_current_curve_callback(hooke, command, argument, value):
54 return current_curve_callback(
55 hooke=hooke, command=command, argument=argument, value=value,
58 CurveArgument = Argument(
59 name='curve', type='curve', callback=current_curve_callback,
61 :class:`hooke.curve.Curve` to act on. Defaults to the current curve
62 of the current playlist.
65 def _name_argument(name, default, help):
68 return Argument(name=name, type='string', default=default, help=help)
70 def block_argument(*args, **kwargs):
73 return _name_argument(*args, **kwargs)
75 def column_argument(*args, **kwargs):
78 return _name_argument(*args, **kwargs)
81 # Define useful command subclasses
83 class CurveCommand (Command):
84 """A :class:`~hooke.command.Command` operating on a
85 :class:`~hooke.curve.Curve`.
87 def __init__(self, **kwargs):
88 if 'arguments' in kwargs:
89 kwargs['arguments'].insert(0, CurveArgument)
91 kwargs['arguments'] = [CurveArgument]
92 super(CurveCommand, self).__init__(**kwargs)
94 def _curve(self, hooke, params):
95 """Get the selected curve.
99 `hooke` is intended to attach the selected curve to the local
100 playlist; the returned curve should not be effected by the
101 state of `hooke`. This is important for reliable
102 :class:`~hooke.command_stack.CommandStack`\s.
104 # HACK? rely on params['curve'] being bound to the local hooke
105 # playlist (i.e. not a copy, as you would get by passing a
106 # curve through the queue). Ugh. Stupid queues. As an
107 # alternative, we could pass lookup information through the
109 return params['curve']
111 def _add_to_command_stack(self, params):
112 """Store the command name and current `params` values in the
113 curve's `.command_stack`.
115 If this would duplicate the command currently on top of the
116 stack, no action is taken. Call early on, or watch out for
117 repeated param processing.
119 Recommended practice is to *not* lock in argument values that
120 are loaded from the plugin's :attr:`.config`.
124 Perhaps we should subclass :meth:`_run` and use :func:`super`,
125 or embed this in :meth:`run` to avoid subclasses calling this
126 method explicitly, with all the tedium and brittality that
127 implies. On the other hand, the current implemtnation allows
128 CurveCommands that don't effect the curve itself
129 (e.g. :class:`GetCommand`) to avoid adding themselves to the
132 if params['stack'] == True:
133 curve = self._curve(hooke=None, params=params)
134 if (len(curve.command_stack) > 0
135 and curve.command_stack[-1].command == self.name
136 and curve.command_stack[-1].arguments == params):
137 pass # no need to place duplicate calls on the stack.
139 curve.command_stack.append(CommandMessage(
143 class BlockCommand (CurveCommand):
144 """A :class:`CurveCommand` operating on a :class:`~hooke.curve.Data` block.
146 def __init__(self, blocks=None, **kwargs):
148 blocks = [('block', None, 'Name of the data block to act on.')]
150 for name,default,help in blocks:
151 block_args.append(block_argument(name, default, help))
152 self._block_arguments = block_args
153 if 'arguments' not in kwargs:
154 kwargs['arguments'] = []
155 kwargs['arguments'] = block_args + kwargs['arguments']
156 super(BlockCommand, self).__init__(**kwargs)
158 def _block_names(self, hooke, params):
159 curve = self._curve(hooke, params)
160 return [b.info['name'] for b in curve.data]
162 def _block_index(self, hooke, params, name=None):
164 name = self._block_arguments[0].name
165 block_name = params[name]
166 if block_name == None:
167 curve = self._curve(hooke=hooke, params=params)
168 if len(curve.data) == 0:
169 raise Failure('no blocks in %s' % curve)
170 block_name = curve.data[0].info['name']
171 names = self._block_names(hooke=hooke, params=params)
173 return names.index(block_name)
174 except ValueError, e:
175 curve = self._curve(hooke, params)
176 raise Failure('no block named %s in %s (%s): %s'
177 % (block_name, curve, names, e))
179 def _block(self, hooke, params, name=None):
180 # HACK? rely on params['block'] being bound to the local hooke
181 # playlist (i.e. not a copy, as you would get by passing a
182 # block through the queue). Ugh. Stupid queues. As an
183 # alternative, we could pass lookup information through the
185 curve = self._curve(hooke, params)
186 index = self._block_index(hooke, params, name)
187 return curve.data[index]
190 class ColumnAccessCommand (BlockCommand):
191 """A :class:`BlockCommand` accessing a :class:`~hooke.curve.Data`
194 def __init__(self, columns=None, **kwargs):
196 columns = [('column', None, 'Name of the data column to act on.')]
198 for name,default,help in columns:
199 column_args.append(column_argument(name, default, help))
200 self._column_arguments = column_args
201 if 'arguments' not in kwargs:
202 kwargs['arguments'] = []
203 kwargs['arguments'] = column_args + kwargs['arguments']
204 super(ColumnAccessCommand, self).__init__(**kwargs)
206 def _get_column(self, hooke, params, block_name=None, column_name=None):
207 if column_name == None:
208 column_name = self._column_arguments[0].name
209 column_name = params[column_name]
210 block = self._block(hooke, params, block_name)
211 columns = block.info['columns']
213 column_index = columns.index(column_name)
214 except ValueError, e:
215 raise Failure('%s not in %s (%s): %s'
216 % (column_name, block.info['name'], columns, e))
217 return block[:,column_index]
220 class ColumnAddingCommand (ColumnAccessCommand):
221 """A :class:`ColumnAccessCommand` that also adds columns.
223 def __init__(self, new_columns=None, **kwargs):
224 if new_columns == None:
227 for name,default,help in new_columns:
228 column_args.append(column_argument(name, default, help))
229 self._new_column_arguments = column_args
230 if 'arguments' not in kwargs:
231 kwargs['arguments'] = []
232 kwargs['arguments'] = column_args + kwargs['arguments']
233 super(ColumnAddingCommand, self).__init__(**kwargs)
235 def _get_column(self, hooke, params, block_name=None, column_name=None):
236 if column_name == None and len(self._column_arguments) == 0:
237 column_name = self._new_column_arguments[0].name
238 return super(ColumnAddingCommand, self)._get_column(
239 hooke=hooke, params=params, block_name=block_name,
240 column_name=column_name)
242 def _set_column(self, hooke, params, block_name=None, column_name=None,
244 if column_name == None:
245 column_name = self._column_arguments[0].name
246 column_name = params[column_name]
247 block = self._block(hooke=hooke, params=params, name=block_name)
248 if column_name not in block.info['columns']:
249 new = Data((block.shape[0], block.shape[1]+1), dtype=block.dtype)
250 new.info = copy.deepcopy(block.info)
252 new.info['columns'].append(column_name)
254 block_index = self._block_index(hooke, params, name=block_name)
255 self._curve(hooke, params).data[block_index] = block
256 column_index = block.info['columns'].index(column_name)
257 block[:,column_index] = values
262 class CurvePlugin (Builtin):
264 super(CurvePlugin, self).__init__(name='curve')
266 GetCommand(self), InfoCommand(self), DeltaCommand(self),
267 ExportCommand(self), DifferenceCommand(self),
268 DerivativeCommand(self), PowerSpectrumCommand(self),
269 ClearStackCommand(self)]
274 class GetCommand (CurveCommand):
275 """Return a :class:`hooke.curve.Curve`.
277 def __init__(self, plugin):
278 super(GetCommand, self).__init__(
279 name='get curve', help=self.__doc__, plugin=plugin)
281 def _run(self, hooke, inqueue, outqueue, params):
282 outqueue.put(self._curve(hooke, params))
285 class InfoCommand (CurveCommand):
286 """Get selected information about a :class:`hooke.curve.Curve`.
288 def __init__(self, plugin):
290 Argument(name='all', type='bool', default=False, count=1,
291 help='Get all curve information.'),
293 self.fields = ['name', 'path', 'experiment', 'driver', 'filetype',
294 'note', 'command stack', 'blocks', 'block sizes']
295 for field in self.fields:
296 args.append(Argument(
297 name=field, type='bool', default=False, count=1,
298 help='Get curve %s' % field))
299 super(InfoCommand, self).__init__(
300 name='curve info', arguments=args,
301 help=self.__doc__, plugin=plugin)
303 def _run(self, hooke, inqueue, outqueue, params):
304 curve = self._curve(hooke, params)
306 for key in self.fields:
307 fields[key] = params[key]
308 if reduce(lambda x,y: x or y, fields.values()) == False:
309 params['all'] = True # No specific fields set, default to 'all'
310 if params['all'] == True:
311 for key in self.fields:
314 for key in self.fields:
315 if fields[key] == True:
316 get = getattr(self, '_get_%s' % key.replace(' ', '_'))
317 lines.append('%s: %s' % (key, get(curve)))
318 outqueue.put('\n'.join(lines))
320 def _get_name(self, curve):
323 def _get_path(self, curve):
326 def _get_experiment(self, curve):
327 return curve.info.get('experiment', None)
329 def _get_driver(self, curve):
332 def _get_filetype(self, curve):
333 return curve.info.get('filetype', None)
335 def _get_note(self, curve):
336 return curve.info.get('note', None)
338 def _get_command_stack(self, curve):
339 return curve.command_stack
341 def _get_blocks(self, curve):
342 return len(curve.data)
344 def _get_block_sizes(self, curve):
345 return [block.shape for block in curve.data]
348 class DeltaCommand (BlockCommand):
349 """Get distance information between two points.
351 With two points A and B, the returned distances are A-B.
353 def __init__(self, plugin):
354 super(DeltaCommand, self).__init__(
357 Argument(name='point', type='point', optional=False, count=2,
359 Indicies of points bounding the selected data.
361 Argument(name='SI', type='bool', default=False,
363 Return distances in SI notation.
366 help=self.__doc__, plugin=plugin)
368 def _run(self, hooke, inqueue, outqueue, params):
369 data = self._block(hooke, params)
370 As = data[params['point'][0],:]
371 Bs = data[params['point'][1],:]
372 ds = [A-B for A,B in zip(As, Bs)]
373 if params['SI'] == False:
374 out = [(name, d) for name,d in zip(data.info['columns'], ds)]
377 for name,d in zip(data.info['columns'], ds):
378 n,units = split_data_label(name)
380 (n, ppSI(value=d, unit=units, decimals=2)))
384 class ExportCommand (BlockCommand):
385 """Export a :class:`hooke.curve.Curve` data block as TAB-delimeted
388 A "#" prefixed header will optionally appear at the beginning of
389 the file naming the columns.
391 def __init__(self, plugin):
392 super(ExportCommand, self).__init__(
395 Argument(name='output', type='file', default='curve.dat',
397 File name for the output data. Defaults to 'curve.dat'
399 Argument(name='header', type='bool', default=True,
401 True if you want the column-naming header line.
404 help=self.__doc__, plugin=plugin)
406 def _run(self, hooke, inqueue, outqueue, params):
407 data = self._block(hooke, params)
409 with open(params['output'], 'w') as f:
410 if params['header'] == True:
411 f.write('# %s \n' % ('\t'.join(data.info['columns'])))
412 numpy.savetxt(f, data, delimiter='\t')
415 class DifferenceCommand (ColumnAddingCommand):
416 """Calculate the difference between two columns of data.
418 The difference is added to block A as a new column.
420 Note that the command will fail if the columns have different
421 lengths, so be careful when differencing columns from different
424 def __init__(self, plugin):
425 super(DifferenceCommand, self).__init__(
429 'Name of block A in A-B. Defaults to the first block'),
431 'Name of block B in A-B. Defaults to matching `block A`.'),
436 Column of data from block A to difference. Defaults to the first column.
440 Column of data from block B to difference. Defaults to matching `column A`.
444 ('output column', None,
446 Name of the new column for storing the difference (without units, defaults to
447 `difference of <block A> <column A> and <block B> <column B>`).
450 help=self.__doc__, plugin=plugin)
452 def _run(self, hooke, inqueue, outqueue, params):
453 self._add_to_command_stack(params)
454 params = self._setup_params(hooke=hooke, params=params)
455 data_A = self._get_column(hooke=hooke, params=params,
456 block_name='block A',
457 column_name='column A')
458 data_B = self._get_column(hooke=hooke, params=params,
459 block_name='block B',
460 column_name='column B')
461 out = data_A - data_B
462 self._set_column(hooke=hooke, params=params,
463 block_name='block A',
464 column_name='output column',
467 def _setup_params(self, hooke, params):
468 curve = self._curve(hooke, params)
469 if params['block A'] == None:
470 params['block A'] = curve.data[0].info['name']
471 if params['block B'] == None:
472 params['block B'] = params['block A']
473 block_A = self._block(hooke, params=params, name='block A')
474 block_B = self._block(hooke, params=params, name='block B')
475 if params['column A'] == None:
476 params['column A'] = block.info['columns'][0]
477 if params['column B'] == None:
478 params['column B'] = params['column A']
479 a_name,a_unit = split_data_label(params['column A'])
480 b_name,b_unit = split_data_label(params['column B'])
482 raise Failure('Unit missmatch: %s != %s' % (a_unit, b_unit))
483 if params['output column'] == None:
484 params['output column'] = join_data_label(
485 'difference of %s %s and %s %s' % (
486 block_A.info['name'], a_name,
487 block_B.info['name'], b_name),
490 params['output column'] = join_data_label(
491 params['output column'], a_unit)
495 class DerivativeCommand (ColumnAddingCommand):
496 """Calculate the derivative (actually, the discrete differentiation)
499 See :func:`hooke.util.calculus.derivative` for implementation
502 def __init__(self, plugin):
503 super(DerivativeCommand, self).__init__(
507 'Column of data block to differentiate with respect to.'),
509 'Column of data block to differentiate.'),
512 ('output column', None,
514 Name of the new column for storing the derivative (without units, defaults to
515 `derivative of <f column> with respect to <x column>`).
519 Argument(name='weights', type='dict', default={-1:-0.5, 1:0.5},
521 Weighting scheme dictionary for finite differencing. Defaults to
522 central differencing.
525 help=self.__doc__, plugin=plugin)
527 def _run(self, hooke, inqueue, outqueue, params):
528 self._add_to_command_stack(params)
529 params = self._setup_params(hooke=hooke, params=params)
530 x_data = self._get_column(hooke=hooke, params=params,
531 column_name='x column')
532 f_data = self._get_column(hooke=hooke, params=params,
533 column_name='f column')
535 x_data=x_data, f_data=f_data, weights=params['weights'])
536 self._set_column(hooke=hooke, params=params,
537 column_name='output column',
540 def _setup_params(self, hooke, params):
541 curve = self._curve(hooke, params)
542 x_name,x_unit = split_data_label(params['x column'])
543 f_name,f_unit = split_data_label(params['f column'])
544 d_unit = '%s/%s' % (f_unit, x_unit)
545 if params['output column'] == None:
546 params['output column'] = join_data_label(
547 'derivative of %s with respect to %s' % (
551 params['output column'] = join_data_label(
552 params['output column'], d_unit)
556 class PowerSpectrumCommand (ColumnAddingCommand):
557 """Calculate the power spectrum of a data column.
559 def __init__(self, plugin):
560 super(PowerSpectrumCommand, self).__init__(
561 name='power spectrum',
563 Argument(name='output block', type='string',
565 Name of the new data block for storing the power spectrum (defaults to
566 `power spectrum of <source block name> <source column name>`).
568 Argument(name='bounds', type='point', optional=True, count=2,
570 Indicies of points bounding the selected data.
572 Argument(name='freq', type='float', default=1.0,
576 Argument(name='freq units', type='string', default='Hz',
578 Units for the sampling frequency.
580 Argument(name='chunk size', type='int', default=2048,
582 Number of samples per chunk. Use a power of two.
584 Argument(name='overlap', type='bool', default=False,
586 If `True`, each chunk overlaps the previous chunk by half its length.
587 Otherwise, the chunks are end-to-end, and not overlapping.
590 help=self.__doc__, plugin=plugin)
592 def _run(self, hooke, inqueue, outqueue, params):
593 self._add_to_command_stack(params)
594 params = self._setup_params(hooke=hooke, params=params)
595 data = self._get_column(hooke=hooke, params=params)
596 bounds = params['bounds']
598 data = data[bounds[0]:bounds[1]]
599 freq_axis,power = unitary_avg_power_spectrum(
600 data, freq=params['freq'],
601 chunk_size=params['chunk size'],
602 overlap=params['overlap'])
603 b = Data((len(freq_axis),2), dtype=data.dtype)
604 b.info['name'] = params['output block']
605 b.info['columns'] = [
606 params['output freq column'],
607 params['output power column'],
609 self._curve(hooke, params).data.append(b)
610 self._set_column(hooke, params, block_name='output block',
611 column_name='output freq column',
613 self._set_column(hooke, params, block_name='output block',
614 column_name='output power column',
618 def _setup_params(self, hooke, params):
619 if params['output block'] in self._block_names(hooke, params):
620 raise Failure('output block %s already exists in %s.'
621 % (params['output block'],
622 self._curve(hooke, params)))
623 data = self._get_column(hooke=hooke, params=params)
624 d_name,d_unit = split_data_label(data.info['name'])
625 if params['output block'] == None:
626 params['output block'] = 'power spectrum of %s %s' % (
627 data.info['name'], params['column'])
628 self.params['output freq column'] = join_data_label(
629 'frequency axis', params['freq units'])
630 self.params['output power column'] = join_data_label(
631 'power density', '%s^2/%s' % (data_units, params['freq units']))
635 class ClearStackCommand (CurveCommand):
636 """Empty a curve's command stack.
638 def __init__(self, plugin):
639 super(ClearStackCommand, self).__init__(
640 name='clear curve command stack',
641 help=self.__doc__, plugin=plugin)
642 i,arg = [(i,arg) for i,arg in enumerate(self.arguments)
643 if arg.name == 'curve'][0]
645 arg.callback = unloaded_current_curve_callback
646 self.arguments[i] = arg
648 def _run(self, hooke, inqueue, outqueue, params):
649 curve = self._curve(hooke, params)
650 curve.command_stack = CommandStack()
653 class OldCruft (object):
655 def do_forcebase(self,args):
659 Measures the difference in force (in pN) between a point and a baseline
660 took as the average between two points.
662 The baseline is fixed once for a given curve and different force measurements,
663 unless the user wants it to be recalculated
665 Syntax: forcebase [rebase]
666 rebase: Forces forcebase to ask again the baseline
667 max: Instead of asking for a point to measure, asks for two points and use
668 the maximum peak in between
670 rebase=False #if true=we select rebase
671 maxpoint=False #if true=we measure the maximum peak
673 plot=self._get_displayed_plot()
674 whatset=1 #fixme: for all sets
675 if 'rebase' in args or (self.basecurrent != self.current.path):
681 print 'Select baseline'
682 self.basepoints=self._measure_N_points(N=2, whatset=whatset)
683 self.basecurrent=self.current.path
686 print 'Select two points'
687 points=self._measure_N_points(N=2, whatset=whatset)
688 boundpoints=[points[0].index, points[1].index]
691 y=min(plot.vectors[whatset][1][boundpoints[0]:boundpoints[1]])
693 print 'Chosen interval not valid. Try picking it again. Did you pick the same point as begin and end of interval?'
695 print 'Select point to measure'
696 points=self._measure_N_points(N=1, whatset=whatset)
697 #whatplot=points[0].dest
698 y=points[0].graph_coords[1]
700 #fixme: code duplication
701 boundaries=[self.basepoints[0].index, self.basepoints[1].index]
703 to_average=plot.vectors[whatset][1][boundaries[0]:boundaries[1]] #y points to average
705 avg=np.mean(to_average)
707 print str(forcebase*(10**12))+' pN'
708 to_dump='forcebase '+self.current.path+' '+str(forcebase*(10**12))+' pN'
709 self.outlet.push(to_dump)
712 def do_slope(self,args):
716 Measures the slope of a delimited chunk on the return trace.
717 The chunk can be delimited either by two manual clicks, or have
718 a fixed width, given as an argument.
720 Syntax: slope [width]
721 The facultative [width] parameter specifies how many
722 points will be considered for the fit. If [width] is
723 specified, only one click will be required.
724 (c) Marco Brucale, Massimo Sandal 2008
727 # Reads the facultative width argument
733 # Decides between the two forms of user input, as per (args)
735 # Gets the Xs of two clicked points as indexes on the current curve vector
736 print 'Click twice to delimit chunk'
737 points=self._measure_N_points(N=2,whatset=1)
739 print 'Click once on the leftmost point of the chunk (i.e.usually the peak)'
740 points=self._measure_N_points(N=1,whatset=1)
742 slope=self._slope(points,fitspan)
744 # Outputs the relevant slope parameter
747 to_dump='slope '+self.current.path+' '+str(slope)
748 self.outlet.push(to_dump)
750 def _slope(self,points,fitspan):
751 # Calls the function linefit_between
752 parameters=[0,0,[],[]]
754 clickedpoints=[points[0].index,points[1].index]
757 clickedpoints=[points[0].index-fitspan,points[0].index]
760 parameters=self.linefit_between(clickedpoints[0],clickedpoints[1])
762 print 'Cannot fit. Did you click twice the same point?'
765 # Outputs the relevant slope parameter
767 print str(parameters[0])
768 to_dump='slope '+self.curve.path+' '+str(parameters[0])
769 self.outlet.push(to_dump)
771 # Makes a vector with the fitted parameters and sends it to the GUI
772 xtoplot=parameters[2]
776 ytoplot.append((x*parameters[0])+parameters[1])
778 clickvector_x, clickvector_y=[], []
780 clickvector_x.append(item.graph_coords[0])
781 clickvector_y.append(item.graph_coords[1])
783 lineplot=self._get_displayed_plot(0) #get topmost displayed plot
785 lineplot.add_set(xtoplot,ytoplot)
786 lineplot.add_set(clickvector_x, clickvector_y)
789 if lineplot.styles==[]:
790 lineplot.styles=[None,None,None,'scatter']
792 lineplot.styles+=[None,'scatter']
793 if lineplot.colors==[]:
794 lineplot.colors=[None,None,'black',None]
796 lineplot.colors+=['black',None]
799 self._send_plot([lineplot])
804 def linefit_between(self,index1,index2,whatset=1):
806 Creates two vectors (xtofit,ytofit) slicing out from the
807 current return trace a portion delimited by the two indexes
809 Then does a least squares linear fit on that slice.
810 Finally returns [0]=the slope, [1]=the intercept of the
811 fitted 1st grade polynomial, and [2,3]=the actual (x,y) vectors
813 (c) Marco Brucale, Massimo Sandal 2008
815 # Translates the indexes into two vectors containing the x,y data to fit
816 xtofit=self.plots[0].vectors[whatset][0][index1:index2]
817 ytofit=self.plots[0].vectors[whatset][1][index1:index2]
819 # Does the actual linear fitting (simple least squares with numpy.polyfit)
821 linefit=np.polyfit(xtofit,ytofit,1)
823 return (linefit[0],linefit[1],xtofit,ytofit)