import numpy as np
import scipy as sp
import numpy.linalg as npalg
import scipy.linalg as spalg

# Better not use matrix type, instead use ndarray
# Q: real orthogonal matrix
# R: upper-triangle matrix
# S: posiive definite symmetric matrix

A  = np.array([[2, 5], 
               [1, 3]])
B  = np.array([[-1], 
               [ 1]])
S  = np.array([[ 2.0, -1.0]
             , [-1.0,  5.0]])
V1 = np.array([1, 1, 2])
V2 = np.array([2, 1, 3])

print("A:\n", A)
print("B:\n", B)
print("S:\n", S)
print("V1:\n", V1)
print("V2:\n", V2)
print("")

# Inner product: Use np.inner, .dot, .matmul for 1D list, ndarray
inner = np.inner(V1, V2)
#inner = np.dot(V1, V2)
print("np.inner V1 dot V2 = {}".format(inner))
print("")

# Outer product
V3 = np.cross(V1, V2)
print("np.cross V1 x V2:")
print(V3)
print("")

# Inverse matrix: linalg.inv for list
A_i = npalg.inv(A)
check = A @ A_i
print("np.inv A^-1:")
print(A_i)
print("A*A.i:")
print(check)
print("")

# Determinant
arr = np.array([[0, 1], [2, 3]])
det = np.linalg.det(arr)
print("determinant for A = {}".format(det))

# Eigen problem
lA, vA = npalg.eig(A)
print("nympy.linalg.eig(A):")
print("Eigen values: ", lA)
print("Eigen vectors:")
print(vA)
print("");

# Solve simultaneous linear equations
print("Solve simultaneous linear equations AX = B:")
X = np.linalg.solve(A, B)
print(X)
print("")

# LU decomposition
print("LU decomposition A = PLU:")
print("A:\n", A)
P, L, U = spalg.lu(A)
check = P @ L @ U
print("P:\n", P)
print("L:\n", L)
print("U:\n", U)
print("P@L@U:\n", check)
print("")

# Cholesky decomosition: S = LL^T
print("Cholesky decompositioin S = LL^T:")
L = npalg.cholesky(S)
check = L @ L.T
print("S:\n", S)
print("L:\n", L)
print("L*L^T:\n", check)
print("")

# QR decomposition: A = QR
# good for anomarly matrix
Q, R = spalg.qr(A)
check = Q @ R
print("QR decomposition:")
print("A:\n", A)
print("Q:\n", Q)
print("R:\n", R)
print("Q@R:\n", check)