OpenWalnut  1.3.1
WLinearAlgebraFunctions_test.h
1 //---------------------------------------------------------------------------
2 //
3 // Project: OpenWalnut ( http://www.openwalnut.org )
4 //
5 // Copyright 2009 OpenWalnut Community, BSV@Uni-Leipzig and CNCF@MPI-CBS
6 // For more information see http://www.openwalnut.org/copying
7 //
8 // This file is part of OpenWalnut.
9 //
10 // OpenWalnut is free software: you can redistribute it and/or modify
11 // it under the terms of the GNU Lesser General Public License as published by
12 // the Free Software Foundation, either version 3 of the License, or
13 // (at your option) any later version.
14 //
15 // OpenWalnut is distributed in the hope that it will be useful,
16 // but WITHOUT ANY WARRANTY; without even the implied warranty of
17 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18 // GNU Lesser General Public License for more details.
19 //
20 // You should have received a copy of the GNU Lesser General Public License
21 // along with OpenWalnut. If not, see <http://www.gnu.org/licenses/>.
22 //
23 //---------------------------------------------------------------------------
24 
25 #ifndef WLINEARALGEBRAFUNCTIONS_TEST_H
26 #define WLINEARALGEBRAFUNCTIONS_TEST_H
27 
28 #include <string>
29 
30 #include <cxxtest/TestSuite.h>
31 
32 #include "../../WException.h"
33 #include "../../WLimits.h"
34 #include "../WLinearAlgebraFunctions.h"
35 #include "../WMatrix.h"
36 #include "../linearAlgebra/WLinearAlgebra.h"
37 #include "WVector3dTraits.h"
38 
39 /**
40  * Tests for WMatrix.
41  */
42 class WLinearAlgebraFunctionsTest : public CxxTest::TestSuite
43 {
44 public:
45  /**
46  * The vector multiplied with the matrix is just a new vector where
47  * each component is the dot product of the corresponding row with the vector.
48  */
50  {
51  WVector3d v( 9, 10, 11 );
52  WMatrix< double > m( 3, 3 );
53  int i = 0;
54  for( size_t r = 0; r < 3; ++r)
55  {
56  for( size_t c = 0; c < 3; ++c, ++i )
57  {
58  m( r, c ) = i;
59  }
60  }
61  WVector3d result = multMatrixWithVector3D( m, v );
62  WVector3d expected( 32, 122, 212 );
63  TS_ASSERT_EQUALS( result, expected );
64  }
65 
66  /**
67  * If the matrix is not singular then an inverse should exist and be definite.
68  */
70  {
71  int i = 0;
72  WMatrix< double > m( 3, 3 );
73  for( size_t r = 0; r < 3; ++r)
74  {
75  for( size_t c = 0; c < 3; ++c, ++i )
76  {
77  m( r, c ) = i;
78  }
79  }
80  m( 2, 1 ) = 8;
81  WMatrix< double > expected( 3, 3 );
82  expected( 0, 0 ) = 1./6 * -8;
83  expected( 0, 1 ) = 1./6 * 8;
84  expected( 0, 2 ) = 1./6 * -3;
85  expected( 1, 0 ) = 1./6 * 6;
86  expected( 1, 1 ) = 1./6 * -12;
87  expected( 1, 2 ) = 1./6 * 6;
88  expected( 2, 0 ) = 1./6 * 0;
89  expected( 2, 1 ) = 1./6 * 6;
90  expected( 2, 2 ) = 1./6 * -3;
91  TS_ASSERT_EQUALS( invertMatrix3x3( m ), expected );
92  }
93 
94  /**
95  * On singular matrices no inverse exists!
96  */
98  {
99  int i = 0;
100  WMatrix< double > m( 3, 3 );
101  for( size_t r = 0; r < 3; ++r)
102  {
103  for( size_t c = 0; c < 3; ++c, ++i )
104  {
105  m( r, c ) = i;
106  }
107  }
108  TS_ASSERT_THROWS( invertMatrix3x3( m ), WException& e );
109  // 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
110  // "Assertion failed: det != 0 (in file /home/lmath/repos/OpenWalnut/src/common/math/WLinearAlgebraFunctions.cpp at line 71),....
111  }
112 
113  /**
114  * Test the inversion of 4x4 matrices
115  */
117  {
118  WMatrix<double> m( 4, 4 );
119 
120  for( size_t r = 0; r < 4; ++r)
121  {
122  for( size_t c = 0; c < 4; ++c )
123  {
124  m( c, r ) = r + c * 4 + 1;
125  }
126  }
127  m( 2, 2 ) = 12;
128  m( 3, 3 ) = 15;
129 
130  WMatrix<double> m_inv = invertMatrix4x4( m );
131 
132  WMatrix<double> id = WMatrix<double>( 4, 4 ).makeIdentity();
133 
134  WMatrix<double> m_m_inv = m * m_inv;
135 
136  TS_ASSERT( m_m_inv == id );
137  }
138 
139  /**
140  * Two vectors are linear independent if the are not parallel
141  */
143  {
144  WVector3d u( 1, 0, 0 );
145  WVector3d v( 0, 1, 0 );
146  TS_ASSERT( linearIndependent( u, v ) );
147  TS_ASSERT( linearIndependent( v, u ) );
148  TS_ASSERT( !linearIndependent( v, v ) );
149  }
150 
151  /**
152  * Two vectors are linear independent if the are not parallel
153  */
155  {
156  WVector3d u( 0, 0, 0 );
157  WVector3d v( 0, 0, 1 );
158  TS_ASSERT( !linearIndependent( u, v ) );
159  TS_ASSERT( !linearIndependent( v, u ) );
160  TS_ASSERT( !linearIndependent( u, u ) );
161  }
162 
163  /**
164  * Small changes should nothing do to correctness
165  */
167  {
170  TS_ASSERT( !linearIndependent( u, v ) );
171  TS_ASSERT( !linearIndependent( v, u ) );
172  TS_ASSERT( !linearIndependent( u, u ) );
173  u[0] = wlimits::DBL_EPS * 2;
174  TS_ASSERT( linearIndependent( u, v ) );
175  TS_ASSERT( linearIndependent( v, u ) );
176  TS_ASSERT( !linearIndependent( u, u ) );
177  }
178 
179  /**
180  * Test SVD calculation
181  */
182  void testComputeSVD( void )
183  {
184  const size_t nbRows = 3, nbCols = 3;
185  const double a = 1.2, b = 2.3, c = 3.4,
186  d = 4.5, e = 5.6, f = 6.7,
187  g = 3.4, h = 1.2, i = 7.0;
188  WMatrix<double> A( nbRows, nbCols );
189  A( 0, 0 ) = a;
190  A( 0, 1 ) = b;
191  A( 0, 2 ) = c;
192  A( 1, 0 ) = d;
193  A( 1, 1 ) = e;
194  A( 1, 2 ) = f;
195  A( 2, 0 ) = g;
196  A( 2, 1 ) = h;
197  A( 2, 2 ) = i;
198  WMatrix<double> U( A.getNbRows(), A.getNbCols() );
199  WMatrix<double> V( A.getNbCols(), A.getNbCols() );
200  WValue<double> Svec( A.getNbCols() );
201  computeSVD( A, U, V, Svec );
202  WMatrix<double> S( Svec.size(), Svec.size() );
203  S.setZero();
204  for( size_t i = 0; i < Svec.size(); ++i )
205  {
206  S( i, i ) = Svec[ i ];
207  }
208 
209  WMatrix<double> A2( U*S*V.transposed() );
210 
211  for( size_t row = 0; row < A.getNbRows(); ++row )
212  {
213  for( size_t col = 0; col < A.getNbCols(); ++col )
214  {
215  TS_ASSERT_DELTA( A( row, col ), A2( row, col ), 0.0001 );
216  }
217  }
218  }
219 
220  /**
221  * Test pseudoInverse calculation
222  */
223  void testPseudoInverse( void )
224  {
225  {
226  const size_t nbRows = 3, nbCols = 3;
227  const double a = 1.2, b = 2.3, c = 3.4,
228  d = 4.5, e = 5.6, f = 6.7,
229  g = 3.4, h = 1.2, i = 7.0;
230  WMatrix<double> A( nbRows, nbCols );
231 
232  A( 0, 0 ) = a;
233  A( 0, 1 ) = b;
234  A( 0, 2 ) = c;
235  A( 1, 0 ) = d;
236  A( 1, 1 ) = e;
237  A( 1, 2 ) = f;
238  A( 2, 0 ) = g;
239  A( 2, 1 ) = h;
240  A( 2, 2 ) = i;
241  WMatrix<double> Ainvers( pseudoInverse( A ) );
242  WMatrix<double> I( A*Ainvers );
243 
244  for( size_t row = 0; row < I.getNbRows(); row++ )
245  {
246  for( size_t col = 0; col < I.getNbCols(); col++ )
247  {
248  if( row == col )
249  {
250  TS_ASSERT_DELTA( I( row, col ), 1.0, 0.0001 );
251  }
252  else
253  {
254  TS_ASSERT_DELTA( I( row, col ), 0.0, 0.0001 );
255  }
256  }
257  }
258  }
259  {
260  WMatrix<double> m( 6, 6 );
261  for( int j = 0; j < 6; ++j )
262  {
263  for( int i = 0; i < 6; ++i )
264  {
265  m( i, j ) = pow( i + 1, j );
266  }
267  }
268  m = m * pseudoInverse( m );
269  for( int j = 0; j < 6; ++j )
270  {
271  for( int i = 0; i < 6; ++i )
272  {
273  TS_ASSERT_DELTA( m( i, j ), i == j ? 1.0 : 0.0, 0.0001 );
274  }
275  }
276  }
277  }
278 };
279 
280 #endif // WLINEARALGEBRAFUNCTIONS_TEST_H