OpenWalnut 1.2.5

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 private:
00469     /**
00470      * A helper function performing a similarity transform using a givens rotation.
00471      *
00472      * \param m The symmetric tensor to transform.
00473      * \param i A row index.
00474      * \param j A column index.
00475      * \param angle The rotation angle (in radians).
00476      *
00477      * \note i must not have the same values as j
00478      *
00479      * \return The new tensor
00480      */
00481     template< std::size_t dim, typename Data_T >
00482     WTensorSym< 2, dim, Data_T > similarity_rotate_givens( WTensorSym< 2, dim, Data_T > const& m,
00483                                                                   std::size_t i,
00484                                                                   std::size_t j,
00485                                                                   double angle )
00486     {
00487         WAssert( i != j, "" );
00488 
00489         double s = sin( angle );
00490         double c = cos( angle );
00491         WTensor< 2, dim, Data_T > r;
00492         for( std::size_t k = 0; k < dim; ++k )
00493         {
00494             if( k == i || k == j )
00495             {
00496                 r( k, k ) = c;
00497             }
00498             else
00499             {
00500                 r( k, k ) = 1.0;
00501             }
00502         }
00503         r( i, j ) = s;
00504         r( j, i ) = -s;
00505         WTensor< 2, dim, Data_T > t( r );
00506         std::swap( t( i, j ), t( j, i ) );
00507 
00508         r = t * m * r;
00509         WTensorSym< 2, dim, Data_T > res;
00510         res( 0, 0 ) = r( 0, 0 );
00511         res( 1, 0 ) = r( 1, 0 );
00512         res( 2, 0 ) = r( 2, 0 );
00513         res( 1, 1 ) = r( 1, 1 );
00514         res( 2, 1 ) = r( 2, 1 );
00515         res( 2, 2 ) = r( 2, 2 );
00516 
00517         return res;
00518     }
00519 
00520     /**
00521      * Test if the given vectors are eigenvectors to the given eigenvalues of a
00522      * symmetric matrix.
00523      *
00524      * \param m A symmetric matrix.
00525      * \param sys The eigen system ( eigenvalues and eigenvectors )
00526      */
00527     template< std::size_t dim, typename Data_T >
00528     void compare_results( WTensorSym< 2, dim, Data_T > const& m, RealEigenSystem const& sys )
00529     {
00530         WTensor< 2, dim, Data_T > t, r;
00531         for( std::size_t i = 0; i < dim; ++i )
00532         {
00533             for( std::size_t j = 0; j < dim; ++j )
00534             {
00535                 t( j, i ) = sys[i].second[j];
00536                 r( j, i ) = sys[i].first * sys[i].second[j];
00537             }
00538         }
00539 
00540         t = m * t;
00541         for( std::size_t i = 0; i < dim; ++i )
00542         {
00543             for( std::size_t j = 0; j < dim; ++j )
00544             {
00545                 TS_ASSERT_DELTA( t( i, j ), r( i, j ), 1e-15 );
00546             }
00547         }
00548     }
00549 };
00550 
00551 /**
00552  * Test class for all tensor operators.
00553  */
00554 class WTensorOperatorsTest : public CxxTest::TestSuite
00555 {
00556 public:
00557     /**
00558      * Test order 2 tensor multiplication.
00559      */
00560     void testMultiplyTensorsOperator()
00561     {
00562         // some special cases, WTensor * WTensor
00563         TS_ASSERT_EQUALS( zero * zero, zero );
00564         TS_ASSERT_EQUALS( zero * one, zero );
00565         TS_ASSERT_EQUALS( one * zero, zero );
00566         TS_ASSERT_EQUALS( one * one, one );
00567         TS_ASSERT_EQUALS( one * rdm1, rdm1 );
00568         TS_ASSERT_EQUALS( one * rdm2, rdm2 );
00569         TS_ASSERT_EQUALS( rdm1 * one, rdm1 );
00570         TS_ASSERT_EQUALS( rdm2 * one, rdm2 );
00571 
00572         // a normal case
00573         TS_ASSERT_EQUALS( rdm1 * rdm2, res1 );
00574 
00575         // some special cases, WTensorSym * WTensorSym
00576         TS_ASSERT_EQUALS( szero * szero, zero );
00577         TS_ASSERT_EQUALS( szero * sone, zero );
00578         TS_ASSERT_EQUALS( sone * szero, zero );
00579         TS_ASSERT_EQUALS( sone * sone, one );
00580         TS_ASSERT_EQUALS( sone * srdm1, res2 );
00581         TS_ASSERT_EQUALS( srdm1 * sone, res2 );
00582         TS_ASSERT_EQUALS( srdm1 * srdm2, res3 );
00583 
00584         // WTensor * WTensorSym
00585         TS_ASSERT_EQUALS( zero * szero, zero );
00586         TS_ASSERT_EQUALS( one * sone, one );
00587         TS_ASSERT_EQUALS( zero * sone, zero );
00588         TS_ASSERT_EQUALS( one * szero, zero );
00589 
00590         // WTensorSym * WTensor
00591         TS_ASSERT_EQUALS( szero * zero, zero );
00592         TS_ASSERT_EQUALS( sone * one, one );
00593         TS_ASSERT_EQUALS( szero * one, zero );
00594         TS_ASSERT_EQUALS( sone * zero, zero );
00595         TS_ASSERT_EQUALS( srdm1 * rdm1, res4 );
00596     }
00597 
00598     /**
00599      * The optimizations for symmetric tensors should not corrupt the result.
00600      */
00601     void testEvaluateSphericalFunction()
00602     {
00603         WTensorSym< 4, 3, double > t;
00604         // the tensor
00605         t( 0, 0, 0, 0 ) = 2.5476;
00606         t( 1, 1, 1, 1 ) = 3.5476;
00607         t( 2, 2, 2, 2 ) = 4.5476;
00608         t( 0, 0, 0, 1 ) = 5.5476;
00609         t( 0, 0, 0, 2 ) = 6.5476;
00610         t( 1, 1, 1, 0 ) = 7.5476;
00611         t( 1, 1, 1, 2 ) = 8.5476;
00612         t( 2, 2, 2, 0 ) = 9.5476;
00613         t( 2, 2, 2, 1 ) = 10.5476;
00614         t( 0, 0, 1, 2 ) = 11.5476;
00615         t( 1, 1, 0, 2 ) = 12.5476;
00616         t( 2, 2, 0, 1 ) = 13.5476;
00617         t( 0, 0, 1, 1 ) = 14.5476;
00618         t( 0, 0, 2, 2 ) = 15.5476;
00619         t( 1, 1, 2, 2 ) = 16.5476;
00620 
00621         // the gradients
00622         std::vector< WVector3d > gradients;
00623         gradients.push_back( WVector3d( 1.0, 0.0, 0.0 ) );
00624         gradients.push_back( WVector3d( 0.0, 1.0, 0.0 ) );
00625         gradients.push_back( normalize( WVector3d( 1.0, 1.0, 0.0 ) ) );
00626         gradients.push_back( normalize( WVector3d( 0.3, 0.4, 0.5 ) ) );
00627         gradients.push_back( normalize( WVector3d( -7.0, 3.0, -1.0 ) ) );
00628 
00629         for( int k = 0; k < 5; ++k )
00630         {
00631             double res = calcTens( t, gradients[ k ] );
00632             TS_ASSERT_DELTA( res, evaluateSphericalFunction( t, gradients[ k ] ), 0.001 );
00633         }
00634     }
00635 
00636 private:
00637     /**
00638      * Initialize a lot of tensors.
00639      */
00640     void setUp()
00641     {
00642         one( 0, 0 ) = one( 1, 1 ) = one( 2, 2 ) = 1.0;
00643         sone( 0, 0 ) = sone( 1, 1 ) = sone( 2, 2 ) = 1.0;
00644 
00645         rdm1( 0, 0 ) = 2;
00646         rdm1( 0, 1 ) = 3;
00647         rdm1( 0, 2 ) = 1;
00648         rdm1( 1, 0 ) = 0;
00649         rdm1( 1, 1 ) = 4;
00650         rdm1( 1, 2 ) = 0;
00651         rdm1( 2, 0 ) = 5;
00652         rdm1( 2, 1 ) = -3;
00653         rdm1( 2, 2 ) = -7;
00654 
00655         rdm2( 0, 0 ) = -1;
00656         rdm2( 0, 1 ) = -1;
00657         rdm2( 0, 2 ) = -1;
00658         rdm2( 1, 0 ) = 3;
00659         rdm2( 1, 1 ) = 0;
00660         rdm2( 1, 2 ) = -2;
00661         rdm2( 2, 0 ) = -1;
00662         rdm2( 2, 1 ) = -2;
00663         rdm2( 2, 2 ) = -3;
00664 
00665         res1( 0, 0 ) = 6;
00666         res1( 0, 1 ) = -4;
00667         res1( 0, 2 ) = -11;
00668         res1( 1, 0 ) = 12;
00669         res1( 1, 1 ) = 0;
00670         res1( 1, 2 ) = -8;
00671         res1( 2, 0 ) = -7;
00672         res1( 2, 1 ) = 9;
00673         res1( 2, 2 ) = 22;
00674 
00675         srdm1( 0, 0 ) = 2;
00676         srdm1( 0, 1 ) = 3;
00677         srdm1( 0, 2 ) = 1;
00678         srdm1( 1, 1 ) = -1;
00679         srdm1( 1, 2 ) = 4;
00680         srdm1( 2, 2 ) = -3;
00681 
00682         srdm2( 0, 0 ) = -2;
00683         srdm2( 0, 1 ) = 5;
00684         srdm2( 0, 2 ) = 2;
00685         srdm2( 1, 1 ) = 0;
00686         srdm2( 1, 2 ) = -3;
00687         srdm2( 2, 2 ) = 1;
00688 
00689         // copy symmetric tensor srdm1 to an asymmetric tensor
00690         res2 = srdm1;
00691 
00692         res3( 0, 0 ) = 13;
00693         res3( 0, 1 ) = 7;
00694         res3( 0, 2 ) = -4;
00695         res3( 1, 0 ) = -3;
00696         res3( 1, 1 ) = 3;
00697         res3( 1, 2 ) = 13;
00698         res3( 2, 0 ) = 12;
00699         res3( 2, 1 ) = 14;
00700         res3( 2, 2 ) = -13;
00701 
00702         res4( 0, 0 ) = 9;
00703         res4( 0, 1 ) = 15;
00704         res4( 0, 2 ) = -5;
00705         res4( 1, 0 ) = 26;
00706         res4( 1, 1 ) = -7;
00707         res4( 1, 2 ) = -25;
00708         res4( 2, 0 ) = -13;
00709         res4( 2, 1 ) = 28;
00710         res4( 2, 2 ) = 22;
00711     }
00712 
00713     /**
00714      * A helper function that implements the simple approach to tensor evaluation.
00715      *
00716      * \param t The tensor.
00717      * \param v The gradient.
00718      *
00719      * \return value
00720      */
00721     double calcTens( WTensorSym< 4, 3, double > const& t, WVector3d const& v )
00722     {
00723         double res = 0.0;
00724         for( int a = 0; a < 3; ++a )
00725         {
00726             for( int b = 0; b < 3; ++b )
00727             {
00728                 for( int c = 0; c < 3; ++c )
00729                 {
00730                     for( int d = 0; d < 3; ++d )
00731                     {
00732                         res += v[ a ] * v[ b ] * v[ c ] * v[ d ] * t( a, b, c, d );
00733                     }
00734                 }
00735             }
00736         }
00737         return res;
00738     }
00739 
00740     //! a test tensor
00741     WTensor< 2, 3, int > one;
00742     //! a test tensor
00743     WTensor< 2, 3, int > zero;
00744     //! a test tensor
00745     WTensor< 2, 3, int > rdm1;
00746     //! a test tensor
00747     WTensor< 2, 3, int > rdm2;
00748     //! a test tensor
00749     WTensor< 2, 3, int > res1;
00750     //! a test tensor
00751     WTensor< 2, 3, int > res2;
00752     //! a test tensor
00753     WTensor< 2, 3, int > res3;
00754     //! a test tensor
00755     WTensor< 2, 3, int > res4;
00756     //! a test tensor
00757     WTensorSym< 2, 3, int > sone;
00758     //! a test tensor
00759     WTensorSym< 2, 3, int > szero;
00760     //! a test tensor
00761     WTensorSym< 2, 3, int > srdm1;
00762     //! a test tensor
00763     WTensorSym< 2, 3, int > srdm2;
00764 };
00765 
00766 #endif  // WTENSORFUNCTIONS_TEST_H
 All Classes Namespaces Functions Variables Typedefs Enumerations Enumerator Friends