[Scipy-svn] r3900 - trunk/scipy/splinalg/eigen/arpack/tests
scipy-svn at scipy.org
scipy-svn at scipy.org
Wed Feb 6 13:30:28 EST 2008
Author: hagberg
Date: 2008-02-06 12:30:22 -0600 (Wed, 06 Feb 2008)
New Revision: 3900
Modified:
trunk/scipy/splinalg/eigen/arpack/tests/test_arpack.py
Log:
Improved arpack tests
Modified: trunk/scipy/splinalg/eigen/arpack/tests/test_arpack.py
===================================================================
--- trunk/scipy/splinalg/eigen/arpack/tests/test_arpack.py 2008-02-06 18:28:08 UTC (rev 3899)
+++ trunk/scipy/splinalg/eigen/arpack/tests/test_arpack.py 2008-02-06 18:30:22 UTC (rev 3900)
@@ -1,345 +1,255 @@
#!/usr/bin/env python
+__usage__ = """
+To run tests locally:
+ python tests/test_arpack.py [-l<int>] [-v<int>]
+"""
+
from scipy.testing import *
-
+from numpy import array,real,imag,finfo,concatenate,\
+ column_stack,argsort,dot,round,conj,sort
from scipy.splinalg.eigen.arpack import eigen_symmetric,eigen
-from scipy.linalg import eig,eigh,norm
-from numpy import array,abs,take,concatenate,dot
+def assert_almost_equal_cc(actual,desired,decimal=7,err_msg='',verbose=True):
+ # almost equal or complex conjugates almost equal
+ try:
+ assert_almost_equal(actual,desired,decimal,err_msg,verbose)
+ except:
+ assert_almost_equal(actual,conj(desired),decimal,err_msg,verbose)
-class TestEigenNonsymmetric(TestCase):
- def get_a1(self,typ):
- mat=array([[-2., -8., 1., 2., -5.],
- [ 6., 6., 0., 2., 1.],
- [ 0., 4., -2., 11., 0.],
- [ 1., 6., 1., 0., -4.],
- [ 2., -6., 4., 9., -3]],typ)
+def assert_array_almost_equal_cc(actual,desired,decimal=7,
+ err_msg='',verbose=True):
+ # almost equal or complex conjugates almost equal
+ try:
+ assert_array_almost_equal(actual,desired,decimal,err_msg,verbose)
+ except:
+ assert_array_almost_equal(actual,conj(desired),decimal,err_msg,verbose)
- w=array([-2.21691+8.59661*1j,-2.21691-8.59661*1j,\
- 4.45961+3.80078*1j, 4.45961-3.80078*1j,\
- -5.48541+0j],typ.upper())
- return mat,w
- def large_magnitude(self,typ,k):
- a,aw = self.get_a1(typ)
- w,v = eigen(a,k,which='LM')
- for i in range(k):
- assert_array_almost_equal(dot(a,v[:,i]),w[i]*v[:,i],decimal=5)
- exact=abs(aw)
- num=abs(w)
- exact.sort()
- num.sort()
- assert_array_almost_equal(num[-k:],exact[-k:],decimal=5)
+# precision for tests
+_ndigits = {'f':5, 'd':12, 'F':5, 'D':12}
- def small_magnitude(self,typ,k):
- a,aw = self.get_a1(typ)
- w,v = eigen(a,k,which='SM')
- for i in range(k):
- assert_array_almost_equal(dot(a,v[:,i]),w[i]*v[:,i],decimal=5)
- exact=abs(aw)
- num=abs(w)
- exact.sort()
- num.sort()
- assert_array_almost_equal(num[:k],exact[:k],decimal=5)
+class TestArpack(TestCase):
+ def setUp(self):
+ self.symmetric=[]
+ self.nonsymmetric=[]
- def large_real(self,typ,k):
- a,aw = self.get_a1(typ)
- w,v = eigen(a,k,which='LR')
- for i in range(k):
- assert_array_almost_equal(dot(a,v[:,i]),w[i]*v[:,i],decimal=5)
- exact=aw.real
- num=w.real
- exact.sort()
- num.sort()
- assert_array_almost_equal(num[-k:],exact[-k:],decimal=5)
+ S1={}
+ S1['mat']=\
+ array([[ 2., 0., 0., -1., 0., -1.],
+ [ 0., 2., 0., -1., 0., -1.],
+ [ 0., 0., 2., -1., 0., -1.],
+ [-1., -1., -1., 4., 0., -1.],
+ [ 0., 0., 0., 0., 1., -1.],
+ [-1., -1., -1., -1., -1., 5.]])
- def small_real(self,typ,k):
- a,aw = self.get_a1(typ)
- w,v = eigen(a,k,which='SR')
- for i in range(k):
- assert_array_almost_equal(dot(a,v[:,i]),w[i]*v[:,i],decimal=5)
- exact=aw.real
- num=w.real
- exact.sort()
- num.sort()
- assert_array_almost_equal(num[:k],exact[:k],decimal=5)
+ S1['eval']=array([0,1,2,2,5,6])
+ self.symmetric.append(S1)
-
- def large_imag(self,typ,k):
- a,aw = self.get_a1(typ)
- w,v = eigen(a,k,which='LI')
- for i in range(k):
- assert_array_almost_equal(dot(a,v[:,i]),w[i]*v[:,i],decimal=5)
- print w
- print aw
- exact=aw.imag
- num=w.imag
- exact.sort()
- num.sort()
- assert_array_almost_equal(num[-k:],exact[-k:],decimal=5)
-
- def small_imag(self,typ,k):
- a,aw = self.get_a1(typ)
- w,v = eigen(a,k,which='SI')
- for i in range(k):
- assert_array_almost_equal(dot(a,v[:,i]),w[i]*v[:,i],decimal=5)
- exact=aw.imag
- num=w.imag
- exact.sort()
- num.sort()
- print num
- assert_array_almost_equal(num[:k],exact[:k],decimal=5)
-
-
- def test_type(self):
- k=2
- for typ in 'fd':
- self.large_magnitude(typ,k)
- self.small_magnitude(typ,k)
- self.large_real(typ,k)
- self.small_real(typ,k)
-# Maybe my understanding of small imaginary and large imaginary
-# isn't too keen. I don't understand why these return
-# different answers than in the complex case (the latter seems correct)
-# self.large_imag(typ,k)
-# self.small_imag(typ,k)
-
-
-
-class TestEigenComplexNonsymmetric(TestCase):
-
- def get_a1(self,typ):
- mat=array([[-2., -8., 1., 2., -5.],
+ N1={}
+ N1['mat']=\
+ array([[-2., -8., 1., 2., -5.],
[ 6., 6., 0., 2., 1.],
[ 0., 4., -2., 11., 0.],
[ 1., 6., 1., 0., -4.],
- [ 2., -6., 4., 9., -3]],typ)
+ [ 2., -6., 4., 9., -3]])
- w=array([-2.21691+8.59661*1j,-2.21691-8.59661*1j,\
- 4.45961+3.80078*1j, 4.45961-3.80078*1j,\
- -5.48541+0j],typ.upper())
- return mat,w
+ N1['eval']=\
+ array([ -5.4854094033782888+0.0j,
+ -2.2169058544873783+8.5966096591588261j,
+ -2.2169058544873783-8.5966096591588261j,
+ 4.4596105561765107+3.8007839204319454j,
+ 4.4596105561765107-3.8007839204319454j],'D')
- def large_magnitude(self,typ,k):
- a,aw = self.get_a1(typ)
- w,v = eigen(a,k,which='LM')
- for i in range(k):
- assert_array_almost_equal(dot(a,v[:,i]),w[i]*v[:,i],decimal=5)
- exact=abs(aw)
- num=abs(w)
- exact.sort()
- num.sort()
- assert_array_almost_equal(num,exact[-k:],decimal=5)
- def small_magnitude(self,typ,k):
- a,aw = self.get_a1(typ)
- w,v = eigen(a,k,which='SM')
- for i in range(k):
- assert_array_almost_equal(dot(a,v[:,i]),w[i]*v[:,i],decimal=5)
- exact=abs(aw)
- num=abs(w)
- exact.sort()
- num.sort()
- assert_array_almost_equal(num,exact[:k],decimal=5)
+ self.nonsymmetric.append(N1)
+
+class TestEigenSymmetric(TestArpack):
- def large_real(self,typ,k):
- a,aw = self.get_a1(typ)
- w,v = eigen(a,k,which='LR')
- for i in range(k):
- assert_array_almost_equal(dot(a,v[:,i]),w[i]*v[:,i],decimal=5)
- exact=aw.real
- num=w.real
- exact.sort()
- num.sort()
- assert_array_almost_equal(num,exact[-k:],decimal=5)
+ def get_exact_eval(self,d,typ,k,which):
+ eval=d['eval'].astype(typ)
+ ind=argsort(eval)
+ eval=eval[ind]
+ if which=='LM':
+ return eval[-k:]
+ if which=='SM':
+ return eval[:k]
+ if which=='BE':
+ # one ev from each end - if k is odd, extra ev on high end
+ l=k/2
+ h=k/2+k%2
+ low=range(len(eval))[:l]
+ high=range(len(eval))[-h:]
+ return eval[low+high]
- def small_real(self,typ,k):
- a,aw = self.get_a1(typ)
- w,v = eigen(a,k,which='SR')
+ def eval_evec(self,d,typ,k,which):
+ a=d['mat'].astype(typ)
+ exact_eval=self.get_exact_eval(d,typ,k,which)
+ eval,evec=eigen_symmetric(a,k,which=which)
+ # check eigenvalues
+ assert_array_almost_equal(eval,exact_eval,decimal=_ndigits[typ])
+ # check eigenvectors A*evec=eval*evec
for i in range(k):
- assert_array_almost_equal(dot(a,v[:,i]),w[i]*v[:,i],decimal=5)
- exact=aw.real
- num=w.real
- exact.sort()
- num.sort()
- assert_array_almost_equal(num,exact[:k],decimal=5)
+ assert_array_almost_equal(dot(a,evec[:,i]),
+ eval[i]*evec[:,i],
+ decimal=_ndigits[typ])
+ def test_symmetric(self):
+ k=2
+ for typ in 'fd':
+ for which in ['LM','SM','BE']:
+ self.eval_evec(self.symmetric[0],typ,k,which)
- def large_imag(self,typ,k):
- a,aw = self.get_a1(typ)
- w,v = eigen(a,k,which='LI')
- for i in range(k):
- assert_array_almost_equal(dot(a,v[:,i]),w[i]*v[:,i],decimal=5)
- exact=aw.imag
- num=w.imag
- exact.sort()
- num.sort()
- assert_array_almost_equal(num,exact[-k:],decimal=5)
+
+class TestEigenComplexSymmetric(TestArpack):
- def small_imag(self,typ,k):
- a,aw = self.get_a1(typ)
- w,v = eigen(a,k,which='SI')
+ def sort_choose(self,eval,typ,k,which):
+ # sort and choose the eigenvalues and eigenvectors
+ # both for the exact answer and that returned from ARPACK
+ reval=round(eval,decimals=_ndigits[typ])
+ ind=argsort(reval)
+ if which=='LM' or which=='LR':
+ return ind[-k:]
+ if which=='SM' or which=='SR':
+ return ind[:k]
+
+ def eval_evec(self,d,typ,k,which):
+ a=d['mat'].astype(typ)
+ # get exact eigenvalues
+ exact_eval=d['eval'].astype(typ)
+ ind=self.sort_choose(exact_eval,typ,k,which)
+ exact_eval=exact_eval[ind]
+ # compute eigenvalues
+ eval,evec=eigen(a,k,which=which)
+ ind=self.sort_choose(eval,typ,k,which)
+ eval=eval[ind]
+ evec=evec[:,ind]
+
+ # check eigenvalues
+ assert_array_almost_equal(eval,exact_eval,decimal=_ndigits[typ])
+ # check eigenvectors A*evec=eval*evec
for i in range(k):
- assert_array_almost_equal(dot(a,v[:,i]),w[i]*v[:,i],decimal=5)
- exact=aw.imag
- num=w.imag
- exact.sort()
- num.sort()
- assert_array_almost_equal(num,exact[:k],decimal=5)
+ assert_array_almost_equal(dot(a,evec[:,i]),
+ eval[i]*evec[:,i],
+ decimal=_ndigits[typ])
-# def test_type(self):
+# def test_complex_symmetric(self):
# k=2
# for typ in 'FD':
-# self.large_magnitude(typ,k)
-# self.small_magnitude(typ,k)
-# self.large_real(typ,k)
-# self.small_real(typ,k)
-# self.large_imag(typ,k)
-# self.small_imag(typ,k)
+# for which in ['LM','SM','LR','SR']:
+# self.eval_evec(self.symmetric[0],typ,k,which)
+
+class TestEigenNonSymmetric(TestArpack):
-class TestEigenSymmetric(TestCase):
+ def sort_choose(self,eval,typ,k,which):
+ reval=round(eval,decimals=_ndigits[typ])
+ if which in ['LR','SR']:
+ ind=argsort(reval.real)
+ elif which in ['LI','SI']:
+ # for LI,SI ARPACK returns largest,smallest abs(imaginary) why?
+ ind=argsort(abs(reval.imag))
+ else:
+ ind=argsort(abs(reval))
- def get_a1(self,typ):
- mat_a1=array([[ 2., 0., 0., -1., 0., -1.],
- [ 0., 2., 0., -1., 0., -1.],
- [ 0., 0., 2., -1., 0., -1.],
- [-1., -1., -1., 4., 0., -1.],
- [ 0., 0., 0., 0., 1., -1.],
- [-1., -1., -1., -1., -1., 5.]],
- typ)
- w = [0,1,2,2,5,6] # eigenvalues of a1
- return mat_a1,w
+ if which in ['LR','LM','LI']:
+ return ind[-k:]
+ if which in ['SR','SM','SI']:
+ return ind[:k]
- def large_eigenvalues(self,typ,k):
- a,aw = self.get_a1(typ)
- w,v = eigen_symmetric(a,k,which='LM',tol=1e-7)
- assert_array_almost_equal(w,aw[-k:])
- def small_eigenvalues(self,typ,k):
- a,aw = self.get_a1(typ)
- w,v = eigen_symmetric(a,k,which='SM')
- assert_array_almost_equal(w,aw[:k])
-
- def end_eigenvalues(self,typ,k):
- a,aw = self.get_a1(typ)
- w,v = eigen_symmetric(a,k,which='BE')
- exact=[aw[0],aw[-1]]
- assert_array_almost_equal(w,exact)
-
- def large_eigenvectors(self,typ,k):
- a,aw = self.get_a1(typ)
- w,v = eigen_symmetric(a,k,which='LM')
- ew,ev = eigh(a)
- ind=ew.argsort()
- assert_array_almost_equal(w,take(ew,ind[-k:]))
+ def eval_evec(self,d,typ,k,which):
+ a=d['mat'].astype(typ)
+ # get exact eigenvalues
+ exact_eval=d['eval'].astype(typ.upper())
+ ind=self.sort_choose(exact_eval,typ,k,which)
+ exact_eval=exact_eval[ind]
+ # compute eigenvalues
+ eval,evec=eigen(a,k,which=which)
+ ind=self.sort_choose(eval,typ,k,which)
+ eval=eval[ind]
+ evec=evec[:,ind]
+ # check eigenvalues
+ # check eigenvectors A*evec=eval*evec
for i in range(k):
- assert_array_almost_equal(dot(a,v[:,i]),w[i]*v[:,i])
+ assert_almost_equal_cc(eval[i],exact_eval[i],decimal=_ndigits[typ])
+ assert_array_almost_equal_cc(dot(a,evec[:,i]),
+ eval[i]*evec[:,i],
+ decimal=_ndigits[typ])
- def small_eigenvectors(self,typ,k):
- a,aw = self.get_a1(typ)
- w,v = eigen_symmetric(a,k,which='SM',tol=1e-7)
- ew,ev = eigh(a)
- ind=ew.argsort()
- assert_array_almost_equal(w,take(ew,ind[:k]))
- for i in range(k):
- assert_array_almost_equal(dot(a,v[:,i]),w[i]*v[:,i])
- def end_eigenvectors(self,typ,k):
- a,aw = self.get_a1(typ)
- w,v = eigen_symmetric(a,k,which='BE')
- ew,ev = eigh(a)
- ind=ew.argsort()
- exact=concatenate(([ind[:k/2],ind[-k/2:]]))
- assert_array_almost_equal(w,take(ew,exact))
- for i in range(k):
- assert_array_almost_equal(dot(a,v[:,i]),w[i]*v[:,i])
-
- def test_eigenvectors(self):
+ def test_nonsymmetric(self):
k=2
for typ in 'fd':
- self.large_eigenvectors(typ,k)
- self.small_eigenvectors(typ,k)
- self.end_eigenvectors(typ,k)
+ for which in ['LI','LR','LM','SM','SR','SI']:
+ for m in self.nonsymmetric:
+ self.eval_evec(m,typ,k,which)
- def test_type(self):
- k=2
- for typ in 'fd':
- self.large_eigenvalues(typ,k)
- self.small_eigenvalues(typ,k)
- self.end_eigenvalues(typ,k)
-class TestEigenComplexSymmetric(TestCase):
- def get_a1(self,typ):
- mat_a1=array([[ 2., 0., 0., -1., 0., -1.],
- [ 0., 2., 0., -1., 0., -1.],
- [ 0., 0., 2., -1., 0., -1.],
- [-1., -1., -1., 4., 0., -1.],
- [ 0., 0., 0., 0., 1., -1.],
- [-1., -1., -1., -1., -1., 5.]],
- typ)
- w = array([0+0j,1+0j,2+0j,2+0j,5+0j,6+0j]) # eigenvalues of a1
- return mat_a1,w
- def large_magnitude(self,typ,k):
- a,aw = self.get_a1(typ)
- w,v = eigen(a,k,which='LM')
- for i in range(k):
- assert_array_almost_equal(dot(a,v[:,i]),w[i]*v[:,i],decimal=5)
- aw.real.sort()
- w.real.sort()
- assert_array_almost_equal(w,aw[-k:])
+class TestEigenComplexNonSymmetric(TestArpack):
+ def sort_choose(self,eval,typ,k,which):
+ eps=finfo(typ).eps
+ reval=round(eval,decimals=_ndigits[typ])
+ if which in ['LR','SR']:
+ ind=argsort(reval)
+ elif which in ['LI','SI']:
+ ind=argsort(reval.imag)
+ else:
+ ind=argsort(abs(reval))
- def small_magnitude(self,typ,k):
- a,aw = self.get_a1(typ)
- w,v = eigen(a,k,which='SM')
- for i in range(k):
- assert_array_almost_equal(dot(a,v[:,i]),w[i]*v[:,i])
- aw.real.sort()
- w.real.sort()
- assert_array_almost_equal(w,aw[:k])
+ if which in ['LR','LI','LM']:
+ return ind[-k:]
+ if which in ['SR','SI','SM']:
+ return ind[:k]
- def large_real(self,typ,k):
- a,aw = self.get_a1(typ)
- w,v = eigen(a,k,which='LR')
- for i in range(k):
- assert_array_almost_equal(dot(a,v[:,i]),w[i]*v[:,i],decimal=5)
- aw.real.sort()
- w.real.sort()
- assert_array_almost_equal(w,aw[-k:],decimal=5)
+ def eval_evec(self,d,typ,k,which):
+ a=d['mat'].astype(typ)
+ # get exact eigenvalues
+ exact_eval=d['eval'].astype(typ.upper())
+ ind=self.sort_choose(exact_eval,typ,k,which)
+ exact_eval=exact_eval[ind]
+ print "exact"
+ print exact_eval
- def small_real(self,typ,k):
- a,aw = self.get_a1(typ)
- w,v = eigen(a,k,which='SR')
+ # compute eigenvalues
+ eval,evec=eigen(a,k,which=which)
+ ind=self.sort_choose(eval,typ,k,which)
+ eval=eval[ind]
+ evec=evec[:,ind]
+ print eval
+ # check eigenvalues
+ # check eigenvectors A*evec=eval*evec
for i in range(k):
- assert_array_almost_equal(dot(a,v[:,i]),w[i]*v[:,i])
- aw.real.sort()
- w.real.sort()
- assert_array_almost_equal(w,aw[:k])
+ assert_almost_equal_cc(eval[i],exact_eval[i],decimal=_ndigits[typ])
+ assert_array_almost_equal_cc(dot(a,evec[:,i]),
+ eval[i]*evec[:,i],
+ decimal=_ndigits[typ])
-# def test_complex_symmetric(self):
+
+# def test_complex_nonsymmetric(self):
# k=2
# for typ in 'FD':
-# self.large_magnitude(typ,k)
-# self.small_magnitude(typ,k)
-# self.large_real(typ,k)
-# self.small_real(typ,k)
+# for which in ['LI','LR','LM','SI','SR','SM']:
+# for m in self.nonsymmetric:
+# self.eval_evec(m,typ,k,which)
+
if __name__ == "__main__":
nose.run(argv=['', __file__])
More information about the Scipy-svn
mailing list