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

# # Backward Stability by Example
# 
# Copyright (C) 2024 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[2]:


import numpy as np
import numpy.linalg as la


# ## A numerical method for solving linear systems
# 
# Here's code for (mostly) straightforward Gauss-Jordan elimination to solve a linear system:

# In[3]:


def mylinsolve(A, b):
    n = len(A)
    Ab = np.concatenate([A, b.reshape(-1, 1)], axis=1)

    # reverse the rows
    Ab = Ab[::-1]

    for col in range(n):
        Ab[col] /= Ab[col, col]

        for row in range(n):
            if col == row:
                continue
            Ab[row] -= Ab[col]*Ab[row, col]

    return Ab[:, -1]


# Test that it solves linear systems:

# In[4]:


Atest = np.random.randn(5, 5)
btest = np.random.randn(5)

xtest = la.solve(Atest, btest)

la.norm(mylinsolve(Atest, btest) - xtest, 2) / la.norm(xtest, 2)


# ## Setting up a "bad" problem
# 
# Make an example of a very poorly conditioned matrix, called the [Vandermonde matrix](https://en.wikipedia.org/wiki/Vandermonde_matrix).

# In[5]:


n = 20
nodes = np.linspace(0, 1, n)
A = nodes.reshape(-1, 1) ** np.arange(n)


# Check the condition number:

# In[6]:


la.cond(A, 2)


# Set up a "true" solution, and a right-hand side for a linear system:

# In[7]:


xtrue = np.random.randn(n)
b = A @ xtrue


# ## Scenario 1: "Backward Stability"
# 
# Solve it computationally, call the solution `xcomp`:

# In[39]:


xcomp = mylinsolve(A, b)


# Compute the **relative forward error** (i.e. the error in the solution):

# In[40]:


la.norm(xcomp - xtrue, 2)/la.norm(xtrue, 2)


# Compute the
# 
# - **relative backward error** (assign to `bw_err`)
# - i.e. the relative error in the right-hand side
# - also known as the **relative residual**.

# In[41]:


bw_err = la.norm(A@xcomp - b, 2)/la.norm(b, 2)
bw_err


# Compare the forward error with the guarantee from the conditioning bound.

# In[42]:


bw_err * la.cond(A, 2)


# -   Repeat this scenario with `la.solve`.    
# -   Why is this not *exactly* the backward error scenario from class (hence the quotes above)?
#     
#     <details>
#     <summary>Solution</summary>
#     Because we are unable to apply the "true" inverse (i.e. matrix-vector product), i.e. a solve undisturbed by numerical error.
#     </details>

# ## Scenario 2: Input Perturbation
# 
# Now *intentionally* perturb the right-hand side and repeat the experiment

# In[28]:


delta_b = np.random.randn(n) * 1e-14 
b_perturbed = b + delta_b
xcomp_perturbed = mylinsolve(A, b_perturbed)


# Find the relative forward perturbation vs `xtrue`:

# In[29]:


la.norm(xcomp_perturbed - xtrue, 2)/la.norm(xtrue, 2)


# Compare with the perturbation bound obtained from conditioning:

# In[30]:


la.norm(delta_b, 2)/la.norm(b, 2) * la.cond(A, 2)


# In[ ]: