hooke/util/calculus.py

   1 # Copyright (C) 2010 W. Trevor King <wking@drexel.edu>
   2 #
   3 # This file is part of Hooke.
   4 #
   5 # Hooke is free software: you can redistribute it and/or modify it
   6 # under the terms of the GNU Lesser General Public License as
   7 # published by the Free Software Foundation, either version 3 of the
   8 # License, or (at your option) any later version.
   9 #
  10 # Hooke is distributed in the hope that it will be useful, but WITHOUT
  11 # ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  12 # or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General
  13 # Public License for more details.
  14 #
  15 # You should have received a copy of the GNU Lesser General Public
  16 # License along with Hooke.  If not, see
  17 # <http://www.gnu.org/licenses/>.
  18
  19 """The `calculus` module provides functions for calculating
  20 derivatives and integrals of discrete data.
  21 """
  22
  23 import copy
  24
  25 import numpy
  26
  27 from ..curve import Data
  28
  29
  30 def derivative(x_data, f_data, weights={-1:-0.5, 1:0.5}):
  31     """Calculate the discrete derivative (finite difference) of
  32     f_data with respect to x_data.
  33
  34     Examples
  35     --------
  36
  37     >>> d = Data((5,2), dtype=numpy.float,
  38     ...          info={'columns':['x', 'x**2']})
  39     >>> for i in range(5):
  40     ...     d[i,0] = i
  41     ...     d[i,1] = i**2
  42     >>> d
  43     Data([[  0.,   0.],
  44            [  1.,   1.],
  45            [  2.,   4.],
  46            [  3.,   9.],
  47            [  4.,  16.]])
  48     >>> dd = derivative(x_data=d[:,0], f_data=d[:,1])
  49     >>> dd
  50     Data([ 1.,  2.,  4.,  6.,  7.])
  51
  52     Notes
  53     -----
  54
  55     *Weights*
  56
  57     The returned :class:`Data` block shares its x vector with the
  58     input data.  The `i`\th df/dx value in the returned data is
  59     caclulated with::
  60
  61         (df/dx)[i] = (SUM_j w[j] f[i+j]) / h
  62
  63     where ``h = x[i+1]-x[i]`` is the x coordinate spacing (assumed
  64     constant) and ``j`` ranges over the keys of `weights`.
  65
  66     There standard schemes translate as follows:
  67
  68     ========  ======================  ===================
  69     scheme    formula                 weights
  70     ========  ======================  ===================
  71     forward   ``(f[i+1]-f[i])/h``     ``{0:-1,1:1}``
  72     backward  ``(f[i]-f[i-1])/h``     ``{0:1,-1:-1}``
  73     central   ``(f[i+1]-f[i-1])/2h``  ``{-1:-0.5,1:0.5}``
  74     ========  ======================  ===================
  75
  76     The default scheme is central differencing.
  77
  78     *Boundary conditions*
  79
  80     The boundary conditions could be configurable in principle.  The
  81     current scheme just extrapolates virtual points out to negative
  82     `i` following::
  83
  84         f[i<0] = 2*f[0] - f[-i]
  85
  86     With analogous treatment for `i > data.shape[0]`.  This ensures that
  87     `f[i]-f[0]` is odd about `i=0`, which keeps derivatives smooth.::
  88
  89         f[i] - f[0] = f[0] - f[-i] == -(f[-i] - f[0])
  90     """
  91     output = Data(f_data.shape, dtype=f_data.dtype)
  92     h = x_data[1] - x_data[0]
  93     chunks = []
  94     for i,w in weights.items():
  95         chunk = numpy.roll(w*f_data, -i)
  96         if i > 0: # chunk shifted down, replace the high `i`s
  97             zero = len(chunk) - 1 - i
  98             for j in range(1,i+1):
  99                 chunk[zero+j] = 2*chunk[zero] - chunk[zero-j]
 100         elif i < 0: # chunk shifted up, replace the low `i`s
 101             zero = -i
 102             for j in range(1,zero+1):
 103                 chunk[zero-j] = 2*chunk[zero] - chunk[zero+j]
 104         chunks.append(chunk)
 105     return sum(chunks)