|
#! /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):
30 ↛ exitline 30 didn't exit the module, because the condition on line 30 was never falseif __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)
470 ↛ 471line 470 didn't jump to line 471, because the condition on line 470 was never true 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
|
