Update DifferenceCommand to use flexible column names.
[hooke.git] / hooke / plugin / curve.py
index 1e0ad061687914fb1bd6576cee485419b8ca4d73..86ea7d8099203e6c2428f5ad565a7442fbd43abf 100644 (file)
@@ -1,17 +1,19 @@
-# Copyright (C) 2010 Fibrin's Benedetti
-#                    W. Trevor King <wking@drexel.edu>
+# Copyright (C) 2008-2010 Alberto Gomez-Casado
+#                         Fabrizio Benedetti
+#                         Massimo Sandal <devicerandom@gmail.com>
+#                         W. Trevor King <wking@drexel.edu>
 #
 # This file is part of Hooke.
 #
-# Hooke is free software: you can redistribute it and/or
-# modify it under the terms of the GNU Lesser General Public
-# License as published by the Free Software Foundation, either
-# version 3 of the License, or (at your option) any later version.
+# Hooke is free software: you can redistribute it and/or modify it
+# under the terms of the GNU Lesser General Public License as
+# published by the Free Software Foundation, either version 3 of the
+# License, or (at your option) any later version.
 #
-# Hooke is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-# GNU Lesser General Public License for more details.
+# Hooke is distributed in the hope that it will be useful, but WITHOUT
+# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
+# or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General
+# Public License for more details.
 #
 # You should have received a copy of the GNU Lesser General Public
 # License along with Hooke.  If not, see
@@ -22,17 +24,27 @@ associated :class:`hooke.command.Command`\s for handling
 :mod:`hooke.curve` classes.
 """
 
+import copy
+
+import numpy
+
 from ..command import Command, Argument, Failure
+from ..curve import Data
 from ..plugin import Builtin
 from ..plugin.playlist import current_playlist_callback
+from ..util.calculus import derivative
+from ..util.fft import unitary_avg_power_spectrum
+from ..util.si import ppSI, join_data_label, split_data_label
+
 
 
 class CurvePlugin (Builtin):
     def __init__(self):
         super(CurvePlugin, self).__init__(name='curve')
-
-    def commands(self):
-        return [InfoCommand(), ]
+        self._commands = [
+            GetCommand(self), InfoCommand(self), DeltaCommand(self),
+            ExportCommand(self), DifferenceCommand(self),
+            DerivativeCommand(self), PowerSpectrumCommand(self)]
 
 
 # Define common or complicated arguments
@@ -56,23 +68,35 @@ of the current playlist.
 
 # Define commands
 
+class GetCommand (Command):
+    """Return a :class:`hooke.curve.Curve`.
+    """
+    def __init__(self, plugin):
+        super(GetCommand, self).__init__(
+            name='get curve', arguments=[CurveArgument],
+            help=self.__doc__, plugin=plugin)
+
+    def _run(self, hooke, inqueue, outqueue, params):
+        outqueue.put(params['curve'])
+
 class InfoCommand (Command):
-    """Print selected information about a :class:`hooke.curve.Curve`.
+    """Get selected information about a :class:`hooke.curve.Curve`.
     """
-    def __init__(self):
+    def __init__(self, plugin):
         args = [
             CurveArgument,                    
             Argument(name='all', type='bool', default=False, count=1,
-                     help='Print all curve information.'),
+                     help='Get all curve information.'),
             ]
         self.fields = ['name', 'path', 'experiment', 'driver', 'filetype', 'note',
                        'blocks', 'block sizes']
         for field in self.fields:
             args.append(Argument(
                     name=field, type='bool', default=False, count=1,
-                    help='Print curve %s' % field))
+                    help='Get curve %s' % field))
         super(InfoCommand, self).__init__(
-            name='curve info', arguments=args, help=self.__doc__)
+            name='curve info', arguments=args,
+            help=self.__doc__, plugin=plugin)
 
     def _run(self, hooke, inqueue, outqueue, params):
         fields = {}
@@ -115,35 +139,294 @@ class InfoCommand (Command):
         return [block.shape for block in curve.data]
 
 
+class DeltaCommand (Command):
+    """Get distance information between two points.
+
+    With two points A and B, the returned distances are A-B.
+    """
+    def __init__(self, plugin):
+        super(DeltaCommand, self).__init__(
+            name='delta',
+            arguments=[
+                CurveArgument,
+                Argument(name='block', type='int', default=0,
+                    help="""
+Data block that points are selected from.  For an approach/retract
+force curve, `0` selects the approaching curve and `1` selects the
+retracting curve.
+""".strip()),
+                Argument(name='point', type='point', optional=False, count=2,
+                         help="""
+Indicies of points bounding the selected data.
+""".strip()),
+                Argument(name='SI', type='bool', default=False,
+                         help="""
+Return distances in SI notation.
+""".strip())
+                ],
+            help=self.__doc__, plugin=plugin)
+
+    def _run(self, hooke, inqueue, outqueue, params):
+        data = params['curve'].data[params['block']]
+        As = data[params['point'][0],:]
+        Bs = data[params['point'][1],:]
+        ds = [A-B for A,B in zip(As, Bs)]
+        if params['SI'] == False:
+            out = [(name, d) for name,d in zip(data.info['columns'], ds)]
+        else:
+            out = []
+            for name,d in zip(data.info['columns'], ds):
+                n,units = split_data_label(name)
+                out.append(
+                  (n, ppSI(value=d, unit=units, decimals=2)))
+        outqueue.put(out)
+
+
 class ExportCommand (Command):
     """Export a :class:`hooke.curve.Curve` data block as TAB-delimeted
     ASCII text.
+
+    A "#" prefixed header will optionally appear at the beginning of
+    the file naming the columns.
     """
-    def __init__(self):
-        self.fields = ['name', 'path', 'experiment', 'driver', 'filetype', 'note',
-                       'blocks', 'block sizes']
-        for field in self.fields:
-            args.append(Argument(
-                    name=field, type='bool', default=False, count=1,
-                    help='Print curve %s' % field))
-        super(InfoCommand, self).__init__(
-            name='curve info',
+    def __init__(self, plugin):
+        super(ExportCommand, self).__init__(
+            name='export block',
             arguments=[
                 CurveArgument,
-                Argument(name='block', aliases=['set'], type='int', default=0,
-                    help="""
+                Argument(name='block', type='int', default=0,
+                         help="""
 Data block to save.  For an approach/retract force curve, `0` selects
-the approacing curve and `1` selects the retracting curve.
+the approaching curve and `1` selects the retracting curve.
 """.strip()),
                 Argument(name='output', type='file', default='curve.dat',
                          help="""
 File name for the output data.  Defaults to 'curve.dat'
+""".strip()),
+                Argument(name='header', type='bool', default=True,
+                         help="""
+True if you want the column-naming header line.
 """.strip()),
                 ],
-            help=self.__doc__)
+            help=self.__doc__, plugin=plugin)
 
     def _run(self, hooke, inqueue, outqueue, params):
-        data = params['curve'].data[params['index']]
+        data = params['curve'].data[params['block']]
+
         f = open(params['output'], 'w')
-        data.tofile(f, sep='\t')
+        if params['header'] == True:
+            f.write('# %s \n' % ('\t'.join(data.info['columns'])))
+        numpy.savetxt(f, data, delimiter='\t')
         f.close()
+
+class DifferenceCommand (Command):
+    """Calculate the difference between two columns of data.
+
+    The difference is added to block A as a new column.
+
+    Note that the command will fail if the columns have different
+    lengths, so be careful when differencing columns from different
+    blocks.
+    """
+    def __init__(self, plugin):
+        super(DifferenceCommand, self).__init__(
+            name='difference',
+            arguments=[
+                CurveArgument,
+                Argument(name='block A', type='int',
+                         help="""
+Block A in A-B.  For an approach/retract force curve, `0` selects the
+approaching curve and `1` selects the retracting curve.  Defaults to
+the first block.
+""".strip()),
+                Argument(name='block B', type='int',
+                         help="""
+Block B in A-B.  Defaults to matching `block A`.
+""".strip()),
+                Argument(name='column A', type='string',
+                         help="""
+Column of data from block A to difference.  Defaults to the first column.
+""".strip()),
+                Argument(name='column B', type='string', default=1,
+                         help="""
+Column of data from block B to difference.  Defaults to matching `column A`.
+""".strip()),
+                Argument(name='output column name', type='string',
+                         help="""
+Name of the new column for storing the difference (without units, defaults to
+`difference of <block A> <column A> and <block B> <column B>`).
+""".strip()),
+                ],
+            help=self.__doc__, plugin=plugin)
+
+    def _run(self, hooke, inqueue, outqueue, params):
+        data_A = params['curve'].data[params['block A']]
+        data_B = params['curve'].data[params['block B']]
+        # HACK? rely on params['curve'] being bound to the local hooke
+        # playlist (i.e. not a copy, as you would get by passing a
+        # curve through the queue).  Ugh.  Stupid queues.  As an
+        # alternative, we could pass lookup information through the
+        # queue...
+        new = Data((data_A.shape[0], data_A.shape[1]+1), dtype=data_A.dtype)
+        new.info = copy.deepcopy(data.info)
+        new[:,:-1] = data_A
+
+        a_col = data_A.info['columns'].index(params['column A'])
+        b_col = data_A.info['columns'].index(params['column A'])
+        out = data_A[:,a_col] - data_B[:,b_col]
+
+        a_name,a_units = split_data_label(params['column A'])
+        b_name,b_units = split_data_label(params['column B'])
+        assert a_units == b_units, (
+            'Unit missmatch: %s != %s' % (a_units, b_units))
+        if params['output column name'] == None:
+            params['output column name'] = (
+                'difference of %s %s and %s %s' % (
+                    block_A.info['name'], params['column A'],
+                    block_B.info['name'], params['column B']))
+        new.info['columns'].append(
+            join_data_label(params['output distance column'], a_units)
+        new[:,-1] = out
+        params['curve'].data[params['block A']] = new
+
+
+class DerivativeCommand (Command):
+    """Calculate the derivative (actually, the discrete differentiation)
+    of a curve data block.
+
+    See :func:`hooke.util.calculus.derivative` for implementation
+    details.
+    """
+    def __init__(self, plugin):
+        super(DerivativeCommand, self).__init__(
+            name='derivative',
+            arguments=[
+                CurveArgument,
+                Argument(name='block', type='int', default=0,
+                         help="""
+Data block to differentiate.  For an approach/retract force curve, `0`
+selects the approaching curve and `1` selects the retracting curve.
+""".strip()),
+                Argument(name='x column', type='string',
+                         help="""
+Column of data block to differentiate with respect to.
+""".strip()),
+                Argument(name='f column', type='string',
+                         help="""
+Column of data block to differentiate.
+""".strip()),
+                Argument(name='weights', type='dict', default={-1:-0.5, 1:0.5},
+                         help="""
+Weighting scheme dictionary for finite differencing.  Defaults to
+central differencing.
+""".strip()),
+                Argument(name='output column name', type='string',
+                         help="""
+Name of the new column for storing the derivative (without units, defaults to
+`derivative of <f column name> with respect to <x column name>`).
+""".strip()),
+                ],
+            help=self.__doc__, plugin=plugin)
+
+    def _run(self, hooke, inqueue, outqueue, params):
+        data = params['curve'].data[params['block']]
+        # HACK? rely on params['curve'] being bound to the local hooke
+        # playlist (i.e. not a copy, as you would get by passing a
+        # curve through the queue).  Ugh.  Stupid queues.  As an
+        # alternative, we could pass lookup information through the
+        # queue...
+        new = Data((data.shape[0], data.shape[1]+1), dtype=data.dtype)
+        new.info = copy.deepcopy(data.info)
+        new[:,:-1] = data
+
+        x_col = data.info['columns'].index(params['x column'])
+        f_col = data.info['columns'].index(params['f column'])
+        d = derivative(
+            block, x_col=x_col, f_col=f_col, weights=params['weights'])
+
+        x_name,x_units = split_data_label(params['x column'])
+        f_name,f_units = split_data_label(params['f column'])
+        if params['output column name'] == None:
+            params['output column name'] = (
+                'derivative of %s with respect to %s' % (
+                    params['f column'], params['x column']))
+
+        new.info['columns'].append(
+            join_data_label(params['output distance column'],
+                            '%s/%s' % (f_units/x_units)))
+        new[:,-1] = d[:,1]
+        params['curve'].data[params['block']] = new
+
+
+class PowerSpectrumCommand (Command):
+    """Calculate the power spectrum of a data block.
+    """
+    def __init__(self, plugin):
+        super(PowerSpectrumCommand, self).__init__(
+            name='power spectrum',
+            arguments=[
+                CurveArgument,
+                Argument(name='block', type='int', default=0,
+                         help="""
+Data block to act on.  For an approach/retract force curve, `0`
+selects the approaching curve and `1` selects the retracting curve.
+""".strip()),
+                Argument(name='column', type='string', optional=False,
+                         help="""
+Name of the data block column containing to-be-transformed data.
+""".strip()),
+                Argument(name='bounds', type='point', optional=True, count=2,
+                         help="""
+Indicies of points bounding the selected data.
+""".strip()),
+                Argument(name='freq', type='float', default=1.0,
+                         help="""
+Sampling frequency.
+""".strip()),
+                Argument(name='freq units', type='string', default='Hz',
+                         help="""
+Units for the sampling frequency.
+""".strip()),
+                Argument(name='chunk size', type='int', default=2048,
+                         help="""
+Number of samples per chunk.  Use a power of two.
+""".strip()),
+                Argument(name='overlap', type='bool', default=False,
+                         help="""
+If `True`, each chunk overlaps the previous chunk by half its length.
+Otherwise, the chunks are end-to-end, and not overlapping.
+""".strip()),
+                Argument(name='output block name', type='string',
+                         help="""
+Name of the new data block for storing the power spectrum (defaults to
+`power spectrum of <source block name> <source column name>`).
+""".strip()),
+                ],
+            help=self.__doc__, plugin=plugin)
+
+    def _run(self, hooke, inqueue, outqueue, params):
+        data = params['curve'].data[params['block']]
+        col = data.info['columns'].index(params['column'])
+        d = data[:,col]
+        if bounds != None:
+            d = d[params['bounds'][0]:params['bounds'][1]]
+        freq_axis,power = unitary_avg_power_spectrum(
+            d, freq=params['freq'],
+            chunk_size=params['chunk size'],
+            overlap=params['overlap'])
+
+        name,data_units = split_data_label(params['column'])
+        b = Data((len(freq_axis),2), dtype=data.dtype)
+        if params['output block name'] == None:
+            params['output block name'] = 'power spectrum of %s %s' % (
+              params['output block name'], data.info['name'], params['column'])
+        b.info['name'] = params['output block name']
+        b.info['columns'] = [
+            join_data_label('frequency axis', params['freq units']),
+            join_data_label('power density',
+                            '%s^2/%s' % (data_units, params['freq units'])),
+            ]
+        b[:,0] = freq_axis
+        b[:,1] = power
+        params['curve'].data.append(b)
+        outqueue.put(b)