hooke/util/calculus.py

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