## 2019年5月8日水曜日

### Python - Matrices - multiplying, identity, transposing, inverting(determinant, submatrix, minor, cofactor)

The Ray Tracer Challenge: A Test-Driven Guide to Your First 3D Renderer (Jamis Buck(著)、Pragmatic Bookshelf)、Chapter 3(Matrices)のPut It Together(42)を取り組んでみる。

コード

Python 3

matrices_test.py

```#!/usr/bin/env python3
from unittest import TestCase, main
from matrices import Matrix
from tuples import Tuple

class MatrixTest(TestCase):
def setUp(self):
cols = []
for i in range(4):
row = []
for j in range(4):
if i == j:
row.append(1)
else:
row.append(0)
cols.append(row)
self.identity = Matrix(cols)

def tearDown(self):
pass

def test_4_4_construct_and_inspect(self):
m = Matrix(((1, 2, 3, 4),
(5.5, 6.5, 7.5, 8.5),
(9, 10, 11, 12),
(13.5, 14.5, 15.5, 16.6)))
pairs = zip([0, 0, 1, 1, 2, 3, 3],
[0, 3, 0, 2, 2, 0, 2])
nums = [1, 4, 5.5, 7.5, 11, 13.5, 15.5]
for i, (row, col) in enumerate(pairs):
self.assertEqual(m[row][col], nums[i])

def test_2_2_construct_and_inspect(self):
m = Matrix(((-3, 5),
(1, -2)))
pairs = zip([0, 0, 1, 1],
[0, 1, 0, 1])
nums = [-3, 5, 1, -2]
for i, (row, col) in enumerate(pairs):
self.assertEqual(m[row][col], nums[i])

def test_3_3_construct_and_inspect(self):
m = Matrix([[-3, 5, 0],
[1, -2, -7],
[0, 1, 1]])
pairs = zip([0, 1, 2], [0, 1, 2])
nums = [-3, -2, 1]
for i, (row, col) in enumerate(pairs):
self.assertEqual(m[row][col], nums[i])

def test_eq(self):
a = Matrix(((1, 2, 3, 4),
(5, 6, 7., 8),
(9, 8, 7, 6),
(5, 4, 3, 2)))
b = Matrix(((1, 2, 3, 4),
(5, 6, 7., 8),
(9, 8, 7, 6),
(5, 4, 3, 2)))
self.assertEqual(a, b)

def test_ne(self):
a = Matrix(((1, 2, 3, 4),
(5, 6, 7, 8),
(9, 8, 7, 6),
(5, 4, 3, 2)))
b = Matrix(((2, 3, 4, 5),
(6, 7, 8, 9),
(8, 7, 6, 5),
(4, 3, 2, 1)))
self.assertNotEqual(a, b)

def test_mul(self):
a = Matrix(((1, 2, 3, 4),
(5, 6, 7, 8),
(9, 8, 7, 6),
(5, 4, 3, 2)))
b = Matrix(((-2, 1, 2, 3),
(3, 2, 1, -1),
(4, 3, 6, 5),
(1, 2, 7, 8)))
self.assertEqual(a * b, Matrix(((20, 22, 50, 48),
(44, 54, 114, 108),
(40, 58, 110, 102),
(16, 26, 46, 42))))

def test_mul_tuple(self):
A = Matrix(((1, 2, 3, 4),
(2, 4, 4, 2),
(8, 6, 4, 1),
(0, 0, 0, 1)))
b = Tuple(1, 2, 3, 1)
self.assertEqual(A * b, Tuple(18, 24, 33, 1))

def test_mul_identity(self):
A = Matrix([[0, 1, 2, 4],
[1, 2, 3, 8],
[2, 4, 8, 16],
[4, 8, 16, 32]])
self.assertEqual(A * self.identity, A)

a = Tuple(1, 2, 3, 4)
self.assertEqual(self.identity * a, a)

def test_transpose(self):
a = Matrix([[0, 9, 3, 0],
[9, 8, 0, 8],
[1, 8, 5, 3],
[0, 0, 5, 8]])
self.assertEqual(a.transpose(), Matrix([[0, 9, 1, 0],
[9, 8, 8, 0],
[3, 0, 5, 5],
[0, 8, 3, 8]]))
self.assertEqual(self.identity.transpose(), self.identity)

def test_determinant_2_2(self):
A = Matrix([[1, 5],
[-3, 2]])
self.assertEqual(A.determinant(), 17)

def test_submarix_3_3_2_2(self):
tests = [(Matrix(((1, 5, 0),
(-3, 2, 7),
(0, 6, -3))),
Matrix([[-3, 2],
[0, 6]]), 0, 2),
(Matrix([[-6, 1, 1, 6],
[-8, 5, 8, 6],
[-1, 0, 8, 2],
[-7, 1, -1, 1]]),
Matrix([[-6, 1, 6],
[-8, 8, 6],
[-7, -1, 1]]), 2, 1)]
for a, b, row, col in tests:
self.assertEqual(a.submatrix(row, col), b)

def test_minor_3_3(self):
A = Matrix([[3, 5, 0],
[2, -1, -7],
[6, -1, 5]])
B = A.submatrix(1, 0)
self.assertEqual(B.determinant(), 25)
self.assertEqual(A.minor(1, 0), 25)

def test_cofactor(self):
A = Matrix([[3, 5, 0],
[2, -1, -7],
[6, -1, 5]])
self.assertEqual(A.minor(0, 0), -12)
self.assertEqual(A.cofactor(0, 0), -12)
self.assertEqual(A.minor(1, 0), 25)
self.assertEqual(A.cofactor(1, 0), -25)

def test_determinant_3_3(self):
A = Matrix([[1, 2, 6],
[-5, 8, -4],
[2, 6, 4]])
tests = [(A.cofactor(0, 0), 56),
(A.cofactor(0, 1), 12),
(A.cofactor(0, 2), -46),
(A.determinant(), -196)]
for a, b in tests:
self.assertEqual(a, b)

def test_determinant_4_4(self):
A = Matrix([[-2, -8, 3, 5],
[-3, 1, 7, 3],
[1, 2, -9, 6],
[-6, 7, 7, - 9]])
tests = [(A.cofactor(0, 0), 690),
(A.cofactor(0, 1), 447),
(A.cofactor(0, 2), 210),
(A.cofactor(0, 3), 51),
(A.determinant(), -4071)]
for a, b in tests:
self.assertEqual(a, b)

def test_is_invertible(self):
A = Matrix([[6, 4, 4, 4],
[5, 5, 7, 6],
[4, -9, 3, -7],
[9, 1, 7, -6]])
self.assertTrue(A.is_invertible())
A = Matrix([[-4, 2, -2, -3],
[9, 6, 2, 6],
[0, -5, 1, -5],
[0, 0, 0, 0]])
self.assertFalse(A.is_invertible())

def test_inverse(self):
A = Matrix([[-5, 2, 6, -8],
[1, -5, 1, 8],
[7, 7, -6, -7],
[1, - 3, 7, 4]])
B = A.inverse()
self.assertEqual(A.determinant(), 532)
self.assertEqual(A.cofactor(2, 3), -160)
self.assertEqual(B[3][2], -160/532)
self.assertEqual(A.cofactor(3, 2), 105)
self.assertEqual(B[2][3], 105 / 532)
self.assertEqual(B, Matrix([[0.21805, 0.45113, 0.24060, -0.04511],
[-0.80827, -1.45677, -0.44361, 0.52068],
[-0.07895, -0.22368, -0.05263, 0.19737],
[-0.52256, -0.81391, -0.30075, 0.30639]]))

def test_mul_by_inverse(self):
A = Matrix([[3, -9, 7, 3],
[3, -8, 2, -9],
[-4, 4, 4, 1],
[-6, 5, 1, 1]])
B = Matrix([[8, 2, 2, 2],
[3, -1, 7, 0],
[7, 0, 5, 4],
[6, -2, 0, 5]])
C = A * B
self.assertEqual(C * B.inverse(), A)

if __name__ == '__main__':
main()
```

matrices.py

```from tuples import is_equal, Tuple

class Matrix:
def __init__(self, matrix: tuple):
self.matrix = matrix
self.rows = len(matrix)
self.cols = len(matrix[0])

def __repr__(self):
return f'Matrix({self.matrix})'

def __getitem__(self, y):
return self.matrix[y]

def __eq__(self, other):
m = self.rows
n = self.cols
if self.rows != other.rows or self.cols != other.cols:
return False
for row in range(self.rows):
for col in range(self.cols):
if not is_equal(self[row][col], other[row][col]):
return False
return True

def __mul__(self, other):
if isinstance(other, Tuple):
M = Matrix(((other.x,),
(other.y,),
(other.z,),
(other.w,)))
M = self * M
return Tuple(*[M[k][0] for k in range(4)])
return Matrix([[sum([self[i][k] * other[k][j]
for k in range(self.cols)])
for j in range(other.cols)]
for i in range(self.rows)])

def transpose(self):
return Matrix([[self[j][i] for j in range(self.rows)]
for i in range(self.cols)])

def determinant(self):
if self.rows == 1:
return self[0][0]
if self.rows == 2:
return self[0][0] * self[1][1] - self[0][1] * self[1][0]
return sum([self[0][j] * self.cofactor(0, j) for j in range(self.cols)])

def submatrix(self, row, col):
rows = []
for i in range(self.rows):
if i == row:
continue
cols = []
for j in range(self.cols):
if j == col:
continue
cols.append(self[i][j])
rows.append(cols)
return Matrix(rows)

def minor(self, row, col):
return self.submatrix(row, col).determinant()

def cofactor(self, row, col):
return (-1) ** (row + col) * self.minor(row, col)

def is_invertible(self):
return self.determinant() != 0

def inverse(self):
if not self.is_invertible():
raise ValueError('Matrix det == 0; not ivertible')
d = self.determinant()
return Matrix([[self.cofactor(j, i) / d for j in range(self.cols)]
for i in range(self.rows)])
```

sample1.py

```#!/usr/bin/env python3
from matrices import Matrix, Tuple

def f(i, j):
if i == j:
return 1
return 0

if __name__ == '__main__':
import random

print('1.')
identity = Matrix([[f(i, j) for j in range(4)]
for i in range(4)])
for o in [identity, identity.inverse()]:
print(o)

print('2.')
M = Matrix([[random.randrange(10) for _ in range(4)]
for _ in range(4)])
for o in [M, M * M.inverse(), M.inverse() * M]:
print(o)

print('3.')
A = M.transpose().inverse()
B = M.inverse().transpose()
for o in [M, A, B, A == B]:
print(o)

print('4.')
t = Tuple(1, 1, 1, 1)

def g(i, j, n):
if i == j:
if i == 2:
return n
return 1
return 0

print(t)
for n in range(1, 6):
identity = Matrix([[g(i, j, n) for j in range(4)]
for i in range(4)])
for o in [identity, identity * t]:
print(o)
```

```C:\Users\...>py matrices_test.py
..................
----------------------------------------------------------------------
Ran 18 tests in 0.003s

OK

C:\Users\...>py matrices_test.py -v
test_2_2_construct_and_inspect (__main__.MatrixTest) ... ok
test_3_3_construct_and_inspect (__main__.MatrixTest) ... ok
test_4_4_construct_and_inspect (__main__.MatrixTest) ... ok
test_cofactor (__main__.MatrixTest) ... ok
test_determinant_2_2 (__main__.MatrixTest) ... ok
test_determinant_3_3 (__main__.MatrixTest) ... ok
test_determinant_4_4 (__main__.MatrixTest) ... ok
test_eq (__main__.MatrixTest) ... ok
test_inverse (__main__.MatrixTest) ... ok
test_is_invertible (__main__.MatrixTest) ... ok
test_minor_3_3 (__main__.MatrixTest) ... ok
test_mul (__main__.MatrixTest) ... ok
test_mul_by_inverse (__main__.MatrixTest) ... ok
test_mul_identity (__main__.MatrixTest) ... ok
test_mul_tuple (__main__.MatrixTest) ... ok
test_ne (__main__.MatrixTest) ... ok
test_submarix_3_3_2_2 (__main__.MatrixTest) ... ok
test_transpose (__main__.MatrixTest) ... ok

----------------------------------------------------------------------
Ran 18 tests in 0.003s

OK

C:\Users\...>py sample1.py
1.
Matrix([[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1]])
Matrix([[1.0, 0.0, 0.0, 0.0], [0.0, 1.0, 0.0, 0.0], [0.0, 0.0, 1.0, 0.0], [0.0, 0.0, 0.0, 1.0]])
2.
Matrix([[9, 9, 9, 7], [4, 8, 4, 7], [3, 6, 8, 0], [7, 1, 2, 7]])
Matrix([[1.0, 4.440892098500626e-16, 8.881784197001252e-16, -8.881784197001252e-16], [-8.881784197001252e-16, 1.0, 0.0, 0.0], [0.0, -2.220446049250313e-16, 1.0, 0.0], [8.881784197001252e-16, 0.0, 0.0, 0.9999999999999996]])
Matrix([[0.9999999999999996, 3.3306690738754696e-16, -2.220446049250313e-16, 0.0], [0.0, 0.9999999999999991, -3.3306690738754696e-16, 0.0], [4.440892098500626e-16, -2.7755575615628914e-16, 1.0000000000000004, 0.0], [0.0, 2.7755575615628914e-16, 1.1102230246251565e-16, 0.9999999999999996]])
3.
Matrix([[9, 9, 9, 7], [4, 8, 4, 7], [3, 6, 8, 0], [7, 1, 2, 7]])
Matrix([[0.8, 0.5454545454545454, -0.7090909090909091, -0.6753246753246753], [-0.4, -0.09090909090909091, 0.21818181818181817, 0.35064935064935066], [-0.6, -0.45454545454545453, 0.6909090909090909, 0.4675324675324675], [-0.4, -0.45454545454545453, 0.4909090909090909, 0.4675324675324675]])
Matrix([[0.8, 0.5454545454545454, -0.7090909090909091, -0.6753246753246753], [-0.4, -0.09090909090909091, 0.21818181818181817, 0.35064935064935066], [-0.6, -0.45454545454545453, 0.6909090909090909, 0.4675324675324675], [-0.4, -0.45454545454545453, 0.4909090909090909, 0.4675324675324675]])
True
4.
Tuple(1,1,1,1)
Matrix([[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1]])
Tuple(1,1,1,1)
Matrix([[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 2, 0], [0, 0, 0, 1]])
Tuple(1,1,2,1)
Matrix([[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 3, 0], [0, 0, 0, 1]])
Tuple(1,1,3,1)
Matrix([[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 4, 0], [0, 0, 0, 1]])
Tuple(1,1,4,1)
Matrix([[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 5, 0], [0, 0, 0, 1]])
Tuple(1,1,5,1)

C:\Users\...>
```