Create a gist now

Instantly share code, notes, and snippets.

anonymous /matrix.py
Created Mar 20, 2009

What would you like to do?
An implementation of Matrix and Vector classes for Python from Ruby's library "matrix.rb"
#!/usr/bin/env python
#
#-- matrix.py v 1.00
# This library is converted from Ruby to Python by CanI.
#
#*-- Original Version from Ruby version
# v 1.13 2001/12/09 14:22:23 by Keiju ISHITSUKA
#**-- Original Version from Smalltalk-80 version
# on July 23, 1985 at 8:37:17 am by Keiju ISHITSUKA
#
# This library's documentations are almost from Ruby source.
# If you'd want more detail usage, please see Ruby's documentation.
'''An implementation of Matrix and Vector classes'''
import math
class Matrix :
# class methods
'''
Creates a matrix where +rows+ is an array of arrays, each of which is a row
to the matrix. If the optional argument +copy+ is false, use the given
arrays as the internal structure of the matrix without copying.
>>>Matrix.rows([[25, 93], [-1, 66]])
=> 25 93
-1 66
'''
@classmethod
def rows(cls, rows, copy = True) :
return cls(rows, copy)
'''
Creates a matrix using +columns+ as an array of column vectors.
>>>Matrix.columns([[25, 93], [-1, 66]])
=> 25 -1
93 66
'''
@classmethod
def columns(cls, columns) :
rows = [[columns[j][i] for j in range(len(columns))] for i in range(len(columns[0]))]
return cls.rows(rows, False)
'''
Creates a matrix where the diagonal elements are composed of +values+.
>>>Matrix.diagonal(9, 5, -3)
=> 9 0 0
0 5 0
0 0 -3
'''
@classmethod
def diagonal(cls, *values) :
size = len(values)
rows = []
for j in range(0, size) :
row = [0 for _ in range(size)]
row[j] = values[j]
rows.append(row)
return cls.rows(rows, False)
'''
Creates an +n+ by +n+ diagonal matrix where each diagonal element is
+value+.
>>>Matrix.scalar(2, 5)
=> 5 0
0 5
'''
@classmethod
def scalar(cls, n, value) :
return cls.diagonal(*[value for _ in range(n)])
'''
Creates an +n+ by +n+ identity matrix.
>>>Matrix.identity(2)
=> 1 0
0 1
'''
@classmethod
def identity(cls, n) :
return cls.scalar(n, 1)
#alias
@classmethod
def unit(cls, n) :
return cls.identity(n)
@classmethod
def I(cls, n) :
return cls.identity(n)
'''
Creates an +n+ by +n+ zero matrix.
>>>Matrix.zero(2)
=> 0 0
0 0
'''
@classmethod
def zero(cls, n) :
return cls.scalar(n, 0)
'''
Creates a single-row matrix where the values of that row are as given in
+row+.
>>>Matrix.row_vector([4,5,6])
=> 4 5 6
'''
@classmethod
def row_vector(cls, row) :
if isinstance(row, Vector) :
return cls.rows([row.to_a()], False)
elif isinstance(row, list) :
return cls.rows([row])
else :
return cls.rows([[row]], False)
'''
Creates a single-column matrix where the values of that column are as given
in +column+.
>>>Matrix.column_vector([4,5,6])
=> 4
5
6
'''
@classmethod
def column_vector(cls, column) :
if isinstance(column, Vector) :
return cls.columns([list(column)])
elif isinstance(column, list) :
return cls.columns([column])
else :
return cls.columns([[column]])
# instance methods
'''
Initializes a matrix where +rows+ is an array of arrays, each of which is a row
to the matrix. If the optional argument +copy+ is false, use the given
arrays as the internal structure of the matrix without copying.
>>>Matrix([[25, 93], [-1, 66]])
=> 25 93
-1 66
'''
def __init__(self, rows, copy = True) :
if copy :
self.rows = [row for row in rows]
else :
self.rows = rows
'''
Returns element (+i+,+j+) of the matrix. That is: row +i+, column +j+.
>>>Matrix([[25, 93], [-1, 66]])[0,1]
=> 93
'''
def __getitem__(self, key) :
return self.rows[key[0]][key[1]]
#alias
def element(self, key) :
return self.__getitem__(key)
def component(self, key) :
return self.__getitem__(key)
'''
Sets element (+i+,+j+) of the matrix. That is: row +i+, column +j+.
>>>Matrix([[25, 93], [-1, 66]])[0,1] = 53
=> 25 53
-1 66
'''
def __setitem__(self, key, value) :
self.rows[key[0]][key[1]] = value
#alias
def __set_element(self, key, value) :
self.__setitem__(key, value)
def __set_component(self, key, value) :
self.__setitem__(key, value)
'''
Returns the number of rows.
>>>Matrix([[25, 93], [-1, 66]]).row_size()
=> 2
'''
def row_size(self) :
return len(self.rows)
'''
Returns the number of columns. Note that it is possible to construct a
matrix with uneven columns (e.g. Matrix[ [1,2,3], [4,5] ]), but this is
mathematically unsound. This method uses the first row to determine the
result.
>>>Matrix([[25, 93], [-1, 66]]).column_size()
=> 2
'''
def column_size(self) :
return len(self.rows[0])
'''
Returns row vector number +i+ of the matrix as a Vector (starting at 0 like
an array) Note that it is changed from Ruby, even if you'd want to use it for iterater,
it can't iterated.
>>>Matrix([[25, 93], [-1, 66]]).row(0)
=> Vector[25, 93]
'''
def row(self, i) :
return Vector.elements(self.rows[i])
'''
Returns column vector number +j+ of the matrix as a Vector (starting at 0
like an array) Note that it is changed from Ruby, too. Please see comment above one.
'''
def column(self, j) :
col = [self.rows[i][j] for i in range(self.row_size())]
return Vector.elements(col, False)
'''
Returns a matrix that is the result of function of the given argument over all
elements of the matrix. Note that it is changed from Ruby, too. In Python, can't give
block, so give function. If you give no argument, it returns itself.
>>>Matrix([ [1,2], [3,4] ]).collect(lambda e : e**2)
=> 1 4
9 16
'''
def collect(self, func = lambda e : e) :
rows = [[func(e) for e in row] for row in self.rows]
return rows
#alias
def map(self, func = lambda e: e) :
yield self.collect(func)
'''
Returns a section of the matrix. The parameters are either:
* start_row, nrows, start_col, ncols; OR
* col_range, row_range
>>>Matrix.diagonal(9, 5, -3).minor(range(2), range(3))
=> 9 0 0
0 5 0
'''
def minor(self, *param) :
if len(param) == 2 :
from_row = param[0][0]
size_row = param[0][-1] + 1
from_col = param[1][0]
size_col = param[1][-1] + 1
elif len(param) == 4 :
from_row = param[0]
size_row = param[1]
from_col = param[2]
size_col = param[3]
else :
raise TypeError()
rows = [row[from_col:size_col] for row in self.rows[from_row:size_row]]
return Matrix.rows(rows, False)
'''
Returns +true+ if this is a regular matrix.
Note that it is changed from Ruby, too. Changed method's name.
'''
def isregular(self) :
return self.issquare() and self.rank == self.column_size()
'''
Returns +true+ is this is a singular (i.e. non-regular) matrix.
Note that it is changed from Ruby, too. Please see comment above one.
'''
def issingular(self) :
return not self.isregular()
'''
Returns +true+ is this is a square matrix. See note in column_size about this
being unreliable, though. Note that it is changed from Ruby, too. Please see comment
above one.
'''
def issquare(self) :
return self.column_size() == self.row_size()
'''
Returns +true+ if and only if the two matrices contain equal elements.
'''
def __eq__(self, other) :
if not isinstance(other, Matrix) :
return False
return other.compare_by_row_vectors(self.rows)
#alias Note that it is changed from Ruby, too. Please see comment above two.
def iseql(self, other) :
return self.__eq__(other)
def __ne__(self, other) :
return not self.__eq__(other)
'''
Not really intended for general consumption.
'''
def compare_by_row_vectors(self, rows) :
if not len(self.rows) == len(rows) :
return False
for i in range(len(self.rows)) :
if not self.rows[i] == rows[i] :
return False
return True
'''
Returns a clone of the matrix, so that the contents of each do not reference
identical objects.
'''
def clone(self) :
return Matrix.rows(self.rows)
'''
Returns a hash-code for the matrix.
'''
def hash(self) :
value = 0
for row in self.rows :
for e in row :
value ^= hash(e)
return value
'''
Matrix multiplication.
>>>Matrix([[2,4], [6,8]]) * Matrix.identity(2)
=> 2 4
6 8
'''
def __mul__(self, other) :
if isinstance(other, int) or isinstance(other, long) or isinstance(other, float) :
rows = [[e * other for e in row] for row in self.rows]
return Matrix.rows(rows, False)
elif isinstance(other, Vector) :
other = Matrix.column_vector(other)
r = self * other
return r.column(0)
elif isinstance(other, Matrix) :
if self.column_size() != other.row_size() :
raise ValueError()
rows = []
for i in range(len(self.rows)) :
columns = []
for j in range(len(other.rows[0])) :
vij = 0
for k in range(len(self.rows[0])) :
vij += self[i, k] * other[k, j]
columns.append(vij)
rows.append(columns)
return Matrix.rows(rows, False)
else :
raise TypeError()
'''
Matrix addition.
>>>Matrix.scalar(2,5) + Matrix([[1,0], [-4,7]])
=> 6 0
-4 12
'''
def __add__(self, other) :
if isinstance(other, Vector) :
other = Matrix.column_vector(other)
elif isinstance(other, Matrix) :
pass
else :
raise TypeError()
if self.row_size() != other.row_size() or self.column_size() != other.column_size() :
raise ValueError()
rows = [[self[i,j] + other[i,j] for j in range(self.column_size())] for i in range(self.row_size())]
return Matrix.rows(rows, False)
'''
Matrix subtraction.
>>>Matrix([[1,5], [4,2]]) - Matrix([[9,3], [-4,1]])
=> -8 2
8 1
'''
def __sub__(self, other) :
if isinstance(other, Vector) :
other = Matrix.column_vector(other)
elif isinstance(other, Matrix) :
pass
else :
raise TypeError()
if self.row_size() != other.row_size() or self.column_size() != other.column_size() :
raise ValueError()
rows = [[self[i,j] - other[i,j] for j in range(self.column_size())] for i in range(self.row_size())]
return Matrix.rows(rows, False)
'''
Matrix division (multiplication by the inverse).
>>>Matrix([[7,6], [3,9]]) / Matrix([[2,9], [3,1]])
=> -7 1
-3 -6
'''
def __div__(self, other) :
if isinstance(other, int) or isinstance(other, long) or isinstance(other, float) :
rows = [[e / other for e in row] for row in self.rows]
return Matrix.rows(rows, False)
elif isinstance(other, Matrix) :
return self * other.inverse()
else :
raise TypeError()
'''
Returns the inverse of the matrix.
>>>Matrix([[1, 2], [2, 1]]).inverse()
=> -1 1
0 -1
'''
def inverse(self) :
if not self.issquare() :
raise ValueError()
return Matrix.I(self.row_size()).inverse_from(self)
#alias
def inv(self) :
return self.inverse()
def inverse_from(self, src) :
size = self.row_size()
a = src.to_a()
for k in range(size) :
i = k
akk = abs(a[k][k])
for j in range(k+1, size) :
v = abs(a[j][k])
if v > akk :
i = j
akk = v
if akk == 0 :
raise ValueError()
if i != k :
a[i], a[k] = a[k], a[i]
self.rows[i], self.rows[k] = self.rows[k], self.rows[i]
akk = a[k][k]
for i in range(size) :
if i == k :
continue
q = a[i][k]/akk
a[i][k] = 0
for j in range(k+1, size) :
a[i][j] -= a[k][j] * q
for j in range(size) :
self.rows[i][j] -= self.rows[k][j] * q
for j in range(k+1, size) :
a[k][j] = a[k][j]/akk
for j in range(size) :
self.rows[k][j] = self.rows[k][j]/akk
return self
'''
Matrix exponentiation. Defined for integer powers only. Equivalent to
multiplying the matrix by itself N times.
>>>Matrix([[7,6], [3,9]]) ** 2
=> 67 96
48 99
'''
def __pow__ (self, other) :
if isinstance(other, int) or isinstance(other, long) :
x = self
if other <= 0 :
x = self.inverse()
if other == 0 :
return Matrix.identity(self.column_size())
other = -other
z = x
n = other - 1
while n != 0 :
while True:
div, mod = divmod(n, 2)
if mod != 0 :
break
x = x * x
n = div
z *= x
n -= 1
return z
else :
raise TypeError()
'''
Returns the determinant of the matrix. If the matrix is not square, the
result is 0. This method's algorism is Gaussian elimination method
and using Numeric#quo(). Beware that using Float values, with their
usual lack of precision, can affect the value returned by this method. Use
Rational values or Matrix#det_e instead if this is important to you.
>>>Matrix([[7,6], [3,9]]).determinant()
=> 63
'''
def determinant(self) :
if not self.issquare() :
return 0
size = self.row_size() - 1
a = self.to_a()
det = 1
k = 0
while True :
akk = a[k][k]
if akk == 0 :
i = k
while True :
i += 1
if i > size :
return 0
if a[i][k] != 0 :
break
a[i], a[k] = a[k], a[i]
akk = a[k][k]
det *= -1
for i in range(k+1, size+1) :
q = a[i][k]/akk
for j in range(k+1, size+1) :
a[i][j] -= a[k][j] * q
det *= akk
k += 1
if not k <= size :
break
return det
#alias
def det(self) :
return self.determinant()
'''
Returns the determinant of the matrix. If the matrix is not square, the
result is 0. This method's algorism is Gaussian elimination method.
This method uses Euclidean algorism. If all elements are integer,
really exact value. But, if an element is a float, can't return
exact value.
>>>Matrix([[7,6], [3,9]]).determinant_e()
=> 63
'''
def determinant_e(self) :
if not self.issquare() :
return 0
size = self.row_size() - 1
a = self.to_a()
det = 1
k = 0
while True :
if a[k][k] == 0 :
i = k
while True :
i += 1
if i > size :
return 0
if a[i][k] != 0 :
break
a[i], a[k] = a[k], a[i]
det *= -1
tmp = range(k+1, size+1)
for i in tmp :
q = a[i][k]/a[k][k]
for j in range(k, size+1) :
a[i][j] -= a[k][j] * q
if a[i][k] != 0 :
a[i], a[k] = a[k], a[i]
det *= -1
tmp.insert(i+1, i)
det *= a[k][k]
k += 1
if k > size :
break
return det
#alias
def det_e(self) :
return self.determinant_e()
'''
Returns the rank of the matrix. Beware that using Float values,
probably return faild value. Use Rational values or Matrix#rank_e
for getting exact result.
>>>Matrix([[7,6], [3,9]]).rank()
=> 2
'''
def rank(self) :
if self.column_size() > self.row_size() :
a = list(self.transpose)
a_column_size = self.row_size()
a_row_size = self.column_size()
else :
a = self.to_a()
a_column_size = self.column_size()
a_row_size = self.row_size()
rank = 0
k = 0
while True :
akk = a[k][k]
if akk == 0 :
i = k
exists = True
while True :
i += 1
if i > a_column_size - 1 :
exists = False
break
if a[i][k] != 0 :
break
if exists :
a[i], a[k] = a[k], a[i]
akk = a[k][k]
else :
i = k
exists = True
while True :
i += 1
if i > a_row_size - 1 :
exists = False
break
if a[k][i] != 0 :
break
if exists :
for j in range(k, a_column_size) :
a[j][k], a[j][i] = a[j][i], a[j][k]
akk = a[k][k]
else :
continue
for i in range(k+1, a_row_size) :
q = a[i][k]/akk
for j in range(k+1, a_column_size) :
a[i][j] -= a[k][j] * q
rank += 1
k += 1
if k > a_column_size - 1 :
break
return rank
'''
Returns the rank of the matrix. This method uses Euclidean
algorism. If all elements are integer, really exact value. But, if
an element is a float, can't return exact value.
>>>Matrix([[7,6], [3,9]]).rank_e()
=> 2
'''
def rank_e(self) :
a = self.to_a()
a_column_size = self.column_size()
a_row_size = self.row_size()
pi = 0
for j in range(a_column_size) :
i = [i0 for i0 in range(pi, a_row_size) if a[i0][j] != 0 ][0]
if i != [] :
if i != pi :
a[pi], a[i] = a[i], a[pi]
tmp = range(pi+1, a_row_size)
for k in tmp :
q = a[k][j]/a[pi][j]
for j0 in range(pi, a_column_size) :
a[k][j0] -= q * a[pi][j0]
if k > pi and a[k][j] != 0 :
a[k], a[pi] = a[pi], a[k]
tmp.insert(k+1, k)
pi += 1
return pi
'''
Returns the trace (sum of diagonal elements) of the matrix.
>>>Matrix([[7,6], [3,9]]).trace()
=> 16
'''
def trace(self) :
tr = 0
for i in range(0, self.column_size()) :
tr += self.rows[i][i]
return tr
#alias
def tr(self) :
return self.trace()
'''
Returns the transpose of the matrix.
>>>Matrix([[1,2], [3,4], [5,6]])
=> 1 2
3 4
5 6
>>>Matrix([[1,2], [3,4], [5,6]]).transpose()
=> 1 3 5
2 4 6
'''
def transpose(self) :
return Matrix.columns(self.rows)
#alias
def t(self) :
return self.transpose()
def coerce(self, other) :
if isinstance(other, int) or isinstance(other, long) or isinstance(other, float) :
return Matrix.Scalar(other), self
else :
raise TypeError()
'''
Returns an array of the row vectors of the matrix. See Vector.
'''
def row_vectors(self) :
rows = [self.row(i) for i in range(self.row_size())]
return rows
'''
Returns an array of the column vectors of the matrix. See Vector.
'''
def column_vectors(self) :
columns = [self.column(i) for i in range(self.column_size())]
return columns
'''
Returns an array of arrays that describe the rows of the matrix.
'''
def to_a(self) :
return [[e for e in row] for row in self.rows]
def elements_to_f(self) :
return [float(e) for e in self.collect()]
def elements_to_i(self) :
return [int(e) for e in self.collect()]
def __str__(self) :
return "Matrix[" + ", ".join(["["+", ".join([str(e) for e in row])+"]" for row in self.rows])+"]"
def __repr__(self) :
return "Matrix"+repr(self.rows)
class Scalar (long) :
'''Private Class'''
def __init__(self, value) :
self.value = value
def __add__(self, other) :
if isinstance(other, int) or isinstance(other, long) or isinstance(other, float) :
return Matrix.Scalar(self.value + other)
elif isinstance(other, Matrix) :
return Matrix.Scalar(self.value + other.value)
else :
raise TypeError()
def __sub__(self, other) :
if isinstance(other, int) or isinstance(other, long) or isinstance(other, float) :
return self.Scalar(self.value - other)
elif isinstance(other, Matrix.Scalar) :
return Matrix.Scalar(self.value - other.value)
else :
raise TypeError()
def __mul__(self, other) :
if isinstance(other, int) or isinstance(other, long) or isinstance(other, float) :
return Matrix.Scalar(self.value * other)
elif isinstance(other, Vector) or isinstance(other, Matrix) :
return [self.value * e for e in other.collect()]
else :
raise TypeError()
def __div__(self, other) :
if isinstance(other, int) or isinstance(other, long) or isinstance(other, float) :
return Matrix.Scalar(self.value / other)
elif isinstance(other, Matrix) :
return self * other.inverse()
else :
raise TypeError()
def __pow__ (self, other) :
if isinstance(other, int) or isinstance(other, long) or isinstance(other, float) :
return Matrix.Scalar(self.value ** other)
else :
raise TypeError()
class Vector :
# class methods
'''
Creates a vector from an Array. The optional second argument specifies
whether the array itself or a copy is used internally.
'''
@classmethod
def elements(cls, array, copy = True) :
return cls(array, copy)
# instance methods
def __init__(self, array, copy = True) :
self.init_elements(array, copy)
def init_elements(self, array, copy) :
if copy :
self.elements = [val for val in array]
else :
self.elements = array
'''
Returns element number +i+ (starting at zero) of the vector.
'''
def __getitem__(self, key) :
return self.elements[key]
#alias
def element(self, key) :
return self.__getitem__(key)
def component(self, key) :
return self.__getitem__(key)
def __setitem__(self, key, value) :
self.elements[key]= value
#alias
def __set_element(self, key, value) :
self.__setitem__(key, value)
def __set_component(self, key, value) :
self.__setitem__(key, value)
'''
Returns the number of elements in the vector.
'''
def size(self) :
return len(self.elements)
#alias
def __len__(self) :
return self.size()
'''
Iterate over the elements of this vector and +v+ in conjunction.
'''
def each2(self, v) :
if self.size() != len(v) :
raise ValueError()
for i in range(self.size()) :
yield self.elements[i], v[i]
'''
Collects over the elements of this vector and +v+ in conjunction.
'''
def collect2(self, v) :
if self.size() != len(v) :
raise ValueError()
for i in range(self.size()) :
yield self.elements[i], v[i]
'''
Returns +true+ iff the two vectors have the same elements in the same order.
'''
def __eq__(self, other) :
if not isinstance(other, Vector) :
return False
return other.compare_by(self.elements)
#alias
def iseql(self, other) :
return self.__eq__(other)
def compare_by(self, elements) :
return self.elements == elements
'''
Return a copy of the vector.
'''
def clone(self) :
return Vector.elements(self.elements)
'''
Return a hash-code for the vector.
'''
def hash(self) :
return hash(self.elements)
'''
Multiplies the vector by +x+, where +x+ is a number or another vector.
'''
def __mul__(self, other) :
if isinstance(other, int) or isinstance(other, long) or isinstance(other, float) :
els = [e * other for e in self.elements]
return Vector.elements(els, False)
elif isinstance(other, Matrix) :
return Matrix.column_vector(self) * x
else :
raise TypeError()
'''
Vector addition.
'''
def __add__(self, other) :
if isinstance(other, Vector) :
if self.size() != other.size() :
raise ValueError()
els = [v1+v2 for v1, v2 in self.collect2(other)]
return Vector.elements(els, False)
elif isinstance(other, Matrix) :
return Matrix.column_vector(self) + other
else :
raise TypeError()
'''
Vector subtraction.
'''
def __sub__(self, other) :
if isinstance(other, Vector) :
if self.size() != other.size() :
raise ValueError()
els = [v1-v2 for v1, v2 in self.collect2(other)]
return Vector.elements(els, False)
elif isinstance(other, Matrix) :
return Matrix.column_vector(self) - other
else :
raise TypeError()
'''
Returns the inner product of this vector with the other.
>>>Vector([4,7]).inner_product(Vector([10,1]))
=> 47
'''
def inner_product(self, v) :
if self.size() != v.size() :
raise ValueError()
p = 0
for v1, v2 in self.collect2(v) :
p += v1 * v2
return p
'''
Like Ruby's Array#collect.
Note that it is changed from Ruby, too. This change is like Matrix#collect.
'''
def collect(self, func = lambda v : v) :
return [func(v) for v in self.elements]
#alias
def map(self) :
yield self.collect()
def map2(self, v) :
for v1, v2 in self.collect2(v) :
yield v1, v2
'''
Returns the modulus (Pythagorean distance) of the vector.
>>>Vector([5,8,2]).r()
=> 9.643650761
'''
def r(self) :
v = 0
for e in self.elements :
v += e*e
return math.sqrt(v)
'''
Creates a single-row matrix from this vector.
'''
def covector(self) :
return Matrix.row_vector(self)
'''
Returns the elements of the vector in an array.
'''
def to_a(self) :
return [val for val in self.elements]
def elements_to_f(self) :
return [float(e) for e in self.collect()]
def elements_to_i(self) :
return [int(e) for e in self.collect()]
def coerce(self, other) :
if isinstance(other, int) or isinstance(other, long) or isinstance(other, float) :
return Matrix.Scalar(other), self
else :
raise TypeError()
def __str__(self) :
return "Vector[" + ", ".join(self.elements)+ "]"
def __repr__(self) :
return "Vector" + repr(self.elements)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment