Demonstration for eand.py package (Easy Algebraic Numerical Differentiation for Python)

We present here the problem of numerical differentiation for a signal characterized by the following features: a. multi-dimensionality: the signal depends on more than one variable b. irregular grid: the signal's sampling points are unorganized (e.g. random) c. noisy environment: samples are subject to uncertainties

eand.py is an open-source project released under the terms of GNU General Public License version 3

In [2]:
'''
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]

In [3]:
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

In [6]:
ax = Axes3D(plt.figure())
ax.scatter3D(tVec[0],tVec[1],signal)
plt.show()
# Image(filename='01-initial_signal.png') 
Out[6]:

We perform numerical differentiation using minimal estimation parameters with a sliding integration window of 0.50 width along each axis

In [8]:
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

In [9]:
multiDiffRep = MultiDiffRep(paramVecSeq,tVec)

Once all differentiators are constructed, method differentiate() computes successive derivative estimates

In [10]:
(tPostVecSeq,dPostSeq) = multiDiffRep.differentiate(signal)

Numerical differentiation is performed by partitioning the sampling space using the Delaunay triangulation method in the 2D case

In [22]:
multiDiffRep.differentiators[0].plotPartition()
# Image(filename='02-initial_partition.png') 
Out[22]: