hooke/util/calculus.py

   1 # Copyright
   2
   3 """The `calculus` module provides functions for calculating
   4 derivatives and integrals of discrete data.
   5 """
   6
   7 import copy
   8
   9 import numpy
  10
  11 from ..curve import Data
  12
  13
  14 def derivative(data, x_col=0, f_col=1, weights={-1:-0.5, 1:0.5}):
  15     """Calculate the discrete derivative (finite difference) of
  16     data[:,f_col] with respect to data[:,x_col].
  17
  18     Examples
  19     --------
  20
  21     >>> import pprint
  22     >>> d = Data((5,2), dtype=numpy.float,
  23     ...          info={'columns':['x', 'x**2']})
  24     >>> for i in range(5):
  25     ...     d[i,0] = i
  26     ...     d[i,1] = i**2
  27     >>> d
  28     Data([[  0.,   0.],
  29            [  1.,   1.],
  30            [  2.,   4.],
  31            [  3.,   9.],
  32            [  4.,  16.]])
  33     >>> dd = derivative(d)
  34     >>> dd
  35     Data([[ 0.,  1.],
  36            [ 1.,  2.],
  37            [ 2.,  4.],
  38            [ 3.,  6.],
  39            [ 4.,  7.]])
  40     >>> pprint.pprint(dd.info)
  41     {'columns': ['x', 'deriv x**2 with respect to x']}
  42
  43     Notes
  44     -----
  45
  46     *Weights*
  47
  48     The returned :class:`Data` block shares its x vector with the
  49     input data.  The ith df/dx value in the returned data is
  50     caclulated with::
  51
  52         (df/dx)[i] = (SUM_j w[j] f[i+j]) / h
  53
  54     where ``h = x[i+1]-x[i]`` is the x coordinate spacing (assumed
  55     constant) and ``j`` ranges over the keys of `weights`.
  56
  57     There standard schemes translate as follows:
  58
  59     ========  ======================  ===================
  60     scheme    formula                 weights
  61     ========  ======================  ===================
  62     forward   ``(f[i+1]-f[i])/h``     ``{0:-1,1:1}``
  63     backward  ``(f[i]-f[i-1])/h``     ``{0:1,-1:-1}``
  64     central   ``(f[i+1]-f[i-1])/2h``  ``{-1:-0.5,1:0.5}``
  65     ========  ======================  ===================
  66
  67     The default scheme is central differencing.
  68
  69     *Boundary conditions*
  70
  71     The boundary conditions could be configurable in principle.  The
  72     current scheme just extrapolates virtual points out to negative
  73     `i` following::
  74
  75         f[i<0] = 2*f[0] - f[-i]
  76
  77     With analogous treatment for `i > data.shape[0]`.  This ensures that
  78     `f[i]-f[0]` is odd about `i=0`, which keeps derivatives smooth.::
  79
  80         f[i] - f[0] = f[0] - f[-i] == -(f[-i] - f[0])
  81     """
  82     output = Data((data.shape[0],2), dtype=data.dtype)
  83     output.info = copy.copy(data.info)
  84     output.info['columns'] = [
  85         data.info['columns'][x_col],
  86         'deriv %s with respect to %s' \
  87         % (data.info['columns'][f_col], data.info['columns'][x_col]),
  88         ]
  89     h = data[1,x_col] - data[0,x_col]
  90     chunks = []
  91     for i,w in weights.items():
  92         chunk = numpy.roll(w*data[:,f_col], -i)
  93         if i > 0: # chunk shifted down, replace the high `i`s
  94             zero = len(chunk) - 1 - i
  95             for j in range(1,i+1):
  96                 chunk[zero+j] = 2*chunk[zero] - chunk[zero-j]
  97         elif i < 0: # chunk shifted up, replace the low `i`s
  98             zero = -i
  99             for j in range(1,zero+1):
 100                 chunk[zero-j] = 2*chunk[zero] - chunk[zero+j]
 101         chunks.append(chunk)
 102     print chunks
 103     output[:,0] = data[:,x_col]
 104     output[:,1] = sum(chunks)
 105     return output