from __future__ import division
import numpy as np
import sympy as sp
from scipy.linalg import *

"""
Linear Algebra (Matrix) Operations in Python
"""

print('Define a 2x3 Matrix [A] \n')
A = np.array(([1,2,3],      	# Define 2x3 matrix [A]
              [4,5,6]))
print('A=',A)
print('\n\n')

print('Define another 2x3 Matrix [B]')
B = np.array(([1,0,1],      	# Define 2x3 matrix [B]
              [0,1,0]))
print('B=',B)
print()

print('Add [A] and [B] and store as [C]')
C = A+B                     	# Defines C equal to addition of A & B
print('C=A+B')
print('C=',C)
print()

D = A-B                     	# A-B
print('D=A-B')
print('D=',D)
print()

E = B-A                     	# B-A
print('E=B-A')
print('E=',E)
print()

s=2                         	# Scalar

print('Multiply by a Scalar (s)')
F = s*A                     	# Multiplies scalar (s) * matrix [A]
print('F=s*A')
print('F=',F)
print()

B = np.array(([1,4,7],
              [2,5,8],
              [3,6,9]))     	# Redefines B to a 3 x 3 matrix
print('B=',B)
print()

C = A.dot(B)                	# Multiplication of [A] and [B]
print('C=A*B')               	# Note: B*A not possible as dimensions don't match
print('C=',C)
print()

print('\n\nTranspose Matrix [A]')
AT = A.transpose()          	# Transposes matrix A
print('A=',A)
print('AT=',AT)
print()

O = (A.dot(B)).transpose() - (B.transpose()).dot(A.transpose())
print('O=',O)
print()                       	# Transpose equivalency

print('\n','Identity Matrix [I]')
I = np.identity(3)          	# 3x3 Identity matrix
print('I=',I)
print

print('\n\n\nSolve System of Equations','\n')
print('Define Coefficients\n')
A = np.array(([3,2,1],
              [2,4,3],
              [1,3,5]))     	# Redefine A to system of equation coefficients
print('A=',A)
print()

print('Define Right Hand Side\n')
B = np.array(([10],
              [19],
              [22]))        	# Create right-hand-side column vector 
print('B=',B)
print()

print('Solve System of Equations using Inversion (Computationally Expensive)\n')
X = np.linalg.inv(A).dot(B)    	# Solves linear system by inversion (expensive)
print('X=',X)
print('\n')

print('Solve System of Equations using LU decomposition (Computationally Efficient)\n')
x = solve(A,B)              	# Solves linear system A*x=B
print('x=',x)
print('\n\n')


print('Define an Orthogonal Matrix [Q]')
Q = np.array(([0.8,-0.6,0],
              [0.6,0.8,0],
              [0,0,1]))     	# Orthogonal matrix
print('\nQ=',Q)
print()

print('Difference between [Q]^-1 and [Q]^T\n')
null = np.linalg.inv(Q)-Q.transpose()   # Orthogonal equivalency 
print('null=',null,'\n')
print('Note Scientfic Notation: Residual from how numbers are stored\n\n')


print('Using a for loop to populate an array:')
print('Preallocated array of just zeros')
A = np.zeros((3,3))           	# Pre-allocation of 3x3 matrix space    
print('A=',A,'\n\n')             


print('Nested for loop where:\ni = Row Number\nj = Col Number\n')
for i in range(3):
    for j in range(3):
        A[i][j] = (i+1)+2*(j+1)
        
print('A=',A)                   	# Construction of A matrix using for loops