#! /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 :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 ``_. - *Date:* Tuesday, April 07, 14:05:29 2015. - *Author:* `Lilian Besson `_ for the `CS101 course `_ at `Mahindra Ecole Centrale `_, 2015, - *Licence:* `MIT Licence `_. """ #from __future__ import division from matrix import Matrix from matrix import Matrix as M from matrix import Fraction, 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 f = lambda x: x**2 + x + 1 S = A.map(f) 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])] # FIXME 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])]) # FIXME 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 xrange(A.n): for j in xrange(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 xrange(A.n): for j in xrange(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