#-----------------------------------------------------------
# Lecture 12 -- Example -- Linear Equations
#-----------------------------------------------------------
#
#-----------------------------------------------------------
# A. Matrix Operations
#-----------------------------------------------------------
#
# We can do a number of matrix operations in Python using the functions in
# numpy.
#
# List of functions:
# * addition/subtraction: <var A> + <var B>, <var A> - <var B>
# * scalar multiplication: <scalar> * <var A>
# * matrix/array multiplication: np.dot(<var A>, <var B>)
# * transpose: <var A>.T
# * norm: np.linalg.norm(<var A>)
# * inverse: np.linalg.inv(<var A>)
#
# Additional information at the link:
#   https://docs.scipy.org/doc/numpy/reference/routines.linalg.html

#%%
import numpy as np

A = np.array([[-2, 4],[5, -1]])
B = np.array([[3, -4],[7, 2]])
s = 4
b = np.array([1, 3])

print('A = \n', A)
print('B = \n', B)
print('s = ', s)
print('b = ', b)

#%% Addition and subtraction

C = A+B
print('\nA + B = \n', C)

C = A-B
print('A - B = \n', C)

#%% Scalar multiplication

C = s*A
print('\ns * A = \n', C)

#%% Matrix/array multiplication

C = np.dot(A, B)
print('\nA . B = \n', C)

x = np.dot(A, b)
print('A . b = ', x)

#%% Norm
Anorm = np.linalg.norm(A)
print('\nnorm(A) = ', Anorm)

#%% Transpose
print('\nA transpose = \n', A.T)

#%%
#-----------------------------------------------------------
# B. Solution to Linear Equations
#-----------------------------------------------------------
#
# To solve Ax = b:
# * define matrix A and array b
# * Use either:
#       x = np.linalg.solve( A, b )
#   or you can also use
#       Ainv = np.linalg.inv(A)
#       x = np.dot(Ainv, b)
#
# * np.lingalg.solve is better than np.lingalg.inv() for large matrices. This 
#   is because solve() can use iterative methods, but inv() is stuck 
#   doing Gauss Elimination!

#%% Solve

A = np.array([[-2, 4],[5, -1]])
b = np.array([1, 3])

x = np.linalg.solve(A, b)

print('\nsolution, x = ', x)
print('b = : ', b)
print('check answer: ', np.dot(A, x))

#%% Inverse

Ainv = np.linalg.inv(A)
x = np.dot(Ainv, b)

print('\nsolution, x = ', x)
print('b = : ', b)
print('check answer: ', np.dot(A, x))

#%% Practice 

# Solve the linear system, (A+B).x = C.b
# 
# where A = 
# [  0  1  0  0  0 ]
# [  1  0  1  0  0 ]
# [  0  1  0  1  0 ]
# [  0  0  1  0  1 ]
# [  0  0  0  1  0 ]
##
# where B = 
# [ -2  0  0  0  0 ]
# [  0 -2  0  0  0 ]
# [  0  0 -2  0  0 ]
# [  0  0  0 -2  0 ]
# [  0  0  0  0 -2 ]
#
# where C = 
# [  1  0  0  0  5 ]
# [  0  2  0  0  0 ]
# [  0  3  1  0  0 ]
# [  0  0  0  4  0 ]
# [  0  4  0  0  5 ]
#
# and b = 
# [   2  ]
# [  1/2 ]
# [ -1/2 ]
# [  1/4 ]
# [ -1/5 ]