############################################################################
## CALCULUS - BIOINFORMATICS, 2018-19 T3
## Authors: Toni Guillamon, Alex Haro
###########################################################################

###Interpolation
#https://docs.scipy.org/doc/scipy/reference/tutorial/interpolate.html
# Take a look at the file Lab0010-ExponentialDecay.py to learn how to save the plots into files. It can be very useful to keep the different plots of this practice.

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 (see Lecture04-Interpolation)

#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(35)
p3(35)

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