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

# In[2]:


from numpy import array
A = array([[3, 4, 5], [2, -1, 6]])


# In[4]:


A.shape


# In[5]:


A[1, 0]


# In[8]:


C = A.transpose()


# In[11]:


D = array([[5], [2.1], [6]])


# In[12]:


D.shape


# In[13]:


A @ D #matrix multiplication


# In[14]:


from numpy import dot
dot(A, D)


# In[15]:


2 * A


# In[17]:


E = array([[-5], [2.1], [6]])


# In[18]:


D * E #element-by-element multiplication


# In[19]:


A[0, :] #select a row of matrix


# In[20]:


A[0, 0:2] #range of columns


# In[21]:


A - 1 #subtraction of a scalar


# In[22]:


D + E #element-by-element addition


# In[25]:


from numpy import diag
diag(C) #main diagonal


# In[26]:


C.diagonal() #also works


# In[28]:


from numpy import sum
sum(C == 3)


# In[29]:


sum(D)


# In[30]:


C.max() #largest element in array


# In[31]:


C.min()


# In[33]:


D.size


# In[11]:


def gauss(A, b):
    #Gauss elimination without pivoting
    n = b.size
    for i in range(n):
        for j in range(i+1, n):
            M = A[j, i] / A[i, i]
            A[j, :] = A[j, :] - M*A[i, :]
            b[j] = b[j] - M*b[i]
    return A, b #A is now upper triangular


# In[50]:


A = array([[-1, 2, 2], [1, 1, 1], 
           [1, 3, 2]])
b = array([[8], [1], [4]])


# In[51]:


U, c = gauss(A.copy(), b.copy())


# In[52]:


U, c


# In[48]:


b.size


# In[53]:


A


# In[ ]:





# In[9]:


def bsub(U, b):
 #solve Ux = b with U upper triangular
 n = b.size
 x = b.copy() #new array to store x   
 for i in range(n, 0, -1):
    terms = U[i-1, i:n] @ x[i:n]
    x[i-1] = (b[i-1] - terms) / U[i-1, i-1]
 return x
     


# In[103]:


x = bsub(U, c)
x


# In[100]:


U @ x, c


# In[101]:


x


# In[105]:


def bsub2(U, b):
 #solve Ux = b with U upper triangular
 #with double loop (should be same)
 n = b.size
 x = b.copy() #new array to store x   
 for i in range(n, 0, -1):
    terms = 0 
    for j in range(i, n):
        terms = terms + U[i-1, j] * x[j]
    x[i-1] = (b[i-1] - terms) / U[i-1, i-1]
 return x


# In[106]:


bsub2(U, c)


# In[12]:


def mysolve(A, b):
    #solve a linear system Ax=b
    #Gauss elimination without pivoting
    U, c = gauss(A.copy(), b.copy())
    x = bsub(U, c)
    return x


# In[109]:


x = mysolve(A, b)
x


# In[28]:


def mylu(A):
    #LU factorization without pivoting
    n = A.shape[0]
    L = A.copy()
    L[:] = 0
    for i in range(n):
        L[i, i] = 1
        for j in range(i+1, n):
            M = A[j, i] / A[i, i]
            A[j, :] = A[j, :] - M*A[i, :]
            L[j, i] = M
    return A, L #A is now upper triangular


# In[119]:


from numpy import double
U, L = mylu(double(A.copy()))
U, L


# In[120]:


L @ U, A


# In[121]:


(L @ U) -  A


# In[128]:


def fsub(L, b):
 #solve Lx = b with L lower triangular
 n = b.size
 x = b.copy() #new array to store x   
 for i in range(n):
    terms = L[i, 0:i] @ x[0:i]
    x[i] = (b[i] - terms) / L[i, i]
 return x


# In[132]:


def lusolve(L, U, b):
    #solve Ax = b
    #if A = LU
    y = fsub(L, b)
    x = bsub(U, y)
    return x


# In[133]:


U, L = mylu(double(A.copy()))
x = lusolve(L, U, b)
x


# In[142]:


def myinv(A):
    #compute matrix inverse column by column
    U, L = mylu(double(A.copy()))
    X = U.copy()
    X[:] = 0
    n = A.shape[0]
    for i in range(n):
        b = zeros((n, 1))
        b[i] = 1
        x = lusolve(L, U, b)
        X[:, i] = x.flatten()
    return X


# In[143]:


X = myinv(A)


# In[144]:


X


# In[145]:


A @ X


# In[146]:


X @ A


# In[149]:


b


# In[150]:


from numpy.linalg import solve
A = array([[-1, 2, 2], [1, 1, 1], 
           [1, 3, 2]])
b = array([[8], [1], [4]])
x = solve(A, b) #uses Gauss elimination with pivoting
x


# In[151]:


from numpy.linalg import inv
inv(A) #matrix inverse


# In[152]:


from scipy.linalg import lu
P, L, U = lu(A) #LU decomposition with pivoting, L*U = P*A


# In[153]:


P, L, U


# In[154]:


L@U - P@A


# In[155]:


#a singular matrix example
S = array([[5, 3, 2], [1, -1, -4], [6, 2, -2]])


# In[156]:


inv(S)


# In[157]:


from numpy.linalg import cond
cond(S) #condition number


# In[158]:


cond(A)


# In[6]:


from numpy import argmax
def gauss_p(A, b):
    #Gauss elimination with pivoting
    n = b.size
    for i in range(n):
        p = argmax(abs(A[i:n, i])) + i
        A[[i, p]] = A[[p, i]] #interchange rows
        b[[i, p]] = b[[p, i]]
        for j in range(i+1, n):
            M = A[j, i] / A[i, i]
            A[j, :] = A[j, :] - M*A[i, :]
            b[j] = b[j] - M*b[i]
    return A, b #A is now upper triangular


# In[7]:


def mysolve_p(A, b):
    #solve a linear system Ax=b
    #Gauss elimination with pivoting
    U, c = gauss_p(A.copy(), b.copy())
    x = bsub(U, c)
    return x


# In[23]:


from numpy import array
A = array([[-1, 2, 2], [1, 1, 1], 
           [1, 3, 2]])


# In[24]:


i = 0; p = 1
A[[i, p], 0:2] = A[[p, i], 0:2]


# In[25]:


A


# In[19]:


from numpy import double
from numpy.linalg import solve
A = array([[-1, 2, 2], [1, 1, 1], 
           [1, 3, 2]])
b = array([[8], [1], [4]])
x = solve(A, b) 
x1 = mysolve(double(A.copy()), double(b.copy()))
x2 = mysolve_p(double(A.copy()), double(b.copy()))


# In[20]:


x, x1, x2


# In[26]:


def mylu_p(A):
    #LU factorization with pivoting
    n = A.shape[0]
    L = A.copy()
    L[:] = 0
    P = eye(n)
    for i in range(n):
        L[i, i] = 1
        p = argmax(abs(A[i:n, i])) + i
        A[[i, p]] = A[[p, i]]
        L[[i, p], 0:i] = L[[p, i], 0:i]
        P[[i, p]] = P[[p, i]]
        for j in range(i+1, n):
            M = A[j, i] / A[i, i]
            A[j, :] = A[j, :] - M*A[i, :]
            L[j, i] = M
    return A, L, P #A is now upper triangular


# In[21]:


from numpy import eye
P = eye(3)


# In[22]:


P


# In[30]:


A = array([[-1, 2, 2], [1, 1, 1], 
           [1, 3, 2]])
from scipy.linalg import lu
P, L, U = lu(A)
U1, L1 = mylu(double(A.copy()))
U2, L2, P2 = mylu_p(double(A.copy()))


# In[31]:


P, L, U


# In[32]:


L1, U1


# In[33]:


P2, L2, U2


# In[36]:


from numpy import sqrt, sum
def mychol(A):
    #Cholesky factorization of a symmetric matrix
    n = A.shape[0]
    L = A.copy()
    L[:] = 0   
    for i in range(n):
        L[i, i] = sqrt(A[i, i] - sum(L[i, 0:i] ** 2))
        for j in range(i+1, n):
            L[j, i] = (A[j, i] - 
                       sum(L[j, 0:i]*L[i, 0:i])) / L[i, i]
    return L
    


# In[37]:


S = [[1, 2, 3], [2, 5, 7], 
     [3, 7, 14]]
L1 = mychol(double(S.copy()))


# In[38]:


L1


# In[41]:


A.shape[0]


# In[45]:


from scipy.linalg import cholesky
L = cholesky(S, lower=True)


# In[46]:


L


# In[47]:


s = [3, 7, 14, 5]
argmax(s)


# In[ ]: