%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import seaborn as sns

"""
Explanation: Interpolation
Learning Objective: Learn to interpolate 1d and 2d datasets of structured and unstructured points using SciPy.
End of explanation
"""


x = np.linspace(0,4*np.pi,10)
x

"""
Explanation: Overview
We have already seen how to evaluate a Python function at a set of numerical points:
$$ f(x) \rightarrow f_i = f(x_i) $$
Here is an array of points:
End of explanation
"""


f = np.sin(x)
print(f)

c = np.cos(x)
print (c)

plt.plot(x, f, marker='o')
plt.xlabel('x')
plt.ylabel('f(x)');

"""
Explanation: This creates a new array of points that are the values of $\sin(x_i)$ at each point $x_i$:
End of explanation
"""


from scipy.interpolate import interp1d

"""
Explanation: This plot shows that the points in this numerical array are an approximation to the actual function as they don't have the function's value at all possible points. In this case we know the actual function ($\sin(x)$). What if we only know the value of the function at a limited set of points, and don't know the analytical form of the function itself? This is common when the data points come from a set of measurements.
Interpolation is a numerical technique that enables you to construct an approximation of the actual function from a set of points:
$$ {x_i,f_i} \rightarrow f(x) $$
It is important to note that unlike curve fitting or regression, interpolation doesn't not allow you to incorporate a statistical model into the approximation. Because of this, interpolation has limitations:

It cannot accurately construct the function's approximation outside the limits of the original points.
It cannot tell you the analytical form of the underlying function.

Once you have performed interpolation you can:

Evaluate the function at other points not in the original dataset.
Use the function in other calculations that require an actual function.
Compute numerical derivatives or integrals.
Plot the approximate function on a finer grid that the original dataset.

Warning:
The different functions in SciPy work with a range of different 1d and 2d arrays. To help you keep all of that straight, I will use lowercase variables for 1d arrays (x, y) and uppercase variables (X,Y) for 2d arrays. 
1d data
We begin with a 1d interpolation example with regularly spaced data. The function we will use it interp1d:
End of explanation
"""


x = np.linspace(0,4*np.pi,10) # only use 10 points to emphasize this is an approx
f = np.sin(x)

"""
Explanation: Let's create the numerical data we will use to build our interpolation.
End of explanation
"""


sin_approx = interp1d(x, f, kind='cubic')

"""
Explanation: To create our approximate function, we call interp1d as follows, with the numerical data. Options for the kind argument includes:

linear: draw a straight line between initial points.
nearest: return the value of the function of the nearest point.
slinear, quadratic, cubic: use a spline (particular kinds of piecewise polynomial of a given order.

The most common case you will want to use is cubic spline (try other options):
End of explanation
"""


newx = np.linspace(0,4*np.pi,100)
newf = sin_approx(newx)

"""
Explanation: The sin_approx variabl that interp1d returns is a callable object that can be used to compute the approximate function at other points. Compute the approximate function on a fine grid:
End of explanation
"""


plt.plot(x, f, marker='o', linestyle='', label='original data')
plt.plot(newx, newf, marker='.', label='interpolated');
plt.legend();
plt.xlabel('x')
plt.ylabel('f(x)');

"""
Explanation: Plot the original data points, along with the approximate interpolated values. It is quite amazing to see how the interpolation has done a good job of reconstructing the actual function with relatively few points.
End of explanation
"""


plt.plot(newx, np.abs(np.sin(newx)-sin_approx(newx)))
plt.xlabel('x')
plt.ylabel('Absolute error');

"""
Explanation: Let's look at the absolute error between the actual function and the approximate interpolated function:
End of explanation
"""


x = 4*np.pi*np.random.rand(15)
f = np.sin(x)

sin_approx = interp1d(x, f, kind='cubic')

# We have to be careful about not interpolating outside the range
newx = np.linspace(np.min(x), np.max(x),100)
newf = sin_approx(newx)

plt.plot(x, f, marker='o', linestyle='', label='original data')
plt.plot(newx, newf, marker='.', label='interpolated');
plt.legend();
plt.xlabel('x')
plt.ylabel('f(x)');

plt.plot(newx, np.abs(np.sin(newx)-sin_approx(newx)))
plt.xlabel('x')
plt.ylabel('Absolute error');

"""
Explanation: 1d non-regular data
It is also possible to use interp1d when the x data is not regularly spaced. To show this, let's repeat the above analysis with randomly distributed data in the range $[0,4\pi]$. Everything else is the same.
End of explanation
"""


from scipy.interpolate import interp2d

"""
Explanation: Notice how the absolute error is larger in the intervals where there are no points.
2d structured
For the 2d case we want to construct a scalar function of two variables, given
$$ {x_i, y_i, f_i} \rightarrow f(x,y) $$
For now, we will assume that the points ${x_i,y_i}$ are on a structured grid of points. This case is covered by the interp2d function:
End of explanation
"""


def wave2d(x, y):
    return np.sin(2*np.pi*x)*np.sin(3*np.pi*y)

"""
Explanation: Here is the actual function we will use the generate our original dataset:
End of explanation
"""


x = np.linspace(0.0, 1.0, 10)
y = np.linspace(0.0, 1.0, 10)

"""
Explanation: Build 1d arrays to use as the structured grid:
End of explanation
"""


X, Y = np.meshgrid(x, y)
Z = wave2d(X, Y)

"""
Explanation: Build 2d arrays to use in computing the function on the grid points:
End of explanation
"""


plt.pcolor(X, Y, Z)
plt.colorbar();
plt.scatter(X, Y);
plt.xlim(0,1)
plt.ylim(0,1)
plt.xlabel('x')
plt.ylabel('y');

"""
Explanation: Here is a scatter plot of the points overlayed with the value of the function at those points:
End of explanation
"""


wave2d_approx = interp2d(X, Y, Z, kind='cubic')

"""
Explanation: You can see in this plot that the function is not smooth as we don't have its value on a fine grid.
Now let's compute the interpolated function using interp2d. Notice how we are passing 2d arrays to this function:
End of explanation
"""


xnew = np.linspace(0.0, 1.0, 40)
ynew = np.linspace(0.0, 1.0, 40)
Xnew, Ynew = np.meshgrid(xnew, ynew) # We will use these in the scatter plot below
Fnew = wave2d_approx(xnew, ynew) # The interpolating function automatically creates the meshgrid!

Fnew.shape

"""
Explanation: Compute the interpolated function on a fine grid:
End of explanation
"""


plt.pcolor(xnew, ynew, Fnew);
plt.colorbar();
plt.scatter(X, Y, label='original points')
plt.scatter(Xnew, Ynew, marker='.', color='green', label='interpolated points')
plt.xlim(0,1)
plt.ylim(0,1)
plt.xlabel('x')
plt.ylabel('y');
plt.legend(bbox_to_anchor=(1.2, 1), loc=2, borderaxespad=0.);

"""
Explanation: Plot the original course grid of points, along with the interpolated function values on a fine grid:
End of explanation
"""


from scipy.interpolate import griddata

"""
Explanation: Notice how the interpolated values (green points) are now smooth and continuous. The amazing thing is that the interpolation algorithm doesn't know anything about the actual function. It creates this nice approximation using only the original course grid (blue points).
2d unstructured
It is also possible to perform interpolation when the original data is not on a regular grid. For this, we will use the griddata function:
End of explanation
"""


x = np.random.rand(100)
y = np.random.rand(100)

"""
Explanation: There is an important difference between griddata and the interp1d/interp2d:

interp1d and interp2d return callable Python objects (functions).
griddata returns the interpolated function evaluated on a finer grid.

This means that you have to pass griddata an array that has the finer grid points to be used. Here is the course unstructured grid we will use:
End of explanation
"""


f = wave2d(x, y)

"""
Explanation: Notice how we pass these 1d arrays to our function and don't use meshgrid:
End of explanation
"""


plt.scatter(x, y);
plt.xlim(0,1)
plt.ylim(0,1)
plt.xlabel('x')
plt.ylabel('y');

"""
Explanation: It is clear that our grid is very unstructured:
End of explanation
"""


xnew = np.linspace(x.min(), x.max(), 40)
ynew = np.linspace(y.min(), y.max(), 40)
Xnew, Ynew = np.meshgrid(xnew, ynew)

Xnew.shape, Ynew.shape

Fnew = griddata((x,y), f, (Xnew, Ynew), method='cubic', fill_value=0.0)

Fnew.shape

plt.pcolor(Xnew, Ynew, Fnew, label="points")
plt.colorbar()
plt.scatter(x, y, label='original points')
plt.scatter(Xnew, Ynew, marker='.', color='green', label='interpolated points')
plt.xlim(0,1)
plt.ylim(0,1)
plt.xlabel('x')
plt.ylabel('y');
plt.legend(bbox_to_anchor=(1.2, 1), loc=2, borderaxespad=0.);

"""
Explanation: To use griddata we need to compute the final (strcutured) grid we want to compute the interpolated function on:
End of explanation
"""


