'''
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 

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]))]

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

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]

multiDiffRep = MultiDiffRep(paramVecSeq,tVec)

(tPostVecSeq,dPostSeq) = multiDiffRep.differentiate(signal)

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

# 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') 

# 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') 

# 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')