# coding: utf-8

# # LU and upper echelon form

# In[5]:

import numpy as np
import scipy.linalg as la


# ## Part I: The Problem

# ---------------
# Here's a matrix:

# In[6]:

A = np.array([[0.,  0., -1., -1.,  0.,  0.,  0., -1.,  0.],
         [0.,  0.,  0.,  0.,  0., -1.,  0.,  1.,  0.],
         [0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
         [-1.,  0.,  1.,  0.,  0.,  0., -1.,  0.,  0.],
         [1., -1.,  0.,  1., -1.,  1.,  0.,  0., -1.],
         [0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  1.],
         [0.,  0.,  0.,  0.,  0.,  0.,  1.,  0.,  0.],
         [0.,  1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
         [0.,  0.,  0.,  0.,  1.,  0.,  0.,  0.,  0.]])


# In[12]:

P, L, U = la.lu(A)


# Check that it is actually a factorization:

# In[13]:

la.norm(A-P.dot(L).dot(U))


# Now look at `U`:

# In[10]:

U


# * What do you notice about the last two rows?
# * Would this be allowed if `U` were upper echelon?
# ----

# # Part II: Getting to Echelon Form

# In[19]:

def swap_rows(mat, i, j):
    temp = mat[i].copy()
    mat[i] = mat[j]
    mat[j] = temp


# In[36]:

def m_echelon(A):
    m, n = A.shape
    M = np.eye(m)
    U = A.copy()

    row = 0 
    for col in range(min(m, n)):
        piv_row = row + np.argmax(np.abs(U[row:, col]))
                 
        if abs(U[piv_row, col]) == 0:
            # column is exactly zero
            continue
                 
        swap_rows(U, row, piv_row)
        swap_rows(M, row, piv_row)
                 
        for el_row in range(row + 1, m):
            fac = -U[el_row, col]/U[row, col]
            U[el_row] += fac*U[row]
            M[el_row] += fac*M[row]

        row += 1 

    return M, U


# Compute M and U, and check that $MA=U$:

# In[29]:

M, U = m_echelon(A)

diff = M.dot(A)-U

print(la.norm(diff))


# Let's see if $U$ is actually in echelon form:

# In[25]:

U


# And what does $M$ look like?

# In[27]:

M


# ... not much structure here.

# --------------
# But we can still have *something* a little like LU:

# In[34]:

A - la.inv(M).dot(U)


# In[ ]: