import numpy as np
import scipy.linalg

print("Section 1.3: #30")

coeff = (
    [[ 1, 1, 1 ],
     [ 1, 2, 2 ],
     [ 2, 3, -4]]
    )
print(
np.linalg.solve(
    coeff,
    [6, 11, 3])
)

print(
np.linalg.solve(
    coeff,
    [7, 10, 3])
)

print("Section 1.3: #32")

dimension = 3
tot = [0]*3
ct = 10000
for x in range(ct):
    P, L, U = scipy.linalg.lu(np.random.rand(dimension, dimension))
    for i in range(dimension):
        tot[i] += U[i][i]
print([x/ct for x in tot]);

print("Section 1.4: #21")

A = [[.5, .5], [.5, .5]]
B = [[1, 0], [0, -1]]
C = np.matmul(A, B)
print("A^2")
print(np.linalg.matrix_power(A, 2))
print("A^3")
print(np.linalg.matrix_power(A, 3))
print("A^k = A");
print("B^2")
print(np.linalg.matrix_power(B, 2))
print("B^3")
print(np.linalg.matrix_power(B, 3))
print("B^{2k} = B^2. B^{2k+1} = B.")
print("C^2")
print(np.linalg.matrix_power(C, 2))
print("C^3")
print(np.linalg.matrix_power(C, 3))
print("C^k = 0")

print("Section 1.4: #59")

print("A * v = [3, 4, 5]")
print("v' * v = 50")
print(f"v * A = {np.matmul([3,4,5], np.identity(3))}")

print("Section 1.4: #60")

print("Av = [4. 4. 4. 4.]")
print("Bw = [10. 10. 10. 10.]")

print("Section 1.6: #12")

print("Section 1.6: #32")

print("the inverse of the first matrix is")
print(np.linalg.inv(5*np.identity(4) - np.ones((4,4))))
print("(meaning a = .4, b = .2)")
print("ones(4)*ones(4) =")
print(np.matmul(np.ones((4,4)),np.ones((4,4))))
print("so we can solve this like the following equation: (a*eye +\n\
b*ones)(c*eye + d*ones) = ac*eye + (ad+bc+dimension*bd)*ones = eye")
print("c = 1/a. d = -b/a/(a+dimension*b), giving for the next matrix")
print("a = 1/6, d = 6/(1/6+5*-1)") # WRONG

print("Section 1.6: #47")

print("I'm not sure. This is Python. It says (in either case):")
print("numpy.linalg.LinAlgError: Singular matrix")
'''
print(np.linalg.solve(np.ones((4, 4)), np.random.rand(4, 1)))
print(np.linalg.solve(np.ones((4, 4)), np.ones((4, 1))))
'''

''' #68
dimension = 500
northwest = np.random.rand(dimension, dimension)
for x in range(0, dimension):
    for y in range(x+1, dimension):
        northwest[y][x] = 0
print(northwest)
'''
print("Section 1.6: #69")


from time import perf_counter
matrix = np.random.rand(500, 500)
start = perf_counter()
np.linalg.inv(matrix)
first_time = perf_counter() - start
print("500x500 inverse:", str(first_time) + "s")
matrix = np.random.rand(1000, 1000)
start = perf_counter()
np.linalg.inv(matrix)
second_time = perf_counter() - start
print("1000x1000 inverse:", str(second_time) + "s")
print("factor:", second_time / first_time)

print("Section 1.6: #70")

dimension = 1000
I = np.identity(dimension)
A = np.random.rand(dimension, dimension)
U = np.random.rand(dimension, dimension)
for x in range(0, dimension):
    for y in range(x+1, dimension):
        U[y][x] = 0

start = perf_counter()
np.linalg.inv(U)
print("inverse of U:", str(perf_counter() - start) + "s")
start = perf_counter()
np.linalg.solve(U, I)
print("U\I:", str(perf_counter() - start) + "s")
start = perf_counter()
np.linalg.inv(A)
print("inverse of A:", str(perf_counter() - start) + "s")
start = perf_counter()
np.linalg.solve(A, I)
print("A\I:", str(perf_counter() - start) + "s")

print("Section 1.6: #71")

def proof(dim):
    L = np.identity(dim) - np.diag([x/(1+x) for x in range(1,dim)], -1)
    print(L)
    Linv = np.linalg.inv(L)
    print("L^{-1}:")
    print(Linv)
    right = True
    for j in range(1,dim+1):
        for i in range(1,j+1):
            if not np.isclose(Linv[j-1][i-1], i/j):
                print("the assertion in the book is wrong at",i,j)
                right = False
    if right:
        print("The book is right")

proof(4)
proof(5)

print("Section 2.2: #33")

print("For the first matrix:")
A = [[1, 3, 3], [2, 6, 9], [-1, -3, 3]]
b = [1, 5, 5]
print("null space:")
print(scipy.linalg.null_space(A))
print("particular solution:")
print(np.linalg.lstsq(A, b, rcond=None)[0])
print("For the second matrix:")
A = [[1, 3, 1, 2], [2, 6, 4, 8], [0, 0, 2, 4]]
b = [1, 3, 1]
print("null space:")
print(scipy.linalg.null_space(A))
print("particular solution:")
print(np.linalg.lstsq(A, b, rcond=None)[0])

print("Section 2.2: #35")
print("TODO")
print("For the first system:")
A = [[1, 2], [2, 4], [2, 5], [3, 9]]

print("Section 2.2: #36")

print("(a)")
A = np.array([[1, 2, 1], [2, 6, 3], [0, 2, 5]])
print("Basis of the column space:")
print(scipy.linalg.orth(A))
print("Multiplying the rows by this gives zero (for this matrix, \
equivalent to [0 0 0]^T)")
print(scipy.linalg.null_space(A.transpose()))

print("(b)")
A = np.array([[1, 1, 1], [1, 2, 4], [2, 4, 8]])
print("Basis of the column space:")
print(scipy.linalg.orth(A))
print("Multiplying the rows by this gives zero (for this matrix, \
equivalent to [0 2 -1]^T)")
print(scipy.linalg.null_space(A.transpose()))

print("Section 2.3: #2")

print("Largest possible number of independent vectors is the dimension \
of the space spanned by these vectors, or the number of vectors in the \
basis, which in this case is:")
v = [[ 1,  1,  1,  0,  0,  0],
     [-1,  0,  0,  1,  1,  0],
     [ 0, -1,  0, -1,  0,  1],
     [ 0,  0, -1,  0, -1, -1]]
print(len(scipy.linalg.orth(v)))

print("Section 2.3: #5")

print("(a)")

v = [[1, 2, 3],
     [3, 1, 2],
     [2, 3, 1]]
if len(scipy.linalg.orth(v)) == 3: print("they are independent")
else: print("they are not independent")

print("(b)")

v = [[ 1,  2, -3],
     [-3,  1,  2],
     [ 2, -3,  1]]
if len(scipy.linalg.orth(v)) == 3: print("they are independent")
else: print("they are not independent")

print("Section 2.3: #13")
# probably out of spirit of the matlab hw
print("This one can be done easily without a calculator because the \
        echelon matrix has already been computed.")
print("The dimensions of all of these is 2, and the row spaces are the \
        same.")

print("Section 2.3: #16")
print("TODO")

print("Section 2.3: #18")
print("TODO")

print("Section 2.3: #24")