OpenWalnut  1.4.0
WLinearAlgebraFunctions_test.h
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