OpenWalnut  1.4.0
WTensorFunctions_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 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