MOAB: Mesh Oriented datABase  (version 5.4.0)
TMetricTest.cpp
Go to the documentation of this file.
00001 /* *****************************************************************
00002     MESQUITE -- The Mesh Quality Improvement Toolkit
00003 
00004     Copyright 2006 Sandia National Laboratories.  Developed at the
00005     University of Wisconsin--Madison under SNL contract number
00006     624796.  The U.S. Government and the University of Wisconsin
00007     retain certain rights to this software.
00008 
00009     This library is free software; you can redistribute it and/or
00010     modify it under the terms of the GNU Lesser General Public
00011     License as published by the Free Software Foundation; either
00012     version 2.1 of the License, or (at your option) any later version.
00013 
00014     This library is distributed in the hope that it will be useful,
00015     but WITHOUT ANY WARRANTY; without even the implied warranty of
00016     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00017     Lesser General Public License for more details.
00018 
00019     You should have received a copy of the GNU Lesser General Public License
00020     (lgpl.txt) along with this library; if not, write to the Free Software
00021     Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
00022 
00023     (2010) kraftche@cae.wisc.edu
00024 
00025   ***************************************************************** */
00026 
00027 /** \file TMetricTest.cpp
00028  *  \brief Unit tests TMetric base class functionality
00029  *  \author Jason Kraftcheck
00030  */
00031 
00032 #include "Mesquite.hpp"
00033 #include "TMetric.hpp"
00034 #include "UnitUtil.hpp"
00035 #include "MsqError.hpp"
00036 
00037 using namespace MBMesquite;
00038 
00039 class TMetricTest : public CppUnit::TestFixture
00040 {
00041     CPPUNIT_TEST_SUITE( TMetricTest );
00042     CPPUNIT_TEST( test_numerical_gradient_2D );
00043     CPPUNIT_TEST( test_numerical_hessian_2D );
00044     CPPUNIT_TEST( test_numerical_gradient_3D );
00045     CPPUNIT_TEST( test_numerical_hessian_3D );
00046     CPPUNIT_TEST_SUITE_END();
00047 
00048   public:
00049     void test_numerical_gradient_2D();
00050     void test_numerical_hessian_2D();
00051     void test_numerical_gradient_3D();
00052     void test_numerical_hessian_3D();
00053 };
00054 
00055 CPPUNIT_TEST_SUITE_NAMED_REGISTRATION( TMetricTest, "Unit" );
00056 CPPUNIT_TEST_SUITE_NAMED_REGISTRATION( TMetricTest, "TMetricTest" );
00057 
00058 // implement metric such that the expected derivatives
00059 // at each location r,c in dm/dT are 3r+c+1
00060 class GradTestMetricRel : public TMetric
00061 {
00062   public:
00063     std::string get_name() const
00064     {
00065         return "GradTest";
00066     }
00067 
00068     static double grad( int r, int c )
00069     {
00070         return 3 * r + c + 1;
00071     }
00072 
00073     bool evaluate( const MsqMatrix< 2, 2 >& T, double& result, MsqError& )
00074     {
00075         result = 0;
00076         for( int r = 0; r < 2; ++r )
00077             for( int c = 0; c < 2; ++c )
00078                 result += grad( r, c ) * T( r, c );
00079         return true;
00080     }
00081 
00082     bool evaluate( const MsqMatrix< 3, 3 >& T, double& result, MsqError& )
00083     {
00084         result = 0;
00085         for( int r = 0; r < 3; ++r )
00086             for( int c = 0; c < 3; ++c )
00087                 result += grad( r, c ) * T( r, c );
00088         return true;
00089     }
00090 };
00091 
00092 // implement metric that is the |2T - T^t|^2
00093 // such that the Hessian is the constant:
00094 //  _          _
00095 // | 2  0  0  0 |
00096 // | 0 10 -8  0 |
00097 // | 0 -8 10  0 |
00098 // |_0  0  0  2_|
00099 //  _                         _
00100 // | 2  0  0  0  0  0  0  0  0 |
00101 // | 0 10  0 -8  0  0  0  0  0 |
00102 // | 0  0 10  0  0  0 -8  0  0 |
00103 // | 0 -8  0 10  0  0  0  0  0 |
00104 // | 0  0  0  0  2  0  0  0  0 |
00105 // | 0  0  0  0  0 10  0 -8  0 |
00106 // | 0  0 -8  0  0  0 10  0  0 |
00107 // | 0  0  0  0  0 -8  0 10  0 |
00108 // |_0  0  0  0  0  0  0  0  2_|
00109 class HessTestMetricRel : public TMetric
00110 {
00111   public:
00112     std::string get_name() const
00113     {
00114         return "HessTest";
00115     }
00116 
00117     bool evaluate( const MsqMatrix< 2, 2 >& T, double& result, MsqError& )
00118     {
00119         result = sqr_Frobenius( 2 * T - transpose( T ) );
00120         return true;
00121     }
00122     bool evaluate( const MsqMatrix< 3, 3 >& T, double& result, MsqError& )
00123     {
00124         result = sqr_Frobenius( 2 * T - transpose( T ) );
00125         return true;
00126     }
00127     bool evaluate_with_grad( const MsqMatrix< 2, 2 >& T, double& result, MsqMatrix< 2, 2 >& wrt_T, MsqError& )
00128     {
00129         result = sqr_Frobenius( 2 * T - transpose( T ) );
00130         wrt_T  = 10 * T - 8 * transpose( T );
00131         return true;
00132     }
00133     bool evaluate_with_grad( const MsqMatrix< 3, 3 >& T, double& result, MsqMatrix< 3, 3 >& wrt_T, MsqError& )
00134     {
00135         result = sqr_Frobenius( 2 * T - transpose( T ) );
00136         wrt_T  = 10 * T - 8 * transpose( T );
00137         return true;
00138     }
00139 };
00140 
00141 /** Simple target metric for testing second partial derivatives.
00142  *  \f$\mu(T) = |T|\f$
00143  *  \f$\frac{\partial\mu}{\partial T} = \frac{1}{|T|} T \f$
00144  *  \f$\frac{\partial^{2}\mu}{\partial t_{i,i}^2} = \frac{1}{|T|} - \frac{t_{i,i}^2}{|T|^3}\f$
00145  *  \f$\frac{\partial^{2}\mu}{\partial t_{i,j} \partial t_{k,l} (i \ne k or j \ne l)} =
00146  * -\frac{t_{i,j} a_{k,l}}{|T|^3}\f$
00147  */
00148 class HessTestMetricRel_2 : public TMetric
00149 {
00150   public:
00151     std::string get_name() const
00152     {
00153         return "HessTest2";
00154     }
00155 
00156     bool evaluate( const MsqMatrix< 2, 2 >& T, double& result, MsqError& err )
00157     {
00158         result = Frobenius( T );
00159         return true;
00160     }
00161 
00162     bool evaluate_with_grad( const MsqMatrix< 2, 2 >& T, double& result, MsqMatrix< 2, 2 >& d, MsqError& err )
00163     {
00164         result = Frobenius( T );
00165         d      = T / result;
00166         return true;
00167     }
00168 
00169     bool evaluate_with_hess( const MsqMatrix< 2, 2 >& T,
00170                              double& result,
00171                              MsqMatrix< 2, 2 >& d,
00172                              MsqMatrix< 2, 2 > d2[3],
00173                              MsqError& err )
00174     {
00175         result = Frobenius( T );
00176         d      = T / result;
00177         int h  = 0;
00178         for( int r = 0; r < 2; ++r )
00179         {
00180             int i = h;
00181             for( int c = r; c < 2; ++c )
00182                 d2[h++] = transpose( T.row( r ) ) * T.row( c ) / -( result * result * result );
00183             d2[i] += MsqMatrix< 2, 2 >( 1.0 / result );
00184         }
00185         return true;
00186     }
00187 
00188     bool evaluate( const MsqMatrix< 3, 3 >& T, double& result, MsqError& err )
00189     {
00190         result = Frobenius( T );
00191         return true;
00192     }
00193 
00194     bool evaluate_with_grad( const MsqMatrix< 3, 3 >& T, double& result, MsqMatrix< 3, 3 >& d, MsqError& err )
00195     {
00196         result = Frobenius( T );
00197         d      = T / result;
00198         return true;
00199     }
00200 
00201     bool evaluate_with_hess( const MsqMatrix< 3, 3 >& T,
00202                              double& result,
00203                              MsqMatrix< 3, 3 >& d,
00204                              MsqMatrix< 3, 3 > d2[6],
00205                              MsqError& err )
00206     {
00207         result = Frobenius( T );
00208         d      = T / result;
00209         int h  = 0;
00210         for( int r = 0; r < 3; ++r )
00211         {
00212             int i = h;
00213             for( int c = r; c < 3; ++c )
00214                 d2[h++] = transpose( T.row( r ) ) * T.row( c ) / -( result * result * result );
00215             d2[i] += MsqMatrix< 3, 3 >( 1.0 / result );
00216         }
00217         return true;
00218     }
00219 };
00220 
00221 void TMetricTest::test_numerical_gradient_2D()
00222 {
00223     GradTestMetricRel metric;
00224     HessTestMetricRel_2 metric2;
00225     const double Avals[] = { 1, 2, 2, 5 };
00226     const double Bvals[] = { -0.1, -0.15, -0.25, -0.8 };
00227     const MsqMatrix< 2, 2 > A( Avals );
00228     const MsqMatrix< 2, 2 > B( Bvals );
00229 
00230     MsqError err;
00231     MsqMatrix< 2, 2 > d;
00232     bool valid;
00233     double val, gval;
00234 
00235     MsqMatrix< 2, 2 > expected;
00236     for( int r = 0; r < 2; ++r )
00237         for( int c = 0; c < 2; ++c )
00238             expected( r, c ) = metric.grad( r, c );
00239 
00240     valid = metric.evaluate( A, val, err );
00241     ASSERT_NO_ERROR( err );
00242     CPPUNIT_ASSERT( valid );
00243     valid = metric.evaluate_with_grad( A, gval, d, err );
00244     ASSERT_NO_ERROR( err );
00245     CPPUNIT_ASSERT( valid );
00246     CPPUNIT_ASSERT_DOUBLES_EQUAL( val, gval, 1e-6 );
00247     ASSERT_MATRICES_EQUAL( expected, d, 1e-6 );
00248 
00249     valid = metric.evaluate( B, val, err );
00250     ASSERT_NO_ERROR( err );
00251     CPPUNIT_ASSERT( valid );
00252     valid = metric.evaluate_with_grad( B, gval, d, err );
00253     ASSERT_NO_ERROR( err );
00254     CPPUNIT_ASSERT( valid );
00255     CPPUNIT_ASSERT_DOUBLES_EQUAL( val, gval, 1e-6 );
00256     ASSERT_MATRICES_EQUAL( expected, d, 1e-6 );
00257 
00258     valid = metric.evaluate( inverse( A ), val, err );
00259     ASSERT_NO_ERROR( err );
00260     CPPUNIT_ASSERT( valid );
00261     valid = metric.evaluate_with_grad( inverse( A ), gval, d, err );
00262     ASSERT_NO_ERROR( err );
00263     CPPUNIT_ASSERT( valid );
00264     ASSERT_MATRICES_EQUAL( expected, d, 1e-6 );
00265 
00266     valid = metric.evaluate( inverse( B ), val, err );
00267     ASSERT_NO_ERROR( err );
00268     CPPUNIT_ASSERT( valid );
00269     valid = metric.evaluate_with_grad( inverse( B ), gval, d, err );
00270     ASSERT_NO_ERROR( err );
00271     CPPUNIT_ASSERT( valid );
00272     CPPUNIT_ASSERT_DOUBLES_EQUAL( val, gval, 1e-6 );
00273     ASSERT_MATRICES_EQUAL( expected, d, 1e-6 );
00274 
00275     valid = metric.evaluate( A * inverse( B ), val, err );
00276     ASSERT_NO_ERROR( err );
00277     CPPUNIT_ASSERT( valid );
00278     valid = metric.evaluate_with_grad( A * inverse( B ), gval, d, err );
00279     ASSERT_NO_ERROR( err );
00280     CPPUNIT_ASSERT( valid );
00281     CPPUNIT_ASSERT_DOUBLES_EQUAL( val, gval, 1e-6 );
00282     ASSERT_MATRICES_EQUAL( expected, d, 1e-6 );
00283 
00284     valid = metric.evaluate( B * inverse( A ), val, err );
00285     ASSERT_NO_ERROR( err );
00286     CPPUNIT_ASSERT( valid );
00287     valid = metric.evaluate_with_grad( B * inverse( A ), gval, d, err );
00288     ASSERT_NO_ERROR( err );
00289     CPPUNIT_ASSERT( valid );
00290     CPPUNIT_ASSERT_DOUBLES_EQUAL( val, gval, 1e-6 );
00291     ASSERT_MATRICES_EQUAL( expected, d, 1e-6 );
00292 
00293     MsqMatrix< 2, 2 > da;
00294     valid = metric2.evaluate_with_grad( A, val, da, err );
00295     ASSERT_NO_ERROR( err );
00296     CPPUNIT_ASSERT( valid );
00297     valid = metric2.TMetric::evaluate_with_grad( A, gval, d, err );
00298     ASSERT_NO_ERROR( err );
00299     CPPUNIT_ASSERT( valid );
00300     CPPUNIT_ASSERT_DOUBLES_EQUAL( val, gval, 1e-6 );
00301     ASSERT_MATRICES_EQUAL( da, d, 1e-6 );
00302 
00303     valid = metric2.evaluate_with_grad( B, val, da, err );
00304     ASSERT_NO_ERROR( err );
00305     CPPUNIT_ASSERT( valid );
00306     valid = metric2.TMetric::evaluate_with_grad( B, gval, d, err );
00307     ASSERT_NO_ERROR( err );
00308     CPPUNIT_ASSERT( valid );
00309     CPPUNIT_ASSERT_DOUBLES_EQUAL( val, gval, 1e-6 );
00310     ASSERT_MATRICES_EQUAL( da, d, 1e-6 );
00311 }
00312 
00313 void TMetricTest::test_numerical_gradient_3D()
00314 {
00315     GradTestMetricRel metric;
00316     HessTestMetricRel_2 metric2;
00317     const double Avals[] = { 1, 2, 3, 4, 1, 4, 3, 2, 1 };
00318     const double Bvals[] = { 0.1, 0.15, 0.05, 0.2, -0.1, -0.15, -0.05, -0.2, 2 };
00319     const MsqMatrix< 3, 3 > A( Avals );
00320     const MsqMatrix< 3, 3 > B( Bvals );
00321 
00322     MsqError err;
00323     MsqMatrix< 3, 3 > d;
00324     bool valid;
00325     double val, gval;
00326 
00327     MsqMatrix< 3, 3 > expected;
00328     for( int r = 0; r < 3; ++r )
00329         for( int c = 0; c < 3; ++c )
00330             expected( r, c ) = metric.grad( r, c );
00331 
00332     valid = metric.evaluate( A, val, err );
00333     ASSERT_NO_ERROR( err );
00334     CPPUNIT_ASSERT( valid );
00335     valid = metric.evaluate_with_grad( A, gval, d, err );
00336     ASSERT_NO_ERROR( err );
00337     CPPUNIT_ASSERT( valid );
00338     CPPUNIT_ASSERT_DOUBLES_EQUAL( val, gval, 1e-6 );
00339     ASSERT_MATRICES_EQUAL( expected, d, 1e-6 );
00340 
00341     valid = metric.evaluate( B, val, err );
00342     ASSERT_NO_ERROR( err );
00343     CPPUNIT_ASSERT( valid );
00344     valid = metric.evaluate_with_grad( B, gval, d, err );
00345     ASSERT_NO_ERROR( err );
00346     CPPUNIT_ASSERT( valid );
00347     CPPUNIT_ASSERT_DOUBLES_EQUAL( val, gval, 1e-6 );
00348     ASSERT_MATRICES_EQUAL( expected, d, 1e-6 );
00349 
00350     valid = metric.evaluate( inverse( A ), val, err );
00351     ASSERT_NO_ERROR( err );
00352     CPPUNIT_ASSERT( valid );
00353     valid = metric.evaluate_with_grad( inverse( A ), gval, d, err );
00354     ASSERT_NO_ERROR( err );
00355     CPPUNIT_ASSERT( valid );
00356     ASSERT_MATRICES_EQUAL( expected, d, 1e-6 );
00357 
00358     valid = metric.evaluate( inverse( B ), val, err );
00359     ASSERT_NO_ERROR( err );
00360     CPPUNIT_ASSERT( valid );
00361     valid = metric.evaluate_with_grad( inverse( B ), gval, d, err );
00362     ASSERT_NO_ERROR( err );
00363     CPPUNIT_ASSERT( valid );
00364     CPPUNIT_ASSERT_DOUBLES_EQUAL( val, gval, 1e-6 );
00365     ASSERT_MATRICES_EQUAL( expected, d, 1e-6 );
00366 
00367     valid = metric.evaluate( A * inverse( B ), val, err );
00368     ASSERT_NO_ERROR( err );
00369     CPPUNIT_ASSERT( valid );
00370     valid = metric.evaluate_with_grad( A * inverse( B ), gval, d, err );
00371     ASSERT_NO_ERROR( err );
00372     CPPUNIT_ASSERT( valid );
00373     CPPUNIT_ASSERT_DOUBLES_EQUAL( val, gval, 1e-6 );
00374     ASSERT_MATRICES_EQUAL( expected, d, 1e-6 );
00375 
00376     valid = metric.evaluate( B * inverse( A ), val, err );
00377     ASSERT_NO_ERROR( err );
00378     CPPUNIT_ASSERT( valid );
00379     valid = metric.evaluate_with_grad( B * inverse( A ), gval, d, err );
00380     ASSERT_NO_ERROR( err );
00381     CPPUNIT_ASSERT( valid );
00382     CPPUNIT_ASSERT_DOUBLES_EQUAL( val, gval, 1e-6 );
00383     ASSERT_MATRICES_EQUAL( expected, d, 1e-6 );
00384 
00385     MsqMatrix< 3, 3 > da;
00386     valid = metric2.evaluate_with_grad( A, val, da, err );
00387     ASSERT_NO_ERROR( err );
00388     CPPUNIT_ASSERT( valid );
00389     valid = metric2.TMetric::evaluate_with_grad( A, gval, d, err );
00390     ASSERT_NO_ERROR( err );
00391     CPPUNIT_ASSERT( valid );
00392     CPPUNIT_ASSERT_DOUBLES_EQUAL( val, gval, 1e-6 );
00393     ASSERT_MATRICES_EQUAL( da, d, 1e-6 );
00394 
00395     valid = metric2.evaluate_with_grad( B, val, da, err );
00396     ASSERT_NO_ERROR( err );
00397     CPPUNIT_ASSERT( valid );
00398     valid = metric2.TMetric::evaluate_with_grad( B, gval, d, err );
00399     ASSERT_NO_ERROR( err );
00400     CPPUNIT_ASSERT( valid );
00401     CPPUNIT_ASSERT_DOUBLES_EQUAL( val, gval, 1e-6 );
00402     ASSERT_MATRICES_EQUAL( da, d, 1e-6 );
00403 }
00404 
00405 void TMetricTest::test_numerical_hessian_2D()
00406 {
00407     HessTestMetricRel metric;
00408     HessTestMetricRel_2 metric2;
00409     const double Avals[] = { 1, 2, 2, 5 };
00410     const double Bvals[] = { -0.1, -0.15, -0.25, -0.8 };
00411     const MsqMatrix< 2, 2 > A( Avals );
00412     const MsqMatrix< 2, 2 > B( Bvals );
00413 
00414     MsqError err;
00415     MsqMatrix< 2, 2 > g, gh;
00416     MsqMatrix< 2, 2 > h[3];
00417     bool valid;
00418     double val, hval;
00419 
00420     const double h_00[] = { 2, 0, 0, 10 };
00421     const double h_01[] = { 0, 0, -8, 0 };
00422     const double h_11[] = { 10, 0, 0, 2 };
00423     MsqMatrix< 2, 2 > h00( h_00 ), h01( h_01 ), h11( h_11 );
00424 
00425     valid = metric.evaluate_with_grad( A, val, g, err );
00426     ASSERT_NO_ERROR( err );
00427     CPPUNIT_ASSERT( valid );
00428     valid = metric.evaluate_with_hess( A, hval, gh, h, err );
00429     ASSERT_NO_ERROR( err );
00430     CPPUNIT_ASSERT( valid );
00431     CPPUNIT_ASSERT_DOUBLES_EQUAL( val, hval, 1e-6 );
00432     ASSERT_MATRICES_EQUAL( g, gh, 1e-6 );
00433     ASSERT_MATRICES_EQUAL( h00, h[0], 1e-6 );
00434     ASSERT_MATRICES_EQUAL( h01, h[1], 1e-6 );
00435     ASSERT_MATRICES_EQUAL( h11, h[2], 1e-6 );
00436 
00437     valid = metric.evaluate_with_grad( B, val, g, err );
00438     ASSERT_NO_ERROR( err );
00439     CPPUNIT_ASSERT( valid );
00440     valid = metric.evaluate_with_hess( B, hval, gh, h, err );
00441     ASSERT_NO_ERROR( err );
00442     CPPUNIT_ASSERT( valid );
00443     CPPUNIT_ASSERT_DOUBLES_EQUAL( val, hval, 1e-6 );
00444     ASSERT_MATRICES_EQUAL( g, gh, 1e-6 );
00445     ASSERT_MATRICES_EQUAL( h00, h[0], 1e-6 );
00446     ASSERT_MATRICES_EQUAL( h01, h[1], 1e-6 );
00447     ASSERT_MATRICES_EQUAL( h11, h[2], 1e-6 );
00448 
00449     valid = metric.evaluate_with_grad( inverse( A ), val, g, err );
00450     ASSERT_NO_ERROR( err );
00451     CPPUNIT_ASSERT( valid );
00452     valid = metric.evaluate_with_hess( inverse( A ), hval, gh, h, err );
00453     ASSERT_NO_ERROR( err );
00454     CPPUNIT_ASSERT( valid );
00455     CPPUNIT_ASSERT_DOUBLES_EQUAL( val, hval, 1e-6 );
00456     ASSERT_MATRICES_EQUAL( g, gh, 1e-6 );
00457     ASSERT_MATRICES_EQUAL( h00, h[0], 1e-6 );
00458     ASSERT_MATRICES_EQUAL( h01, h[1], 1e-6 );
00459     ASSERT_MATRICES_EQUAL( h11, h[2], 1e-6 );
00460 
00461     valid = metric.evaluate_with_grad( inverse( B ), val, g, err );
00462     ASSERT_NO_ERROR( err );
00463     CPPUNIT_ASSERT( valid );
00464     valid = metric.evaluate_with_hess( inverse( B ), hval, gh, h, err );
00465     ASSERT_NO_ERROR( err );
00466     CPPUNIT_ASSERT( valid );
00467     CPPUNIT_ASSERT_DOUBLES_EQUAL( val, hval, 1e-6 );
00468     ASSERT_MATRICES_EQUAL( g, gh, 1e-6 );
00469     ASSERT_MATRICES_EQUAL( h00, h[0], 1e-6 );
00470     ASSERT_MATRICES_EQUAL( h01, h[1], 1e-6 );
00471     ASSERT_MATRICES_EQUAL( h11, h[2], 1e-6 );
00472 
00473     valid = metric.evaluate_with_grad( A * inverse( B ), val, g, err );
00474     ASSERT_NO_ERROR( err );
00475     CPPUNIT_ASSERT( valid );
00476     valid = metric.evaluate_with_hess( A * inverse( B ), hval, gh, h, err );
00477     ASSERT_NO_ERROR( err );
00478     CPPUNIT_ASSERT( valid );
00479     CPPUNIT_ASSERT_DOUBLES_EQUAL( val, hval, 1e-6 );
00480     ASSERT_MATRICES_EQUAL( g, gh, 1e-6 );
00481     ASSERT_MATRICES_EQUAL( h00, h[0], 1e-6 );
00482     ASSERT_MATRICES_EQUAL( h01, h[1], 1e-6 );
00483     ASSERT_MATRICES_EQUAL( h11, h[2], 1e-6 );
00484 
00485     valid = metric.evaluate_with_grad( B * inverse( A ), val, g, err );
00486     ASSERT_NO_ERROR( err );
00487     CPPUNIT_ASSERT( valid );
00488     valid = metric.evaluate_with_hess( B * inverse( A ), hval, gh, h, err );
00489     ASSERT_NO_ERROR( err );
00490     CPPUNIT_ASSERT( valid );
00491     CPPUNIT_ASSERT_DOUBLES_EQUAL( val, hval, 1e-6 );
00492     ASSERT_MATRICES_EQUAL( g, gh, 1e-6 );
00493     ASSERT_MATRICES_EQUAL( h00, h[0], 1e-6 );
00494     ASSERT_MATRICES_EQUAL( h01, h[1], 1e-6 );
00495     ASSERT_MATRICES_EQUAL( h11, h[2], 1e-6 );
00496 
00497     MsqMatrix< 2, 2 > ah[3];
00498     valid = metric2.evaluate_with_hess( A, val, g, ah, err );
00499     ASSERT_NO_ERROR( err );
00500     CPPUNIT_ASSERT( valid );
00501     valid = metric2.TMetric::evaluate_with_hess( A, hval, gh, h, err );
00502     ASSERT_NO_ERROR( err );
00503     CPPUNIT_ASSERT( valid );
00504     CPPUNIT_ASSERT_DOUBLES_EQUAL( val, hval, 1e-6 );
00505     ASSERT_MATRICES_EQUAL( g, gh, 1e-6 );
00506     ASSERT_MATRICES_EQUAL( ah[0], h[0], 1e-6 );
00507     ASSERT_MATRICES_EQUAL( ah[1], h[1], 1e-6 );
00508     ASSERT_MATRICES_EQUAL( ah[2], h[2], 1e-6 );
00509 
00510     valid = metric2.evaluate_with_hess( B, val, g, ah, err );
00511     ASSERT_NO_ERROR( err );
00512     CPPUNIT_ASSERT( valid );
00513     valid = metric2.TMetric::evaluate_with_hess( B, hval, gh, h, err );
00514     ASSERT_NO_ERROR( err );
00515     CPPUNIT_ASSERT( valid );
00516     CPPUNIT_ASSERT_DOUBLES_EQUAL( val, hval, 1e-6 );
00517     ASSERT_MATRICES_EQUAL( g, gh, 1e-6 );
00518     ASSERT_MATRICES_EQUAL( ah[0], h[0], 1e-6 );
00519     ASSERT_MATRICES_EQUAL( ah[1], h[1], 1e-6 );
00520     ASSERT_MATRICES_EQUAL( ah[2], h[2], 1e-6 );
00521 }
00522 
00523 void TMetricTest::test_numerical_hessian_3D()
00524 {
00525     HessTestMetricRel metric;
00526     HessTestMetricRel_2 metric2;
00527     const double Avals[] = { 1, 2, 3, 4, 1, 4, 3, 2, 1 };
00528     const double Bvals[] = { 0.1, 0.15, 0.05, 0.2, -0.1, -0.15, -0.05, -0.2, 2 };
00529     const MsqMatrix< 3, 3 > A( Avals );
00530     const MsqMatrix< 3, 3 > B( Bvals );
00531 
00532     MsqError err;
00533     MsqMatrix< 3, 3 > g, gh;
00534     MsqMatrix< 3, 3 > h[6];
00535     bool valid;
00536     double val, hval;
00537 
00538     const double h_00[] = { 2, 0, 0, 0, 10, 0, 0, 0, 10 };
00539     const double h_01[] = { 0, 0, 0, -8, 0, 0, 0, 0, 0 };
00540     const double h_02[] = { 0, 0, 0, 0, 0, 0, -8, 0, 0 };
00541     const double h_11[] = { 10, 0, 0, 0, 2, 0, 0, 0, 10 };
00542     const double h_12[] = { 0, 0, 0, 0, 0, 0, 0, -8, 0 };
00543     const double h_22[] = { 10, 0, 0, 0, 10, 0, 0, 0, 2 };
00544     MsqMatrix< 3, 3 > h00( h_00 ), h01( h_01 ), h02( h_02 ), h11( h_11 ), h12( h_12 ), h22( h_22 );
00545 
00546     valid = metric.evaluate_with_grad( A, val, g, err );
00547     ASSERT_NO_ERROR( err );
00548     CPPUNIT_ASSERT( valid );
00549     valid = metric.evaluate_with_hess( A, hval, gh, h, err );
00550     ASSERT_NO_ERROR( err );
00551     CPPUNIT_ASSERT( valid );
00552     CPPUNIT_ASSERT_DOUBLES_EQUAL( val, hval, 1e-6 );
00553     ASSERT_MATRICES_EQUAL( g, gh, 1e-6 );
00554     ASSERT_MATRICES_EQUAL( h00, h[0], 1e-6 );
00555     ASSERT_MATRICES_EQUAL( h01, h[1], 1e-6 );
00556     ASSERT_MATRICES_EQUAL( h02, h[2], 1e-6 );
00557     ASSERT_MATRICES_EQUAL( h11, h[3], 1e-6 );
00558     ASSERT_MATRICES_EQUAL( h12, h[4], 1e-6 );
00559     ASSERT_MATRICES_EQUAL( h22, h[5], 1e-6 );
00560 
00561     valid = metric.evaluate_with_grad( B, val, g, err );
00562     ASSERT_NO_ERROR( err );
00563     CPPUNIT_ASSERT( valid );
00564     valid = metric.evaluate_with_hess( B, hval, gh, h, err );
00565     ASSERT_NO_ERROR( err );
00566     CPPUNIT_ASSERT( valid );
00567     CPPUNIT_ASSERT_DOUBLES_EQUAL( val, hval, 1e-6 );
00568     ASSERT_MATRICES_EQUAL( g, gh, 1e-6 );
00569     ASSERT_MATRICES_EQUAL( h00, h[0], 1e-6 );
00570     ASSERT_MATRICES_EQUAL( h01, h[1], 1e-6 );
00571     ASSERT_MATRICES_EQUAL( h02, h[2], 1e-6 );
00572     ASSERT_MATRICES_EQUAL( h11, h[3], 1e-6 );
00573     ASSERT_MATRICES_EQUAL( h12, h[4], 1e-6 );
00574     ASSERT_MATRICES_EQUAL( h22, h[5], 1e-6 );
00575 
00576     valid = metric.evaluate_with_grad( inverse( A ), val, g, err );
00577     ASSERT_NO_ERROR( err );
00578     CPPUNIT_ASSERT( valid );
00579     valid = metric.evaluate_with_hess( inverse( A ), hval, gh, h, err );
00580     ASSERT_NO_ERROR( err );
00581     CPPUNIT_ASSERT( valid );
00582     CPPUNIT_ASSERT_DOUBLES_EQUAL( val, hval, 1e-6 );
00583     ASSERT_MATRICES_EQUAL( g, gh, 1e-6 );
00584     ASSERT_MATRICES_EQUAL( h00, h[0], 1e-6 );
00585     ASSERT_MATRICES_EQUAL( h01, h[1], 1e-6 );
00586     ASSERT_MATRICES_EQUAL( h02, h[2], 1e-6 );
00587     ASSERT_MATRICES_EQUAL( h11, h[3], 1e-6 );
00588     ASSERT_MATRICES_EQUAL( h12, h[4], 1e-6 );
00589     ASSERT_MATRICES_EQUAL( h22, h[5], 1e-6 );
00590 
00591     valid = metric.evaluate_with_grad( inverse( B ), val, g, err );
00592     ASSERT_NO_ERROR( err );
00593     CPPUNIT_ASSERT( valid );
00594     valid = metric.evaluate_with_hess( inverse( B ), hval, gh, h, err );
00595     ASSERT_NO_ERROR( err );
00596     CPPUNIT_ASSERT( valid );
00597     CPPUNIT_ASSERT_DOUBLES_EQUAL( val, hval, 1e-6 );
00598     ASSERT_MATRICES_EQUAL( g, gh, 1e-6 );
00599     ASSERT_MATRICES_EQUAL( h00, h[0], 1e-6 );
00600     ASSERT_MATRICES_EQUAL( h01, h[1], 1e-6 );
00601     ASSERT_MATRICES_EQUAL( h02, h[2], 1e-6 );
00602     ASSERT_MATRICES_EQUAL( h11, h[3], 1e-6 );
00603     ASSERT_MATRICES_EQUAL( h12, h[4], 1e-6 );
00604     ASSERT_MATRICES_EQUAL( h22, h[5], 1e-6 );
00605 
00606     valid = metric.evaluate_with_grad( A * inverse( B ), val, g, err );
00607     ASSERT_NO_ERROR( err );
00608     CPPUNIT_ASSERT( valid );
00609     valid = metric.evaluate_with_hess( A * inverse( B ), hval, gh, h, err );
00610     ASSERT_NO_ERROR( err );
00611     CPPUNIT_ASSERT( valid );
00612     CPPUNIT_ASSERT_DOUBLES_EQUAL( val, hval, 1e-6 );
00613     ASSERT_MATRICES_EQUAL( g, gh, 1e-6 );
00614     ASSERT_MATRICES_EQUAL( h00, h[0], 1e-6 );
00615     ASSERT_MATRICES_EQUAL( h01, h[1], 1e-6 );
00616     ASSERT_MATRICES_EQUAL( h02, h[2], 1e-6 );
00617     ASSERT_MATRICES_EQUAL( h11, h[3], 1e-6 );
00618     ASSERT_MATRICES_EQUAL( h12, h[4], 1e-6 );
00619     ASSERT_MATRICES_EQUAL( h22, h[5], 1e-6 );
00620 
00621     valid = metric.evaluate_with_grad( B * inverse( A ), val, g, err );
00622     ASSERT_NO_ERROR( err );
00623     CPPUNIT_ASSERT( valid );
00624     valid = metric.evaluate_with_hess( B * inverse( A ), hval, gh, h, err );
00625     ASSERT_NO_ERROR( err );
00626     CPPUNIT_ASSERT( valid );
00627     CPPUNIT_ASSERT_DOUBLES_EQUAL( val, hval, 1e-6 );
00628     ASSERT_MATRICES_EQUAL( g, gh, 1e-6 );
00629     ASSERT_MATRICES_EQUAL( h00, h[0], 1e-6 );
00630     ASSERT_MATRICES_EQUAL( h01, h[1], 1e-6 );
00631     ASSERT_MATRICES_EQUAL( h02, h[2], 1e-6 );
00632     ASSERT_MATRICES_EQUAL( h11, h[3], 1e-6 );
00633     ASSERT_MATRICES_EQUAL( h12, h[4], 1e-6 );
00634     ASSERT_MATRICES_EQUAL( h22, h[5], 1e-6 );
00635 
00636     MsqMatrix< 3, 3 > ah[6];
00637     valid = metric2.evaluate_with_hess( A, val, g, ah, err );
00638     ASSERT_NO_ERROR( err );
00639     CPPUNIT_ASSERT( valid );
00640     valid = metric2.TMetric::evaluate_with_hess( A, hval, gh, h, err );
00641     ASSERT_NO_ERROR( err );
00642     CPPUNIT_ASSERT( valid );
00643     CPPUNIT_ASSERT_DOUBLES_EQUAL( val, hval, 1e-6 );
00644     ASSERT_MATRICES_EQUAL( g, gh, 1e-6 );
00645     ASSERT_MATRICES_EQUAL( ah[0], h[0], 1e-6 );
00646     ASSERT_MATRICES_EQUAL( ah[1], h[1], 1e-6 );
00647     ASSERT_MATRICES_EQUAL( ah[2], h[2], 1e-6 );
00648     ASSERT_MATRICES_EQUAL( ah[3], h[3], 1e-6 );
00649     ASSERT_MATRICES_EQUAL( ah[4], h[4], 1e-6 );
00650     ASSERT_MATRICES_EQUAL( ah[5], h[5], 1e-6 );
00651 
00652     valid = metric2.evaluate_with_hess( B, val, g, ah, err );
00653     ASSERT_NO_ERROR( err );
00654     CPPUNIT_ASSERT( valid );
00655     valid = metric2.TMetric::evaluate_with_hess( B, hval, gh, h, err );
00656     ASSERT_NO_ERROR( err );
00657     CPPUNIT_ASSERT( valid );
00658     CPPUNIT_ASSERT_DOUBLES_EQUAL( val, hval, 1e-6 );
00659     ASSERT_MATRICES_EQUAL( g, gh, 1e-6 );
00660     ASSERT_MATRICES_EQUAL( ah[0], h[0], 1e-6 );
00661     ASSERT_MATRICES_EQUAL( ah[1], h[1], 1e-6 );
00662     ASSERT_MATRICES_EQUAL( ah[2], h[2], 1e-6 );
00663     ASSERT_MATRICES_EQUAL( ah[3], h[3], 1e-6 );
00664     ASSERT_MATRICES_EQUAL( ah[4], h[4], 1e-6 );
00665     ASSERT_MATRICES_EQUAL( ah[5], h[5], 1e-6 );
00666 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines