MOAB: Mesh Oriented datABase
(version 5.2.1)
|
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 }