############################################################################
## SESSIÓ PROGRAMARIA EPSEB 20190603
## INTERPOLACIÓ
## Toni Guillamon
###########################################################################

###Interpolation
#https://docs.scipy.org/doc/scipy/reference/tutorial/interpolate.html
###

import numpy as np
import matplotlib.pyplot as plt
from scipy import interpolate
from scipy.interpolate import interp1d
#from pylab import figure, title, show, savefig

###EXAMPLE GLYCERIN

#We introduce the full table first and then specific subsets that we will work separately.
xp=[0,10,20,30,40,50,60,70,80,90,100]
yp=[0,-1.6,-4.8,-9.5,-15.4,-21.9,-33.6,-37.8,-19.1,-1.6,17]
xp1=[30,40,50,60]
yp1=[-9.5,-15.4,-21.9,-33.6]

#Refined vectors for the plots
#xf = np.arange(0.,100.,0.5)
xf = np.linspace(0, 100, num=201, endpoint=True)
#xf1 = np.arange(30.,60.,0.5)
xf1 = np.linspace(30, 60, num=61, endpoint=True)

#Linear interpolation
fL = interp1d(xp, yp)
#Cubic (spline) interpolation
fS = interp1d(xp, yp, kind='cubic')

#We plot together the data, the linear and the cubic-spline interpolation:
#1. Overlapping representation:
plt.plot(xp, yp, 'o', xf, fL(xf), '-', xf, fS(xf), '--')
plt.legend(['data', 'linear', 'cubic'], loc='best')
plt.show()

#2. In separate subplots (not so striking, just to introduce subplots):
#plt.subplot(311)
#plt.plot(xp, yp, 'o')
#plt.legend(['data'], loc='best')
#plt.subplot(312)
#plt.plot(xf, f(xf), '-')
#plt.legend(['linear'], loc='best')
#plt.subplot(313)
#plt.plot(xf, f2(xf), '--')
#plt.legend(['cubic'], loc='best')
#plt.show()

#We proceed with the polynomial interpolation.
#The command "polyfit" can both interpolate (if the number of points equals the degree plus one) and approximate (when the number of points is bigger than the degree plus one).

#We first interpolate the 4 data points given in the interval [30,60], and compare with the spline computed above:
z3 = np.polyfit(xp1, yp1, 3)
p3 = np.poly1d(z3)

#To evaluate the interpolating function at some point (for instance, x=35), just do...
fS(45)
p3(45)

#plt.plot(xp,yp,'o',xf,fS(xf),'--',xf,p3(xf),'r')
plt.plot(xp1,yp1,'o',xf1,fS(xf1),'--',xf1,p3(xf1),'r')
plt.legend(['data', 'cubic spline', 'polynomial'], loc='best')
plt.show()

#We repeat the exercise with the 11 data points given in the interval [0,100], and compare with the spline computed above. Note the strong (unlike) oscillations of the polynomial of degree 10:
z10 = np.polyfit(xp, yp, 10)
p10 = np.poly1d(z10)

plt.plot(xp,yp,'o',xf,fS(xf),'--',xf,p10(xf),'r')
plt.legend(['data', 'cubic spline', 'polynomial'], loc='best')
plt.show()

###NEWTON DIVIDED DIFFERENCES ALGORITHM: self-made function
def PNewt(xp,yp,xa):
    n=len(yp)
    m=len(xa)
    ypLoc=np.zeros(n)
    ya=np.zeros(m)
    for i in np.arange(0,n):
        ypLoc[i]=yp[i]
    for i in np.arange(1,n):
        for j in np.arange(n-1,i-1,-1):
            ypLoc[j]=(ypLoc[j]-ypLoc[j-1])/(xp[j]-xp[j-i])
            print(xp[j]-xp[j-i],i,j)
#   print(ypLoc)
    for jm in np.arange(0,m):
        p=ypLoc[0]
        for j in np.arange(1,n):
            z=1
            for k in np.arange(0,j):
                z=(xa[jm]-xp[k])*z
            p=p+z*ypLoc[j]
        ya[jm]=p
    return ya

QNewt=PNewt(xp,yp,xf)
QNewt

plt.plot(xp,yp,'o',xf,QNewt,'r')
plt.legend(['data','polynomial'], loc='best')
plt.show()