#!/usr/bin/env python
# <<BEGIN-copyright>>
# <<END-copyright>>
from __future__ import print_function
from numpy import *
from numpy import linalg as LA
from fudge.core.math.linearAlgebra import *
import unittest
# --------- base class for all matrix tests -----------
[docs]class MatrixTests( unittest.TestCase ):
[docs] def assertMatrixEqual( self, a, b ):
"""The numpy "==" operator is a universal function, so it operates on each element in the matrices.
The result is another matrix whose elements are all a[i,j] == b[i,j].
You need to call the .all() function to check that all elements evaluated to "True"
Or, could just use array_equal()..."""
#return self.assertTrue( ( a == b ).all() )
return self.assertTrue( array_equal( a, b ) )
[docs] def assertMatrixAlmostEqual( self, a, b, rtol=1e-05, atol=1e-08 ):
return self.assertTrue( allclose( a, b, rtol=rtol, atol=atol ) )
[docs] def assertMatrixTrue( self, a ):
"""You need to call the .all() function to check that all elements evaluated to "True" """
return self.assertTrue( a.all() )
# --------- basic tests -----------
[docs]class BasicLinearAlgebraTests( MatrixTests ):
"""Very basic matrix tests"""
[docs] def test_transpose( self ):
"""Transpose"""
x = mat([[1, 0, 0],[1, 0, 0],[0, 0, 0]])
self.assertMatrixEqual( x.T, mat([[1, 1, 0],[0, 0, 0],[0, 0, 0]]) )
[docs] def test_inner_products(self):
"""Inner product test"""
x = mat( [1,0,0])
# * test
self.assertEqual( x*x.T, [[1]] )
# dot test
self.assertEqual( dot(x,x.T), [[1]] )
[docs] def test_outer_product(self):
"""Outer product test"""
x = mat( [1,0,0])
# outer product test
if False: print( outer(x,x.T) == mat([[1, 0, 0],[0, 0, 0],[0, 0, 0]]) )
self.assertMatrixEqual( outer(x,x.T), mat([[1, 0, 0],[0, 0, 0],[0, 0, 0]]) )
[docs] def test_matrix_addition( self ):
"""Addition test"""
x = mat( [[1,0,0],[0,0,0],[0,0,0]])
y = mat( [[0,0,0],[0,1,0],[0,0,0]])
z = zeros_like( x )
z += x
z += y
self.assertMatrixEqual( x+y, mat([[1,0,0],[0,1,0],[0,0,0]]) )
self.assertMatrixEqual( z, mat([[1,0,0],[0,1,0],[0,0,0]]) )
# --------- matrix composition tests -----------
[docs]class MatrixCompositionTests( MatrixTests ):
"""Very basic matrix tests"""
[docs] def setUp( self ):
self.a = mat([[ 1. , 0.1, 0.1], [ 0.1, 1. , 0.1], [ 0.1, 0.1, 1. ]])
self.b = zeros( ( 3,6 ) )
self.c = None
self.d = mat([[1, 2], [3, 4]])
self.e = mat([[5, 6], [7, 8]])
[docs] def test_stackH( self ):
"""Testing horizontal stacking"""
self.assertMatrixEqual( stackHorizontal( [ self.a, self.b, self.c ] ),
mat([[ 1.0, 0.1, 0.1, 0., 0., 0., 0., 0., 0.],
[ 0.1, 1.0, 0.1, 0., 0., 0., 0., 0., 0.],
[ 0.1, 0.1, 1.0, 0., 0., 0., 0., 0., 0.]]) )
[docs] def test_stackV( self ):
"""Testing vertical stacking"""
self.assertMatrixEqual( stackVertical( [ self.d, self.e, self.c, self.d, self.e ] ), mat( [[ 1., 2.], [ 3., 4.], [ 5., 6.], [ 7., 8.], [ 1., 2.], [ 3., 4.], [ 5., 6.], [ 7., 8.]] ) )
[docs] def test_stackD( self ):
"""Testing diagonal stacking"""
self.assertMatrixEqual( stackDiagonal( [ self.d, self.e, self.a ] ), mat( [[ 1., 2., 0., 0., 0., 0., 0. ], [ 3., 4., 0., 0., 0., 0., 0. ], [ 0., 0., 5., 6., 0., 0., 0. ], [ 0., 0., 7., 8., 0., 0., 0. ], [ 0., 0., 0., 0., 1., 0.1, 0.1], [ 0., 0., 0., 0., 0.1, 1., 0.1], [ 0., 0., 0., 0., 0.1, 0.1, 1. ]]) )
# --------- CGLSQR tests -----------
[docs]class CGLSQRTests( MatrixTests ):
[docs] def setUp( self ):
self.answer = numpy.matrix([[ 1.34883721,-0.69767442, 0.34883721, 0.1, 42.0]])
self.kernel = numpy.matrix( [ [ 1.0, 2.0, 3.0, 0.0, 0.0 ], [ 2.0, 3.0, 4.0, 0.0, 0.0 ], [ 4.5, 5.4, 2.0, 0.0, 0.0 ] ] ) # Note: lower two subspaces map to 0
self.data = numpy.mat([ [ 1.1, 1.89, 3.05 ] ] )
self.dataCov = numpy.mat([ [ 1.0, 0.1, 0.1 ], [ 0.1, 1.0, 0.1 ] , [ 0.1, 0.1, 1.0 ] ] )
self.constraintVector = numpy.matrix( [[ 42.1 ]] )
self.constraintMatrix = numpy.matrix( [[ 0.0, 0.0, 0.0, 1.0, 1.0 ]] )
self.prior = numpy.matrix([[ 1.3, -0.7, 0.3, 0.11, 42.2 ]])
self.priorCov = numpy.matrix(
[[ 0.07 , 0. , 0. , 0. , 0. ],
[ 0. , 0.5 , 0. , 0. , 0. ],
[ 0. , 0. , 0.4 , 0. , 0. ],
[ 0. , 0. , 0. , 0.1 , 0. ],
[ 0. , 0. , 0. , 0. , 20. ]])
[docs] def test_data_only( self ):
ans, ansCov, ansResid, ansChi2 = cglsqrSolve( self.data, dataUnc = None, dataCov = self.dataCov, kernel = self.kernel, prior = None, priorCov = None, constraintVector = None, constraintMatrix = None )
if False:
print('\ndata_only')
print( 'model\n', ans )
print( 'data\n', self.data )
print( 'residual\n', self.data.T - self.kernel*ans )
print( 'error\n', self.answer.T - ans)
self.assertMatrixAlmostEqual( ( self.answer - ans ).T, numpy.matrix( [[ 0.68651163],[ -0.64302326],[ 0.16651163],[ 0.1 ],[ 42. ]]) )
[docs] def test_data_prior( self ):
ans, ansCov, ansResid, ansChi2 = cglsqrSolve( self.data, dataUnc = None, dataCov = self.dataCov, kernel = self.kernel, prior = self.prior, priorCov = self.priorCov, constraintVector = None, constraintMatrix = None )
if False:
print('\ndata_prior')
print( 'model\n', ans )
print( 'data\n', self.data )
print( 'residual\n', self.data.T - self.kernel*ans )
print( 'error\n', self.answer.T - ans)
self.assertMatrixAlmostEqual( ( self.answer - ans ).T, numpy.matrix([[ 0.04524624],[-0.05105358],[ 0.02455487],[-0.01 ],[-0.2 ]]) )
[docs] def test_data_constraint( self ):
ans, ansCov, ansResid, ansChi2 = cglsqrSolve( self.data, dataUnc = None, dataCov = self.dataCov, kernel = self.kernel, prior = None, priorCov = None, constraintVector = self.constraintVector, constraintMatrix = self.constraintMatrix )
if False:
print('\ndata_constraint')
print( 'model\n', ans )
print( 'data\n', self.data )
print( 'kernel\n', self.kernel )
#print( 'kernel*ans\n', self.kernel*ans.T )
print( 'residual\n', self.data.T - self.kernel*ans.T )
print( 'error\n', self.answer - ans)
print( 'disobey constraint:\n', self.constraintVector - self.constraintMatrix*ans.T)
self.assertMatrixAlmostEqual( ( self.answer - ans ).T, numpy.matrix([[ 0.68651163],[ -0.64302326],[ 0.16651163],[-20.95 ],[ 20.95 ]]) )
[docs] def test_data_constraint_prior( self ):
ans, ansCov, ansResid, ansChi2 = cglsqrSolve( self.data, dataUnc = None, dataCov = self.dataCov, kernel = self.kernel, prior = self.prior, priorCov = self.priorCov, constraintVector = self.constraintVector, constraintMatrix = self.constraintMatrix )
if False:
print('\ndata_constraint_prior')
print( 'model\n', ans )
print( 'data\n', self.data )
print( 'residual\n', self.data.T - self.kernel*ans.T )
print( 'error\n', self.answer - ans)
print( 'disobey constraint:\n', self.constraintVector - self.constraintMatrix*ans.T)
self.assertMatrixAlmostEqual( ( self.answer - ans ).T, numpy.matrix( [[ 0.04524624], [-0.05105358],[ 0.02455487],[-0.00895522],[ 0.00895522]] ) )
# ------------- eigendecomposition tests -----------------
[docs]class Eigendecomposition_base( MatrixTests ):
[docs] def setUpMtx( self ):
raise Exception("Must be overridden by derived classes")
[docs] def setUp( self ):
# The starting matrix, it is symmetric and real, but maybe not positive definite
self.setUpMtx()
self.ndim = self.A.shape[0]
# The eigenvalue decomposition
self.e, self.O = LA.eig( self.A )
self.v = []
for i in range( self.ndim ): self.v.append( self.O.T[i] )
[docs]class EigendecompositionTests:
""" Define some tests that should be run for DERIVED classes only.
This doesn't inherit from unittest.TestCase since the tests aren't meant to be run inside this class.
Derived classes need to inherit from both this and from Eigendecomposition_base
"""
[docs] def test_eigenmodes_are_really_eigenmodes( self ):
"""Check that A * v[i] = e[i] * v[i]"""
for i in range( self.ndim ):
self.assertMatrixAlmostEqual( dot( self.A, self.v[i].T ), self.e[i] * self.v[i].T )
[docs] def test_eigenvector_orthogonality( self ):
"""Eigenvector orthogonality test"""
for i in range( self.ndim ):
for j in range( self.ndim ):
if i == j: continue
self.assertAlmostEqual( ( self.v[i] * self.v[j].T )[0,0], 0.0 )
[docs] def test_eigenvector_normalization( self ):
"""Eigenvector normalization test:"""
for i in range( self.ndim ): self.assertAlmostEqual( ( self.v[i] * self.v[i].T )[0,0], 1.0 )
[docs] def test_eigenvector_orthonormality( self ):
"""Check that eigenvectors are orthonormal, better get back identity matrix here"""
self.assertMatrixAlmostEqual( self.O.T * self.O, identity( self.ndim ) )
self.assertMatrixAlmostEqual( self.O * self.O.T, identity( self.ndim ) )
I = zeros_like( self.A )
for i in range( self.ndim ): I += outer( self.v[i], self.v[i].T )
self.assertMatrixAlmostEqual( I, identity( self.ndim ) )
[docs] def test_reconstruct_matrix_from_eigendecomposition( self ):
"""Try to reconstruct the matrix using the eigenvalue decomposition"""
B = self.O * diag( self.e ) * self.O.T
self.assertMatrixAlmostEqual( self.A, B )
[docs] def test_construct_matrixinverse_from_eigendecomposition( self ):
"""Try to construct the matrix inverse using the eigenvalue decomposition"""
B = self.O * diag( 1.0/self.e ) * self.O.T
self.assertMatrixAlmostEqual( self.A * B, identity( self.ndim ) )
self.assertMatrixAlmostEqual( B * self.A, identity( self.ndim ) )
[docs] @unittest.expectedFailure
def test_construct_matrixinverse_by_pruning( self ):
"""Try to construct the matrix inverse using the PRUNED eigenvalue decomposition"""
with numpy.errstate(divide='ignore'):
B = pruned_matrix_inverse( self.A )
self.assertMatrixAlmostEqual( B * self.A, identity( self.ndim ), rtol=1e-05, atol=1e-08 )
[docs] def test_construct_matrix_by_pruning( self ):
"""Try to reconstruct the matix using the PRUNED eigenvalue decomposition"""
B = pruned_matrix( self.A )
self.assertMatrixAlmostEqual( self.A, B, rtol=1e-04, atol=1e-06 )
[docs]class EigendecompositionTests_allShouldFail( EigendecompositionTests ):
""" For some input matrices, all tests are expected to fail. They should inherit from
this class instead of EigendecompositionTests """
[docs] @unittest.expectedFailure
def test_eigenmodes_are_really_eigenmodes( self ):
super(self).test_eigenmodes_are_really_eigenmodes()
[docs] @unittest.expectedFailure
def test_eigenvector_orthogonality( self ):
super(self).test_eigenvector_orthogonality()
[docs] @unittest.expectedFailure
def test_eigenvector_normalization( self ):
super(self).test_eigenvector_normalization()
[docs] @unittest.expectedFailure
def test_eigenvector_orthonormality( self ):
super(self).test_eigenvector_orthonormality()
[docs] @unittest.expectedFailure
def test_reconstruct_matrix_from_eigendecomposition( self ):
super(self).test_reconstruct_matrix_from_eigendecomposition()
[docs] @unittest.expectedFailure
def test_construct_matrixinverse_from_eigendecomposition( self ):
super(self).test_construct_matrixinverse_from_eigendecomposition()
[docs] @unittest.expectedFailure
def test_construct_matrixinverse_by_pruning( self ):
super(self).test_construct_matrixinverse_by_pruning()
[docs] @unittest.expectedFailure
def test_construct_matrix_by_pruning( self ):
super(self).test_construct_matrix_by_pruning()
# --------- Specific Covariance Test Cases -----------
[docs]class BoringCovarianceTests( Eigendecomposition_base, EigendecompositionTests ):
"""Boring, a diagonal matrix"""
[docs] def setUpMtx( self ): self.A = mat( [[1.,0.,0.],[0.,2.,0.],[0.,0.,3.]] )
[docs]class SomeOffDiagonalCovarianceTests( Eigendecomposition_base, EigendecompositionTests ):
"""Some off-diagonal-ness"""
[docs] def setUpMtx( self ): self.A = mat( [[1.,1.,0.],[1.,2.,0.],[0.,0.,3.]] )
[docs]class LotsOffDiagonalCovarianceTests( Eigendecomposition_base, EigendecompositionTests ):
"""Lots of off-diagonal"""
[docs] def setUpMtx( self ): self.A = mat( [[1.,1.,1.],[1.,2.,1.],[1.,1.,3.]] )
[docs]class ZeroSubspaceCovarianceTests( Eigendecomposition_base, EigendecompositionTests_allShouldFail ):
"""Zeros? it still works!, but has negative eigenvalues"""
[docs] def setUpMtx( self ): self.A = mat( [[0.,0.,1.],[0.,2.,1.],[1.,1.,3.]] )
[docs]class OnDiagonalBarelyPathologicalCovarianceTests( Eigendecomposition_base, EigendecompositionTests_allShouldFail ):
"""Pathological, all correlated, barely illegal and very much not invertible"""
[docs] def setUpMtx( self ): self.A = mat( [[1.,1.,1.],[1.,1.,1.],[1.,1.,0.9999]] )
[docs]class OffDiagonalBarelyPathologicalCovarianceTests( Eigendecomposition_base, EigendecompositionTests_allShouldFail ):
"""Pathological, all correlated, barely illegal and very much not invertible, and slightly off-diagonal"""
[docs] def setUpMtx( self ): self.A = mat( [[1.,1.,1.],[1.,1.,1.00001],[1.,1.,0.9999]] )
[docs]class IllegalCovarianceTest( Eigendecomposition_base, EigendecompositionTests_allShouldFail ):
"""
A bad matrix for testing:
Eigensystem for A. It is OK using numpy.linalg.eig, but not (apparently) with numpy.linalg.eigh, I checked against mathematica:
>>> m = {{1., 2., 3.}, {1., 2., 1.}, {3., 2., 1.}}
>>> {e, v} = Eigensystem[ m ]
>>> {{5.23607, -2., 0.763932}, {{-0.647936, -0.400447, -0.647936}, {-0.707107, 1.08315*10^-16, 0.707107}, {0.465341, -0.752938, 0.465341}}}
numpy.linalg.eig gives:
>>> evals:
>>> [-2.14644241 0.78156156 5.36488085]
>>> evecs:
>>> [[-0.64135318 0.52747554 -0.55716753]
>>> [-0.20230251 -0.81675359 -0.54035845]
>>> [ 0.74009445 0.23384422 -0.63053714]]
However, note that this A is not a valid covariance as the off-diagonal elements are bigger than the on-diagonal ones.
"""
[docs] def setUpMtx( self ): self.A = mat( [[1.,2.,3.],[1.,2.,1.],[3.,2.,1.]] )
[docs] @unittest.expectedFailure
def test_construct_matrix_by_pruning( self ):
super(self).test_construct_matrix_by_pruning()
[docs] @unittest.expectedFailure
def test_construct_matrixinverse_from_eigendecomposition( self ):
super(self).test_construct_matrixinverse_from_eigendecomposition()
[docs]class RealCovarianceTests( Eigendecomposition_base, EigendecompositionTests ):
"""Real covariance for testing, H(n,tot) cs cov."""
[docs] def setUpMtx( self ):
self.A = mat([
[ 8.77542100e-06, 1.38848800e-05, 2.79801400e-05, 3.88661200e-05,
5.25228700e-05, 6.77747100e-05, 7.46361700e-05, 7.52854500e-05,
7.06024400e-05, 6.12296700e-05, 4.77549800e-05, 3.07353700e-05,
1.07019000e-05],
[ 1.38848800e-05, 2.21044000e-05, 4.46560200e-05, 6.20403900e-05,
8.38345200e-05, 1.08155700e-04, 1.19094300e-04, 1.20119000e-04,
1.12637300e-04, 9.76801600e-05, 7.61837900e-05, 4.90328600e-05,
1.70763000e-05],
[ 2.79801400e-05, 4.46560200e-05, 9.03092300e-05, 1.25482100e-04,
1.69564800e-04, 2.18748600e-04, 2.40870000e-04, 2.42936900e-04,
2.27793700e-04, 1.97555100e-04, 1.54076600e-04, 9.91667500e-05,
3.45412500e-05],
[ 3.88661200e-05, 6.20403900e-05, 1.25482100e-04, 1.74359700e-04,
2.35628000e-04, 3.03993100e-04, 3.34742600e-04, 3.37618300e-04,
3.16593500e-04, 2.74561300e-04, 2.14138000e-04, 1.37822000e-04,
4.80018200e-05],
[ 5.25228700e-05, 8.38345200e-05, 1.69564800e-04, 2.35628000e-04,
3.18449800e-04, 4.10888600e-04, 4.52472000e-04, 4.56392600e-04,
4.27988300e-04, 3.71158600e-04, 2.89474100e-04, 1.86312900e-04,
6.48883500e-05],
[ 6.77747100e-05, 1.08155700e-04, 2.18748600e-04, 3.03993100e-04,
4.10888600e-04, 5.30240600e-04, 5.83960200e-04, 5.89030300e-04,
5.52377700e-04, 4.79052800e-04, 3.73630300e-04, 2.40474700e-04,
8.37559700e-05],
[ 7.46361700e-05, 1.19094300e-04, 2.40870000e-04, 3.34742600e-04,
4.52472000e-04, 5.83960200e-04, 6.43158400e-04, 6.48730300e-04,
6.08404400e-04, 5.27657300e-04, 4.11540800e-04, 2.64883200e-04,
9.22628000e-05],
[ 7.52854500e-05, 1.20119000e-04, 2.42936900e-04, 3.37618300e-04,
4.56392600e-04, 5.89030300e-04, 6.48730300e-04, 6.54432800e-04,
6.13741100e-04, 5.32291300e-04, 4.15170400e-04, 2.67219300e-04,
9.30924600e-05],
[ 7.06024400e-05, 1.12637300e-04, 2.27793700e-04, 3.16593500e-04,
4.27988300e-04, 5.52377700e-04, 6.08404400e-04, 6.13741100e-04,
5.75582300e-04, 4.99212100e-04, 3.89369700e-04, 2.50637600e-04,
8.73319300e-05],
[ 6.12296700e-05, 9.76801600e-05, 1.97555100e-04, 2.74561300e-04,
3.71158600e-04, 4.79052800e-04, 5.27657300e-04, 5.32291300e-04,
4.99212100e-04, 4.32971200e-04, 3.37730200e-04, 2.17401500e-04,
7.57803500e-05],
[ 4.77549800e-05, 7.61837900e-05, 1.54076600e-04, 2.14138000e-04,
2.89474100e-04, 3.73630300e-04, 4.11540800e-04, 4.15170400e-04,
3.89369700e-04, 3.37730200e-04, 2.63446700e-04, 1.69609800e-04,
5.91586900e-05],
[ 3.07353700e-05, 4.90328600e-05, 9.91667500e-05, 1.37822000e-04,
1.86312900e-04, 2.40474700e-04, 2.64883200e-04, 2.67219300e-04,
2.50637600e-04, 2.17401500e-04, 1.69609800e-04, 1.09231900e-04,
3.81553500e-05],
[ 1.07019000e-05, 1.70763000e-05, 3.45412500e-05, 4.80018200e-05,
6.48883500e-05, 8.37559700e-05, 9.22628000e-05, 9.30924600e-05,
8.73319300e-05, 7.57803500e-05, 5.91586900e-05, 3.81553500e-05,
1.34232900e-05]] )
if __name__ == "__main__":
unittest.main(verbosity=0)
#In numpy, have around( obj[, decimals[, out] ), _round() too
#numpy.testing.assert_approx_equal(actual, desired, significant=7, err_msg='', verbose=True)