# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

import numpy as np
import matplotlib.pyplot as plt

# # Ortographic and oblique projections
#
# In this notebook you play with ortographic and oblique projections using non-orthogonal basis. 
# In the first part you will find a best approximation of a Gaussian function $x$ using symmetric trygonometric polynomials (powers of $\cos(t)$) via ortographic projection.
#
# In the second part you will calculate an oblique projection using different two different bases. 
#
# ## Orographic projection
#
# We generated for you both the function and the basis $\Phi$ of the subspace of powers of $\cos(t)$ of dimension `dimension_phi` that you can change and play with. Run the code below to visualise them.
#
# You can also change $x$ to different functions.

t = np.linspace(-np.pi, np.pi, 900)
step_size = 2 / t.size
x = np.exp(-3*t**2)
plt.plot(t, x, label='$x(t)$')
plt.legend()
plt.show()

# +
dimension_phi = 4
Phi = np.stack([np.cos(t)**k for k in range(dimension_phi)], axis=-1)

for k in range(dimension_phi):
    plt.plot(t, Phi[:,k], label=f'$\\phi_{k}$')
plt.gca().axhline(c="k", linewidth=1)
plt.legend()
plt.show()
# -

# To calculate the best approximation, we will find the ortogonal projection of $x$ into the subspace.
# First, fill in the code bellow to calculate the Gramian. What is the `step_size` for? Are the elements of $\Phi$ orthogonal?

# +
Gram = step_size * ...

plt.imshow(Gram)
plt.colorbar()
plt.show()
# -

# Now you can use the theorem about Orthogonal Projection Onto a Subspace to from $\Phi$ an orthogonal projection $P_s$.
#
# Fill in the code below to calculate and plot the projection of $x$ onto trygonometric polynomials $P_S x$. Calculate and plot the coefficents $a_k$ representing $\Phi x$ in the basis $\Phi$.

# +
basis_coeffients = ...
projection = ...

plt.plot(t, x, label='$x$')
plt.plot(t, x - projection, label='$x-P_S x$')
plt.plot(t, projection, label='$P_S x$')
plt.legend()
plt.gca().axhline(c="k", linewidth=1)
plt.show()
# -

# You can plot basis coefficients $a_k$ or can also plot each basis component of $x$ separately:

# +
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 4))

for k in range(dimension_phi):
    ax1.stem([k], [basis_coeffients[k]], label=f"$a_{k}$", markerfmt=f"C{k}o", linefmt=f"C{k}-")
ax1.legend()
ax1.axhline(c="k", linewidth=1)

for k in range(dimension_phi):
    ax2.plot(t, basis_coeffients[k]*Phi[:, k], label=f"$a_{k}\\Phi_{k}$")
ax2.legend()
ax2.axhline(c="k", linewidth=1)
plt.show()
# -

# ## Oblique projection
# Now we will consider a situation when we have two different, non orthogonal basis. 
#
# Lets keep basis $\Phi$ for the analysis, but project into space of shifted Gaussians of dimension `dimension_psi`. The basis $\Psi$ of this space is generated for you below:

# +
dimension_psi = 6
shifts = np.linspace(-np.pi/2, np.pi/2, 10)
Psi = np.exp(-(t[:,None]-shifts[None,:])**2/0.4)

for k in range(dimension_psi):
    plt.plot(t, Psi[:,k], label=f'$\\psi_{k}$')
plt.gca().axhline(c="k", linewidth=1)
plt.legend()
plt.show()
# -

# First, calculate generalised Gramian, $\Phi^* \Psi$. Are the of $\Phi$ and $\Psi$ biorthogonal?

# +
Gram_gen = ...

plt.imshow(Gram_gen)
plt.colorbar()
plt.show()
# -

# Now you can use the theorem about General Oblique Projections to from projection on to the space of shifted Gaussians.

# +
...

plt.plot(t, x, label='$x$')
plt.plot(t, x - oblique_projection, label='$x-P_S x$')
plt.plot(t, oblique_projection, label='$P_S x$')
plt.legend()
plt.gca().axhline(c="k", linewidth=1)
plt.show()

# +
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 4))

for k in range(dimension_psi):
    ax1.stem([k], [oblique_basis_coefficients[k]], label=f"$a_{k}$", markerfmt=f"C{k}o", linefmt=f"C{k}-")
ax1.legend()
ax1.axhline(c="k", linewidth=1)

for k in range(dimension_psi):
    ax2.plot(t, oblique_basis_coefficients[k]*Psi[:, k], label=f"$a_{k}\\Psi_{k}$")
ax2.legend()
ax2.axhline(c="k", linewidth=1)
plt.show()

StarcoderdataJupyterScriptsDedupFiltered