[Scipy-svn] r3611 - in trunk/scipy/io: . tests
scipy-svn at scipy.org
scipy-svn at scipy.org
Mon Dec 3 16:15:10 EST 2007
Author: wnbell
Date: 2007-12-03 15:14:51 -0600 (Mon, 03 Dec 2007)
New Revision: 3611
Modified:
trunk/scipy/io/mmio.py
trunk/scipy/io/tests/test_mmio.py
Log:
made sparse matrix market reading faster
added more comprehensive unittests for sparse format
Modified: trunk/scipy/io/mmio.py
===================================================================
--- trunk/scipy/io/mmio.py 2007-12-03 03:31:42 UTC (rev 3610)
+++ trunk/scipy/io/mmio.py 2007-12-03 21:14:51 UTC (rev 3611)
@@ -1,5 +1,3 @@
-## Automatically adapted for scipy Oct 19, 2005 by convertcode.py
-
"""
Matrix Market I/O in Python.
"""
@@ -10,10 +8,10 @@
# References:
# http://math.nist.gov/MatrixMarket/
#
-# TODO: support for sparse matrices, need spmatrix.tocoo().
import os
-from numpy import asarray, real, imag, conj, zeros, ndarray
+from numpy import asarray, real, imag, conj, zeros, ndarray, \
+ empty, concatenate, ones
from itertools import izip
__all__ = ['mminfo','mmread','mmwrite']
@@ -184,38 +182,46 @@
assert k==entries,`k,entries`
elif rep=='coordinate':
- k = 0
- data,row,col = [],[],[]
- row_append = row.append
- col_append = col.append
- data_append = data.append
- line = '%'
- while line:
- if not line.startswith('%'):
- l = line.split()
- i = int(l[0])-1
- j = int(l[1])-1
- if is_pattern:
- aij = 1.0 #use 1.0 for pattern matrices
- elif is_complex:
- aij = complex(*map(float,l[2:]))
- else:
- aij = float(l[2])
- row_append(i)
- col_append(j)
- data_append(aij)
- if has_symmetry and i!=j:
- if is_skew:
- aij = -aij
- elif is_herm:
- aij = conj(aij)
- row_append(j)
- col_append(i)
- data_append(aij)
- k += 1
- line = source.readline()
- assert k==entries,`k,entries`
- a = coo_matrix((data, (row, col)), dims=(rows, cols), dtype=dtype)
+ from numpy import fromfile
+ flat_data = fromfile(source,sep=' ')
+ if is_pattern:
+ flat_data = flat_data.reshape(-1,2)
+ I = flat_data[:,0].astype('i')
+ J = flat_data[:,1].astype('i')
+ V = ones(len(I))
+ elif is_complex:
+ flat_data = flat_data.reshape(-1,4)
+ I = flat_data[:,0].astype('i')
+ J = flat_data[:,1].astype('i')
+ V = empty(len(I),dtype='complex')
+ V.real = flat_data[:,2]
+ V.imag = flat_data[:,3]
+ else:
+ flat_data = flat_data.reshape(-1,3)
+ I = flat_data[:,0].astype('i')
+ J = flat_data[:,1].astype('i')
+ V = flat_data[:,2].copy()
+
+ I -= 1 #adjust indices (base 1 -> base 0)
+ J -= 1
+
+ if has_symmetry:
+ mask = (I != J) #off diagonal mask
+ od_I = I[mask]
+ od_J = J[mask]
+ od_V = V[mask]
+
+ I = concatenate((I,od_J))
+ J = concatenate((J,od_I))
+
+ if is_skew:
+ od_V *= -1
+ elif is_herm:
+ od_V = od_V.conjugate()
+
+ V = concatenate((V,od_V))
+
+ a = coo_matrix((V, (I, J)), dims=(rows, cols), dtype=dtype)
else:
raise NotImplementedError,`rep`
Modified: trunk/scipy/io/tests/test_mmio.py
===================================================================
--- trunk/scipy/io/tests/test_mmio.py 2007-12-03 03:31:42 UTC (rev 3610)
+++ trunk/scipy/io/tests/test_mmio.py 2007-12-03 21:14:51 UTC (rev 3611)
@@ -102,7 +102,7 @@
b = mmread(fn)
assert_array_almost_equal(a,b)
-_exmpl_mtx = '''\
+_general_example = '''\
%%MatrixMarket matrix coordinate real general
%=================================================================================
%
@@ -136,12 +136,60 @@
5 5 1.200e+01
'''
+_hermitian_example = '''\
+%%MatrixMarket matrix coordinate complex hermitian
+ 5 5 7
+ 1 1 1.0 0
+ 2 2 10.5 0
+ 4 2 250.5 22.22
+ 3 3 1.5e-2 0
+ 4 4 -2.8e2 0
+ 5 5 12. 0
+ 5 4 0 33.32
+'''
+
+_skew_example = '''\
+%%MatrixMarket matrix coordinate real skew-symmetric
+ 5 5 7
+ 1 1 1.0
+ 2 2 10.5
+ 4 2 250.5
+ 3 3 1.5e-2
+ 4 4 -2.8e2
+ 5 5 12.
+ 5 4 0
+'''
+
+_symmetric_example = '''\
+%%MatrixMarket matrix coordinate real symmetric
+ 5 5 7
+ 1 1 1.0
+ 2 2 10.5
+ 4 2 250.5
+ 3 3 1.5e-2
+ 4 4 -2.8e2
+ 5 5 12.
+ 5 4 8
+'''
+
+_symmetric_pattern_example = '''\
+%%MatrixMarket matrix coordinate pattern symmetric
+ 5 5 7
+ 1 1
+ 2 2
+ 4 2
+ 3 3
+ 4 4
+ 5 5
+ 5 4
+'''
+
class TestMMIOCoordinate(NumpyTestCase):
-
- def check_simple_todense(self):
+ def check_read_geneal(self):
+ """read a general matrix"""
fn = mktemp()
f = open(fn,'w')
- f.write(_exmpl_mtx)
+ f.write(_general_example)
f.close()
assert_equal(mminfo(fn),(5,5,8,'coordinate','real','general'))
a = [[1, 0, 0, 6, 0],
@@ -152,6 +200,67 @@
b = mmread(fn).todense()
assert_array_almost_equal(a,b)
+ def check_read_hermitian(self):
+ """read a hermitian matrix"""
+ fn = mktemp()
+ f = open(fn,'w')
+ f.write(_hermitian_example)
+ f.close()
+ assert_equal(mminfo(fn),(5,5,7,'coordinate','complex','hermitian'))
+ a = [[1, 0, 0, 0, 0],
+ [0, 10.5, 0, 250.5 - 22.22j, 0],
+ [0, 0, .015, 0, 0],
+ [0, 250.5 + 22.22j, 0, -280, -33.32j],
+ [0, 0, 0, 33.32j, 12]]
+ b = mmread(fn).todense()
+ assert_array_almost_equal(a,b)
+
+ def check_read_skew(self):
+ """read a skew-symmetric matrix"""
+ fn = mktemp()
+ f = open(fn,'w')
+ f.write(_skew_example)
+ f.close()
+ assert_equal(mminfo(fn),(5,5,7,'coordinate','real','skew-symmetric'))
+ a = [[1, 0, 0, 0, 0],
+ [0, 10.5, 0, -250.5, 0],
+ [0, 0, .015, 0, 0],
+ [0, 250.5, 0, -280, 0],
+ [0, 0, 0, 0, 12]]
+ b = mmread(fn).todense()
+ assert_array_almost_equal(a,b)
+
+ def check_read_symmetric(self):
+ """read a symmetric matrix"""
+ fn = mktemp()
+ f = open(fn,'w')
+ f.write(_symmetric_example)
+ f.close()
+ assert_equal(mminfo(fn),(5,5,7,'coordinate','real','symmetric'))
+ a = [[1, 0, 0, 0, 0],
+ [0, 10.5, 0, 250.5, 0],
+ [0, 0, .015, 0, 0],
+ [0, 250.5, 0, -280, 8],
+ [0, 0, 0, 8, 12]]
+ b = mmread(fn).todense()
+ assert_array_almost_equal(a,b)
+
+ def check_read_symmetric(self):
+ """read a symmetric pattern matrix"""
+ fn = mktemp()
+ f = open(fn,'w')
+ f.write(_symmetric_pattern_example)
+ f.close()
+ assert_equal(mminfo(fn),(5,5,7,'coordinate','pattern','symmetric'))
+ a = [[1, 0, 0, 0, 0],
+ [0, 1, 0, 1, 0],
+ [0, 0, 1, 0, 0],
+ [0, 1, 0, 1, 1],
+ [0, 0, 0, 1, 1]]
+ b = mmread(fn).todense()
+ assert_array_almost_equal(a,b)
+
+
def check_real_write_read(self):
I = array([0, 0, 1, 2, 3, 3, 3, 4])
J = array([0, 3, 1, 2, 1, 3, 4, 4])
More information about the Scipy-svn
mailing list