OpenWalnut
1.4.0
|
00001 //--------------------------------------------------------------------------- 00002 // 00003 // Project: OpenWalnut ( http://www.openwalnut.org ) 00004 // 00005 // Copyright 2009 OpenWalnut Community, BSV@Uni-Leipzig and CNCF@MPI-CBS 00006 // For more information see http://www.openwalnut.org/copying 00007 // 00008 // This file is part of OpenWalnut. 00009 // 00010 // OpenWalnut is free software: you can redistribute it and/or modify 00011 // it under the terms of the GNU Lesser General Public License as published by 00012 // the Free Software Foundation, either version 3 of the License, or 00013 // (at your option) any later version. 00014 // 00015 // OpenWalnut is distributed in the hope that it will be useful, 00016 // but WITHOUT ANY WARRANTY; without even the implied warranty of 00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00018 // GNU Lesser General Public License for more details. 00019 // 00020 // You should have received a copy of the GNU Lesser General Public License 00021 // along with OpenWalnut. If not, see <http://www.gnu.org/licenses/>. 00022 // 00023 //--------------------------------------------------------------------------- 00024 00025 #ifndef WLINEARALGEBRAFUNCTIONS_TEST_H 00026 #define WLINEARALGEBRAFUNCTIONS_TEST_H 00027 00028 #include <string> 00029 00030 #include <cxxtest/TestSuite.h> 00031 00032 #include "../../WException.h" 00033 #include "../../WLimits.h" 00034 #include "../WLinearAlgebraFunctions.h" 00035 #include "../WMatrix.h" 00036 #include "../linearAlgebra/WVectorFixed.h" 00037 #include "WVector3dTraits.h" 00038 00039 /** 00040 * Tests for WMatrix. 00041 */ 00042 class WLinearAlgebraFunctionsTest : public CxxTest::TestSuite 00043 { 00044 public: 00045 /** 00046 * The vector multiplied with the matrix is just a new vector where 00047 * each component is the dot product of the corresponding row with the vector. 00048 */ 00049 void testMatrixVectorMultiply( void ) 00050 { 00051 WVector3d v( 9, 10, 11 ); 00052 WMatrix< double > m( 3, 3 ); 00053 int i = 0; 00054 for( size_t r = 0; r < 3; ++r) 00055 { 00056 for( size_t c = 0; c < 3; ++c, ++i ) 00057 { 00058 m( r, c ) = i; 00059 } 00060 } 00061 WVector3d result = multMatrixWithVector3D( m, v ); 00062 WVector3d expected( 32, 122, 212 ); 00063 TS_ASSERT_EQUALS( result, expected ); 00064 } 00065 00066 /** 00067 * If the matrix is not singular then an inverse should exist and be definite. 00068 */ 00069 void test3x3MatrixInversion( void ) 00070 { 00071 int i = 0; 00072 WMatrix< double > m( 3, 3 ); 00073 for( size_t r = 0; r < 3; ++r) 00074 { 00075 for( size_t c = 0; c < 3; ++c, ++i ) 00076 { 00077 m( r, c ) = i; 00078 } 00079 } 00080 m( 2, 1 ) = 8; 00081 WMatrix< double > expected( 3, 3 ); 00082 expected( 0, 0 ) = 1./6 * -8; 00083 expected( 0, 1 ) = 1./6 * 8; 00084 expected( 0, 2 ) = 1./6 * -3; 00085 expected( 1, 0 ) = 1./6 * 6; 00086 expected( 1, 1 ) = 1./6 * -12; 00087 expected( 1, 2 ) = 1./6 * 6; 00088 expected( 2, 0 ) = 1./6 * 0; 00089 expected( 2, 1 ) = 1./6 * 6; 00090 expected( 2, 2 ) = 1./6 * -3; 00091 TS_ASSERT_EQUALS( invertMatrix3x3( m ), expected ); 00092 } 00093 00094 /** 00095 * On singular matrices no inverse exists! 00096 */ 00097 void test3x3MatrixInversionOnSingularMatrix( void ) 00098 { 00099 int i = 0; 00100 WMatrix< double > m( 3, 3 ); 00101 for( size_t r = 0; r < 3; ++r) 00102 { 00103 for( size_t c = 0; c < 3; ++c, ++i ) 00104 { 00105 m( r, c ) = i; 00106 } 00107 } 00108 TS_ASSERT_THROWS( invertMatrix3x3( m ), WException& e ); 00109 // ATTENTION we can't compare the messages from WAssert since there is now a path string in side which is different on any developers machine 00110 // "Assertion failed: det != 0 (in file /home/lmath/repos/OpenWalnut/src/common/math/WLinearAlgebraFunctions.cpp at line 71),.... 00111 } 00112 00113 /** 00114 * Test the inversion of 4x4 matrices 00115 */ 00116 void test4x4Inverse() 00117 { 00118 WMatrix<double> m( 4, 4 ); 00119 00120 for( size_t r = 0; r < 4; ++r) 00121 { 00122 for( size_t c = 0; c < 4; ++c ) 00123 { 00124 m( c, r ) = r + c * 4 + 1; 00125 } 00126 } 00127 m( 2, 2 ) = 12; 00128 m( 3, 3 ) = 15; 00129 00130 WMatrix<double> m_inv = invertMatrix4x4( m ); 00131 00132 WMatrix<double> id = WMatrix<double>( 4, 4 ).makeIdentity(); 00133 00134 WMatrix<double> m_m_inv = m * m_inv; 00135 00136 TS_ASSERT( m_m_inv == id ); 00137 } 00138 00139 /** 00140 * Two vectors are linear independent if the are not parallel 00141 */ 00142 void testLinearIndependeceOfTwoVectors( void ) 00143 { 00144 WVector3d u( 1, 0, 0 ); 00145 WVector3d v( 0, 1, 0 ); 00146 TS_ASSERT( linearIndependent( u, v ) ); 00147 TS_ASSERT( linearIndependent( v, u ) ); 00148 TS_ASSERT( !linearIndependent( v, v ) ); 00149 } 00150 00151 /** 00152 * Two vectors are linear independent if the are not parallel 00153 */ 00154 void testLinearIndependeceOfTheNullVector( void ) 00155 { 00156 WVector3d u( 0, 0, 0 ); 00157 WVector3d v( 0, 0, 1 ); 00158 TS_ASSERT( !linearIndependent( u, v ) ); 00159 TS_ASSERT( !linearIndependent( v, u ) ); 00160 TS_ASSERT( !linearIndependent( u, u ) ); 00161 } 00162 00163 /** 00164 * Small changes should nothing do to correctness 00165 */ 00166 void testLinearIndependenceOnNumericalStability( void ) 00167 { 00168 WVector3d u( wlimits::DBL_EPS, wlimits::DBL_EPS, wlimits::DBL_EPS ); 00169 WVector3d v( wlimits::DBL_EPS, wlimits::DBL_EPS, 1 ); 00170 TS_ASSERT( !linearIndependent( u, v ) ); 00171 TS_ASSERT( !linearIndependent( v, u ) ); 00172 TS_ASSERT( !linearIndependent( u, u ) ); 00173 u[0] = wlimits::DBL_EPS * 2; 00174 TS_ASSERT( linearIndependent( u, v ) ); 00175 TS_ASSERT( linearIndependent( v, u ) ); 00176 TS_ASSERT( !linearIndependent( u, u ) ); 00177 } 00178 00179 /** 00180 * Test SVD calculation 00181 */ 00182 void testComputeSVD( void ) 00183 { 00184 const size_t nbRows = 3, nbCols = 3; 00185 const double a = 1.2, b = 2.3, c = 3.4, 00186 d = 4.5, e = 5.6, f = 6.7, 00187 g = 3.4, h = 1.2, i = 7.0; 00188 WMatrix<double> A( nbRows, nbCols ); 00189 A( 0, 0 ) = a; 00190 A( 0, 1 ) = b; 00191 A( 0, 2 ) = c; 00192 A( 1, 0 ) = d; 00193 A( 1, 1 ) = e; 00194 A( 1, 2 ) = f; 00195 A( 2, 0 ) = g; 00196 A( 2, 1 ) = h; 00197 A( 2, 2 ) = i; 00198 WMatrix<double> U( A.getNbRows(), A.getNbCols() ); 00199 WMatrix<double> V( A.getNbCols(), A.getNbCols() ); 00200 WValue<double> Svec( A.getNbCols() ); 00201 computeSVD( A, U, V, Svec ); 00202 WMatrix<double> S( Svec.size(), Svec.size() ); 00203 S.setZero(); 00204 for( size_t i = 0; i < Svec.size(); ++i ) 00205 { 00206 S( i, i ) = Svec[ i ]; 00207 } 00208 00209 WMatrix<double> A2( U*S*V.transposed() ); 00210 00211 for( size_t row = 0; row < A.getNbRows(); ++row ) 00212 { 00213 for( size_t col = 0; col < A.getNbCols(); ++col ) 00214 { 00215 TS_ASSERT_DELTA( A( row, col ), A2( row, col ), 0.0001 ); 00216 } 00217 } 00218 } 00219 00220 /** 00221 * Test pseudoInverse calculation 00222 */ 00223 void testPseudoInverse( void ) 00224 { 00225 { 00226 const size_t nbRows = 3, nbCols = 3; 00227 const double a = 1.2, b = 2.3, c = 3.4, 00228 d = 4.5, e = 5.6, f = 6.7, 00229 g = 3.4, h = 1.2, i = 7.0; 00230 WMatrix<double> A( nbRows, nbCols ); 00231 00232 A( 0, 0 ) = a; 00233 A( 0, 1 ) = b; 00234 A( 0, 2 ) = c; 00235 A( 1, 0 ) = d; 00236 A( 1, 1 ) = e; 00237 A( 1, 2 ) = f; 00238 A( 2, 0 ) = g; 00239 A( 2, 1 ) = h; 00240 A( 2, 2 ) = i; 00241 WMatrix<double> Ainvers( pseudoInverse( A ) ); 00242 WMatrix<double> I( A*Ainvers ); 00243 00244 for( size_t row = 0; row < I.getNbRows(); row++ ) 00245 { 00246 for( size_t col = 0; col < I.getNbCols(); col++ ) 00247 { 00248 if( row == col ) 00249 { 00250 TS_ASSERT_DELTA( I( row, col ), 1.0, 0.0001 ); 00251 } 00252 else 00253 { 00254 TS_ASSERT_DELTA( I( row, col ), 0.0, 0.0001 ); 00255 } 00256 } 00257 } 00258 } 00259 { 00260 WMatrix<double> m( 6, 6 ); 00261 for( int j = 0; j < 6; ++j ) 00262 { 00263 for( int i = 0; i < 6; ++i ) 00264 { 00265 m( i, j ) = pow( i + 1, j ); 00266 } 00267 } 00268 m = m * pseudoInverse( m ); 00269 for( int j = 0; j < 6; ++j ) 00270 { 00271 for( int i = 0; i < 6; ++i ) 00272 { 00273 TS_ASSERT_DELTA( m( i, j ), i == j ? 1.0 : 0.0, 0.0001 ); 00274 } 00275 } 00276 } 00277 } 00278 }; 00279 00280 #endif // WLINEARALGEBRAFUNCTIONS_TEST_H