#! /usr/bin/env python
# -*- coding: utf-8 -*-
""" Complete solution for the CS101 Programming Project about matrices.

This file uses the module :mod:`matrix`, and its class :py:class:`matrix.Matrix`, to do many examples of matrices and linear operations.

Examples and naming conventions for all these functions and methods are mainly inspired of `<http://docs.sympy.org/dev/modules/matrices/matrices.html>`_.


- *Date:* Saturday 18 juin 2016, 10:31:25.
- *Author:* `Lilian Besson <https://bitbucket.org/lbesson/>`_ for the `CS101 course <http://perso.crans.org/besson/cs101/>`_ at `Mahindra Ecole Centrale <http://www.mahindraecolecentrale.edu.in/>`_,  2015,
- *Licence:* `MIT Licence <http://lbesson.mit-license.org>`_.
"""

from __future__ import print_function, absolute_import  # Python 2/3 compatibility

from matrix import Matrix
# from matrix import Matrix as M
M = Matrix
from matrix import Fraction
# from matrix import Decimal
from matrix import det, rank, norm, trace, exp, inv
from matrix import eye, ones, zeros, diag
from matrix import rand_matrix, rand_matrix_float, mat_from_f
from matrix import minor, cofactor, adjugate
from matrix import PLUdecomposition


# Do some quick tests (it should be in the 'tests.py' file):
if __name__ == '__main__':
    print("\nSome examples of matrices operations:\n")
    A = Matrix([[1, 2], [3, 4]])
    print("A is:", A)
    B = Matrix([[5, 6], [7, 8]])
    print("B is:", B)

    # New style class
    print("\ntype(A) is:", type(A))
    print("A.__class__ is:", A.__class__)

    # Reading
    print("\nA[0, 1] is:", A[0, 1])
    print("B[1, 1] is:", B[1, 1])

    # Updating
    print("\nA[0, 1] becomes 0.")
    A[0, 1] = 0
    print("A[0, 1] is now:", A[0, 1])
    print("A is now:", A)

    # Comparing (using __eq__)
    print("\nA == A is:", A == A)
    print("A == B is:", A == B)

    # Sum of two matrices
    C = A + B
    print("\nA + B is:", C)
    D = B + A
    print("B + A is:", D)

    # Sum of one matrix and a constant
    E = A + 3
    print("\nA + 3 is:", E)
    F = 1 + B
    print("1 + B is:", F)

    # Operation in place
    print("\nA += 1 can be done")
    A += 1
    print("A is now:", A)

    # Difference of two matrices
    G = A - B
    print("\nA - B is:", G)
    H = B - A
    print("B - A is:", H)

    # Negation of one matrix
    G2 = -A
    print("\n-A is:", G2)
    H2 = B + (-A)
    print("B + (-A) is:", H2)

    # Positive of one matrix
    G3 = +A
    print("\n+A is:", G3)

    # Product of two matrices
    L = A * B
    print("\nA * B is:", L)
    N = B * A
    print("B * A is:", N)

    # Product of one matrix and a constant
    P = A * 3
    print("\nA * 3 is:", P)
    Q = 2 * B
    print("2 * B is:", Q)

    # Power of a matrix
    print("\nA is:", A)
    P = A ** 3
    print("\nA ** 3 is:", P)
    Q = B ** 91
    print("B ** 91 is:", Q)

    # Exponential of a matrix
    THETA = diag([3, -1])
    print("\nTHETA is:", THETA)
    P = THETA.exp()
    print("\nTHETA.exp() is:", P)
    from math import exp as mathexp
    P2 = diag([mathexp(3), mathexp(-1)])
    print("Conversely diag([exp(3), exp(-1)]) is:", P2)

    print("\nexp(diag([3, -1])) ~= diag([exp(3), exp(-1)]):", P.almosteq(P2))

    # Almost equal with a high precision
    print("\nWith epsilon = 1e-16:")
    print("exp(diag([3, -1])) ~= diag([exp(3), exp(-1)]):", P.almosteq(P2, epsilon=1e-16))

    # Remove rounding errors, then checking
    P2 = P2.map(lambda x: round(x, 10))
    P = P.map(lambda x: round(x, 10))
    print("exp(diag([3, -1])) == diag([exp(3), exp(-1)]):", P == P2)

    print("\nA is:", A)
    P = exp(A)
    print("\nexp(A) is:", P)

    # Difference of one matrix and a constant
    J = A - 3
    print("\nA - 3 is:", J)
    K = 1 - B
    print("1 - B is:", K)

    # Product element-wise of two matrices
    print("\nA is:", A, "and B is:", B)
    print("A.multiply_elementwise(B) is:", A.multiply_elementwise(B))
    print("B.multiply_elementwise(A) is:", B.multiply_elementwise(A))

    # Count and contains
    print("\nA.count(1) is:", A.count(1))
    print("eye(4).count(1) is:", eye(4).count(1))

    print("\n1 in A is:", 1 in A)
    print("0 not in B is:", 0 not in B)

    # Utility functions
    # Identity
    I = eye(3)
    print("\nI = eye(3) is:", I)
    print("I.n is:", I.n)
    print("I.m is:", I.m)
    print("len(I) is:", len(I))
    print("len(I) is:", I.shape)
    print("I.count(0) is:", I.count(0))

    # Diag
    D = diag([1, 2, 3, 4])
    print("\ndiag([1, 2, 3, 4]) is:", D)

    # Ones
    O = ones(3, 4)
    print("\nO = ones(3, 4) is:", O)
    print("I * O is:", I * O)

    # Diag with complex numbers
    D = diag([1 + 1j, 2 + 2j, 3 + 3j, 4 + 4j])
    print("\ndiag([1+1j, 2+2j, 3+3j, 4+4j]) is:", D)
    print("D.real is:", D.real)
    print("D.imag is:", D.imag)
    print("D.conjugate() is:", D.conjugate())

    # Transposition
    print("\nA is:", A)
    print("A.T is:", A.T)
    assert A.T.T == A

    # Row and column
    print("\nA is:", A)
    print("A.col(0) is:", A.col(0))
    print("A.col(1) is:", A.col(1))
    print("A.row(0) is:", A.row(0))
    print("A.row(1) is:", A.row(1))

    # map
    print("\nA is:", A)
    S = A.map(lambda x: x**2 + x + 1)
    print("For f: x -> x² + x + 1, A.map(f) is:", S)

    # iterator
    print("\nA is:", A)
    print("Iterating over values of a matrix: \nfor v in A: print v")
    for v in A:
        print(v)

    # Converting to a list or a tuple
    print("\nB is:", B)
    print("list(B) is:", list(B))
    print("tuple(B) is:", tuple(B))

    # Polynomial computation are also easy now
    U = Matrix([[2, 3], [7, -1]])
    print("\nU is:", U)
    V = 2 * (U**2) + 4 * U + eye(U.n)
    print("2×U² + 4×U + I_n =", V)

    # Random matrix
    M1 = rand_matrix(3, 2, 5)
    print("\nA random matrix of size (3, 2) and integer coefficients between -5 and 5 is:", M1)

    Mf = rand_matrix_float(2, 4, 2)
    print("\nA random matrix of size (2, 4) and float coefficients between -2 and 2 is:", Mf)

    # Object identity
    Aminus = -A
    print("\nA == -(-A) is:", A == -Aminus)
    print("A is -(-A) is:", A is -Aminus)

    # Surprises: boolean operators are working (without the need of anything else)
    A1 = Matrix([[1, 1], [1, 0]])
    print("\nA1 is:", A1)
    B1 = Matrix([[0, 0], [1, 0]])
    print("B1 is:", B1)

    print("A1 and B1 is:", A1 and B1)
    print("A1 or B1 is:", A1 or B1)
    print("not A1 is:", not A1, " ==> WARNING is not the negation element-wise!")

    print("\nWARNING: 'or' and 'and' work weirdly with numbers")
    A2 = Matrix([[1, 1], [1, 0]]) * 7
    print("\nA2 is:", A2)
    B2 = Matrix([[0, 0], [1, 0]]) * (-6)
    print("B2 is:", B2)

    print("A2 and B2 is:", A2 and B2)
    print("A2 or B2 is:", A2 or B2)
    # Example: 4 or 6 is 4, 3 and 2 is 2

    # Modulus and divmod
    print("\nA is:", A)
    print("A % 3 is:", A % 3)
    print("A % 2 is:", A % 2)

    print("\nB is:", B)
    print("B % 4 is:", B % 4)
    print("B % 2 is:", B % 2)

    # ones and zeros
    print("\nones(4) is:", ones(4))
    print("ones(2, 3) is:", ones(2, 3))

    print("\nzeros(4) is:", zeros(4))
    print("zeros(2, 3) is:", zeros(2, 3))

    # Dot product
    Z = rand_matrix(4, 2, 10)
    v = [1] * Z.m
    print("\nZ is:", Z)
    print("Vector v is:", v)
    print("Z.dot(v) is:", Z.dot(v))

    # Norm
    print("\nA is:", A)
    print("A.norm() is:", A.norm())
    print("A.norm(p=1) is:", A.norm(1), "which is like sum(A):", sum(A))
    print("A.norm(3) is:", A.norm(3))
    print("norm(A) is:", norm(A))
    print("norm(A, 42) is:", norm(A, 42))

    # Absolute value (norm max)
    U = rand_matrix(5, 5, 20)
    print("\nRandom U is:", U)
    print("abs(U) is:", abs(U))

    # Trace
    print("\neye(12) is:", eye(12))
    print("eye(12).trace() is:", eye(12).trace())
    print("trace(ones(5)) is:", trace(ones(5)))

    # Symetric or anti-symetric, diagonal or not diagonal, hermitian or not
    print("\nA is:", A)
    print("A.is_symetric is:", A.is_symetric)
    print("A.is_anti_symetric is:", A.is_anti_symetric)
    print("A.is_diagonal is:", A.is_diagonal)
    print("A.is_hermitian is:", A.is_hermitian)

    V = eye(4) * 1j
    print("\nV is:", V)
    print("V.is_square is:", V.is_square)
    print("V.is_symetric is:", V.is_symetric)
    print("V.is_anti_symetric is:", V.is_anti_symetric)
    print("V.is_diagonal is:", V.is_diagonal)
    print("V.is_hermitian is:", V.is_hermitian)

    X = Matrix([[3, 0], [1, -9]])
    print("\nX is:", X)
    print("X.trace() is:", X.trace())
    print("X.is_upper is:", X.is_upper)
    print("X.is_lower is:", X.is_lower)
    print("X.T.is_upper is:", X.T.is_upper)
    print("X.T.is_lower is:", X.T.is_lower)
    print("X.is_zero is:", X.is_zero)

    # Creating a matrix from a function of i, j
    B = mat_from_f(lambda i, j: '{},{}'.format(i, j), 13)
    print("\nB = mat_from_f(lambda i, j: '{},{}'.format(i, j), 13) is:", B)
    print("B[0, 0] is:", B[0, 0])

    # Slice reading
    print("\nSlicing with the first index: row sub-vectors.")
    print("B[1:8:2, 0] is:", B[1:8:2, 0])
    print("B[0::2, 0] is:", B[0::2, 0])
    print("B[:2, 0] is:", B[:2, 0])
    print("B[:10:3, 0] is:", B[:10:3, 0])
    print("B[::3, 0] is:", B[::3, 0])
    print("B[:, 0] is:", B[:, 0])
    print("B[3:, 0] is:", B[3:, 0])

    print("\nSlicing with the second index: column sub-vectors.")
    print("B[0, 1:8:2] is:", B[0, 1:8:2])
    print("B[0, 0::2] is:", B[0, 0::2])
    print("B[0, :2] is:", B[0, :2])
    print("B[0, :10:3] is:", B[0, :10:3])
    print("B[0, ::3] is:", B[0, ::3])
    print("B[0, :] is:", B[0, :])
    print("B[0, 3:] is:", B[0, 3:])

    print("\nSlicing with the two indeces: sub-matrices.")
    print("B[:4, 1:8:2] is:", B[:4, 1:8:2])
    print("B[::4, 0::2] is:", B[::4, 0::2])
    print("B[5:, :2] is:", B[5:, :2])
    print("B[1::2, :10:3] is:", B[1::2, :10:3])
    print("B[1:5, ::3] is:", B[1:5, ::3])
    print("B[:, :] is:", B[:, :])
    print("B[6:, 3:] is:", B[6:, 3:])

    # Slice affectation
    print("\nModifying a slice with the first index: row sub-vectors.")

    print("\nModifying a slice (with the first index) with a constant value:", 0)
    print("B[1:8:2, 0] is:", B[1:8:2, 0])
    B[1:8:2, 0] = 0
    print("B[1:8:2, 0] is:", B[1:8:2, 0])

    print("\nModifying a slice (with the first index) with a list:", [1] * len(B[0::2, 0]))
    print("B[0::2, 0] is:", B[0::2, 0])
    B[0::2, 0] = [1] * len(B[0::2, 0])
    print("B[0::2, 0] is:", B[0::2, 0])

    print("\nModifying a slice (with the first index) with a row vector:", [[6] * len(B[:10:3, 0])])
    print("B[:10:3, 0] is:", B[:10:3, 0])
    B[0::2, 0] = [[6] * len(B[:10:3, 0])]
    # WARNING: Modifying a slice with a row vector (list of 1 list) fails.
    print("B[:10:3, 0] is:", B[:10:3, 0])
    print("WARNING: Modifying a slice with a row vector (list of 1 list) fails.")
    B[0::2, 0] = 6
    print("B[:10:3, 0] is fixed:", B[:10:3, 0])

    print("\nModifying a slice (with the first index) with a column vector:", [[-3j]] * len(B[:, 0]))
    print("B[:, 2] is:", B[:, 2])
    B[:, 2] = [[-3j]] * len(B[:, 0])
    print("B[:, 2] is:", B[:, 2])

    print("\nModifying a slice (with the first index) with a row vector (as a Matrix object):", Matrix([[-9] * len(B[:10:3, 0])]))
    print("B[:10:3, 0] is:", B[:10:3, 0])
    B[0::2, 0] = Matrix([[-9] * len(B[:10:3, 0])])
    # WARNING: Modifying a slice with a row vector (list of 1 list) fails.
    print("B[:10:3, 0] is:", B[:10:3, 0])
    print("WARNING: Modifying a slice with a row vector (list of 1 list) fails.")
    B[0::2, 0] = -9
    print("B[:10:3, 0] is fixed:", B[:10:3, 0])

    print("\nModifying a slice (with the first index) with a column vector (as a Matrix object):", Matrix([[4]] * len(B[:, 0])))
    print("B[:, 2] is:", B[:, 2])
    B[:, 2] = Matrix([[4]] * len(B[:, 0]))
    print("B[:, 2] is:", B[:, 2])
#
#    print("\nModifying a slice with the second index: column sub-vectors.")
#    print("B[0, 1:8:2] is:", B[0, 1:8:2])
#    print("B[0, 0::2] is:", B[0, 0::2])
#    print("B[0, :2] is:", B[0, :2])
#    print("B[0, :10:3] is:", B[0, :10:3])
#    print("B[0, ::3] is:", B[0, ::3])
#    print("B[0, :] is:", B[0, :])
#    print("B[0, 3:] is:", B[0, 3:])
#
#    print("\nModifying a slice with the two indeces: sub-matrices.")
#    print("B[:4, 1:8:2] is:", B[:4, 1:8:2])
#    print("B[::4, 0::2] is:", B[::4, 0::2])
#    print("B[5:, :2] is:", B[5:, :2])
#    print("B[1::2, :10:3] is:", B[1::2, :10:3])
#    print("B[1:5, ::3] is:", B[1:5, ::3])
#    print("B[:, :] is:", B[:, :])
#    print("B[6:, 3:] is:", B[6:, 3:])

    # Gauss elimination
    print("\nGauss elimination examples:")

    print("\nA is:", A)
    row_echelon_A, det_A = A.gauss(det=True)
    print("A.gauss() row echelon form of A is:", row_echelon_A)
    print("det_A is:", det_A)
    print("A.det is:", A.det)
    print("A.rank is:", A.rank)
    print("det(A) is:", det(A))
    print("rank(A) is:", rank(A))

    print("\nChanging A with: A[1, :] = A[0, :]")
    A[1, :] = A[0, :]
    print("Now A is:", A)
    print("A.det is:", A.det)
    print("A.rank is:", A.rank)

    # Another test
    Z = zeros(3)
    print("\nZ = zeros(3) is:", Z)
    row_echelon_Z = Z.gauss()
    _, det_Z = Z.gauss(det=True)
    print("Z.det is:", Z.det)
    print("Z.rank is:", Z.rank)

    # Another test, with a matrix that should not fool Gauss
    Z2 = zeros(3)
    Z2[1, 2] = 1
    print("\n\nZ2 is:", Z2)
    row_echelon_Z2 = Z2.gauss(verb=True)
    _, det_Z2 = Z2.gauss(det=True)
    print("Z2.det is:", Z2.det)
    print("Z2.rank is:", Z2.rank)

    # An example from https://fr.wikipedia.org/wiki/%C3%89limination_de_Gauss-Jordan#Pseudocode
    # A0 = Matrix([[2, -1, 0], [-1, 2, -1], [0, -1, 2]])
    # print("\n\nWith this matrix A0:", A0, "we try the Gauss-Jordan algorithm.")
    # row_echelon_A0 = A0.gauss_jordan(verb=True, inv=False)
    # print("A0.gauss_jordan() row echelon form of A0 is:", row_echelon_A0)
    # print("det(A0) is:", det(A0))

    # Same example from https://fr.wikipedia.org/wiki/%C3%89limination_de_Gauss-Jordan#Pseudocode
    A0 = Matrix([[2, -1, 0], [-1, 2, -1], [0, -1, 2]])
    print("\n\nWith this matrix A0:", A0, "we try the extended Gauss-Jordan algorithm.")

    # Experimental
    # from decimal import Decimal, getcontext
    # getcontext().prec = 100
    # A0 = A0.map(Decimal)
    # print("A0.map(Decimal) is:", A0)

    A0 = A0.map(Fraction)
    print("A0.map(Fraction) is:", A0)

    row_echelon_A0, inv_A0 = A0.gauss_jordan(verb=False, inv=True, mode='f')
    print("A0.gauss_jordan(inv=True): row echelon form of A0 is:", row_echelon_A0)
    print("A0.gauss_jordan(inv=True): inverse of A0 is:", inv_A0)
    print("det(A0) is:", det(A0))

    print("\nA0 * inv_A0 is:", A0 * inv_A0)
    assert eye(A0.n) == A0 * inv_A0
    print("inv_A0 * A0 is:", inv_A0 * A0)
    assert eye(A0.n) == inv_A0 * A0
    assert 1 == det(A0) * det(inv_A0)

    assert inv_A0 == A0.inv()

    # Random matrix
    print("\n\nLooking for a non-singular random matrix R of size (4, 4):")
    R = rand_matrix(4, 4, 10)
    print("Trying", R)
    while R.det == 0:
        R = rand_matrix(4, 4, 10)
        print("Still looking... Is this R good ?", R)
    R = R.map(Fraction)
    print("R.map(Fraction) is:", R)
    print("Note: we use Fraction in order to be exact and not numerically approximative !")

    row_echelon_R, inv_R = R.gauss_jordan(verb=False, inv=True)
    print("R.gauss_jordan(inv=True): row echelon form of R is:", row_echelon_R)
    print("R.gauss_jordan(inv=True): inverse of R is:", inv_R)
    print("det(R) is:", det(R))

    print("\nR * inv_R is:", R * inv_R)
    assert eye(R.n) == R * inv_R
    print("inv_R * R is:", inv_R * R)
    assert eye(R.n) == inv_R * R

    assert inv_R == R.inv()

    assert 1 == det(R) * det(inv_R)

    # Minors, cofactors, and co-matrix
    A = Matrix([[1, 2], [3, 4]])
    A = A.map(Fraction)  # Trick to avoid rounding mistakes
    print("\nFor A:", A)
    for i in range(A.n):
        for j in range(A.m):
            print("\nFor i =", i, "and j =", j)
            print("  - The (i,j) minor is", minor(A, i, j))
            print("  - The (i,j) cofactor is", cofactor(A, i, j))
    print("So the co-matrix of A is", adjugate(A))

    # Checking A**(-1) = 1/det(A) * adjugate(A).T
    assert inv(A) == adjugate(A).T / float(det(A))

    # Second example
    A = Matrix([[-1, 4], [4, 1]])
    A = A.map(Fraction)  # Trick to avoid rounding mistakes
    print("\nFor A:", A)
    for i in range(A.n):
        for j in range(A.m):
            print("\nFor i =", i, "and j =", j)
            print("  - The (i,j) minor is", A.minor(i, j))
            print("  - The (i,j) cofactor is", A.cofactor(i, j))
    print("So the co-matrix of A is", A.adjugate())

    # Checking A**(-1) = 1/det(A) * adjugate(A).T
    assert A.inv() == A.adjugate().T.map(lambda x: Fraction(x, A.det))
#    assert A.inv() == A.adjugate().T / float(A.det)

    # Example of PLUdecomposition
    A = Matrix([[1, 2, -1, 0], [2, 4, -2, -1], [-3, -5, 6, 1], [-1, 2, 8, -2]])
    A = A.map(Fraction)
    print("\nFor A:", A)
    print(A.type())
    print("We compute the permuted LU decomposition...")
    P, L, U = PLUdecomposition(A)
    assert P.count(1) == P.n
    assert P.count(0) == P.n ** 2 - P.n
    print("\nWe have computed the permutation matrix P:", P)
    print(P.type())
    print("\nAnd the lower triangular matrix L:", L)
    print("L.type():\n", L.type())
    print("L.is_lower:", L.is_lower)  # False: NOOOOO!
    print("L.is_upper:", L.is_upper)  # False: OK
    print("\nAnd the upper triangular matrix  U:", U)
    print("U.type():\n", U.type())
    print("U.is_lower:", U.is_lower)  # False: OK
    print("U.is_upper:", U.is_upper)  # True: YES!
    print("\nSo P * A is:", P * A)
    print("So L * U is:", L * U)
    print("P * A == L * U:", P * A == L * U)

    # Done
    print("\n\nExamples are done. Are you satisfied?")
    # from sympy.matrices import Matrix as sMatrix
    # If we need to check some operations

# End of tests.py