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):
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',
289 'note', 'command stack', 'blocks', 'block sizes']
290 for field in self.fields:
291 args.append(Argument(
292 name=field, type='bool', default=False, count=1,
293 help='Get curve %s' % field))
294 super(InfoCommand, self).__init__(
295 name='curve info', arguments=args,
296 help=self.__doc__, plugin=plugin)
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 or y, fields.values()) == False:
304 params['all'] = True # No specific fields set, default to 'all'
305 if params['all'] == True:
306 for key in self.fields:
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_command_stack(self, curve):
334 return curve.command_stack
336 def _get_blocks(self, curve):
337 return len(curve.data)
339 def _get_block_sizes(self, curve):
340 return [block.shape for block in curve.data]
343 class DeltaCommand (BlockCommand):
344 """Get distance information between two points.
346 With two points A and B, the returned distances are A-B.
348 def __init__(self, plugin):
349 super(DeltaCommand, self).__init__(
352 Argument(name='point', type='point', optional=False, count=2,
354 Indicies of points bounding the selected data.
356 Argument(name='SI', type='bool', default=False,
358 Return distances in SI notation.
361 help=self.__doc__, plugin=plugin)
363 def _run(self, hooke, inqueue, outqueue, params):
364 data = self._block(hooke, params)
365 As = data[params['point'][0],:]
366 Bs = data[params['point'][1],:]
367 ds = [A-B for A,B in zip(As, Bs)]
368 if params['SI'] == False:
369 out = [(name, d) for name,d in zip(data.info['columns'], ds)]
372 for name,d in zip(data.info['columns'], ds):
373 n,units = split_data_label(name)
375 (n, ppSI(value=d, unit=units, decimals=2)))
379 class ExportCommand (BlockCommand):
380 """Export a :class:`hooke.curve.Curve` data block as TAB-delimeted
383 A "#" prefixed header will optionally appear at the beginning of
384 the file naming the columns.
386 def __init__(self, plugin):
387 super(ExportCommand, self).__init__(
390 Argument(name='output', type='file', default='curve.dat',
392 File name for the output data. Defaults to 'curve.dat'
394 Argument(name='header', type='bool', default=True,
396 True if you want the column-naming header line.
399 help=self.__doc__, plugin=plugin)
401 def _run(self, hooke, inqueue, outqueue, params):
402 data = self._block(hooke, params)
404 with open(params['output'], 'w') as f:
405 if params['header'] == True:
406 f.write('# %s \n' % ('\t'.join(data.info['columns'])))
407 numpy.savetxt(f, data, delimiter='\t')
410 class DifferenceCommand (ColumnAddingCommand):
411 """Calculate the difference between two columns of data.
413 The difference is added to block A as a new column.
415 Note that the command will fail if the columns have different
416 lengths, so be careful when differencing columns from different
419 def __init__(self, plugin):
420 super(DifferenceCommand, self).__init__(
424 'Name of block A in A-B. Defaults to the first block'),
426 'Name of block B in A-B. Defaults to matching `block A`.'),
431 Column of data from block A to difference. Defaults to the first column.
435 Column of data from block B to difference. Defaults to matching `column A`.
439 ('output column', None,
441 Name of the new column for storing the difference (without units, defaults to
442 `difference of <block A> <column A> and <block B> <column B>`).
445 help=self.__doc__, plugin=plugin)
447 def _run(self, hooke, inqueue, outqueue, params):
448 self._add_to_command_stack(params)
449 params = self.__setup_params(hooke=hooke, params=params)
450 data_A = self._get_column(hooke=hooke, params=params,
451 block_name='block A',
452 column_name='column A')
453 data_B = self._get_column(hooke=hooke, params=params,
454 block_name='block B',
455 column_name='column B')
456 out = data_A - data_B
457 self._set_column(hooke=hooke, params=params,
458 block_name='block A',
459 column_name='output column',
462 def __setup_params(self, hooke, params):
463 curve = self._curve(hooke, params)
464 if params['block A'] == None:
465 params['block A'] = curve.data[0].info['name']
466 if params['block B'] == None:
467 params['block B'] = params['block A']
468 block_A = self._block(hooke, params=params, name='block A')
469 block_B = self._block(hooke, params=params, name='block B')
470 if params['column A'] == None:
471 params['column A'] = block.info['columns'][0]
472 if params['column B'] == None:
473 params['column B'] = params['column A']
474 a_name,a_unit = split_data_label(params['column A'])
475 b_name,b_unit = split_data_label(params['column B'])
477 raise Failure('Unit missmatch: %s != %s' % (a_unit, b_unit))
478 if params['output column'] == None:
479 params['output column'] = join_data_label(
480 'difference of %s %s and %s %s' % (
481 block_A.info['name'], a_name,
482 block_B.info['name'], b_name),
485 params['output column'] = join_data_label(
486 params['output column'], a_unit)
490 class DerivativeCommand (ColumnAddingCommand):
491 """Calculate the derivative (actually, the discrete differentiation)
494 See :func:`hooke.util.calculus.derivative` for implementation
497 def __init__(self, plugin):
498 super(DerivativeCommand, self).__init__(
502 'Column of data block to differentiate with respect to.'),
504 'Column of data block to differentiate.'),
507 ('output column', None,
509 Name of the new column for storing the derivative (without units, defaults to
510 `derivative of <f column> with respect to <x column>`).
514 Argument(name='weights', type='dict', default={-1:-0.5, 1:0.5},
516 Weighting scheme dictionary for finite differencing. Defaults to
517 central differencing.
520 help=self.__doc__, plugin=plugin)
522 def _run(self, hooke, inqueue, outqueue, params):
523 self._add_to_command_stack(params)
524 params = self.__setup_params(hooke=hooke, params=params)
525 x_data = self._get_column(hooke=hooke, params=params,
526 column_name='x column')
527 f_data = self._get_column(hooke=hooke, params=params,
528 column_name='f column')
530 x_data=x_data, f_data=f_data, weights=params['weights'])
531 self._set_column(hooke=hooke, params=params,
532 column_name='output column',
535 def __setup_params(self, hooke, params):
536 curve = self._curve(hooke, params)
537 x_name,x_unit = split_data_label(params['x column'])
538 f_name,f_unit = split_data_label(params['f column'])
539 d_unit = '%s/%s' % (f_unit, x_unit)
540 if params['output column'] == None:
541 params['output column'] = join_data_label(
542 'derivative of %s with respect to %s' % (
546 params['output column'] = join_data_label(
547 params['output column'], d_unit)
551 class PowerSpectrumCommand (ColumnAddingCommand):
552 """Calculate the power spectrum of a data column.
554 def __init__(self, plugin):
555 super(PowerSpectrumCommand, self).__init__(
556 name='power spectrum',
558 Argument(name='output block', type='string',
560 Name of the new data block for storing the power spectrum (defaults to
561 `power spectrum of <source block name> <source column name>`).
563 Argument(name='bounds', type='point', optional=True, count=2,
565 Indicies of points bounding the selected data.
567 Argument(name='freq', type='float', default=1.0,
571 Argument(name='freq units', type='string', default='Hz',
573 Units for the sampling frequency.
575 Argument(name='chunk size', type='int', default=2048,
577 Number of samples per chunk. Use a power of two.
579 Argument(name='overlap', type='bool', default=False,
581 If `True`, each chunk overlaps the previous chunk by half its length.
582 Otherwise, the chunks are end-to-end, and not overlapping.
585 help=self.__doc__, plugin=plugin)
587 def _run(self, hooke, inqueue, outqueue, params):
588 self._add_to_command_stack(params)
589 params = self.__setup_params(hooke=hooke, params=params)
590 data = self._get_column(hooke=hooke, params=params)
591 bounds = params['bounds']
593 data = data[bounds[0]:bounds[1]]
594 freq_axis,power = unitary_avg_power_spectrum(
595 data, freq=params['freq'],
596 chunk_size=params['chunk size'],
597 overlap=params['overlap'])
598 b = Data((len(freq_axis),2), dtype=data.dtype)
599 b.info['name'] = params['output block']
600 b.info['columns'] = [
601 params['output freq column'],
602 params['output power column'],
604 self._curve(hooke, params).data.append(b)
605 self._set_column(hooke, params, block_name='output block',
606 column_name='output freq column',
608 self._set_column(hooke, params, block_name='output block',
609 column_name='output power column',
613 def __setup_params(self, hooke, params):
614 if params['output block'] in self._block_names(hooke, params):
615 raise Failure('output block %s already exists in %s.'
616 % (params['output block'],
617 self._curve(hooke, params)))
618 data = self._get_column(hooke=hooke, params=params)
619 d_name,d_unit = split_data_label(data.info['name'])
620 if params['output block'] == None:
621 params['output block'] = 'power spectrum of %s %s' % (
622 data.info['name'], params['column'])
623 self.params['output freq column'] = join_data_label(
624 'frequency axis', params['freq units'])
625 self.params['output power column'] = join_data_label(
626 'power density', '%s^2/%s' % (data_units, params['freq units']))
630 class ClearStackCommand (CurveCommand):
631 """Empty a curve's command stack.
633 def __init__(self, plugin):
634 super(ClearStackCommand, self).__init__(
635 name='clear curve command stack',
636 help=self.__doc__, plugin=plugin)
638 def _run(self, hooke, inqueue, outqueue, params):
639 curve = self._curve(hooke, params)
640 curve.command_stack = CommandStack()
643 class OldCruft (object):
645 def do_forcebase(self,args):
649 Measures the difference in force (in pN) between a point and a baseline
650 took as the average between two points.
652 The baseline is fixed once for a given curve and different force measurements,
653 unless the user wants it to be recalculated
655 Syntax: forcebase [rebase]
656 rebase: Forces forcebase to ask again the baseline
657 max: Instead of asking for a point to measure, asks for two points and use
658 the maximum peak in between
660 rebase=False #if true=we select rebase
661 maxpoint=False #if true=we measure the maximum peak
663 plot=self._get_displayed_plot()
664 whatset=1 #fixme: for all sets
665 if 'rebase' in args or (self.basecurrent != self.current.path):
671 print 'Select baseline'
672 self.basepoints=self._measure_N_points(N=2, whatset=whatset)
673 self.basecurrent=self.current.path
676 print 'Select two points'
677 points=self._measure_N_points(N=2, whatset=whatset)
678 boundpoints=[points[0].index, points[1].index]
681 y=min(plot.vectors[whatset][1][boundpoints[0]:boundpoints[1]])
683 print 'Chosen interval not valid. Try picking it again. Did you pick the same point as begin and end of interval?'
685 print 'Select point to measure'
686 points=self._measure_N_points(N=1, whatset=whatset)
687 #whatplot=points[0].dest
688 y=points[0].graph_coords[1]
690 #fixme: code duplication
691 boundaries=[self.basepoints[0].index, self.basepoints[1].index]
693 to_average=plot.vectors[whatset][1][boundaries[0]:boundaries[1]] #y points to average
695 avg=np.mean(to_average)
697 print str(forcebase*(10**12))+' pN'
698 to_dump='forcebase '+self.current.path+' '+str(forcebase*(10**12))+' pN'
699 self.outlet.push(to_dump)
702 def do_slope(self,args):
706 Measures the slope of a delimited chunk on the return trace.
707 The chunk can be delimited either by two manual clicks, or have
708 a fixed width, given as an argument.
710 Syntax: slope [width]
711 The facultative [width] parameter specifies how many
712 points will be considered for the fit. If [width] is
713 specified, only one click will be required.
714 (c) Marco Brucale, Massimo Sandal 2008
717 # Reads the facultative width argument
723 # Decides between the two forms of user input, as per (args)
725 # Gets the Xs of two clicked points as indexes on the current curve vector
726 print 'Click twice to delimit chunk'
727 points=self._measure_N_points(N=2,whatset=1)
729 print 'Click once on the leftmost point of the chunk (i.e.usually the peak)'
730 points=self._measure_N_points(N=1,whatset=1)
732 slope=self._slope(points,fitspan)
734 # Outputs the relevant slope parameter
737 to_dump='slope '+self.current.path+' '+str(slope)
738 self.outlet.push(to_dump)
740 def _slope(self,points,fitspan):
741 # Calls the function linefit_between
742 parameters=[0,0,[],[]]
744 clickedpoints=[points[0].index,points[1].index]
747 clickedpoints=[points[0].index-fitspan,points[0].index]
750 parameters=self.linefit_between(clickedpoints[0],clickedpoints[1])
752 print 'Cannot fit. Did you click twice the same point?'
755 # Outputs the relevant slope parameter
757 print str(parameters[0])
758 to_dump='slope '+self.curve.path+' '+str(parameters[0])
759 self.outlet.push(to_dump)
761 # Makes a vector with the fitted parameters and sends it to the GUI
762 xtoplot=parameters[2]
766 ytoplot.append((x*parameters[0])+parameters[1])
768 clickvector_x, clickvector_y=[], []
770 clickvector_x.append(item.graph_coords[0])
771 clickvector_y.append(item.graph_coords[1])
773 lineplot=self._get_displayed_plot(0) #get topmost displayed plot
775 lineplot.add_set(xtoplot,ytoplot)
776 lineplot.add_set(clickvector_x, clickvector_y)
779 if lineplot.styles==[]:
780 lineplot.styles=[None,None,None,'scatter']
782 lineplot.styles+=[None,'scatter']
783 if lineplot.colors==[]:
784 lineplot.colors=[None,None,'black',None]
786 lineplot.colors+=['black',None]
789 self._send_plot([lineplot])
794 def linefit_between(self,index1,index2,whatset=1):
796 Creates two vectors (xtofit,ytofit) slicing out from the
797 current return trace a portion delimited by the two indexes
799 Then does a least squares linear fit on that slice.
800 Finally returns [0]=the slope, [1]=the intercept of the
801 fitted 1st grade polynomial, and [2,3]=the actual (x,y) vectors
803 (c) Marco Brucale, Massimo Sandal 2008
805 # Translates the indexes into two vectors containing the x,y data to fit
806 xtofit=self.plots[0].vectors[whatset][0][index1:index2]
807 ytofit=self.plots[0].vectors[whatset][1][index1:index2]
809 # Does the actual linear fitting (simple least squares with numpy.polyfit)
811 linefit=np.polyfit(xtofit,ytofit,1)
813 return (linefit[0],linefit[1],xtofit,ytofit)