eand.py is an open-source project released under the terms of GNU General Public License version 3
'''
eand.py package (Easy Algebraic Numerical Differentiation for Python)
Copyright (C) 2013 Tu-Hoa Pham
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 3 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program. If not, see <http://www.gnu.org/licenses/>.
'''
import matplotlib.pyplot as plt
from math import cos,sin
import numpy as np
from eand.mddig.multidiffrep import MultiDiffRep
from mpl_toolkits.mplot3d import Axes3D
from IPython.core.display import Image
In this example, we study a function of two variables: signal(x,y) = cos(2(x+y)) The numerical differentiation is based on 4000 samples drawn randomly on domain D = [-1,1]x[-1,1]
Ns = 4000
t1Min = -1.0
t1Max = 1.0
t2Min = -1.0
t2Max = 1.0
t1 = np.array([0. for _ in range(Ns)])
t2 = np.array([0. for _ in range(Ns)])
for i in range(Ns):
t1[i] = t1Min+(t1Max-t1Min)*np.random.random()
t2[i] = t2Min+(t2Max-t2Min)*np.random.random()
tVec = [t1,t2]
signal = [cos(2*sum([t[i] for t in tVec])) for i in range(len(tVec[0]))]
We may represent the initial signal in the usual 3D reference frame
ax = Axes3D(plt.figure())
ax.scatter3D(tVec[0],tVec[1],signal)
plt.show()
# Image(filename='01-initial_signal.png')
We perform numerical differentiation using minimal estimation parameters with a sliding integration window of 0.50 width along each axis
n1 = 1
alpha1 = 1
beta1 = 1
n2 = 0
alpha2 = 0
beta2 = 0
T1 = 0.25 # half time window on first axis
T2 = 0.25 # half time window on second axis
paramVec1 = [[n1,alpha1,beta1,T1],[n2,alpha2,beta2,T2]] # differentiate first along the first axis
paramVec2 = [[n2,alpha2,beta2,T2],[n1,alpha1,beta1,T1]] # then along the second coordinate
paramVecSeq = [paramVec1,paramVec2]
Class MultiDiff of package eand.mddig (MultiDimensional Differentiation on Irregular Grids) performs direct numerical differentiation from initial signal while class MultiDiffRep implements consecutive redifferentiation
multiDiffRep = MultiDiffRep(paramVecSeq,tVec)
Once all differentiators are constructed, method differentiate() computes successive derivative estimates
(tPostVecSeq,dPostSeq) = multiDiffRep.differentiate(signal)
Numerical differentiation is performed by partitioning the sampling space using the Delaunay triangulation method in the 2D case
multiDiffRep.differentiators[0].plotPartition()
# Image(filename='02-initial_partition.png')
The resulting estimates can then be represented as scatter plots
# First order derivative scatter plot
multiDiffRep.differentiators[0].plotScatter(tPostVecSeq[0], dPostSeq[0])
# Image(filename='03-first_order_scatter.png')
# Second order derivative scatter plot
multiDiffRep.differentiators[1].plotScatter(tPostVecSeq[0], dPostSeq[0])
# Image(filename='04-second_order_scatter.png')
The same estimates can be represented as surfaces based on Delaunay partitioning
# First order derivative surface plot
multiDiffRep.differentiators[0].plotSurface(tPostVecSeq[0], dPostSeq[0])
# Image(filename='05-first_order_surface.png')
# Second order derivative surface plot
multiDiffRep.differentiators[1].plotSurface(tPostVecSeq[1], dPostSeq[1])
# Image(filename='06-second_order_surface.png')
It is also possible to represent just a slice of the 3D estimate plot specifying an axis number (e.g. 1 for first axis), an associated coordinate (e.g. 0.) and an intersection threshold (e.g 0.1)
# First order estimates for |x - 0.| <= 0.1
tSlice1,dSlice1 = multiDiffRep.differentiators[0].plotSlice(tPostVecSeq[0],dPostSeq[0],1,0.,0.1,0,plotNow=0)
# First order derivative refence
dRef1 = [-2*sin(2*i) for i in tSlice1]
plt.plot(tSlice1,dRef1,'r.')
# Image(filename='07-first_order_slice.png')
# Second order estimates for |x - 0.| <= 0.1
tSlice2,dSlice2 = multiDiffRep.differentiators[0].plotSlice(tPostVecSeq[0],dPostSeq[0],1,0.,0.1,0,plotNow=0)
# Second order derivative refence
dRef2 = [-4*cos(2*i) for i in tSlice2]
plt.plot(tSlice2,dRef2,'r.')
plt.show()
# Image(filename='08-second_order_slice.png')