import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
import scipy as sp
from scipy import interpolate
import ipywidgets as widgets
import matplotlib as mpl
mpl.rcParams['font.size'] = 14
mpl.rcParams['axes.labelsize'] = 20
mpl.rcParams['xtick.labelsize'] = 14
mpl.rcParams['ytick.labelsize'] = 14
from scipy.interpolate import CubicSpline
M=8

"""
Explanation: <center>
    <h1> ILI285 - Computación Científica I  / INF285 - Computación Científica </h1>
    <h2> Interpolation: Splines </h2>
    <h2> <a href="#acknowledgements"> [S]cientific [C]omputing [T]eam </a> </h2>
    <h2> Version: 1.22</h2>
</center>
Table of Contents

Introduction
Splines
Properties
Solving
The additional Property
Exercises
Acknowledgements
End of explanation
"""


# Code based on Example from: https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.CubicSpline.html#scipy.interpolate.CubicSpline

# The data
x = np.linspace(0,2*np.pi,12)
y = np.sin(x)*x

# Building interpolation object
cs = CubicSpline(x, y)

# Defining a finer mesh to plot the function
xx = np.linspace(0,2*np.pi,1000)
yyo = np.sin(xx)*xx
yyo1 = np.cos(xx)*xx+np.sin(xx)
yyo2 = -np.sin(xx)*xx+2*np.cos(xx)
yyo3 = -np.cos(xx)*xx-3*np.sin(xx)
yyo4 = np.sin(xx)*xx-4*np.cos(xx)

#Interpolating the date with the spline
yy = cs(xx)
yy1 = cs(xx, 1)
yy2 = cs(xx, 2)
yy3 = cs(xx, 3)
yy4 = cs(xx, 4)

# Plotting the splines and its derivatives
plt.figure(figsize=(M,M))
plt.plot(x,y,'k.',markersize=20,label=r'Data Points')
plt.plot(xx,yy, linewidth=4, label=r'S$(x)$')
plt.plot(xx,yy1, linewidth=4, label=r'$\frac{d}{dx}$S$(x)$')
plt.plot(xx,yy2, linewidth=4, label=r'$\frac{d^2}{dx^2}$S$(x)$')
plt.plot(xx,yy3, linewidth=4, label=r'$\frac{d^3}{dx^3}$S$(x)$')
plt.plot(xx,yy4, linewidth=4, label=r'$\frac{d^4}{dx^4}$S$(x)$')

plt.plot(xx,yyo4,'k--',linewidth=4, label='test', alpha=0.4)

plt.plot(x,y,'k.',markersize=20)

plt.title(r'Cubic Spline is defined as S$(x)$')
plt.axis('tight')
plt.xlabel(r'$x$')
plt.ylabel(r'$y$')
plt.grid(True)
plt.legend(loc='center left', bbox_to_anchor=(1, 0.5))
plt.show()


plt.figure()
plt.semilogy(np.abs(yy-yyo))
plt.grid(True)

"""
Explanation: <div id='intro' />

Introduction
Previously in our jupyter notebooks, we learn about interpolation. Methods like Newton's Divided Difference, Lagrange, among others. Other alternative for interpolate a set of data points is using  Cubic Splines.
This technique, avoids the Runge's Phenomenon and creates a 3-degree polynomial easily.
<div id='sp' />

Splines
The most common spline is the linear spline. Given a set of points $(x_{1},y_{1}), (x_{2},y_{2}),...,(x_{n},y_{n}) $, this spline connects each point creating a non-smooth curve. However, this polynomial haves a problem. It's no smooth curve! For to avoid this problem, the cubic splines creates a set of 3-degree polynomial (specifically n-1 polynomials)... a much better curve.
End of explanation
"""


def cubic_spline(x, y, end=None, k1=0, k2=0, p1=0, p2=0):
    #x: x-coordinates of points
    #y: y-coordinates of points
    #end: Natural, Adjusted, Clamped, Parabolically, NaK
    
    n = len(x)
    A = np.zeros((3*n-3, 3*n-3))
    b = np.zeros(3*n-3)
    
    delta_x=np.diff(x)
       
    #Building the linear system of equations
    
    #1st property
    for i in np.arange(n-1):
        b[i]= y[i+1]-y[i]
        A[i,3*i:3*(i+1)] = [delta_x[i],delta_x[i]**2,delta_x[i]**3]
    #2nd property
    for i in np.arange(n-2):
        A[(n-1)+i,3*i:3*(i+1)+1]=[1, 2*delta_x[i], 3*delta_x[i]**2, -1]
    #3rd property
    for i in np.arange(n-2):
        A[(n-1)+(n-2)+i,3*i:3*(i+1)+2] = [0, 2, 6*delta_x[i], 0, -2]
    
    #Ending conditions (4th property)
    if end =='Natural':
        A[-2,1]= 2
        A[-1,-2] = 2
        A[-1,-1] = 6*delta_x[-1]

    elif end == 'Adjusted':
        A[-2,1]= 2
        A[-1,-2] = 2
        A[-1,-1] = 6*delta_x[-1]
        b[-2:] = [k1,k2]
        print('Adjusted',b[-2:])

    elif end == 'Clamped':
        A[-2,0]=1
        A[-1,-3:] = [1,2*delta_x[-1],3*delta_x[-1]**2]
        b[-2:] = [p1,p2]

    elif end == 'Parabolically':
        A[-2,2]=1
        A[-1,-1]=1

    elif end == 'NaK':
        A[-2,2:6]=[6,0,0,-6]
        A[-1,-4:]=[6,0,0,-6]
    
    #Solving the system
    sol = np.linalg.solve(A,b)
    S = {'b':sol[::3],
         'c':sol[1::3],
         'd':sol[2::3],
         'x':x,
         'y':y
        }
    return S

# 'der' computes the 'der'-derivative of the Spline, 
# but it has not been implemented. Can you do it? Please do it!
def cubic_spline_eval(xx,S,der=0):
    x=S['x']
    y=S['y']
    b=S['b']
    c=S['c']
    d=S['d']
    n=len(x)
    yy=np.zeros_like(xx)
    for i in np.arange(n-1):
        jj = np.where(np.logical_and(x[i]<=xx,xx<=x[i+1]))
        yy[jj]=y[i]+b[i]*(xx[jj]-x[i])+c[i]*(xx[jj]-x[i])**2+d[i]*(xx[jj]-x[i])**3
    return yy

"""
Explanation: The orange curve is generated with cubic splines (using the scipy implementation). The other colors are the derivatives of the Cubic Spline as indicated in the legend.
However, if we thinks about this curve, we can say that exists an infinitely many quantity of polynomials such that meets all the points. Our goal is to create an unique polynomial. Given this condition, there is 4 properties that defines the cubic spline we are looking for. 
<div id='pr'/>

Properties of Splines
When we want creates a spline of n data points, we obtains a set of n-1 3-degree polynomials. For example:
Given a set of points $(x_{1},y_{1}), (x_{2},y_{2}),...,(x_{n},y_{n})$, the splines is:
\begin{equation} S_{1}(x) = y_{1} + b_{1}(x-x_{1}) + c_{1}(x-x_{1})^{2} + d_{1}(x-x_{1})^{3} \ 
                 S_{2}(x) = y_{2} + b_{2}(x-x_{2}) + c_{2}(x-x_{2})^{2} + d_{2}(x-x_{2})^{3} \
                 ... \
                 ... \
                 ... \
                 S_{n-1}(x) = y_{n-1} + b_{n-1}(x-x_{n-1}) + c_{n-1}(x-x_{n-1})^{2} + d_{n-1}(x-x_{n-1})^{3}
\end{equation}
Thus, our goal is obtains the $y, b, c$ and $d$ coefficients. With this values, we are creating the spline $S(x)$ that meets all the data points. This spline have the next properties:
Property 1 (Are the points connected?)
The first property checks if each x-coordinate reachs the correspondant y-coordinate in the spline $S(x)$. Thus, checks if the spline meets the data points.
$$S_{i}(x_{i}) = y_{i}$$ $$ S_{i}(x_{i+1}) = y_{i+1}$$ 
$$i \in [1,n-1]$$
Property 2 (Slope Continuity)
The second property ensures that the slopes of the splines at the points' neighborhood to be equals. This guarantees the smoothly of $S(x)$
$$S'{i-1}(x{i}) = S'{i}(x{i})$$
$$i \in [2,n-1]$$
Property 3 (Curvature Continuity)
This property ensure that the curvature between differents polynomials will be equal at the neighborhood of points, avoiding abrupt changes of the curve at the data points. 
$$S''{i-1}(x{i}) = S''{i}(x{i})$$
$$i \in [2,n-1]$$
<div id='so' />

Solving the system
If we have n points, we know that our splines will be composed of n-1 curves $S_{i}(x)$. We have too, (3n-3) unknowns variables ($b_{i}, c_{i}, d_{i}$ for each spline). However, we can build a system of equations for find this variables. How can i do this? Easy.. Using the previous properties!
Using the previously defined splines for n points:
\begin{equation} S_{1}(x) = y_{1} + b_{1}(x-x_{1}) + c_{1}(x-x_{1})^{2} + d_{1}(x-x_{1})^{3} \ 
                 S_{2}(x) = y_{2} + b_{2}(x-x_{2}) + c_{2}(x-x_{2})^{2} + d_{2}(x-x_{2})^{3} \
                 \vdots \
                 S_{n-1}(x) = y_{n-1} + b_{n-1}(x-x_{n-1}) + c_{n-1}(x-x_{n-1})^{2} + d_{n-1}(x-x_{n-1})^{3}
\end{equation}
We need too, the first derivatives of this curves:
\begin{equation} S'{1}(x) = b{1} + 2c_{1}(x-x_{1}) + 3d_{1}(x-x_{1})^{2} \ 
                 S'{1}(x) = b{2} + 2c_{2}(x-x_{2}) + 3d_{2}(x-x_{2})^{2} \
                 \vdots \
                 S'{n-1}(x) = b{n-1} + 2c_{n-1}(x-x_{n-1}) + 3d_{n-1}(x-x_{n-1})^{2} \
\end{equation}
And its second derivatives:
\begin{equation} S''{1}(x) = 2c{1} + 6d_{1}(x-x_{1}) \ 
                 S''{2}(x) = 2c{2} + 6d_{2}(x-x_{2}) \ 
                 \vdots \
                 S''{n-1}(x) = 2c{n-1} + 6d_{n-1}(x-x_{n-1}) \ 
\end{equation}
Using the first property, we get (n-1) equations:
\begin{equation} b_{1}(x_{2}-x_{1}) + c_{1}(x_{2}-x_{1})^2 + d_{1}(x_{2}-x_{1})^3 = y_{2} - y_{1} \hspace{1cm}(1)\
                 b_{2}(x_{3}-x_{2}) + c_{2}(x_{3}-x_{2})^2 + d_{2}(x_{3}-x_{2})^3 = y_{3} - y_{2} \hspace{1cm}(2)\ 
                 \vdots\
                 b_{n-1}(x_{n}-x_{n-1}) + c_{n-1}(x_{n}-x_{n-1})^2 + d_{n-1}(x_{n}-x_{n-1})^3 = y_{n} - y_{n-1} \hspace{1cm}(n-1) 
\end{equation}
Using the second property, we get (n-2) equations:
\begin{equation} b_{1}+2c_{1}(x_{2}-x_{1}) + 3d_{1}(x_{2}-x_{1})^2 - b_{2}= 0 \hspace{1cm}(1)\ 
                 b_{2}+2c_{2}(x_{3}-x_{2}) + 3d_{2}(x_{3}-x_{2})^2 - b_{3}=  0 \hspace{1cm}(2)\ 
                 \vdots\
                 b_{n-2}+2c_{n-2}(x_{n-1}-x_{n-2}) + 3d_{n-2}(x_{n-1}-x_{n-2})^2 -b_{n-1}=0  \hspace{1cm}(n-2)\ 
\end{equation}
Using the third property, we get (n-2) equations:
\begin{equation} 2c_{1}+6d_{1}(x_{2}-x_{1}) - 2c_{2} = 0 \hspace{1cm}(1)\ 
                 2c_{2}+6d_{2}(x_{3}-x_{2}) - 2c_{3}=0 \hspace{1cm}(2)\ 
                 \vdots\
                 2c_{n-2}+6d_{n-2}(x_{n-1}-x_{n-2}) - 2c_{n-1} = 0 \hspace{1cm}(n-2)\ 
\end{equation}
If we adds all our equations, we obtains (3n-5) equations. Clearly, the matrix in that system is not square (we need 2 equations more). For this, we have another property, that defines the edges conditions of the splines
<div id='ad'/>

Splines Ending options
For this special property, we have the following 5 properties:
Natural Spline:
This property create a spline with zero curvature, thus:
\begin{align}
S''{1}(x{1}) &= 2c_{1} = 0\
S''{n-1}(x{n}) &= 2c_{n-1}+6d_{n-1}(x_{n}-x_{n-1}) = 0
\end{align}
Adjusted curvature:
This property create a spline which curvature is equal to a parameter previously defined, not necessarily zero.
\begin{align}
S''{1}(x{1}) &= 2c_{1} = \kappa_{1}\
S''{n-1}(x{n}) &= 2c_{n-1}+6d_{n-1}(x_{n}-x_{n-1}) = \kappa_{2}
\end{align}
Clamped cubic spline:
This property adjust the slopes at the edges of splines to a value previously defined
\begin{align}
S'{1}(x{1}) & = b_{1} = p_{1} \
S'{n-1}(x{n}) & = b_{n-1}+2c_{n-1}(x_{n}-x_{n-1}) + 3d_{n-1}(x_{n}-x_{n-1})^2  = p_{2}
\end{align}
Ended Parabolically
With this property, the edges of the splines are 2-degree polynomials. Hence the coefficients:
\begin{align}
d_{1} &= 0 \
d_{n-1} &= 0
\end{align}

Why does this property does not work for $n\le 3$?

Not-a-Knot cubic spline
This condition, checks the continuity at the edges, for the third derivative:
\begin{align}
S'''{1}(x{2}) &= S'''{2}(x{2})\
6d_{1}&=6d_{2}\
S'''{n-2}(x{n-1}) &= S'''{n-1}(x{n-1})\
6d_{n-2}&=6d_{n-1}
\end{align}

Why does this property does not work for $n\le 4$?

Each property give us the 2 equations needed. Thanks to this, we have (3n-3) unknowns and equations
Finally, for to find the coefficients of the spline, we''ll build the system of equations.
Example:

If we have the following 3 points: $(x_1,y_1),(x_1,y_1),(x_3,y_3)$
We will get 6 unknowns
So, to build a Natural spline, we need to solve the following linear system of equation:

$$ \begin{bmatrix} (x_2-x_1) & (x_2-x_1)^2 & (x_2-x_1)^3 & 0 & 0 & 0 \ 
                  0 & 0 & 0 &(x_3-x_2) & (x_3-x_2)^2 & (x_3-x_2)^3  \
                  1 & 2(x_2-x_1) & 3(x_2-x_1)^2 & -1 & 0 & 0 \
                  0 & 2 & 6(x_2-x_1) & 0 & -2 & 0 \
                  0 & 2 & 0 & 0 & 0 & 0 \
                  0 & 0 & 0 & 0 & 2 & 6(x_3-x_2) \
                  \end{bmatrix} 
                  \left[ \begin{array}{c} b_1 \ c_1 \ d_1 \ b_2 \ c_2 \ d_2 \end{array} \right] = 
                  \left[ \begin{array}{c} y_2-y_1 \ y_3-y_2 \ 0 \ 0 \ 0 \ 0  \end{array} \right] 
$$
Finally, find the coefficients of splines is reduced to solve an equation system, and we already know this from previous notebooks!
Now the code:
End of explanation
"""


x = np.array([1,2,4,5])
y = np.array([2,1,4,3])
S = cubic_spline(x,y,end='Natural')

x1 = np.linspace(1,2,200)
x2 = np.linspace(2,4,200)
x3 = np.linspace(4,5,200)

S1 = y[0]+S['b'][0]*(x1-x[0])+S['c'][0]*(x1-x[0])**2+S['d'][0]*(x1-x[0])**3
S2 = y[1]+S['b'][1]*(x2-x[1])+S['c'][1]*(x2-x[1])**2+S['d'][1]*(x2-x[1])**3
S3 = y[2]+S['b'][2]*(x3-x[2])+S['c'][2]*(x3-x[2])**2+S['d'][2]*(x3-x[2])**3

plt.figure(figsize=(M,M))
plt.plot(x,y,'k.',markersize=20,label='Data Points')
plt.plot(x1,S1,'b',linewidth=5,label=r'S$1(x)$')
plt.plot(x2,S2,'g',linewidth=5,label=r'S$2(x)$')
plt.plot(x3,S3,'r',linewidth=5,label=r'S$2(x)$')
plt.xlabel('$x$')
plt.ylabel('$y$')
plt.legend(loc='center left', bbox_to_anchor=(1, 0.5))
plt.grid(True)
plt.show()

"""
Explanation: Example 1 - Hand made interpolation
End of explanation
"""


def show_spline(type_ending='Natural',k1=0, k2=0, p1=0, p2=0):
    x = np.array([1,2,4,5,7,9])
    y = np.array([2,1,4,3,3,4])
    xx=np.linspace(np.min(x),np.max(x),1000)
    S = cubic_spline(x,y,end=type_ending, k1=k1, k2=k2, p1=p1, p2=p2)
    plt.figure(figsize=(M,M))
    plt.plot(xx,cubic_spline_eval(xx,S),'-',linewidth=5,label=r'S$(x)$')
    plt.plot(x,y,'k.',markersize=20,label='Data Points')
    plt.xlabel('$x$')
    plt.ylabel('$y$')
    plt.legend(loc='center left', bbox_to_anchor=(1, 0.5))
    plt.grid(True)
    plt.show()
widgets.interact(show_spline, type_ending=['Natural','Adjusted','Clamped','Parabolically','NaK'],
                k1=(-20,20,1),k2=(-20,20,1),p1=(-2,2,0.2),p2=(-2,2,0.2))

"""
Explanation: Example 2
End of explanation
"""


