Permalink
Browse files

Making tests for dokmatrix work

  • Loading branch information...
1 parent ff430bf commit 39adceffe68ce1bcd3c4386edfae7060340d795e @sherjilozair committed Aug 14, 2011
View
@@ -1,3 +1,4 @@
from densematrix import DenseMatrix
from lilmatrix import LILMatrix
from dokmatrix import DOKMatrix
+from matrix_ import Matrix_
View
@@ -60,6 +60,8 @@ def __init__(self, *args, **kwargs):
rows = cols = 0
mat = []
+ for i in mat:
+ mat[i] = S(mat[i])
self.rows = rows
self.cols = cols
self.mat = mat
@@ -71,7 +73,7 @@ def __getitem__(self, key):
def __setitem__(self, key, value):
"Key cannot be a slice"
if value != 0:
- self.mat[key] == value
+ self.mat[key] = value
else:
del self.mat[key]
@@ -120,27 +122,28 @@ def __mul__(self, other):
Multiplication a matrix with another matrix or a scalar.
"""
if isinstance(other, DOKMatrix):
- rows1 = _lil_row_major(A)
- rows2 = _lil_row_major(B)
+ from dokmatrix_tools import _lil_row_major
+ rows1 = _lil_row_major(self)
+ rows2 = _lil_row_major(other)
Cdict = {}
- for k in xrange(A.rows):
+ for k in xrange(self.rows):
for _, j in rows1[k]:
for _, n in rows2[j]:
- temp = A[k, j] * B[j, n]
+ temp = self[k, j] * other[j, n]
if (k, n) in Cdict:
Cdict[k, n] += temp
else:
Cdict[k, n] = temp
- C = DOKMatrix(A.rows, B.cols, Cdict)
+ C = DOKMatrix(self.rows, other.cols, Cdict)
if C.shape == (1, 1):
return C[0, 0]
return C
else:
"other is a scalar"
- C = DOKMatrix(matrix.rows, matrix.cols, {})
- if scalar != 0:
- for i in matrix.mat:
- C.mat[i] = scalar * matrix.mat[i]
+ C = DOKMatrix(self.rows, self.cols, {})
+ if B != 0:
+ for i in self.mat:
+ C.mat[i] = B * self.mat[i]
return C
def __neg__(self):
@@ -153,7 +156,7 @@ def __rmul__(self, other):
return self.__mul__(other)
def __eq__(self, other):
- if not isinstance(other, (Matrix, DOKMatrix)):
+ if not isinstance(other, DOKMatrix):
return False
if self.rows != other.rows or self.cols != other.cols:
return False
@@ -213,11 +216,12 @@ def to_dokmatrix(self):
return self
def to_lilmatrix(self):
- from sympy.matrices.lilmatrix import LILMatrix
+ from sympy.linalg.lilmatrix import LILMatrix
return LILMatrix._from_dict(self.rows, self.cols, self.mat)
def to_densematrix(self):
- return Matrix(self.rows, self.cols, lambda i, j: self[i, j])
+ from sympy.linalg.densematrix import DenseMatrix
+ return DenseMatrix(self.rows, self.cols, lambda i, j: self[i, j])
def copy(self):
return DOKMatrix(self.rows, self.cols, self.mat)
@@ -1,3 +1,5 @@
+from __future__ import division
+from sympy.linalg.dokmatrix import DOKMatrix
def rowdecomp(self, num):
nmax = len(self)
if not (0 <= num < nmax) and not (0 <= -num < nmax):
@@ -130,7 +132,7 @@ def liu(self):
return parent
def liupc(self):
- R = self._lower_row_nonzero_structure()
+ R = _lower_row_nonzero_structure(self)
parent = [None] * self.rows
virtual = [None] * self.rows
for j in xrange(self.rows):
@@ -146,7 +148,7 @@ def liupc(self):
def row_structure_symbolic_cholesky(self):
from copy import deepcopy
- R, parent = self.liupc()
+ R, parent = liupc(self)
Lrow = deepcopy(R)
for k in xrange(self.rows):
for _, j in R[k]:
@@ -169,7 +171,7 @@ def elementary_symbolic_cholesky(self):
def fast_symbolic_cholesky(self):
""" implement algorithm 1.3 from 10.1.1.163.7506 """
- C = self._lower_columnar_nonzero_structure()
+ C = _lower_columnar_nonzero_structure(self)
p=[]
for k, col in enumerate(C):
if len(C[k]) > 1:
@@ -185,11 +187,7 @@ def fast_symbolic_cholesky(self):
def _cholesky(self):
- cached_result = self._get_cache('CH')
- if cached_result and self._cached:
- return cached_result
-
- CHstruc = self.fast_symbolic_cholesky()
+ CHstruc = fast_symbolic_cholesky(self)
C = DOKMatrix(self.rows, self.cols, {})
for col in CHstruc:
for i, j in col:
@@ -200,15 +198,11 @@ def _cholesky(self):
C[j, j] = (self[j, j] - sum(C[j, k] ** 2
for k in xrange(j))) ** (0.5)
- self._set_cache('CH', C)
+
return C
def _cholesky_sparse(self):
- cached_result = self._get_cache('CH')
- if cached_result and self._cached:
- return cached_result
-
Crowstruc = self.row_structure_symbolic_cholesky()
C = DOKMatrix.zeros(self.rows)
for row in Crowstruc:
@@ -239,18 +233,14 @@ def _cholesky_sparse(self):
C[j, j] -= summ
C[j, j] = C[j, j] ** 0.5
- self._set_cache('CH', C)
+
return C
def _LDL(self):
- cached_result = self._get_cache('LDL')
- if cached_result and self._cached:
- return cached_result
-
- CHstruc = self.fast_symbolic_cholesky()
- L = DOKMatrix.eye(self.rows)
+ CHstruc = fast_symbolic_cholesky(self)
+ L = eye(self.rows)
D = DOKMatrix(self.rows, self.cols, {})
for col in CHstruc:
@@ -261,19 +251,13 @@ def _LDL(self):
elif i == j:
D[i, i] = self[i, i] - sum(L[i, k]**2 * D[k, k]
for k in range(i))
-
- self._set_cache('LDL', (L, D))
return L, D
def _LDL_sparse(self):
- cached_result = self._get_cache('LDL')
- if cached_result and self._cached:
- return cached_result
-
- Lrowstruc = self.row_structure_symbolic_cholesky()
- L = DOKMatrix.eye(self.rows)
+ Lrowstruc = row_structure_symbolic_cholesky(self)
+ L = eye(self.rows)
D = DOKMatrix(self.rows, self.cols, {})
for row in Lrowstruc:
@@ -303,12 +287,11 @@ def _LDL_sparse(self):
break
D[i, i] -= summ
- self._set_cache('LDL', (L, D))
return L, D
def _lower_triangular_solve(self, rhs):
- rows = self._lil_row_major()
+ rows = _lil_row_major(self)
X = DOKMatrix(rhs.rows, 1, rhs.mat)
for i in xrange(self.rows):
for key in rows[i]:
@@ -321,7 +304,7 @@ def _lower_triangular_solve(self, rhs):
def _upper_triangular_solve(self, rhs):
""" Helper function of function upper_triangular_solve, without the error checks, to be used privately. """
- rows = self._lil_row_major()
+ rows = _lil_row_major(self)
X = DOKMatrix(rhs.rows, 1, rhs.mat)
for i in reversed(xrange(self.rows)):
for key in reversed(rows[i]):
@@ -336,32 +319,22 @@ def _diagonal_solve(self, rhs):
return DOKMatrix(self.rows, 1, lambda i, j: rhs[i, 0] / self[i, i])
def _cholesky_solve(self, rhs):
- L = self._cholesky()
- Y = L._lower_triangular_solve(rhs)
- X = L.T._upper_triangular_solve(Y)
- return X
-
-def LDL_solve(self, rhs):
- if not self.is_square():
- raise Exception("Make a square matrix exception")
- if not self.rows == rhs.rows:
- raise Exception
if not self.is_symmetric():
rhs = self.T * rhs
self = self.T * self
- L, D = self._LDL_sparse()
- z = L._lower_triangular_solve(rhs)
- y = D._diagonal_solve(z)
- x = L.T._upper_triangular_solve(y)
- if x.has(nan) or x.has(oo): # import this
- raise Exception
- return x
+ L = _cholesky(self)
+ Y = _lower_triangular_solve(L, rhs)
+ X = _upper_triangular_solve(L.T, Y)
+ return X
def _LDLsolve(self, rhs):
- L, D = self._LDL_sparse()
- z = L._lower_triangular_solve(rhs)
- y = D._diagonal_solve(z)
- x = L.T._upper_triangular_solve(y)
+ if not self.is_symmetric():
+ rhs = self.T * rhs
+ self = self.T * self
+ L, D = _LDL_sparse(self)
+ z = _lower_triangular_solve(L, rhs)
+ y = _diagonal_solve(D, z)
+ x = _upper_triangular_solve(L.T, y)
return x
def _lower_columnar_nonzero_structure(self):
@@ -416,12 +389,25 @@ def copyin_list(self, key, value):
def inverse_solver(self, solver=None):
from matrixutils import vecs2matrix
if not solver or solver == 'CH':
- solver = self._cholesky_solve
+ solver = _cholesky_solve
elif solver == 'LDL':
- solver = self.LDL_solve
+ solver = LDL_solve
else:
raise Exception('solver method not recognized')
- I = Matrix.eye(self.rows)
- return vecs2matrix([solver(I[:, i])
+ I = eye(self.rows)
+ return vecs2matrix([solver(self, I[:, i])
for i in xrange(self.cols)])
+def zeros(*args):
+ from dokmatrix import DOKMatrix
+ if len(args) == 1:
+ m = n = args[0]
+ elif len(args) == 2:
+ m, n = args
+ return DOKMatrix(m, n, lambda i, j: 0)
@vks

vks Aug 16, 2011

This is horribly inefficient. It should be O(1) and not O(n**2). (At least for sparse matrices.)

@sherjilozair

sherjilozair Aug 16, 2011

Owner

Horrible code ! I'll edit it.

+
+def eye(n):
+ matrix = zeros(n)
+ for i in xrange(n):
+ matrix[i, i] = 1
+ return matrix
@vks

vks Aug 16, 2011

This code is not specific to DOKMatrix. Why put it here? Do you intend to implement a more specialized variant?

Also, this is not mentioned in the commit message. Ideally it should be in another commit. At least it should be mentioned that you implemented a new feature.

@@ -81,6 +81,32 @@ def _lilrepr_from_lil(rows, cols, lil):
mat[i].append((j, val))
return mat
+def _dokrepr_from_lil(rows, cols, lil):
+ mat = {}
+ for i in xrange(rows):
+ for j in xrange(cols):
+ val = lil[i][j]
+ if val != 0:
+ mat[i, j] = val
+ return mat
+
+def _dokrepr_from_dict(rows, cols, di):
+ mat = {}
+ for i in xrange(rows):
+ for j in xrange(cols):
+ if (i, j) in di:
+ mat[i, j] = di[i, j]
+ return mat
+
+def _dokrepr_from_callable(rows, cols, func):
+ mat = {}
+ for i in xrange(rows):
+ for j in xrange(cols):
+ val = func(i, j)
+ if val != 0:
+ mat[i, j] = val
+ return mat
+
@@ -1,5 +1,5 @@
-from sympy.linalg.dokmatrix import DOKMatrix
from __future__ import division
+from sympy.linalg.dokmatrix import DOKMatrix
def test_getitem():
A = DOKMatrix(((1,2,0),(4,0,6),(7,8,0)))
@@ -21,7 +21,7 @@ def test_copyinmatrix():
def test_transpose():
A = DOKMatrix(((1,0,3),(0,5,6),(7,0,9)))
- assert A.T == DOKMatrix(((1,0,3),(2,0,4),(3,0,0)))
+ assert A.T == DOKMatrix(((1,0,7),(0,5,0),(3,6,9)))
def test_add():
A = DOKMatrix(((1,0,3),(0,5,6),(7,0,9)))
@@ -31,7 +31,7 @@ def test_add():
def test_mul():
A = DOKMatrix(((1,2,3),(4,0,1),(9,8,0)))
B = DOKMatrix(((0,0,1),(2,1,0),(0,2,7)))
- C = DOKMatrix(((4,8,27),(0,2,11),(16,8,9)))
+ C = DOKMatrix(((4,8,22),(0,2,11),(16,8,9)))
assert A * B == C
def test_submatrix():
@@ -45,13 +45,6 @@ def test_submatrix():
def test_solve():
A = DOKMatrix(((1,2,3),(4,0,1),(9,8,0)))
- rhs = DOKMatrix((3,4,5)).T
- X = DOKMatrix((41/53),(-13/53),(48/53)).T
+ rhs = DOKMatrix([[3,4,5]]).T
for method in ["LDL", "CH"]:
- assert A.solve(rhs, method=method) == X, method
-
-
-
-
-
-
+ assert A * A.solve(rhs, method=method) == rhs, method

0 comments on commit 39adcef

Please sign in to comment.