1 # Copyright (C) 2010-2012 W. Trevor King <wking@tremily.us>
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/>.
18 """The `calculus` module provides functions for calculating
19 derivatives and integrals of discrete data.
20 """
22 import copy
24 import numpy
26 from ..curve import Data
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.
33     Examples
34     --------
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.])
51     Notes
52     -----
54     *Weights*
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::
60         (df/dx)[i] = (SUM_j w[j] f[i+j]) / h
62     where ``h = x[i+1]-x[i]`` is the x coordinate spacing (assumed
63     constant) and ``j`` ranges over the keys of `weights`.
65     There standard schemes translate as follows:
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     ========  ======================  ===================
75     The default scheme is central differencing.
77     *Boundary conditions*
79     The boundary conditions could be configurable in principle.  The
80     current scheme just extrapolates virtual points out to negative
81     `i` following::
83         f[i<0] = 2*f - f[-i]
85     With analogous treatment for `i > data.shape`.  This ensures that
86     `f[i]-f` is odd about `i=0`, which keeps derivatives smooth.::
88         f[i] - f = f - f[-i] == -(f[-i] - f)
89     """
90     output = Data(f_data.shape, dtype=f_data.dtype)
91     h = x_data - x_data
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)