MOAB: Mesh Oriented datABase  (version 5.2.1)
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, double& result, MsqMatrix< 2, 2 >& d, MsqMatrix< 2, 2 > d2[3],
00170                              MsqError& err )
00171     {
00172         result = Frobenius( T );
00173         d      = T / result;
00174         int h  = 0;
00175         for( int r = 0; r < 2; ++r )
00176         {
00177             int i = h;
00178             for( int c = r; c < 2; ++c )
00179                 d2[h++] = transpose( T.row( r ) ) * T.row( c ) / -( result * result * result );
00180             d2[i] += MsqMatrix< 2, 2 >( 1.0 / result );
00181         }
00182         return true;
00183     }
00184 
00185     bool evaluate( const MsqMatrix< 3, 3 >& T, double& result, MsqError& err )
00186     {
00187         result = Frobenius( T );
00188         return true;
00189     }
00190 
00191     bool evaluate_with_grad( const MsqMatrix< 3, 3 >& T, double& result, MsqMatrix< 3, 3 >& d, MsqError& err )
00192     {
00193         result = Frobenius( T );
00194         d      = T / result;
00195         return true;
00196     }
00197 
00198     bool evaluate_with_hess( const MsqMatrix< 3, 3 >& T, double& result, MsqMatrix< 3, 3 >& d, MsqMatrix< 3, 3 > d2[6],
00199                              MsqError& err )
00200     {
00201         result = Frobenius( T );
00202         d      = T / result;
00203         int h  = 0;
00204         for( int r = 0; r < 3; ++r )
00205         {
00206             int i = h;
00207             for( int c = r; c < 3; ++c )
00208                 d2[h++] = transpose( T.row( r ) ) * T.row( c ) / -( result * result * result );
00209             d2[i] += MsqMatrix< 3, 3 >( 1.0 / result );
00210         }
00211         return true;
00212     }
00213 };
00214 
00215 void TMetricTest::test_numerical_gradient_2D()
00216 {
00217     GradTestMetricRel metric;
00218     HessTestMetricRel_2 metric2;
00219     const double Avals[] = { 1, 2, 2, 5 };
00220     const double Bvals[] = { -0.1, -0.15, -0.25, -0.8 };
00221     const MsqMatrix< 2, 2 > A( Avals );
00222     const MsqMatrix< 2, 2 > B( Bvals );
00223 
00224     MsqError err;
00225     MsqMatrix< 2, 2 > d;
00226     bool valid;
00227     double val, gval;
00228 
00229     MsqMatrix< 2, 2 > expected;
00230     for( int r = 0; r < 2; ++r )
00231         for( int c = 0; c < 2; ++c )
00232             expected( r, c ) = metric.grad( r, c );
00233 
00234     valid = metric.evaluate( A, val, err );
00235     ASSERT_NO_ERROR( err );
00236     CPPUNIT_ASSERT( valid );
00237     valid = metric.evaluate_with_grad( A, gval, d, err );
00238     ASSERT_NO_ERROR( err );
00239     CPPUNIT_ASSERT( valid );
00240     CPPUNIT_ASSERT_DOUBLES_EQUAL( val, gval, 1e-6 );
00241     ASSERT_MATRICES_EQUAL( expected, d, 1e-6 );
00242 
00243     valid = metric.evaluate( B, val, err );
00244     ASSERT_NO_ERROR( err );
00245     CPPUNIT_ASSERT( valid );
00246     valid = metric.evaluate_with_grad( B, gval, d, err );
00247     ASSERT_NO_ERROR( err );
00248     CPPUNIT_ASSERT( valid );
00249     CPPUNIT_ASSERT_DOUBLES_EQUAL( val, gval, 1e-6 );
00250     ASSERT_MATRICES_EQUAL( expected, d, 1e-6 );
00251 
00252     valid = metric.evaluate( inverse( A ), val, err );
00253     ASSERT_NO_ERROR( err );
00254     CPPUNIT_ASSERT( valid );
00255     valid = metric.evaluate_with_grad( inverse( A ), gval, d, err );
00256     ASSERT_NO_ERROR( err );
00257     CPPUNIT_ASSERT( valid );
00258     ASSERT_MATRICES_EQUAL( expected, d, 1e-6 );
00259 
00260     valid = metric.evaluate( inverse( B ), val, err );
00261     ASSERT_NO_ERROR( err );
00262     CPPUNIT_ASSERT( valid );
00263     valid = metric.evaluate_with_grad( inverse( B ), gval, d, err );
00264     ASSERT_NO_ERROR( err );
00265     CPPUNIT_ASSERT( valid );
00266     CPPUNIT_ASSERT_DOUBLES_EQUAL( val, gval, 1e-6 );
00267     ASSERT_MATRICES_EQUAL( expected, d, 1e-6 );
00268 
00269     valid = metric.evaluate( A * inverse( B ), val, err );
00270     ASSERT_NO_ERROR( err );
00271     CPPUNIT_ASSERT( valid );
00272     valid = metric.evaluate_with_grad( A * inverse( B ), gval, d, err );
00273     ASSERT_NO_ERROR( err );
00274     CPPUNIT_ASSERT( valid );
00275     CPPUNIT_ASSERT_DOUBLES_EQUAL( val, gval, 1e-6 );
00276     ASSERT_MATRICES_EQUAL( expected, d, 1e-6 );
00277 
00278     valid = metric.evaluate( B * inverse( A ), val, err );
00279     ASSERT_NO_ERROR( err );
00280     CPPUNIT_ASSERT( valid );
00281     valid = metric.evaluate_with_grad( B * inverse( A ), gval, d, err );
00282     ASSERT_NO_ERROR( err );
00283     CPPUNIT_ASSERT( valid );
00284     CPPUNIT_ASSERT_DOUBLES_EQUAL( val, gval, 1e-6 );
00285     ASSERT_MATRICES_EQUAL( expected, d, 1e-6 );
00286 
00287     MsqMatrix< 2, 2 > da;
00288     valid = metric2.evaluate_with_grad( A, val, da, err );
00289     ASSERT_NO_ERROR( err );
00290     CPPUNIT_ASSERT( valid );
00291     valid = metric2.TMetric::evaluate_with_grad( A, gval, d, err );
00292     ASSERT_NO_ERROR( err );
00293     CPPUNIT_ASSERT( valid );
00294     CPPUNIT_ASSERT_DOUBLES_EQUAL( val, gval, 1e-6 );
00295     ASSERT_MATRICES_EQUAL( da, d, 1e-6 );
00296 
00297     valid = metric2.evaluate_with_grad( B, val, da, err );
00298     ASSERT_NO_ERROR( err );
00299     CPPUNIT_ASSERT( valid );
00300     valid = metric2.TMetric::evaluate_with_grad( B, gval, d, err );
00301     ASSERT_NO_ERROR( err );
00302     CPPUNIT_ASSERT( valid );
00303     CPPUNIT_ASSERT_DOUBLES_EQUAL( val, gval, 1e-6 );
00304     ASSERT_MATRICES_EQUAL( da, d, 1e-6 );
00305 }
00306 
00307 void TMetricTest::test_numerical_gradient_3D()
00308 {
00309     GradTestMetricRel metric;
00310     HessTestMetricRel_2 metric2;
00311     const double Avals[] = { 1, 2, 3, 4, 1, 4, 3, 2, 1 };
00312     const double Bvals[] = { 0.1, 0.15, 0.05, 0.2, -0.1, -0.15, -0.05, -0.2, 2 };
00313     const MsqMatrix< 3, 3 > A( Avals );
00314     const MsqMatrix< 3, 3 > B( Bvals );
00315 
00316     MsqError err;
00317     MsqMatrix< 3, 3 > d;
00318     bool valid;
00319     double val, gval;
00320 
00321     MsqMatrix< 3, 3 > expected;
00322     for( int r = 0; r < 3; ++r )
00323         for( int c = 0; c < 3; ++c )
00324             expected( r, c ) = metric.grad( r, c );
00325 
00326     valid = metric.evaluate( A, val, err );
00327     ASSERT_NO_ERROR( err );
00328     CPPUNIT_ASSERT( valid );
00329     valid = metric.evaluate_with_grad( A, gval, d, err );
00330     ASSERT_NO_ERROR( err );
00331     CPPUNIT_ASSERT( valid );
00332     CPPUNIT_ASSERT_DOUBLES_EQUAL( val, gval, 1e-6 );
00333     ASSERT_MATRICES_EQUAL( expected, d, 1e-6 );
00334 
00335     valid = metric.evaluate( B, val, err );
00336     ASSERT_NO_ERROR( err );
00337     CPPUNIT_ASSERT( valid );
00338     valid = metric.evaluate_with_grad( B, gval, d, err );
00339     ASSERT_NO_ERROR( err );
00340     CPPUNIT_ASSERT( valid );
00341     CPPUNIT_ASSERT_DOUBLES_EQUAL( val, gval, 1e-6 );
00342     ASSERT_MATRICES_EQUAL( expected, d, 1e-6 );
00343 
00344     valid = metric.evaluate( inverse( A ), val, err );
00345     ASSERT_NO_ERROR( err );
00346     CPPUNIT_ASSERT( valid );
00347     valid = metric.evaluate_with_grad( inverse( A ), gval, d, err );
00348     ASSERT_NO_ERROR( err );
00349     CPPUNIT_ASSERT( valid );
00350     ASSERT_MATRICES_EQUAL( expected, d, 1e-6 );
00351 
00352     valid = metric.evaluate( inverse( B ), val, err );
00353     ASSERT_NO_ERROR( err );
00354     CPPUNIT_ASSERT( valid );
00355     valid = metric.evaluate_with_grad( inverse( B ), gval, d, err );
00356     ASSERT_NO_ERROR( err );
00357     CPPUNIT_ASSERT( valid );
00358     CPPUNIT_ASSERT_DOUBLES_EQUAL( val, gval, 1e-6 );
00359     ASSERT_MATRICES_EQUAL( expected, d, 1e-6 );
00360 
00361     valid = metric.evaluate( A * inverse( B ), val, err );
00362     ASSERT_NO_ERROR( err );
00363     CPPUNIT_ASSERT( valid );
00364     valid = metric.evaluate_with_grad( A * inverse( B ), gval, d, err );
00365     ASSERT_NO_ERROR( err );
00366     CPPUNIT_ASSERT( valid );
00367     CPPUNIT_ASSERT_DOUBLES_EQUAL( val, gval, 1e-6 );
00368     ASSERT_MATRICES_EQUAL( expected, d, 1e-6 );
00369 
00370     valid = metric.evaluate( B * inverse( A ), val, err );
00371     ASSERT_NO_ERROR( err );
00372     CPPUNIT_ASSERT( valid );
00373     valid = metric.evaluate_with_grad( B * inverse( A ), gval, d, err );
00374     ASSERT_NO_ERROR( err );
00375     CPPUNIT_ASSERT( valid );
00376     CPPUNIT_ASSERT_DOUBLES_EQUAL( val, gval, 1e-6 );
00377     ASSERT_MATRICES_EQUAL( expected, d, 1e-6 );
00378 
00379     MsqMatrix< 3, 3 > da;
00380     valid = metric2.evaluate_with_grad( A, val, da, err );
00381     ASSERT_NO_ERROR( err );
00382     CPPUNIT_ASSERT( valid );
00383     valid = metric2.TMetric::evaluate_with_grad( A, gval, d, err );
00384     ASSERT_NO_ERROR( err );
00385     CPPUNIT_ASSERT( valid );
00386     CPPUNIT_ASSERT_DOUBLES_EQUAL( val, gval, 1e-6 );
00387     ASSERT_MATRICES_EQUAL( da, d, 1e-6 );
00388 
00389     valid = metric2.evaluate_with_grad( B, val, da, err );
00390     ASSERT_NO_ERROR( err );
00391     CPPUNIT_ASSERT( valid );
00392     valid = metric2.TMetric::evaluate_with_grad( B, gval, d, err );
00393     ASSERT_NO_ERROR( err );
00394     CPPUNIT_ASSERT( valid );
00395     CPPUNIT_ASSERT_DOUBLES_EQUAL( val, gval, 1e-6 );
00396     ASSERT_MATRICES_EQUAL( da, d, 1e-6 );
00397 }
00398 
00399 void TMetricTest::test_numerical_hessian_2D()
00400 {
00401     HessTestMetricRel metric;
00402     HessTestMetricRel_2 metric2;
00403     const double Avals[] = { 1, 2, 2, 5 };
00404     const double Bvals[] = { -0.1, -0.15, -0.25, -0.8 };
00405     const MsqMatrix< 2, 2 > A( Avals );
00406     const MsqMatrix< 2, 2 > B( Bvals );
00407 
00408     MsqError err;
00409     MsqMatrix< 2, 2 > g, gh;
00410     MsqMatrix< 2, 2 > h[3];
00411     bool valid;
00412     double val, hval;
00413 
00414     const double h_00[] = { 2, 0, 0, 10 };
00415     const double h_01[] = { 0, 0, -8, 0 };
00416     const double h_11[] = { 10, 0, 0, 2 };
00417     MsqMatrix< 2, 2 > h00( h_00 ), h01( h_01 ), h11( h_11 );
00418 
00419     valid = metric.evaluate_with_grad( A, val, g, err );
00420     ASSERT_NO_ERROR( err );
00421     CPPUNIT_ASSERT( valid );
00422     valid = metric.evaluate_with_hess( A, hval, gh, h, err );
00423     ASSERT_NO_ERROR( err );
00424     CPPUNIT_ASSERT( valid );
00425     CPPUNIT_ASSERT_DOUBLES_EQUAL( val, hval, 1e-6 );
00426     ASSERT_MATRICES_EQUAL( g, gh, 1e-6 );
00427     ASSERT_MATRICES_EQUAL( h00, h[0], 1e-6 );
00428     ASSERT_MATRICES_EQUAL( h01, h[1], 1e-6 );
00429     ASSERT_MATRICES_EQUAL( h11, h[2], 1e-6 );
00430 
00431     valid = metric.evaluate_with_grad( B, val, g, err );
00432     ASSERT_NO_ERROR( err );
00433     CPPUNIT_ASSERT( valid );
00434     valid = metric.evaluate_with_hess( B, hval, gh, h, err );
00435     ASSERT_NO_ERROR( err );
00436     CPPUNIT_ASSERT( valid );
00437     CPPUNIT_ASSERT_DOUBLES_EQUAL( val, hval, 1e-6 );
00438     ASSERT_MATRICES_EQUAL( g, gh, 1e-6 );
00439     ASSERT_MATRICES_EQUAL( h00, h[0], 1e-6 );
00440     ASSERT_MATRICES_EQUAL( h01, h[1], 1e-6 );
00441     ASSERT_MATRICES_EQUAL( h11, h[2], 1e-6 );
00442 
00443     valid = metric.evaluate_with_grad( inverse( A ), val, g, err );
00444     ASSERT_NO_ERROR( err );
00445     CPPUNIT_ASSERT( valid );
00446     valid = metric.evaluate_with_hess( inverse( A ), hval, gh, h, err );
00447     ASSERT_NO_ERROR( err );
00448     CPPUNIT_ASSERT( valid );
00449     CPPUNIT_ASSERT_DOUBLES_EQUAL( val, hval, 1e-6 );
00450     ASSERT_MATRICES_EQUAL( g, gh, 1e-6 );
00451     ASSERT_MATRICES_EQUAL( h00, h[0], 1e-6 );
00452     ASSERT_MATRICES_EQUAL( h01, h[1], 1e-6 );
00453     ASSERT_MATRICES_EQUAL( h11, h[2], 1e-6 );
00454 
00455     valid = metric.evaluate_with_grad( inverse( B ), val, g, err );
00456     ASSERT_NO_ERROR( err );
00457     CPPUNIT_ASSERT( valid );
00458     valid = metric.evaluate_with_hess( inverse( B ), hval, gh, h, err );
00459     ASSERT_NO_ERROR( err );
00460     CPPUNIT_ASSERT( valid );
00461     CPPUNIT_ASSERT_DOUBLES_EQUAL( val, hval, 1e-6 );
00462     ASSERT_MATRICES_EQUAL( g, gh, 1e-6 );
00463     ASSERT_MATRICES_EQUAL( h00, h[0], 1e-6 );
00464     ASSERT_MATRICES_EQUAL( h01, h[1], 1e-6 );
00465     ASSERT_MATRICES_EQUAL( h11, h[2], 1e-6 );
00466 
00467     valid = metric.evaluate_with_grad( A * inverse( B ), val, g, err );
00468     ASSERT_NO_ERROR( err );
00469     CPPUNIT_ASSERT( valid );
00470     valid = metric.evaluate_with_hess( A * inverse( B ), hval, gh, h, err );
00471     ASSERT_NO_ERROR( err );
00472     CPPUNIT_ASSERT( valid );
00473     CPPUNIT_ASSERT_DOUBLES_EQUAL( val, hval, 1e-6 );
00474     ASSERT_MATRICES_EQUAL( g, gh, 1e-6 );
00475     ASSERT_MATRICES_EQUAL( h00, h[0], 1e-6 );
00476     ASSERT_MATRICES_EQUAL( h01, h[1], 1e-6 );
00477     ASSERT_MATRICES_EQUAL( h11, h[2], 1e-6 );
00478 
00479     valid = metric.evaluate_with_grad( B * inverse( A ), val, g, err );
00480     ASSERT_NO_ERROR( err );
00481     CPPUNIT_ASSERT( valid );
00482     valid = metric.evaluate_with_hess( B * inverse( A ), hval, gh, h, err );
00483     ASSERT_NO_ERROR( err );
00484     CPPUNIT_ASSERT( valid );
00485     CPPUNIT_ASSERT_DOUBLES_EQUAL( val, hval, 1e-6 );
00486     ASSERT_MATRICES_EQUAL( g, gh, 1e-6 );
00487     ASSERT_MATRICES_EQUAL( h00, h[0], 1e-6 );
00488     ASSERT_MATRICES_EQUAL( h01, h[1], 1e-6 );
00489     ASSERT_MATRICES_EQUAL( h11, h[2], 1e-6 );
00490 
00491     MsqMatrix< 2, 2 > ah[3];
00492     valid = metric2.evaluate_with_hess( A, val, g, ah, err );
00493     ASSERT_NO_ERROR( err );
00494     CPPUNIT_ASSERT( valid );
00495     valid = metric2.TMetric::evaluate_with_hess( A, hval, gh, h, err );
00496     ASSERT_NO_ERROR( err );
00497     CPPUNIT_ASSERT( valid );
00498     CPPUNIT_ASSERT_DOUBLES_EQUAL( val, hval, 1e-6 );
00499     ASSERT_MATRICES_EQUAL( g, gh, 1e-6 );
00500     ASSERT_MATRICES_EQUAL( ah[0], h[0], 1e-6 );
00501     ASSERT_MATRICES_EQUAL( ah[1], h[1], 1e-6 );
00502     ASSERT_MATRICES_EQUAL( ah[2], h[2], 1e-6 );
00503 
00504     valid = metric2.evaluate_with_hess( B, val, g, ah, err );
00505     ASSERT_NO_ERROR( err );
00506     CPPUNIT_ASSERT( valid );
00507     valid = metric2.TMetric::evaluate_with_hess( B, hval, gh, h, err );
00508     ASSERT_NO_ERROR( err );
00509     CPPUNIT_ASSERT( valid );
00510     CPPUNIT_ASSERT_DOUBLES_EQUAL( val, hval, 1e-6 );
00511     ASSERT_MATRICES_EQUAL( g, gh, 1e-6 );
00512     ASSERT_MATRICES_EQUAL( ah[0], h[0], 1e-6 );
00513     ASSERT_MATRICES_EQUAL( ah[1], h[1], 1e-6 );
00514     ASSERT_MATRICES_EQUAL( ah[2], h[2], 1e-6 );
00515 }
00516 
00517 void TMetricTest::test_numerical_hessian_3D()
00518 {
00519     HessTestMetricRel metric;
00520     HessTestMetricRel_2 metric2;
00521     const double Avals[] = { 1, 2, 3, 4, 1, 4, 3, 2, 1 };
00522     const double Bvals[] = { 0.1, 0.15, 0.05, 0.2, -0.1, -0.15, -0.05, -0.2, 2 };
00523     const MsqMatrix< 3, 3 > A( Avals );
00524     const MsqMatrix< 3, 3 > B( Bvals );
00525 
00526     MsqError err;
00527     MsqMatrix< 3, 3 > g, gh;
00528     MsqMatrix< 3, 3 > h[6];
00529     bool valid;
00530     double val, hval;
00531 
00532     const double h_00[] = { 2, 0, 0, 0, 10, 0, 0, 0, 10 };
00533     const double h_01[] = { 0, 0, 0, -8, 0, 0, 0, 0, 0 };
00534     const double h_02[] = { 0, 0, 0, 0, 0, 0, -8, 0, 0 };
00535     const double h_11[] = { 10, 0, 0, 0, 2, 0, 0, 0, 10 };
00536     const double h_12[] = { 0, 0, 0, 0, 0, 0, 0, -8, 0 };
00537     const double h_22[] = { 10, 0, 0, 0, 10, 0, 0, 0, 2 };
00538     MsqMatrix< 3, 3 > h00( h_00 ), h01( h_01 ), h02( h_02 ), h11( h_11 ), h12( h_12 ), h22( h_22 );
00539 
00540     valid = metric.evaluate_with_grad( A, val, g, err );
00541     ASSERT_NO_ERROR( err );
00542     CPPUNIT_ASSERT( valid );
00543     valid = metric.evaluate_with_hess( A, hval, gh, h, err );
00544     ASSERT_NO_ERROR( err );
00545     CPPUNIT_ASSERT( valid );
00546     CPPUNIT_ASSERT_DOUBLES_EQUAL( val, hval, 1e-6 );
00547     ASSERT_MATRICES_EQUAL( g, gh, 1e-6 );
00548     ASSERT_MATRICES_EQUAL( h00, h[0], 1e-6 );
00549     ASSERT_MATRICES_EQUAL( h01, h[1], 1e-6 );
00550     ASSERT_MATRICES_EQUAL( h02, h[2], 1e-6 );
00551     ASSERT_MATRICES_EQUAL( h11, h[3], 1e-6 );
00552     ASSERT_MATRICES_EQUAL( h12, h[4], 1e-6 );
00553     ASSERT_MATRICES_EQUAL( h22, h[5], 1e-6 );
00554 
00555     valid = metric.evaluate_with_grad( B, val, g, err );
00556     ASSERT_NO_ERROR( err );
00557     CPPUNIT_ASSERT( valid );
00558     valid = metric.evaluate_with_hess( B, hval, gh, h, err );
00559     ASSERT_NO_ERROR( err );
00560     CPPUNIT_ASSERT( valid );
00561     CPPUNIT_ASSERT_DOUBLES_EQUAL( val, hval, 1e-6 );
00562     ASSERT_MATRICES_EQUAL( g, gh, 1e-6 );
00563     ASSERT_MATRICES_EQUAL( h00, h[0], 1e-6 );
00564     ASSERT_MATRICES_EQUAL( h01, h[1], 1e-6 );
00565     ASSERT_MATRICES_EQUAL( h02, h[2], 1e-6 );
00566     ASSERT_MATRICES_EQUAL( h11, h[3], 1e-6 );
00567     ASSERT_MATRICES_EQUAL( h12, h[4], 1e-6 );
00568     ASSERT_MATRICES_EQUAL( h22, h[5], 1e-6 );
00569 
00570     valid = metric.evaluate_with_grad( inverse( A ), val, g, err );
00571     ASSERT_NO_ERROR( err );
00572     CPPUNIT_ASSERT( valid );
00573     valid = metric.evaluate_with_hess( inverse( A ), hval, gh, h, err );
00574     ASSERT_NO_ERROR( err );
00575     CPPUNIT_ASSERT( valid );
00576     CPPUNIT_ASSERT_DOUBLES_EQUAL( val, hval, 1e-6 );
00577     ASSERT_MATRICES_EQUAL( g, gh, 1e-6 );
00578     ASSERT_MATRICES_EQUAL( h00, h[0], 1e-6 );
00579     ASSERT_MATRICES_EQUAL( h01, h[1], 1e-6 );
00580     ASSERT_MATRICES_EQUAL( h02, h[2], 1e-6 );
00581     ASSERT_MATRICES_EQUAL( h11, h[3], 1e-6 );
00582     ASSERT_MATRICES_EQUAL( h12, h[4], 1e-6 );
00583     ASSERT_MATRICES_EQUAL( h22, h[5], 1e-6 );
00584 
00585     valid = metric.evaluate_with_grad( inverse( B ), val, g, err );
00586     ASSERT_NO_ERROR( err );
00587     CPPUNIT_ASSERT( valid );
00588     valid = metric.evaluate_with_hess( inverse( B ), hval, gh, h, err );
00589     ASSERT_NO_ERROR( err );
00590     CPPUNIT_ASSERT( valid );
00591     CPPUNIT_ASSERT_DOUBLES_EQUAL( val, hval, 1e-6 );
00592     ASSERT_MATRICES_EQUAL( g, gh, 1e-6 );
00593     ASSERT_MATRICES_EQUAL( h00, h[0], 1e-6 );
00594     ASSERT_MATRICES_EQUAL( h01, h[1], 1e-6 );
00595     ASSERT_MATRICES_EQUAL( h02, h[2], 1e-6 );
00596     ASSERT_MATRICES_EQUAL( h11, h[3], 1e-6 );
00597     ASSERT_MATRICES_EQUAL( h12, h[4], 1e-6 );
00598     ASSERT_MATRICES_EQUAL( h22, h[5], 1e-6 );
00599 
00600     valid = metric.evaluate_with_grad( A * inverse( B ), val, g, err );
00601     ASSERT_NO_ERROR( err );
00602     CPPUNIT_ASSERT( valid );
00603     valid = metric.evaluate_with_hess( A * inverse( B ), hval, gh, h, err );
00604     ASSERT_NO_ERROR( err );
00605     CPPUNIT_ASSERT( valid );
00606     CPPUNIT_ASSERT_DOUBLES_EQUAL( val, hval, 1e-6 );
00607     ASSERT_MATRICES_EQUAL( g, gh, 1e-6 );
00608     ASSERT_MATRICES_EQUAL( h00, h[0], 1e-6 );
00609     ASSERT_MATRICES_EQUAL( h01, h[1], 1e-6 );
00610     ASSERT_MATRICES_EQUAL( h02, h[2], 1e-6 );
00611     ASSERT_MATRICES_EQUAL( h11, h[3], 1e-6 );
00612     ASSERT_MATRICES_EQUAL( h12, h[4], 1e-6 );
00613     ASSERT_MATRICES_EQUAL( h22, h[5], 1e-6 );
00614 
00615     valid = metric.evaluate_with_grad( B * inverse( A ), val, g, err );
00616     ASSERT_NO_ERROR( err );
00617     CPPUNIT_ASSERT( valid );
00618     valid = metric.evaluate_with_hess( B * inverse( A ), hval, gh, h, err );
00619     ASSERT_NO_ERROR( err );
00620     CPPUNIT_ASSERT( valid );
00621     CPPUNIT_ASSERT_DOUBLES_EQUAL( val, hval, 1e-6 );
00622     ASSERT_MATRICES_EQUAL( g, gh, 1e-6 );
00623     ASSERT_MATRICES_EQUAL( h00, h[0], 1e-6 );
00624     ASSERT_MATRICES_EQUAL( h01, h[1], 1e-6 );
00625     ASSERT_MATRICES_EQUAL( h02, h[2], 1e-6 );
00626     ASSERT_MATRICES_EQUAL( h11, h[3], 1e-6 );
00627     ASSERT_MATRICES_EQUAL( h12, h[4], 1e-6 );
00628     ASSERT_MATRICES_EQUAL( h22, h[5], 1e-6 );
00629 
00630     MsqMatrix< 3, 3 > ah[6];
00631     valid = metric2.evaluate_with_hess( A, val, g, ah, err );
00632     ASSERT_NO_ERROR( err );
00633     CPPUNIT_ASSERT( valid );
00634     valid = metric2.TMetric::evaluate_with_hess( A, hval, gh, h, err );
00635     ASSERT_NO_ERROR( err );
00636     CPPUNIT_ASSERT( valid );
00637     CPPUNIT_ASSERT_DOUBLES_EQUAL( val, hval, 1e-6 );
00638     ASSERT_MATRICES_EQUAL( g, gh, 1e-6 );
00639     ASSERT_MATRICES_EQUAL( ah[0], h[0], 1e-6 );
00640     ASSERT_MATRICES_EQUAL( ah[1], h[1], 1e-6 );
00641     ASSERT_MATRICES_EQUAL( ah[2], h[2], 1e-6 );
00642     ASSERT_MATRICES_EQUAL( ah[3], h[3], 1e-6 );
00643     ASSERT_MATRICES_EQUAL( ah[4], h[4], 1e-6 );
00644     ASSERT_MATRICES_EQUAL( ah[5], h[5], 1e-6 );
00645 
00646     valid = metric2.evaluate_with_hess( B, val, g, ah, err );
00647     ASSERT_NO_ERROR( err );
00648     CPPUNIT_ASSERT( valid );
00649     valid = metric2.TMetric::evaluate_with_hess( B, hval, gh, h, err );
00650     ASSERT_NO_ERROR( err );
00651     CPPUNIT_ASSERT( valid );
00652     CPPUNIT_ASSERT_DOUBLES_EQUAL( val, hval, 1e-6 );
00653     ASSERT_MATRICES_EQUAL( g, gh, 1e-6 );
00654     ASSERT_MATRICES_EQUAL( ah[0], h[0], 1e-6 );
00655     ASSERT_MATRICES_EQUAL( ah[1], h[1], 1e-6 );
00656     ASSERT_MATRICES_EQUAL( ah[2], h[2], 1e-6 );
00657     ASSERT_MATRICES_EQUAL( ah[3], h[3], 1e-6 );
00658     ASSERT_MATRICES_EQUAL( ah[4], h[4], 1e-6 );
00659     ASSERT_MATRICES_EQUAL( ah[5], h[5], 1e-6 );
00660 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines