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

# ## Wave equation

# + raw_mimetype="text/latex" active=""
# Using our framework, we want to infer the constant c from the Wave Equation. The Wave Equation is given by
#
# \begin{align*}
# \frac{\partial^2 u}{\partial t^2} = c\nabla^2u,
# \end{align*}
#
# where $u = u(x_1,x_2, \dotsc, x_n, t)$ and $c>0$ is some constant \phantomsection\label{\detokenize{index:id17}}{\hyperref[\detokenize{index:lamoureux2006mathematics}]{\sphinxcrossref{{[}17{]}}}}. In one spatial dimension it boils down to:
#
# \begin{align} \label{4}
# \frac{\partial^2 u}{\partial t^2} - c\frac{\partial^2 u}{\partial x^2} = 0.
# \end{align}
#
# We generate the data from a solution of the equation (\ref{4}) corresponding to $c=1$ and get an estimation of $c = 1.0003$. 
# -

# #### Problem Setup

# + raw_mimetype="text/latex" active=""
# \begin{align*}
# u_{tt} - c u_{xx} = 0
# \end{align*}
# -

# The general solution is given by:
# $u(x,t) = F(x-ct) + G(x+ct)$ with F, G some functions.
#
# Take $F(x) = x^2$ and $G(x) = \sin(x)$ and $c=1$.
#
# Thus: $u(x,t) = (x-t)^2 + \sin(x + t)$.
#
# Set $f = 0$.
#
# Consider $u$ to be a Gaussian process:
#
# $u \sim \mathcal{GP}(0, k_{uu}(x_i, x_j; \tilde{\theta}))$ with the hyperparameters $\tilde{\theta} = \{\theta, l_x, l_t\}$.
#
# And the linear operator:
#
# $\mathcal{L}_x^c = \frac{d^2}{dt^2} \cdot - c \frac{d^2}{dx^2} \cdot$
#
# so that
#
# $\mathcal{L}_x^c u = f$
#
# Problem at hand: Estimate $c$ (should be $c = 1$ in the end).
#
#
# #### Step 1: Simulate data

# + nbsphinx="hidden"
import time
import numpy as np
import sympy as sp
from scipy.optimize import minimize
import matplotlib.pyplot as plt
import warnings


# -

# $x \in [0, 1]^n, \; t \in [0,1]^n$

# +
def get_simulated_data(n = 20):
    t = np.random.rand(n)
    x = np.random.rand(n)
    y_u = np.multiply(x-t, x-t) + np.sin(x+t)
    y_f = 0*x
    return(x, t, y_u, y_f)

(x, t, y_u, y_f) = get_simulated_data()
# -

# #### Step 2: Evaluate kernels
#
# 1) $k_{uu}(y_i, y_j; \tilde{\theta}) = \theta exp(-\frac{1}{2l_x}(x_i-x_j)^2 - \frac{1}{2l_t}(t_i-t_j)^2)$, where $y_i = (x_i, t_i)$, $y_j = (x_j, t_j)$.

# + nbsphinx="hidden"
x_i, x_j, t_i, t_j, theta, l_x, l_t, c = sp.symbols('x_i x_j t_i t_j theta l_x l_t c')
kuu_sym = theta*sp.exp(-1/(2*l_x)*((x_i - x_j)**2) - 1/(2*l_t)*((t_i - t_j)**2))
kuu_fn = sp.lambdify((x_i, x_j, t_i, t_j, theta, l_x, l_t), kuu_sym, "numpy")
def kuu(x, t, theta, l_x, l_t):
    k = np.zeros((x.size, x.size))
    for i in range(x.size):
        for j in range(x.size):
            k[i,j] = kuu_fn(x[i], x[j], t[i], t[j], theta, l_x, l_t)
    return k


# -

# 2) $k_{ff}(y_i,y_j; \tilde{\theta}, c) 
# = \mathcal{L}_{y_i}^c \mathcal{L}_{y_j}^c k_{uu}(y_i, y_j; \tilde{\theta}) \\
# = \frac{d^4}{dt_i^2 dt_j^2}k_{uu} - c\frac{d^4}{dt_i^2 dx_j^2}k_{uu} - c\frac{d^4}{dx_i^2 dt_j^2}k_{uu} + c^2\frac{d^4}{dx_i^2 dx_j^2}k_{uu}$

# + nbsphinx="hidden"
kff_sym = sp.diff(kuu_sym, t_i, t_i, t_j, t_j) \
        - c*sp.diff(kuu_sym, t_i, t_i, x_j, x_j) \
        - c*sp.diff(kuu_sym, x_i, x_i, t_j, t_j) \
        + c**2*sp.diff(kuu_sym, x_i, x_i, x_j, x_j)
kff_fn = sp.lambdify((x_i, x_j, t_i, t_j, theta, l_x, l_t, c), kff_sym, "numpy")
def kff(x, t, theta, l_x, l_t, c):
    k = np.zeros((x.size, x.size))
    for i in range(x.size):
        for j in range(x.size):
            k[i,j] = kff_fn(x[i], x[j], t[i], t[j], theta, l_x, l_t, c)
    return k


# -

# 3) $k_{fu}(y_i,y_j;\tilde{\theta}, c)
# = \mathcal{L}_{\tilde{x}_i}^c k_{uu}(y_i, y_j; \tilde{\theta})
# = \frac{d^2}{dt_i^2}k_{uu} - c\frac{d^2}{dx_i^2}k_{uu}$

# + nbsphinx="hidden"
kfu_sym = sp.diff(kuu_sym, t_i, t_i) - c*sp.diff(kuu_sym, x_i, x_i)
kfu_fn = sp.lambdify((x_i, x_j, t_i, t_j, theta, l_x, l_t, c), kfu_sym, "numpy")
def kfu(x, t, theta, l_x, l_t, c):
    k = np.zeros((x.size, x.size))
    for i in range(x.size):
        for j in range(x.size):
            k[i,j] = kfu_fn(x[i], x[j], t[i], t[j], theta, l_x, l_t, c)
    return k


# -

# 4) $k_{uf}(y_i, y_j; \tilde{\theta}, c)$ is given by the transpose of $k_{fu}(y_i, y_j; \tilde{\theta}, c)$.

# + nbsphinx="hidden"
def kuf(x, t, theta, l_x, l_t, c):
    return kfu(x, t, theta, l_x, l_t, c).T


# -

# #### Steps 3 and 4: Compute NLML and optimize the hyperparameters

# + nbsphinx="hidden"
def nlml(params, x, t, y1, y2, s):
    params = np.exp(params)
    K = np.block([
        [kuu(x, t, params[0], params[1], params[2]) + s*np.identity(x.size), kuf(x, t, params[0], params[1], params[2], params[3])],
        [kfu(x, t, params[0], params[1], params[2], params[3]), kff(x, t, params[0], params[1], params[2], params[3]) + s*np.identity(x.size)]
    ])
    y = np.concatenate((y1, y2))
    val = 0.5*(np.log(abs(np.linalg.det(K))) + np.mat(y) * np.linalg.inv(K) * np.mat(y).T)
    return val.item(0)


# + nbsphinx="hidden"
def minimize_restarts(x,t,y_u,y_f,n=10): 
    nlml_wp = lambda params: nlml(params, x, t, y_u, y_f, 1e-7)
    all_results = []
    for it in range(0,n):
        all_results.append(minimize(nlml_wp, np.random.rand(4), method="Nelder-Mead"))
    filtered_results = [m for m in all_results if 0==m.status]
    return min(filtered_results, key = lambda x: x.fun)


# -

m = minimize_restarts(x, t, y_u, y_f, 5)
np.exp(m.x[3]) # This is the optimized value for our parameter c

# #### Step 5: Plotting the behavior for varied parameters

# The logarithms of the optimal hyperparameters are given by (arranged in $[\theta, l_x, l_t, c]$):

m.x 

# We want to plot the behavior of the nlml-function around the minimizer:

# + nbsphinx="hidden"
lin0 = np.linspace(5, 9, 200) # Set to 200
lin1 = np.linspace(0, 5, 200)
lin3 = np.linspace(0, 0.1, 200)
res0 = [nlml((q, m.x[1], m.x[2], m.x[3]), x, t, y_u, y_f, 1e-7) for q in lin0]
res1 = [nlml((m.x[0], q, m.x[2], m.x[3]), x, t, y_u, y_f, 1e-7) for q in lin1]
res2 = [nlml((m.x[0], m.x[1], q, m.x[3]), x, t, y_u, y_f, 1e-7) for q in lin1]
res3 = [nlml((m.x[0], m.x[1], m.x[2], q), x, t, y_u, y_f, 1e-7) for q in lin3]

def show_1(lin0, lin1, lin3, res0, res1, res2, res3):
    f, (ax1, ax2) = plt.subplots(ncols=2, nrows=2, figsize=(13,7))
    f.suptitle("Local behavior of the nlml around the optimum")
    ax1[0].plot(lin0, res0)
    ax1[0].set(xlabel= r"Value of $\theta$", ylabel= "nlml")
    ax1[1].plot(lin1, res1)
    ax1[1].set(xlabel= r"Value of $l_x$", ylabel= "nlml")
    ax2[0].plot(lin1, res2)
    ax2[0].set(xlabel= r"Value of $l_t$", ylabel= "nlml")
    ax2[1].plot(lin3, res3)
    ax2[1].set(xlabel= r"Value of c", ylabel= "nlml")
    plt.show()


# -

show_1(lin0, lin1, lin3, res0, res1, res2, res3);

# + nbsphinx="hidden"
lin = np.linspace(0, 10, 50)
res = [nlml((q, m.x[1], m.x[2], m.x[3]), x, t, y_u, y_f, 1e-7) for q in lin]
plt.plot(lin, res);

# + nbsphinx="hidden"
lin = np.linspace(0, 10, 50)
res = [nlml((m.x[0], q, m.x[2], m.x[3]), x, t, y_u, y_f, 1e-7) for q in lin]
plt.plot(lin, res);

# + nbsphinx="hidden"
lin = np.linspace(0, 10, 50)
res = [nlml((m.x[0], m.x[1], q, m.x[3]), x, t, y_u, y_f, 1e-7) for q in lin]
plt.plot(lin, res);

# + nbsphinx="hidden"
lin = np.linspace(-1, 1, 50)
res = [nlml((m.x[0], m.x[1], m.x[2], q), x, t, y_u, y_f, 1e-7) for q in lin]
plt.plot(lin, res);
# -

# #### Step 6: Analysis of the error

# In this section we want to analyze the error of our algorithm using two different ways and plot its time complexity.

# + nbsphinx="hidden"
res = np.zeros((5,25))
timing = np.zeros((5,25))

# Needed for L2-Norm-calculation: in columns 25 vectors for five runs
X = np.zeros((25, 25, 5))
T = X

warnings.filterwarnings("ignore")
for k in range(5):
    for n in range(25):
        start_time = time.time()
        (x, t, y_u, y_f) = get_simulated_data(n)
        
        # Storing the x and t-values
        for j in range(n):
            X[j][n][k] = x[j]
            T[j][n][k] = t[j]
        
        m = minimize(nlml, np.random.rand(4), args=(x, t, y_u, y_f, 1e-7), method="Nelder-Mead")
        res[k][n] = np.exp(m.x[3])
        timing[k][n] = time.time() - start_time
# -

# **1. Plotting the error in our estimate for c:**

# The error is given by $| c_{estimate} - c_{true} |$.

# + nbsphinx="hidden"
lin = np.linspace(8, res.shape[1] - 1, res.shape[1] - 8)
ones = np.ones(res.shape[1])
est = np.repeat(0.043, len(lin))

def show_2(lin, ones, est):
    f, (ax1, ax2) = plt.subplots(ncols=2, nrows=2, figsize=(13,7))
    to_del = np.linspace(0, 7, 8)
    
    a0 = np.delete(np.abs(res[0,:] - ones), to_del)
    ax1[0].plot(lin, a0)
    ax1[0].plot(lin, est, linestyle='dashed', color='green')
    ax1[0].set(xlabel= r"Number of data points", ylabel= "Error")
    
    a1 = np.delete(np.abs(res[1,:] - ones), to_del)
    ax1[1].plot(lin, a1)
    ax1[1].plot(lin, est, linestyle='dashed', color='green')
    ax1[1].set(xlabel= r"Number of data points", ylabel= "Error")
    
    a2 = np.delete(np.abs(res[2,:] - ones), to_del)
    ax2[0].plot(lin, a2)
    ax2[0].plot(lin, est, linestyle='dashed', color='green')
    ax2[0].set(xlabel= r"Number of data points", ylabel= "Error")
    
    a3 = np.delete(np.abs(res[3,:] - ones), to_del)
    ax2[1].plot(lin, a3)
    ax2[1].plot(lin, est, linestyle='dashed', color='green')
    ax2[1].set(xlabel= r"Number of data points", ylabel= "Error")
    
    plt.show()


# + nbsphinx="hidden"
show_2(lin, ones, est)
# -

# We ran the algorithm five times and plotted the respective outcomes in different colors:

# + nbsphinx="hidden"
lin = np.linspace(8, res.shape[1] - 1, res.shape[1] - 8)

def show_3(lin, ones, res): 
    plt.figure(figsize=(5,3))
    
    for i in range(res.shape[0]):
        to_del = np.linspace(0, 7, 8)
        a_i = np.delete(np.abs(res[i,:] - ones), to_del)
        plt.plot(lin, a_i)
        plt.suptitle('Error in our estimate for c')
        plt.ylabel('Error')
        plt.xlabel('Number of data points')

    est1 = np.repeat(0.041, len(lin))
    plt.plot(lin, est1, color='blue', linestyle='dashed')
    plt.show();


# -

show_3(lin, ones, res)

# We see that for n sufficiently large (in this case $n \geq 10$), we can assume the error to be bounded by 0.041.

# **2. Plotting the error between the solution and the approximative solution:**

# + raw_mimetype="text/latex" active=""
# Another approach of plotting the error is by calculating the difference between the approximative solution and the true solution.
#
# That is: Let $\tilde{c}$ be the parameter, resulting from our algorithm. Set $\Omega := \{(x_i, t_i) \; \vert \; x_i \in x, t_i \in t\} \subseteq [0,1] \times [0,1]$. 
#
# Then we can calculate the solution of the PDE 
# \begin{align} \label{sol}
#     \frac{d^2}{dt^2}\tilde{u}(x,t) - \tilde{c}\frac{d^2}{dx^2}\tilde{u}(x,t) = 0.  
# \end{align}
#
# and set the error to $\lVert \tilde{u}(x,t) - u(x,t) \rVert_{\Omega}$. The norm can be chosen freely.
#
# In our case, finding the solution to a given $\tilde{c}$ is not difficult. It is given by 
# \begin{align}\label{sol2}
#     \tilde{u}(x,t) = u(x,\sqrt{\tilde{c}}t) = (x-\sqrt{\tilde{c}}t)^2 + \sin(x+\sqrt{\tilde{c}}t)
# \end{align}
#
# We thus get:
# \begin{align*}
# \lVert \tilde{u}(x,t) - u(x,t) \rVert_{\Omega} = \lVert (x-\sqrt{\tilde{c}}t)^2 + \sin(x+\sqrt{\tilde{c}}t) - (x-t)^2 - \sin(x+t) \rVert_{\Omega} 
# \end{align*}
#
# With the $L^2$-norm, this is
# \begin{align*}
# (\sum_{(x_i,t_i) \in \Omega} \vert (x_i-\sqrt{\tilde{c}}t_i)^2 + \sin(x_i+\sqrt{\tilde{c}}t_i) - (x_i-t_i)^2 - \sin(x_i+t_i) \vert^2 )^{1/2}
# \end{align*}

# + raw_mimetype="text/latex" active=""
# \textit{Short proof} of $(\ref{sol2})$:
#
# We assume $\tilde{c} \geq 0$ and want to find some $\alpha \in \mathbb{R}$ such that 
# \begin{align*}
#     \frac{d^2}{dt^2}\tilde{u}(x,\alpha t) - \tilde{c}\frac{d^2}{dx^2}\tilde{u}(x,\alpha t) = 0.
# \end{align*}
#
# By setting $\alpha = \sqrt{\tilde{c}}$ we have:
# \begin{align*}
# \frac{d^2}{dt^2}\tilde{u}(x,\alpha t) - \tilde{c}\frac{d^2}{dx^2}\tilde{u}(x,\alpha t) &= \alpha^2 [\frac{d^2}{dt^2}u(x,t)](x, \alpha t) - \tilde{c}[\frac{d^2}{dx^2}u(x,t)](x, \alpha t) \\ &= \tilde{c} \left(\frac{d^2}{dt^2}u(x,t) - \frac{d^2}{dx^2}u(x,t) \right)(x, \alpha t) \stackrel{(\ref{sol})}{=} 0.
# \end{align*}
#

# + nbsphinx="hidden"
lin = np.linspace(8, res.shape[1] - 1, res.shape[1] - 8)
ones = np.ones(res.shape[1])
diff = np.ndarray(res.shape[1])

def show_4(lin, ones, res, diff):
    plt.figure(figsize=(5,3))
    
    to_del = np.linspace(0, 7, 8)
    for i in range(res.shape[0]):
        for j in range(res.shape[1]):
            diff[j] = np.linalg.norm((X[:,j,i] - np.sqrt(res[i,j])*T[:,j,i])**2 + \
                    np.sin(X[:,j,i]+np.sqrt(res[i,j])*T[:,j,i])-(X[:,j,i]-T[:,j,i])**2 - \
                    np.sin(X[:,j,i]+T[:,j,i]))
        diff_i = np.delete(diff, to_del)
        plt.suptitle('$L^2$-error in our estimate for c')
        plt.plot(lin, diff_i)
        plt.ylabel('Error')
        plt.xlabel('Number of data points')

    est = np.repeat(0.015, len(lin))
    plt.plot(lin, est, color='blue', linestyle='dashed')    
    plt.show()


# -

show_4(lin, ones, res, diff)

# The $L^2$-error is in our case bounded by 0.015 for $n \geq 10$. 

# **3. Plotting the execution time:**

# + nbsphinx="hidden"
lin = np.linspace(1, timing.shape[1], timing.shape[1])

for i in range(timing.shape[0]):
    plt.plot(lin, timing[i,:])
    plt.ylabel('Execution time in seconds')
    plt.xlabel('Number of data points')
    plt.show()

# + nbsphinx="hidden"
lin = np.linspace(1, timing.shape[1], timing.shape[1])

def show_5(lin, timing):
    plt.figure(figsize=(5,3))
    
    for i in range(timing.shape[0]):
        plt.suptitle('Execution time of our algorithm')
        plt.plot(lin, timing[i,:])
        plt.ylabel('Seconds')
        plt.xlabel('Number of data points')

    est = lin**(1.33)
    plt.plot(lin, est, color='blue', linestyle='dashed')
    plt.show()


# -

show_5(lin, timing)

# Curiously, the time complexity seems to be around $\mathcal{O}(n^{4/3})$ (blue-dashed line).
#
# Assuming an equal amount of function evaluations in the Nelder-Mead algorithm for different values of n,
# we would have been expecting a time complexity of $\mathcal{O}(n^3)$, due to the computation of the inverse of an $n\times n$-matrix in every evaluation of $\textit{nlml}$. This could probably be seen with larger values of n.

StarcoderdataJupyterScriptsDedupFiltered