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

# In[1]:


from numpy import array
v = array([2, 3, 5])


# In[3]:


v + 1


# In[4]:


v * 2


# In[5]:


w = array([4, 1, 6])


# In[6]:


v + w


# In[7]:


v*w #elementwise multiplication


# In[8]:


A = array([[2, 4, 5], [3, 7, 8], [-9, 0, 1]])


# In[9]:


A


# In[13]:


A + 1


# In[11]:


A * 2


# In[14]:


from scipy.linalg import det
det(A)


# In[15]:


from numpy.linalg import inv
inv(A)


# In[16]:


from numpy.linalg import cond
cond(A) #condition number of a matrix


# In[17]:


A[0, 0]


# In[18]:


A[1, 2] #second row, third column


# In[20]:


A[:, 2] #third column


# In[29]:


from numpy import zeros, ones, eye, empty, nan


# In[22]:


zeros((2, 3))


# In[26]:


empty((4, 3))


# In[30]:


nan #not-a-number (special floating point value)


# In[31]:


eye(4) #4x4 identity matrix


# In[32]:


A@eye(3) #matrix multiplication


# In[34]:


v[0:2] #elements 0 and 1 of v


# In[35]:


len(v)


# In[71]:


#solve the linear system Ux = b where U is upper triangular
def bsub(U, b):
    n = len(b)
    from numpy import empty, sum
    x = empty(n)
    for i in reversed(range(n)):
        x[i] = (b[i] - sum(U[i, (i+1):n] * x[(i+1):n])) / U[i, i]
    return x
    


# In[37]:


U = array([[2, 0, 3], [0, 1, -2], [0, 0, -4]])
b = array([17, -4, -5])


# In[51]:



x = bsub(U.astype(float), b.astype(float))
print(x)
print(U@x)


# In[47]:


U


# In[53]:


#solve the linear system Lx = b where L is lower triangular
def fsub(L, b):
    n = len(b)
    from numpy import empty, sum
    x = empty(n)
    for i in range(n):
        x[i] = (b[i] - sum(L[i, 0:i] * x[0:i])) / L[i, i]
    return x


# In[54]:


L = array([[2, 0, 0], [-1, 1, 0], [1, 2, -4]])
b = array([17, -4, -5])


# In[55]:


x = fsub(L, b)
print(x)
print(L@x)


# In[69]:


#given a linear system Ax = b, return an equivalent system Ux = c 
#such that U is upper triangular using row reduction
def rred(A, b):
    n = len(b)
    U = A.copy()
    c = b.copy()
    for i in range(n-1):
        for j in range(i+1, n):
           m = U[j, i] / U[i, i]
           U[j, (i+1):n] = U[j, (i+1):n] - m*U[i, (i+1):n]
           U[j, i] = 0
           c[j] = c[j] - m*c[i]
    return U, c
    
    


# In[67]:


A = array([[2, 0, 3], [-1, 1, -2], [1, 2, -4]])
b = array([17, -4, -5])
U, c = rred(A.astype(float), b.astype(float))


# In[72]:


x = bsub(U, c)
print(x)
print(A@x)


# In[94]:


#given a linear system Ax = b, return an equivalent system Ux = c 
#such that U is upper triangular using row reduction with pivoting
def rred_p(A, b):
    n = len(b)
    U = A.copy()
    c = b.copy()
    for i in range(n-1):
        emax = abs(U[i, i])
        prow = i
        for j in range(i+1, n):
            ej = abs(U[j, i])
            if ej > emax:
                emax = ej
                prow = j
        if prow != i: #interchange rows
            tmp = U[i, :].copy()
            U[i, :] = U[prow, :]
            U[prow, :] = tmp
            tmp = c[i]
            c[i] = c[prow]
            c[prow] = tmp
        for j in range(i+1, n):
           m = U[j, i] / U[i, i]
           U[j, (i+1):n] = U[j, (i+1):n] - m*U[i, (i+1):n]
           U[j, i] = 0
           c[j] = c[j] - m*c[i]
    return U, c


# In[95]:


Up, cp = rred_p(A.astype(float), b.astype(float))


# In[92]:


Up, cp


# In[93]:


x = bsub(Up, cp)
print(x)
print(A@x)


# In[100]:


#given a square matrix A, return its LU factorization
#L is unit lower triangular, U upper triangular, LU = A
def mylu(A):
    from numpy import shape, zeros, eye
    n = shape(A)[0]
    U = A.copy()
    L = zeros((n, n))
    for i in range(n-1):
        for j in range(i+1, n):
           m = U[j, i] / U[i, i]
           L[j, i] = m
           U[j, (i+1):n] = U[j, (i+1):n] - m*U[i, (i+1):n]
           U[j, i] = 0
           c[j] = c[j] - m*c[i]
    L = L + eye(n)
    return L, U


# In[104]:


L, U = mylu(A.astype(float))


# In[105]:


A - L@U


# In[106]:


L, U


# In[114]:


L, U = mylu(A.astype(float))
b = array([17, -4, -5])
y = fsub(L, b)
x = bsub(U, y)
x


# In[110]:


#given a square matrix A, return its LU factorization with pivoting
#L is unit lower triangular, U upper triangular, 
#P a permutation matrix
#LU = PA
def mylu_p(A):
    from numpy import shape, zeros, eye
    n = shape(A)[0]
    U = A.copy()
    L = zeros((n, n))
    P = eye(n)
    for i in range(n-1):
        emax = abs(U[i, i])
        prow = i
        for j in range(i+1, n):
            ej = abs(U[j, i])
            if ej > emax:
                emax = ej
                prow = j
        if prow != i: #interchange rows
            tmp = U[i, :].copy()
            U[i, :] = U[prow, :]
            U[prow, :] = tmp
            tmp = L[i, :].copy()
            L[i, :] = L[prow, :]
            L[prow, :] = tmp           
            tmp = P[i, :].copy()
            P[i, :] = P[prow, :]
            P[prow, :] = tmp        
        for j in range(i+1, n):
           m = U[j, i] / U[i, i]
           L[j, i] = m
           U[j, (i+1):n] = U[j, (i+1):n] - m*U[i, (i+1):n]
           U[j, i] = 0
           c[j] = c[j] - m*c[i]
    L = L + eye(n)
    return L, U, P


# In[111]:


L, U, P = mylu_p(A.astype(float))


# In[112]:


L, U, P


# In[113]:


L@U - P@A


# In[115]:


#using the LU factorization to solve a linear system
L, U, P = mylu_p(A.astype(float))
b = array([17, -4, -5])
y = fsub(L, P@b)
x = bsub(U, y)
x


# In[116]:


from numpy.linalg import solve
solve(A, b)


# In[117]:


S = array([[2, 3, 4], [1, 0, 1], [3, 3, 5]]) #singular
solve(S, b) #there should be no unique solution


# In[118]:


inv(S) #actually, the inverse does not exist -- roundoff error is giving us the inverse of a nearby matrix


# In[119]:


cond(S) #should be Inf


# In[ ]: