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 WTENSORFUNCTIONS_TEST_H 00026 #define WTENSORFUNCTIONS_TEST_H 00027 00028 #include <string> 00029 #include <vector> 00030 #include <algorithm> 00031 00032 #include <cxxtest/TestSuite.h> 00033 00034 #include "../../WAssert.h" 00035 #include "../../WLimits.h" 00036 #include "../../WStringUtils.h" 00037 #include "../WTensorFunctions.h" 00038 00039 /** 00040 * Test class for some tensor functions. 00041 */ 00042 class WTensorFunctionsTest : public CxxTest::TestSuite 00043 { 00044 public: 00045 /** 00046 * The eigenvalue of the symmetrical matrix: 00047 * 0.000179516, 2.09569e-05, 2.76557e-06, 0.000170189, -5.52619e-07, 0.00015239 00048 * (0.000196397;0.000155074;0.000150625) 00049 */ 00050 void testSpecialSymMatrixEigenvalueTestCaseNumericalStability( void ) 00051 { 00052 WTensorSym< 2, 3 > t; 00053 t( 0, 0 ) = 0.000179516; 00054 t( 0, 1 ) = 2.09569e-05; 00055 t( 0, 2 ) = 2.76557e-06; 00056 t( 1, 1 ) = 0.000170189; 00057 t( 1, 2 ) = -5.52619e-07; 00058 t( 2, 2 ) = 0.00015239; 00059 RealEigenSystem sys; 00060 jacobiEigenvector3D( t, &sys ); 00061 00062 TS_ASSERT_DELTA( sys[0].first, 1.5062467240725114e-04, 1e-9 ); 00063 TS_ASSERT_DELTA( sys[1].first, 1.5507354000104679e-04, 1e-9 ); 00064 TS_ASSERT_DELTA( sys[2].first, 1.9639678759170208e-04, 1e-9 ); 00065 } 00066 00067 /** 00068 * Test the jacobi eigenvector calculation. 00069 */ 00070 void testJacobiEigenvectors() 00071 { 00072 // the test matrix 00073 WTensorSym< 2, 3 > t; 00074 00075 // simple diagonal matrices should be decomposed correctly 00076 // 1 2 3 00077 t( 0, 0 ) = 1.0; 00078 t( 1, 1 ) = 2.0; 00079 t( 2, 2 ) = 3.0; 00080 00081 RealEigenSystem sys; 00082 jacobiEigenvector3D( t, &sys ); 00083 00084 TS_ASSERT_DELTA( sys[0].first, 1.0, 1e-6 ); 00085 TS_ASSERT_DELTA( sys[1].first, 2.0, 1e-6 ); 00086 TS_ASSERT_DELTA( sys[2].first, 3.0, 1e-6 ); 00087 // eigenvectors should be perpendicular 00088 TS_ASSERT_DELTA( dot( sys[0].second, sys[1].second ), 0.0, 1e-9 ); 00089 TS_ASSERT_DELTA( dot( sys[1].second, sys[2].second ), 0.0, 1e-9 ); 00090 TS_ASSERT_DELTA( dot( sys[2].second, sys[0].second ), 0.0, 1e-9 ); 00091 TS_ASSERT_DELTA( length( sys[0].second ), 1.0, 1e-9 ); 00092 TS_ASSERT_DELTA( length( sys[1].second ), 1.0, 1e-9 ); 00093 TS_ASSERT_DELTA( length( sys[2].second ), 1.0, 1e-9 ); 00094 00095 // 1 2 -3 00096 t( 0, 0 ) = 1.0; 00097 t( 1, 1 ) = 2.0; 00098 t( 2, 2 ) = -3.0; 00099 00100 jacobiEigenvector3D( t, &sys ); 00101 00102 TS_ASSERT_DELTA( sys[0].first, -3.0, 1e-6 ); 00103 TS_ASSERT_DELTA( sys[1].first, 1.0, 1e-6 ); 00104 TS_ASSERT_DELTA( sys[2].first, 2.0, 1e-6 ); 00105 TS_ASSERT_DELTA( dot( sys[0].second, sys[1].second ), 0.0, 1e-9 ); 00106 TS_ASSERT_DELTA( dot( sys[1].second, sys[2].second ), 0.0, 1e-9 ); 00107 TS_ASSERT_DELTA( dot( sys[2].second, sys[0].second ), 0.0, 1e-9 ); 00108 TS_ASSERT_DELTA( length( sys[0].second ), 1.0, 1e-9 ); 00109 TS_ASSERT_DELTA( length( sys[1].second ), 1.0, 1e-9 ); 00110 TS_ASSERT_DELTA( length( sys[2].second ), 1.0, 1e-9 ); 00111 00112 // 1 2 2 00113 t( 0, 0 ) = 1.0; 00114 t( 1, 1 ) = 2.0; 00115 t( 2, 2 ) = 2.0; 00116 00117 jacobiEigenvector3D( t, &sys ); 00118 00119 TS_ASSERT_DELTA( sys[0].first, 1.0, 1e-6 ); 00120 TS_ASSERT_DELTA( sys[1].first, 2.0, 1e-6 ); 00121 TS_ASSERT_DELTA( sys[2].first, 2.0, 1e-6 ); 00122 TS_ASSERT_DELTA( dot( sys[0].second, sys[1].second ), 0.0, 1e-9 ); 00123 TS_ASSERT_DELTA( dot( sys[1].second, sys[2].second ), 0.0, 1e-9 ); 00124 TS_ASSERT_DELTA( dot( sys[2].second, sys[0].second ), 0.0, 1e-9 ); 00125 TS_ASSERT_DELTA( length( sys[0].second ), 1.0, 1e-9 ); 00126 TS_ASSERT_DELTA( length( sys[1].second ), 1.0, 1e-9 ); 00127 TS_ASSERT_DELTA( length( sys[2].second ), 1.0, 1e-9 ); 00128 00129 // -1 -1 -1 00130 t( 0, 0 ) = -1.0; 00131 t( 1, 1 ) = -1.0; 00132 t( 2, 2 ) = -1.0; 00133 00134 jacobiEigenvector3D( t, &sys ); 00135 00136 TS_ASSERT_DELTA( sys[0].first, -1.0, 1e-6 ); 00137 TS_ASSERT_DELTA( sys[1].first, -1.0, 1e-6 ); 00138 TS_ASSERT_DELTA( sys[2].first, -1.0, 1e-6 ); 00139 TS_ASSERT_DELTA( dot( sys[0].second, sys[1].second ), 0.0, 1e-9 ); 00140 TS_ASSERT_DELTA( dot( sys[1].second, sys[2].second ), 0.0, 1e-9 ); 00141 TS_ASSERT_DELTA( dot( sys[2].second, sys[0].second ), 0.0, 1e-9 ); 00142 TS_ASSERT_DELTA( length( sys[0].second ), 1.0, 1e-9 ); 00143 TS_ASSERT_DELTA( length( sys[1].second ), 1.0, 1e-9 ); 00144 TS_ASSERT_DELTA( length( sys[2].second ), 1.0, 1e-9 ); 00145 00146 // 1 0 1 00147 t( 0, 0 ) = 1.0; 00148 t( 1, 1 ) = 0.0; 00149 t( 2, 2 ) = 1.0; 00150 00151 jacobiEigenvector3D( t, &sys ); 00152 00153 TS_ASSERT_DELTA( sys[0].first, 0.0, 1e-6 ); 00154 TS_ASSERT_DELTA( sys[1].first, 1.0, 1e-6 ); 00155 TS_ASSERT_DELTA( sys[2].first, 1.0, 1e-6 ); 00156 TS_ASSERT_DELTA( dot( sys[0].second, sys[1].second ), 0.0, 1e-9 ); 00157 TS_ASSERT_DELTA( dot( sys[1].second, sys[2].second ), 0.0, 1e-9 ); 00158 TS_ASSERT_DELTA( dot( sys[2].second, sys[0].second ), 0.0, 1e-9 ); 00159 TS_ASSERT_DELTA( length( sys[0].second ), 1.0, 1e-9 ); 00160 TS_ASSERT_DELTA( length( sys[1].second ), 1.0, 1e-9 ); 00161 TS_ASSERT_DELTA( length( sys[2].second ), 1.0, 1e-9 ); 00162 00163 // 0 0 0 00164 t( 0, 0 ) = 0.0; 00165 t( 1, 1 ) = 0.0; 00166 t( 2, 2 ) = 0.0; 00167 00168 jacobiEigenvector3D( t, &sys ); 00169 00170 TS_ASSERT_DELTA( sys[0].first, 0.0, 1e-6 ); 00171 TS_ASSERT_DELTA( sys[1].first, 0.0, 1e-6 ); 00172 TS_ASSERT_DELTA( sys[2].first, 0.0, 1e-6 ); 00173 TS_ASSERT_DELTA( dot( sys[0].second, sys[1].second ), 0.0, 1e-9 ); 00174 TS_ASSERT_DELTA( dot( sys[1].second, sys[2].second ), 0.0, 1e-9 ); 00175 TS_ASSERT_DELTA( dot( sys[2].second, sys[0].second ), 0.0, 1e-9 ); 00176 TS_ASSERT_DELTA( length( sys[0].second ), 1.0, 1e-9 ); 00177 TS_ASSERT_DELTA( length( sys[1].second ), 1.0, 1e-9 ); 00178 TS_ASSERT_DELTA( length( sys[2].second ), 1.0, 1e-9 ); 00179 00180 // similar eigenvalues 00181 // 2.000001 0.0 1.999998 00182 t( 0, 0 ) = 2.000001; 00183 t( 1, 1 ) = 0.0; 00184 t( 2, 2 ) = 1.999998; 00185 00186 jacobiEigenvector3D( t, &sys ); 00187 00188 TS_ASSERT_DELTA( sys[0].first, 0.0, 1e-6 ); 00189 TS_ASSERT_DELTA( sys[1].first, 1.999998, 1e-6 ); 00190 TS_ASSERT_DELTA( sys[2].first, 2.000001, 1e-6 ); 00191 TS_ASSERT_DELTA( dot( sys[0].second, sys[1].second ), 0.0, 1e-9 ); 00192 TS_ASSERT_DELTA( dot( sys[1].second, sys[2].second ), 0.0, 1e-9 ); 00193 TS_ASSERT_DELTA( dot( sys[2].second, sys[0].second ), 0.0, 1e-9 ); 00194 TS_ASSERT_DELTA( length( sys[0].second ), 1.0, 1e-9 ); 00195 TS_ASSERT_DELTA( length( sys[1].second ), 1.0, 1e-9 ); 00196 TS_ASSERT_DELTA( length( sys[2].second ), 1.0, 1e-9 ); 00197 00198 // very large eigenvalues 00199 // 3.824572321236e1000 1 2 00200 t( 0, 0 ) = 3.824572321236e30; 00201 t( 1, 1 ) = 1.0; 00202 t( 2, 2 ) = 2.0; 00203 00204 jacobiEigenvector3D( t, &sys ); 00205 00206 TS_ASSERT_DELTA( sys[0].first, 1.0, 1e-6 ); 00207 TS_ASSERT_DELTA( sys[1].first, 2.0, 1e-6 ); 00208 TS_ASSERT_DELTA( sys[2].first, 3.824572321236e30, 1e-6 ); 00209 TS_ASSERT_DELTA( dot( sys[0].second, sys[1].second ), 0.0, 1e-9 ); 00210 TS_ASSERT_DELTA( dot( sys[1].second, sys[2].second ), 0.0, 1e-9 ); 00211 TS_ASSERT_DELTA( dot( sys[2].second, sys[0].second ), 0.0, 1e-9 ); 00212 TS_ASSERT_DELTA( length( sys[0].second ), 1.0, 1e-9 ); 00213 TS_ASSERT_DELTA( length( sys[1].second ), 1.0, 1e-9 ); 00214 TS_ASSERT_DELTA( length( sys[2].second ), 1.0, 1e-9 ); 00215 00216 // very small eigenvalues 00217 // 3.824572321236e-1000 1 2 00218 t( 0, 0 ) = 3.824572321236e-30; 00219 t( 1, 1 ) = 1.0; 00220 t( 2, 2 ) = 2.0; 00221 00222 jacobiEigenvector3D( t, &sys ); 00223 00224 TS_ASSERT_DELTA( sys[0].first, 3.824572321236e-30, 1e-6 ); 00225 TS_ASSERT_DELTA( sys[1].first, 1.0, 1e-6 ); 00226 TS_ASSERT_DELTA( sys[2].first, 2.0, 1e-6 ); 00227 TS_ASSERT_DELTA( dot( sys[0].second, sys[1].second ), 0.0, 1e-9 ); 00228 TS_ASSERT_DELTA( dot( sys[1].second, sys[2].second ), 0.0, 1e-9 ); 00229 TS_ASSERT_DELTA( dot( sys[2].second, sys[0].second ), 0.0, 1e-9 ); 00230 TS_ASSERT_DELTA( length( sys[0].second ), 1.0, 1e-9 ); 00231 TS_ASSERT_DELTA( length( sys[1].second ), 1.0, 1e-9 ); 00232 TS_ASSERT_DELTA( length( sys[2].second ), 1.0, 1e-9 ); 00233 00234 // some more sophisticated tests 00235 // (using similarity transformations on diagonal matrices to create test cases) 00236 t( 0, 0 ) = 1; 00237 t( 1, 1 ) = 2; 00238 t( 2, 2 ) = 3; 00239 00240 // note that the jacobi eigenvector functions does not sort the eigenvalues and 00241 // eigenvectors that were found 00242 t = similarity_rotate_givens( t, 0, 2, 2.78 ); 00243 00244 jacobiEigenvector3D( t, &sys ); 00245 00246 TS_ASSERT_DELTA( sys[0].first, 1.0, 1e-6 ); 00247 TS_ASSERT_DELTA( sys[1].first, 2.0, 1e-6 ); 00248 TS_ASSERT_DELTA( sys[2].first, 3.0, 1e-6 ); 00249 TS_ASSERT_DELTA( dot( sys[0].second, sys[1].second ), 0.0, 1e-9 ); 00250 TS_ASSERT_DELTA( dot( sys[1].second, sys[2].second ), 0.0, 1e-9 ); 00251 TS_ASSERT_DELTA( dot( sys[2].second, sys[0].second ), 0.0, 1e-9 ); 00252 TS_ASSERT_DELTA( length( sys[0].second ), 1.0, 1e-9 ); 00253 TS_ASSERT_DELTA( length( sys[1].second ), 1.0, 1e-9 ); 00254 TS_ASSERT_DELTA( length( sys[2].second ), 1.0, 1e-9 ); 00255 compare_results( t, sys ); 00256 00257 t = WTensorSym< 2, 3 >(); 00258 t( 0, 0 ) = 2; 00259 t( 1, 1 ) = 1; 00260 t( 2, 2 ) = 3; 00261 00262 t = similarity_rotate_givens( t, 0, 2, 2.79 ); 00263 t = similarity_rotate_givens( t, 1, 2, -3.44 ); 00264 t = similarity_rotate_givens( t, 1, 0, -0.46 ); 00265 t = similarity_rotate_givens( t, 2, 1, 5.98 ); 00266 00267 jacobiEigenvector3D( t, &sys ); 00268 00269 TS_ASSERT_DELTA( sys[0].first, 1.0, 1e-6 ); 00270 TS_ASSERT_DELTA( sys[1].first, 2.0, 1e-6 ); 00271 TS_ASSERT_DELTA( sys[2].first, 3.0, 1e-6 ); 00272 TS_ASSERT_DELTA( dot( sys[0].second, sys[1].second ), 0.0, 1e-9 ); 00273 TS_ASSERT_DELTA( dot( sys[1].second, sys[2].second ), 0.0, 1e-9 ); 00274 TS_ASSERT_DELTA( dot( sys[2].second, sys[0].second ), 0.0, 1e-9 ); 00275 TS_ASSERT_DELTA( length( sys[0].second ), 1.0, 1e-9 ); 00276 TS_ASSERT_DELTA( length( sys[1].second ), 1.0, 1e-9 ); 00277 TS_ASSERT_DELTA( length( sys[2].second ), 1.0, 1e-9 ); 00278 compare_results( t, sys ); 00279 00280 t = WTensorSym< 2, 3 >(); 00281 t( 0, 0 ) = 2; 00282 t( 1, 1 ) = 2; 00283 t( 2, 2 ) = 2; 00284 00285 t = similarity_rotate_givens( t, 1, 2, -3.44 ); 00286 t = similarity_rotate_givens( t, 2, 1, 5.98 ); 00287 t = similarity_rotate_givens( t, 0, 2, 2.79 ); 00288 t = similarity_rotate_givens( t, 1, 0, -0.46 ); 00289 00290 jacobiEigenvector3D( t, &sys ); 00291 00292 TS_ASSERT_DELTA( sys[0].first, 2.0, 1e-6 ); 00293 TS_ASSERT_DELTA( sys[1].first, 2.0, 1e-6 ); 00294 TS_ASSERT_DELTA( sys[2].first, 2.0, 1e-6 ); 00295 TS_ASSERT_DELTA( dot( sys[0].second, sys[1].second ), 0.0, 1e-9 ); 00296 TS_ASSERT_DELTA( dot( sys[1].second, sys[2].second ), 0.0, 1e-9 ); 00297 TS_ASSERT_DELTA( dot( sys[2].second, sys[0].second ), 0.0, 1e-9 ); 00298 TS_ASSERT_DELTA( length( sys[0].second ), 1.0, 1e-9 ); 00299 TS_ASSERT_DELTA( length( sys[1].second ), 1.0, 1e-9 ); 00300 TS_ASSERT_DELTA( length( sys[2].second ), 1.0, 1e-9 ); 00301 compare_results( t, sys ); 00302 } 00303 00304 /** 00305 * Test the cardano eigenvalue calculation. 00306 */ 00307 void testCardanoEigenvalues() 00308 { 00309 // the test matrix 00310 WTensorSym< 2, 3 > t; 00311 // variables for the output values 00312 std::vector< double > d( 3 ); 00313 00314 // simple diagonal matrices should be decomposed correctly 00315 00316 // 1 2 3 00317 t( 0, 0 ) = 1.0; 00318 t( 1, 1 ) = 2.0; 00319 t( 2, 2 ) = 3.0; 00320 00321 d = getEigenvaluesCardano( t ); 00322 00323 TS_ASSERT_DELTA( d[ 0 ], 3.0, 1e-6 ); 00324 TS_ASSERT_DELTA( d[ 1 ], 2.0, 1e-6 ); 00325 TS_ASSERT_DELTA( d[ 2 ], 1.0, 1e-6 ); 00326 00327 // 1 2 -3 00328 t( 0, 0 ) = 1.0; 00329 t( 1, 1 ) = 2.0; 00330 t( 2, 2 ) = -3.0; 00331 00332 d = getEigenvaluesCardano( t ); 00333 00334 TS_ASSERT_DELTA( d[ 0 ], 2.0, 1e-6 ); 00335 TS_ASSERT_DELTA( d[ 1 ], 1.0, 1e-6 ); 00336 TS_ASSERT_DELTA( d[ 2 ], -3.0, 1e-6 ); 00337 00338 // 1 2 2 00339 t( 0, 0 ) = 1.0; 00340 t( 1, 1 ) = 2.0; 00341 t( 2, 2 ) = 2.0; 00342 00343 d = getEigenvaluesCardano( t ); 00344 00345 TS_ASSERT_DELTA( d[ 0 ], 2.0, 1e-6 ); 00346 TS_ASSERT_DELTA( d[ 1 ], 2.0, 1e-6 ); 00347 TS_ASSERT_DELTA( d[ 2 ], 1.0, 1e-6 ); 00348 00349 // -1 -1 -1 00350 t( 0, 0 ) = -1.0; 00351 t( 1, 1 ) = -1.0; 00352 t( 2, 2 ) = -1.0; 00353 00354 d = getEigenvaluesCardano( t ); 00355 00356 TS_ASSERT_DELTA( d[ 0 ], -1.0, 1e-6 ); 00357 TS_ASSERT_DELTA( d[ 1 ], -1.0, 1e-6 ); 00358 TS_ASSERT_DELTA( d[ 2 ], -1.0, 1e-6 ); 00359 00360 // 1 0 1 00361 t( 0, 0 ) = 1.0; 00362 t( 1, 1 ) = 0.0; 00363 t( 2, 2 ) = 1.0; 00364 00365 d = getEigenvaluesCardano( t ); 00366 00367 TS_ASSERT_DELTA( d[ 0 ], 1.0, 1e-6 ); 00368 TS_ASSERT_DELTA( d[ 1 ], 1.0, 1e-6 ); 00369 TS_ASSERT_DELTA( d[ 2 ], 0.0, 1e-6 ); 00370 00371 // 0 0 0 00372 t( 0, 0 ) = 0.0; 00373 t( 1, 1 ) = 0.0; 00374 t( 2, 2 ) = 0.0; 00375 00376 d = getEigenvaluesCardano( t ); 00377 00378 TS_ASSERT_DELTA( d[ 0 ], 0.0, 1e-6 ); 00379 TS_ASSERT_DELTA( d[ 1 ], 0.0, 1e-6 ); 00380 TS_ASSERT_DELTA( d[ 2 ], 0.0, 1e-6 ); 00381 00382 // similar eigenvalues 00383 // 2.000001 0.0 1.999998 00384 t( 0, 0 ) = 2.000001; 00385 t( 1, 1 ) = 0.0; 00386 t( 2, 2 ) = 1.999998; 00387 00388 d = getEigenvaluesCardano( t ); 00389 00390 TS_ASSERT_DELTA( d[ 0 ], 2.000001, 1e-6 ); 00391 TS_ASSERT_DELTA( d[ 1 ], 1.999998, 1e-6 ); 00392 TS_ASSERT_DELTA( d[ 2 ], 0.0, 1e-6 ); 00393 00394 // very large eigenvalues 00395 // 3.824572321236e10 1 2 00396 t( 0, 0 ) = 3.824572321236e10; 00397 t( 1, 1 ) = 1.0; 00398 t( 2, 2 ) = 2.0; 00399 00400 d = getEigenvaluesCardano( t ); 00401 00402 TS_ASSERT_DELTA( d[ 0 ], 3.824572321236e10, 1e-6 ); 00403 TS_ASSERT_DELTA( d[ 1 ], 2.0, 1e-6 ); 00404 TS_ASSERT_DELTA( d[ 2 ], 1.0, 1e-6 ); 00405 00406 // very small eigenvalues 00407 // 3.824572321236e-10 1 2 00408 t( 0, 0 ) = 3.824572321236e-30; 00409 t( 1, 1 ) = 1.0; 00410 t( 2, 2 ) = 2.0; 00411 00412 d = getEigenvaluesCardano( t ); 00413 00414 TS_ASSERT_DELTA( d[ 0 ], 2.0, 1e-6 ); 00415 TS_ASSERT_DELTA( d[ 1 ], 1.0, 1e-6 ); 00416 TS_ASSERT_DELTA( d[ 2 ], 3.824572321236e-30, 1e-6 ); 00417 00418 // some more sophisticated tests 00419 // (using similarity transformations on diagonal matrices to create test cases) 00420 t( 0, 0 ) = 1; 00421 t( 1, 1 ) = 2; 00422 t( 2, 2 ) = 3; 00423 00424 // note that the jacobi eigenvector functions does not sort the eigenvalues and 00425 // eigenvectors that were found 00426 t = similarity_rotate_givens( t, 0, 2, 2.78 ); 00427 00428 d = getEigenvaluesCardano( t ); 00429 00430 TS_ASSERT_DELTA( d[ 0 ], 3.0, 1e-6 ); 00431 TS_ASSERT_DELTA( d[ 1 ], 2.0, 1e-6 ); 00432 TS_ASSERT_DELTA( d[ 2 ], 1.0, 1e-6 ); 00433 00434 // the eigenvalues calculated for this test are a bit off 00435 //t = WTensorSym< 2, 3 >(); 00436 //t( 0, 0 ) = 2; 00437 //t( 1, 1 ) = 1; 00438 //t( 2, 2 ) = 3; 00439 00440 //t = similarity_rotate_givens( t, 0, 2, 2.79 ); 00441 //t = similarity_rotate_givens( t, 1, 2, -3.44 ); 00442 //t = similarity_rotate_givens( t, 1, 0, -0.46 ); 00443 //t = similarity_rotate_givens( t, 2, 1, 5.98 ); 00444 00445 //d = getEigenvaluesCardano( t ); 00446 00447 //TS_ASSERT_DELTA( d[ 0 ], 3.0, 1e-6 ); 00448 //TS_ASSERT_DELTA( d[ 1 ], 2.0, 1e-6 ); 00449 //TS_ASSERT_DELTA( d[ 2 ], 1.0, 1e-6 ); 00450 00451 t = WTensorSym< 2, 3 >(); 00452 t( 0, 0 ) = 2; 00453 t( 1, 1 ) = 2; 00454 t( 2, 2 ) = 2; 00455 00456 t = similarity_rotate_givens( t, 1, 2, -3.44 ); 00457 t = similarity_rotate_givens( t, 2, 1, 5.98 ); 00458 t = similarity_rotate_givens( t, 0, 2, 2.79 ); 00459 t = similarity_rotate_givens( t, 1, 0, -0.46 ); 00460 00461 d = getEigenvaluesCardano( t ); 00462 00463 TS_ASSERT_DELTA( d[ 0 ], 2.0, 1e-6 ); 00464 TS_ASSERT_DELTA( d[ 1 ], 2.0, 1e-6 ); 00465 TS_ASSERT_DELTA( d[ 2 ], 2.0, 1e-6 ); 00466 } 00467 00468 /** 00469 * Test if tensor log and exp functions behave correctly. 00470 */ 00471 void testLogAndExp() 00472 { 00473 // init some tensor 00474 WTensorSym< 2, 3, double > tens; 00475 tens( 0, 0 ) = 8.0; 00476 tens( 0, 1 ) = 10.0; 00477 tens( 0, 2 ) = 4.0; 00478 tens( 1, 1 ) = 17.0; 00479 tens( 1, 2 ) = 14.0; 00480 tens( 2, 2 ) = 20.0; 00481 00482 WTensorSym< 2, 3, double > s = tensorExp( tensorLog( tens ) ); 00483 00484 TS_ASSERT_DELTA( s( 0, 0 ), tens( 0, 0 ), 1e-4 ); 00485 TS_ASSERT_DELTA( s( 0, 1 ), tens( 0, 1 ), 1e-4 ); 00486 TS_ASSERT_DELTA( s( 0, 2 ), tens( 0, 2 ), 1e-4 ); 00487 TS_ASSERT_DELTA( s( 1, 1 ), tens( 1, 1 ), 1e-4 ); 00488 TS_ASSERT_DELTA( s( 1, 2 ), tens( 1, 2 ), 1e-4 ); 00489 TS_ASSERT_DELTA( s( 2, 2 ), tens( 2, 2 ), 1e-4 ); 00490 } 00491 00492 private: 00493 /** 00494 * A helper function performing a similarity transform using a givens rotation. 00495 * 00496 * \param m The symmetric tensor to transform. 00497 * \param i A row index. 00498 * \param j A column index. 00499 * \param angle The rotation angle (in radians). 00500 * 00501 * \note i must not have the same values as j 00502 * 00503 * \return The new tensor 00504 */ 00505 template< std::size_t dim, typename Data_T > 00506 WTensorSym< 2, dim, Data_T > similarity_rotate_givens( WTensorSym< 2, dim, Data_T > const& m, 00507 std::size_t i, 00508 std::size_t j, 00509 double angle ) 00510 { 00511 WAssert( i != j, "" ); 00512 00513 double s = sin( angle ); 00514 double c = cos( angle ); 00515 WTensor< 2, dim, Data_T > r; 00516 for( std::size_t k = 0; k < dim; ++k ) 00517 { 00518 if( k == i || k == j ) 00519 { 00520 r( k, k ) = c; 00521 } 00522 else 00523 { 00524 r( k, k ) = 1.0; 00525 } 00526 } 00527 r( i, j ) = s; 00528 r( j, i ) = -s; 00529 WTensor< 2, dim, Data_T > t( r ); 00530 std::swap( t( i, j ), t( j, i ) ); 00531 00532 r = t * m * r; 00533 WTensorSym< 2, dim, Data_T > res; 00534 res( 0, 0 ) = r( 0, 0 ); 00535 res( 1, 0 ) = r( 1, 0 ); 00536 res( 2, 0 ) = r( 2, 0 ); 00537 res( 1, 1 ) = r( 1, 1 ); 00538 res( 2, 1 ) = r( 2, 1 ); 00539 res( 2, 2 ) = r( 2, 2 ); 00540 00541 return res; 00542 } 00543 00544 /** 00545 * Test if the given vectors are eigenvectors to the given eigenvalues of a 00546 * symmetric matrix. 00547 * 00548 * \param m A symmetric matrix. 00549 * \param sys The eigen system ( eigenvalues and eigenvectors ) 00550 */ 00551 template< std::size_t dim, typename Data_T > 00552 void compare_results( WTensorSym< 2, dim, Data_T > const& m, RealEigenSystem const& sys ) 00553 { 00554 WTensor< 2, dim, Data_T > t, r; 00555 for( std::size_t i = 0; i < dim; ++i ) 00556 { 00557 for( std::size_t j = 0; j < dim; ++j ) 00558 { 00559 t( j, i ) = sys[i].second[j]; 00560 r( j, i ) = sys[i].first * sys[i].second[j]; 00561 } 00562 } 00563 00564 t = m * t; 00565 for( std::size_t i = 0; i < dim; ++i ) 00566 { 00567 for( std::size_t j = 0; j < dim; ++j ) 00568 { 00569 TS_ASSERT_DELTA( t( i, j ), r( i, j ), 1e-15 ); 00570 } 00571 } 00572 } 00573 }; 00574 00575 /** 00576 * Test class for all tensor operators. 00577 */ 00578 class WTensorOperatorsTest : public CxxTest::TestSuite 00579 { 00580 public: 00581 /** 00582 * Test order 2 tensor multiplication. 00583 */ 00584 void testMultiplyTensorsOperator() 00585 { 00586 // some special cases, WTensor * WTensor 00587 TS_ASSERT_EQUALS( zero * zero, zero ); 00588 TS_ASSERT_EQUALS( zero * one, zero ); 00589 TS_ASSERT_EQUALS( one * zero, zero ); 00590 TS_ASSERT_EQUALS( one * one, one ); 00591 TS_ASSERT_EQUALS( one * rdm1, rdm1 ); 00592 TS_ASSERT_EQUALS( one * rdm2, rdm2 ); 00593 TS_ASSERT_EQUALS( rdm1 * one, rdm1 ); 00594 TS_ASSERT_EQUALS( rdm2 * one, rdm2 ); 00595 00596 // a normal case 00597 TS_ASSERT_EQUALS( rdm1 * rdm2, res1 ); 00598 00599 // some special cases, WTensorSym * WTensorSym 00600 TS_ASSERT_EQUALS( szero * szero, zero ); 00601 TS_ASSERT_EQUALS( szero * sone, zero ); 00602 TS_ASSERT_EQUALS( sone * szero, zero ); 00603 TS_ASSERT_EQUALS( sone * sone, one ); 00604 TS_ASSERT_EQUALS( sone * srdm1, res2 ); 00605 TS_ASSERT_EQUALS( srdm1 * sone, res2 ); 00606 TS_ASSERT_EQUALS( srdm1 * srdm2, res3 ); 00607 00608 // WTensor * WTensorSym 00609 TS_ASSERT_EQUALS( zero * szero, zero ); 00610 TS_ASSERT_EQUALS( one * sone, one ); 00611 TS_ASSERT_EQUALS( zero * sone, zero ); 00612 TS_ASSERT_EQUALS( one * szero, zero ); 00613 00614 // WTensorSym * WTensor 00615 TS_ASSERT_EQUALS( szero * zero, zero ); 00616 TS_ASSERT_EQUALS( sone * one, one ); 00617 TS_ASSERT_EQUALS( szero * one, zero ); 00618 TS_ASSERT_EQUALS( sone * zero, zero ); 00619 TS_ASSERT_EQUALS( srdm1 * rdm1, res4 ); 00620 } 00621 00622 /** 00623 * The optimizations for symmetric tensors should not corrupt the result. 00624 */ 00625 void testEvaluateSphericalFunction() 00626 { 00627 WTensorSym< 4, 3, double > t; 00628 // the tensor 00629 t( 0, 0, 0, 0 ) = 2.5476; 00630 t( 1, 1, 1, 1 ) = 3.5476; 00631 t( 2, 2, 2, 2 ) = 4.5476; 00632 t( 0, 0, 0, 1 ) = 5.5476; 00633 t( 0, 0, 0, 2 ) = 6.5476; 00634 t( 1, 1, 1, 0 ) = 7.5476; 00635 t( 1, 1, 1, 2 ) = 8.5476; 00636 t( 2, 2, 2, 0 ) = 9.5476; 00637 t( 2, 2, 2, 1 ) = 10.5476; 00638 t( 0, 0, 1, 2 ) = 11.5476; 00639 t( 1, 1, 0, 2 ) = 12.5476; 00640 t( 2, 2, 0, 1 ) = 13.5476; 00641 t( 0, 0, 1, 1 ) = 14.5476; 00642 t( 0, 0, 2, 2 ) = 15.5476; 00643 t( 1, 1, 2, 2 ) = 16.5476; 00644 00645 // the gradients 00646 std::vector< WVector3d > gradients; 00647 gradients.push_back( WVector3d( 1.0, 0.0, 0.0 ) ); 00648 gradients.push_back( WVector3d( 0.0, 1.0, 0.0 ) ); 00649 gradients.push_back( normalize( WVector3d( 1.0, 1.0, 0.0 ) ) ); 00650 gradients.push_back( normalize( WVector3d( 0.3, 0.4, 0.5 ) ) ); 00651 gradients.push_back( normalize( WVector3d( -7.0, 3.0, -1.0 ) ) ); 00652 00653 for( int k = 0; k < 5; ++k ) 00654 { 00655 double res = calcTens( t, gradients[ k ] ); 00656 TS_ASSERT_DELTA( res, evaluateSphericalFunction( t, gradients[ k ] ), 0.001 ); 00657 } 00658 } 00659 00660 private: 00661 /** 00662 * Initialize a lot of tensors. 00663 */ 00664 void setUp() 00665 { 00666 one( 0, 0 ) = one( 1, 1 ) = one( 2, 2 ) = 1.0; 00667 sone( 0, 0 ) = sone( 1, 1 ) = sone( 2, 2 ) = 1.0; 00668 00669 rdm1( 0, 0 ) = 2; 00670 rdm1( 0, 1 ) = 3; 00671 rdm1( 0, 2 ) = 1; 00672 rdm1( 1, 0 ) = 0; 00673 rdm1( 1, 1 ) = 4; 00674 rdm1( 1, 2 ) = 0; 00675 rdm1( 2, 0 ) = 5; 00676 rdm1( 2, 1 ) = -3; 00677 rdm1( 2, 2 ) = -7; 00678 00679 rdm2( 0, 0 ) = -1; 00680 rdm2( 0, 1 ) = -1; 00681 rdm2( 0, 2 ) = -1; 00682 rdm2( 1, 0 ) = 3; 00683 rdm2( 1, 1 ) = 0; 00684 rdm2( 1, 2 ) = -2; 00685 rdm2( 2, 0 ) = -1; 00686 rdm2( 2, 1 ) = -2; 00687 rdm2( 2, 2 ) = -3; 00688 00689 res1( 0, 0 ) = 6; 00690 res1( 0, 1 ) = -4; 00691 res1( 0, 2 ) = -11; 00692 res1( 1, 0 ) = 12; 00693 res1( 1, 1 ) = 0; 00694 res1( 1, 2 ) = -8; 00695 res1( 2, 0 ) = -7; 00696 res1( 2, 1 ) = 9; 00697 res1( 2, 2 ) = 22; 00698 00699 srdm1( 0, 0 ) = 2; 00700 srdm1( 0, 1 ) = 3; 00701 srdm1( 0, 2 ) = 1; 00702 srdm1( 1, 1 ) = -1; 00703 srdm1( 1, 2 ) = 4; 00704 srdm1( 2, 2 ) = -3; 00705 00706 srdm2( 0, 0 ) = -2; 00707 srdm2( 0, 1 ) = 5; 00708 srdm2( 0, 2 ) = 2; 00709 srdm2( 1, 1 ) = 0; 00710 srdm2( 1, 2 ) = -3; 00711 srdm2( 2, 2 ) = 1; 00712 00713 // copy symmetric tensor srdm1 to an asymmetric tensor 00714 res2 = srdm1; 00715 00716 res3( 0, 0 ) = 13; 00717 res3( 0, 1 ) = 7; 00718 res3( 0, 2 ) = -4; 00719 res3( 1, 0 ) = -3; 00720 res3( 1, 1 ) = 3; 00721 res3( 1, 2 ) = 13; 00722 res3( 2, 0 ) = 12; 00723 res3( 2, 1 ) = 14; 00724 res3( 2, 2 ) = -13; 00725 00726 res4( 0, 0 ) = 9; 00727 res4( 0, 1 ) = 15; 00728 res4( 0, 2 ) = -5; 00729 res4( 1, 0 ) = 26; 00730 res4( 1, 1 ) = -7; 00731 res4( 1, 2 ) = -25; 00732 res4( 2, 0 ) = -13; 00733 res4( 2, 1 ) = 28; 00734 res4( 2, 2 ) = 22; 00735 } 00736 00737 /** 00738 * A helper function that implements the simple approach to tensor evaluation. 00739 * 00740 * \param t The tensor. 00741 * \param v The gradient. 00742 * 00743 * \return value 00744 */ 00745 double calcTens( WTensorSym< 4, 3, double > const& t, WVector3d const& v ) 00746 { 00747 double res = 0.0; 00748 for( int a = 0; a < 3; ++a ) 00749 { 00750 for( int b = 0; b < 3; ++b ) 00751 { 00752 for( int c = 0; c < 3; ++c ) 00753 { 00754 for( int d = 0; d < 3; ++d ) 00755 { 00756 res += v[ a ] * v[ b ] * v[ c ] * v[ d ] * t( a, b, c, d ); 00757 } 00758 } 00759 } 00760 } 00761 return res; 00762 } 00763 00764 //! a test tensor 00765 WTensor< 2, 3, int > one; 00766 //! a test tensor 00767 WTensor< 2, 3, int > zero; 00768 //! a test tensor 00769 WTensor< 2, 3, int > rdm1; 00770 //! a test tensor 00771 WTensor< 2, 3, int > rdm2; 00772 //! a test tensor 00773 WTensor< 2, 3, int > res1; 00774 //! a test tensor 00775 WTensor< 2, 3, int > res2; 00776 //! a test tensor 00777 WTensor< 2, 3, int > res3; 00778 //! a test tensor 00779 WTensor< 2, 3, int > res4; 00780 //! a test tensor 00781 WTensorSym< 2, 3, int > sone; 00782 //! a test tensor 00783 WTensorSym< 2, 3, int > szero; 00784 //! a test tensor 00785 WTensorSym< 2, 3, int > srdm1; 00786 //! a test tensor 00787 WTensorSym< 2, 3, int > srdm2; 00788 }; 00789 00790 #endif // WTENSORFUNCTIONS_TEST_H