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
19 """The `calculus` module provides functions for calculating
20 derivatives and integrals of discrete data.
21 """
23 import copy
25 import numpy
27 from ..curve import Data
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.
34     Examples
35     --------
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.])
52     Notes
53     -----
55     *Weights*
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::
61         (df/dx)[i] = (SUM_j w[j] f[i+j]) / h
63     where ``h = x[i+1]-x[i]`` is the x coordinate spacing (assumed
64     constant) and ``j`` ranges over the keys of `weights`.
66     There standard schemes translate as follows:
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     ========  ======================  ===================
76     The default scheme is central differencing.
78     *Boundary conditions*
80     The boundary conditions could be configurable in principle.  The
81     current scheme just extrapolates virtual points out to negative
82     `i` following::
84         f[i<0] = 2*f - f[-i]
86     With analogous treatment for `i > data.shape`.  This ensures that
87     `f[i]-f` is odd about `i=0`, which keeps derivatives smooth.::
89         f[i] - f = f - f[-i] == -(f[-i] - f)
90     """
91     output = Data(f_data.shape, dtype=f_data.dtype)
92     h = x_data - x_data
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)