#!/usr/bin/env python
# coding: utf-8

# # Chebyshev polynomials
# 
# Copyright (C) 2020 Andreas Kloeckner
# 
# <details>
# <summary>MIT License</summary>
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
# 
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software.
# 
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
# THE SOFTWARE.
# </details>

# In[18]:


import numpy as np
import numpy.linalg as la
import matplotlib.pyplot as pt


# ## Part I: Plotting the Chebyshev polynomials

# In[2]:


x = np.linspace(-1, 1, 100)

pt.xlim([-1.2, 1.2])
pt.ylim([-1.2, 1.2])

for k in range(10): # crank up
    pt.plot(x, np.cos(k*np.arccos(x)))


# ## Part II: Understanding the Nodes
# 
# What if we interpolate random data?

# In[3]:


n = 50 # crank up


# ### "Extremal" Chebyshev Nodes (or: Chebyshev Nodes of the Second Kind)
# 
# * Most often used for computation
# * Note: Generates $n+1$ nodes -> drop $k$

# In[4]:


k = n-1

i = np.arange(0, k+1)
x = np.linspace(-1, 1, 3000)

def f(x):
    return np.cos(k*np.arccos(x))

nodes = np.cos(i/k*np.pi)

pt.plot(x, f(x))
pt.plot(nodes, f(nodes), "o")


# ### Chebyshev Nodes of the First Kind (Roots)
# 
# * Generates $n$ nodes

# In[5]:


i = np.arange(1, n+1)
x = np.linspace(-1, 1, 3000)

def f(x):
    return np.cos(n*np.arccos(x))

nodes = np.cos((2*i-1)/(2*n)*np.pi)

pt.plot(x, f(x))
pt.plot(nodes, f(nodes), "o")


# ### Observe Spacing

# In[6]:


pt.plot(nodes, 0*nodes, "o")


# ## Part III: Chebyshev Interpolation

# In[15]:


n = 100

i = np.arange(n, dtype=np.float64)
nodes = np.cos((2*(i+1)-1)/(2*n)*np.pi)

V = np.cos(i*np.arccos(nodes.reshape(-1, 1)))
if 1:
    # random data
    data = np.random.randn(n)
else:
    # Runge's example
    data = 1/(1+25*nodes**2)
    
coeffs = la.solve(V, data)


# In[16]:


x = np.linspace(-1, 1, 1000)
Vfull = np.cos(i*np.arccos(x.reshape(-1, 1)))
pt.plot(x, np.dot(Vfull, coeffs))
pt.plot(nodes, data, "o")


# ## Part IV: Conditioning

# In[23]:


n = 10 # crank up

i = np.arange(n, dtype=np.float64)
nodes = np.cos((2*(i+1)-1)/(2*n)*np.pi)
V = np.cos(i*np.arccos(nodes.reshape(-1, 1)))

la.cond(V)


# ## Part V: Error Result
# 
# Plot the product term from the estimate of truncation error in interpolation for the Chebyshev nodes:
# $$\left|\prod_{i=1}^n (x-x_i)\right| $$

# In[45]:


def plot_err_prod(nodes, label):
    eval_pts = np.linspace(-1, 1, 30000)

    product = 1
    for xi in nodes:
        product = product*(eval_pts-xi)
    pt.plot(eval_pts, np.abs(product), label=label)


# In[54]:


n = 10 # crank up

i = np.arange(n, dtype=np.float64)
cheb_nodes = np.cos((2*(i+1)-1)/(2*n)*np.pi)
plot_err_prod(cheb_nodes, label="Chebyshev")

if 0:
    nodes = np.linspace(-1, 1, n)
    plot_err_prod(nodes, label="equispaced")
elif 0:
    nodes = cheb_nodes.copy()
    nodes[3] += 0.1
    plot_err_prod(nodes, label="Perturbed")

pt.legend()


# In[ ]: