MOAB: Mesh Oriented datABase  (version 5.2.1)
ObjectiveFunctionTests.cpp
Go to the documentation of this file.
00001 /* *****************************************************************
00002     MESQUITE -- The Mesh Quality Improvement Toolkit
00003 
00004     Copyright 2010 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 ObjectiveFunctionTests.cpp
00028  *  \brief
00029  *  \author Jason Kraftcheck
00030  */
00031 
00032 #include "ObjectiveFunctionTests.hpp"
00033 
00034 /** get a patchdata to use for testing */
00035 static bool init_pd( PatchData& pd );
00036 PatchData& ObjectiveFunctionTests::patch()
00037 {
00038     static PatchData the_pd;
00039     static bool did_init = init_pd( the_pd );
00040     CPPUNIT_ASSERT( did_init );
00041     return the_pd;
00042 }
00043 static bool init_pd( PatchData& pd )
00044 {
00045     MsqError err;
00046     const double coords[]  = { 0, 0, 0, 1, 1, 1, 2, 2, 2 };
00047     const bool fixed[]     = { false, true, true };
00048     const size_t indices[] = { 0, 1, 2 };
00049     pd.fill( 3, coords, 1, TRIANGLE, indices, fixed, err );
00050     return !err;
00051 }
00052 
00053 /** Internal helper function for test_eval_type */
00054 double ObjectiveFunctionTests::evaluate_internal( ObjectiveFunction::EvalType type, OFTestMode test_mode,
00055                                                   ObjectiveFunction* of )
00056 {
00057     MsqPrintError err( cout );
00058     vector< Vector3D > grad;
00059     vector< SymMatrix3D > diag;
00060     MsqHessian hess;
00061     bool valid = false;
00062     double result;
00063 
00064     switch( test_mode )
00065     {
00066         case EVAL:
00067             valid = of->evaluate( type, patch(), result, OF_FREE_EVALS_ONLY, err );
00068             break;
00069         case GRAD:
00070             valid = of->evaluate_with_gradient( type, patch(), result, grad, err );
00071             break;
00072         case DIAG:
00073             valid = of->evaluate_with_Hessian_diagonal( type, patch(), result, grad, diag, err );
00074             break;
00075         case HESS:
00076             hess.initialize( patch(), err );
00077             ASSERT_NO_ERROR( err );
00078             valid = of->evaluate_with_Hessian( type, patch(), result, grad, hess, err );
00079             break;
00080         default:
00081             CPPUNIT_ASSERT( false );
00082     }
00083 
00084     ASSERT_NO_ERROR( err );
00085     CPPUNIT_ASSERT( valid );
00086     return result;
00087 }
00088 
00089 void ObjectiveFunctionTests::test_eval_type( ObjectiveFunction::EvalType type, OFTestMode test_mode,
00090                                              ObjectiveFunctionTemplate* of )
00091 {
00092     // define two sets of quality metric values
00093     const double vals1[]   = { 1.0, 3.0, 4.0, 8.0 };
00094     const size_t vals1_len = sizeof( vals1 ) / sizeof( vals1[0] );
00095     const double vals2[]   = { 21.5, 11.1, 30.0, 0.5 };
00096     const size_t vals2_len = sizeof( vals2 ) / sizeof( vals2[0] );
00097 
00098     // create a quality metric to use
00099     OFTestQM metric;
00100     of->set_quality_metric( &metric );
00101 
00102     // get some initial values to compare to
00103     of->clear();
00104     metric.set_values( vals1, vals1_len );
00105     const double init1 = evaluate_internal( ObjectiveFunction::CALCULATE, test_mode, of );
00106     of->clear();
00107     metric.set_values( vals2, vals2_len );
00108     const double init2 = evaluate_internal( ObjectiveFunction::CALCULATE, test_mode, of );
00109     of->clear();
00110     metric.append_values( vals1, vals1_len );
00111     const double inits = evaluate_internal( ObjectiveFunction::CALCULATE, test_mode, of );
00112     of->clear();
00113 
00114     double val, expected;
00115     switch( type )
00116     {
00117         case ObjectiveFunction::CALCULATE:
00118 
00119             // first make sure we get back the same values
00120             metric.set_values( vals1, vals1_len );
00121             val = evaluate_internal( ObjectiveFunction::CALCULATE, test_mode, of );
00122             CPPUNIT_ASSERT_DOUBLES_EQUAL( init1, val, 1e-6 );
00123             metric.set_values( vals2, vals2_len );
00124             val = evaluate_internal( ObjectiveFunction::CALCULATE, test_mode, of );
00125             CPPUNIT_ASSERT_DOUBLES_EQUAL( init2, val, 1e-6 );
00126 
00127             // now do something that should modify the accumulated value of the OF
00128             evaluate_internal( ObjectiveFunction::ACCUMULATE, test_mode, of );
00129 
00130             // check that the values are unchanged
00131             metric.set_values( vals1, vals1_len );
00132             val = evaluate_internal( ObjectiveFunction::CALCULATE, test_mode, of );
00133             CPPUNIT_ASSERT_DOUBLES_EQUAL( init1, val, 1e-6 );
00134             metric.set_values( vals2, vals2_len );
00135             val = evaluate_internal( ObjectiveFunction::CALCULATE, test_mode, of );
00136             CPPUNIT_ASSERT_DOUBLES_EQUAL( init2, val, 1e-6 );
00137 
00138             break;
00139 
00140         case ObjectiveFunction::ACCUMULATE:
00141 
00142             // begin with first set
00143             metric.set_values( vals1, vals1_len );
00144             val = evaluate_internal( ObjectiveFunction::ACCUMULATE, test_mode, of );
00145             CPPUNIT_ASSERT_DOUBLES_EQUAL( init1, val, 1e-6 );
00146 
00147             // add in second set
00148             metric.set_values( vals2, vals2_len );
00149             val = evaluate_internal( ObjectiveFunction::ACCUMULATE, test_mode, of );
00150             CPPUNIT_ASSERT_DOUBLES_EQUAL( inits, val, 1e-6 );
00151 
00152             // clear
00153             of->clear();
00154 
00155             // begin with second set
00156             metric.set_values( vals2, vals2_len );
00157             val = evaluate_internal( ObjectiveFunction::ACCUMULATE, test_mode, of );
00158             CPPUNIT_ASSERT_DOUBLES_EQUAL( init2, val, 1e-6 );
00159 
00160             // add in first set
00161             metric.set_values( vals1, vals1_len );
00162             val = evaluate_internal( ObjectiveFunction::ACCUMULATE, test_mode, of );
00163             CPPUNIT_ASSERT_DOUBLES_EQUAL( inits, val, 1e-6 );
00164 
00165             break;
00166 
00167         case ObjectiveFunction::SAVE:
00168 
00169             // calculate value for first set twice
00170             metric.set_values( vals1, vals1_len );
00171             metric.append_values( vals1, vals1_len );
00172             expected = evaluate_internal( ObjectiveFunction::CALCULATE, test_mode, of );
00173 
00174             // begin with the first set
00175             metric.set_values( vals1, vals1_len );
00176             val = evaluate_internal( ObjectiveFunction::ACCUMULATE, test_mode, of );
00177             CPPUNIT_ASSERT_DOUBLES_EQUAL( init1, val, 1e-6 );
00178 
00179             // add the second set
00180             metric.set_values( vals2, vals2_len );
00181             val = evaluate_internal( ObjectiveFunction::ACCUMULATE, test_mode, of );
00182             CPPUNIT_ASSERT_DOUBLES_EQUAL( inits, val, 1e-6 );
00183 
00184             // now save the second set of values - OF value should not change
00185             metric.set_values( vals2, vals2_len );
00186             val = evaluate_internal( ObjectiveFunction::SAVE, test_mode, of );
00187             CPPUNIT_ASSERT_DOUBLES_EQUAL( inits, val, 1e-6 );
00188 
00189             // now replace the second set with the first
00190             metric.set_values( vals1, vals1_len );
00191             val = evaluate_internal( ObjectiveFunction::UPDATE, test_mode, of );
00192             CPPUNIT_ASSERT_DOUBLES_EQUAL( expected, val, 1e-6 );
00193 
00194             // check that saved values are cleared
00195             of->clear();
00196             metric.set_values( vals2, vals2_len );
00197             val = evaluate_internal( ObjectiveFunction::UPDATE, test_mode, of );
00198             CPPUNIT_ASSERT_DOUBLES_EQUAL( init2, val, 1e-6 );
00199 
00200             break;
00201 
00202         case ObjectiveFunction::UPDATE:
00203 
00204             // With no saved data, an update should produce the same
00205             // result as CALCULATE
00206             metric.set_values( vals1, vals1_len );
00207             val = evaluate_internal( ObjectiveFunction::UPDATE, test_mode, of );
00208             CPPUNIT_ASSERT_DOUBLES_EQUAL( init1, val, 1e-6 );
00209 
00210             // Doing an update with a second patch should change
00211             // the global value to that of the second set
00212             metric.set_values( vals2, vals2_len );
00213             val = evaluate_internal( ObjectiveFunction::UPDATE, test_mode, of );
00214             CPPUNIT_ASSERT_DOUBLES_EQUAL( init2, val, 1e-6 );
00215 
00216             // Now add in the second set again
00217             metric.set_values( vals2, vals2_len );
00218             evaluate_internal( ObjectiveFunction::ACCUMULATE, test_mode, of );
00219 
00220             // Now replace one accumulation of the second set with the first
00221             // Should result in the first set + the second set
00222             metric.set_values( vals1, vals1_len );
00223             val = evaluate_internal( ObjectiveFunction::UPDATE, test_mode, of );
00224             CPPUNIT_ASSERT_DOUBLES_EQUAL( inits, val, 1e-6 );
00225 
00226             break;
00227 
00228         case ObjectiveFunction::TEMPORARY:
00229 
00230             // With no saved data, an TEMPORARY should produce the same
00231             // result as CALCULATE
00232             metric.set_values( vals1, vals1_len );
00233             val = evaluate_internal( ObjectiveFunction::TEMPORARY, test_mode, of );
00234             CPPUNIT_ASSERT_DOUBLES_EQUAL( init1, val, 1e-6 );
00235 
00236             // Begin with two instances of the second set in the
00237             // accumulated value, with one instance saved so it
00238             // can be removed later
00239             metric.set_values( vals2, vals2_len );
00240             evaluate_internal( ObjectiveFunction::ACCUMULATE, test_mode, of );
00241             evaluate_internal( ObjectiveFunction::UPDATE, test_mode, of );
00242 
00243             // Now do a temporary eval, replacing one instance of the
00244             // second set wiht the first set
00245             metric.set_values( vals1, vals1_len );
00246             val = evaluate_internal( ObjectiveFunction::TEMPORARY, test_mode, of );
00247 
00248             // TEMPORARY should produce the same value as UPDATE, but without
00249             // modifying any internal state.
00250             expected = evaluate_internal( ObjectiveFunction::UPDATE, test_mode, of );
00251             CPPUNIT_ASSERT_DOUBLES_EQUAL( expected, val, 1e-6 );
00252 
00253             break;
00254 
00255         default:
00256             CPPUNIT_ASSERT_MESSAGE( "No test for specified evaluation type", false );
00257             break;
00258     }
00259 }
00260 
00261 void ObjectiveFunctionTests::test_value( const double* input_values, unsigned num_input_values, double expected_value,
00262                                          OFTestMode test_mode, ObjectiveFunctionTemplate* of )
00263 {
00264     OFTestQM metric( input_values, num_input_values );
00265     of->set_quality_metric( &metric );
00266 
00267     double val = evaluate_internal( ObjectiveFunction::CALCULATE, test_mode, of );
00268     CPPUNIT_ASSERT_DOUBLES_EQUAL( expected_value, val, 1e-6 );
00269 }
00270 
00271 void ObjectiveFunctionTests::test_clone( ObjectiveFunctionTemplate* of )
00272 {
00273     const double some_vals[] = { 1, 2, 3, 4, 5, 6, 0.5 };
00274     const unsigned num_vals  = sizeof( some_vals ) / sizeof( some_vals[0] );
00275     OFTestQM metric( some_vals, num_vals );
00276     of->set_quality_metric( &metric );
00277 
00278     // test that we get the same value from both
00279     unique_ptr< ObjectiveFunction > of2( of->clone() );
00280     double exp_val, val;
00281     exp_val = evaluate_internal( ObjectiveFunction::CALCULATE, EVAL, of );
00282     val     = evaluate_internal( ObjectiveFunction::CALCULATE, EVAL, of2.get() );
00283     CPPUNIT_ASSERT_DOUBLES_EQUAL( exp_val, val, 1e-12 );
00284 
00285     // check if OF supports BCD -- if not then done
00286     MsqError err;
00287     of->evaluate( ObjectiveFunction::UPDATE, patch(), val, false, err );
00288     if( err )
00289     {
00290         err.clear();
00291         return;
00292     }
00293 
00294     // build up some saved state in the objective function
00295     of->clear();
00296     const double vals1[]   = { 1.0, 3.0, 4.0, 8.0 };
00297     const size_t vals1_len = sizeof( vals1 ) / sizeof( vals1[0] );
00298     const double vals2[]   = { 21.5, 11.1, 30.0, 0.5 };
00299     const size_t vals2_len = sizeof( vals2 ) / sizeof( vals2[0] );
00300     metric.set_values( vals1, vals1_len );
00301     evaluate_internal( ObjectiveFunction::SAVE, EVAL, of );
00302     metric.append_values( vals2, vals2_len );
00303     evaluate_internal( ObjectiveFunction::ACCUMULATE, EVAL, of );
00304 
00305     // check that clone has same accumulated data
00306     of2 = unique_ptr< ObjectiveFunction >( of->clone() );
00307     metric.set_values( some_vals, num_vals );
00308     exp_val = evaluate_internal( ObjectiveFunction::UPDATE, EVAL, of );
00309     val     = evaluate_internal( ObjectiveFunction::UPDATE, EVAL, of2.get() );
00310     CPPUNIT_ASSERT_DOUBLES_EQUAL( exp_val, val, 1e-12 );
00311 }
00312 
00313 /** Test correct handling of QM negate flag */
00314 void ObjectiveFunctionTests::test_negate_flag( OFTestMode test_mode, ObjectiveFunctionTemplate* of )
00315 {
00316     const double some_vals[] = { 1, 2, 3, 4, 5, 6, 0.5 };
00317     const unsigned num_vals  = sizeof( some_vals ) / sizeof( some_vals[0] );
00318     OFTestQM metric( some_vals, num_vals );
00319     of->set_quality_metric( &metric );
00320 
00321     MsqPrintError err( cout );
00322     bool rval = false;
00323     double value[2];
00324     vector< Vector3D > grad[2];
00325     vector< SymMatrix3D > diag[2];
00326     MsqHessian hess[2];
00327 
00328     // Do twice, once w/out negate flag set and then once
00329     // with negate flag == -1.
00330     ObjectiveFunction::EvalType type = ObjectiveFunction::CALCULATE;
00331     for( unsigned i = 0; i < 2; ++i )
00332     {
00333         switch( test_mode )
00334         {
00335             case EVAL:
00336                 rval = of->evaluate( type, patch(), value[i], false, err );
00337                 break;
00338             case GRAD:
00339                 rval = of->evaluate_with_gradient( type, patch(), value[i], grad[i], err );
00340                 break;
00341             case DIAG:
00342                 rval = of->evaluate_with_Hessian_diagonal( type, patch(), value[i], grad[i], diag[i], err );
00343                 break;
00344             case HESS:
00345                 hess[i].initialize( patch(), err );
00346                 ASSERT_NO_ERROR( err );
00347                 rval = of->evaluate_with_Hessian( type, patch(), value[i], grad[i], hess[i], err );
00348                 break;
00349             default:
00350                 CPPUNIT_ASSERT_MESSAGE( "Invalid enum value in test code", false );
00351                 break;
00352         }
00353         ASSERT_NO_ERROR( err );
00354         CPPUNIT_ASSERT( rval );
00355         metric.set_negate_flag( -1 );
00356     }
00357 
00358     switch( test_mode )
00359     {
00360         case HESS:
00361             CPPUNIT_ASSERT_EQUAL( hess[0].size(), hess[1].size() );
00362             for( size_t r = 0; r < hess[0].size(); ++r )
00363                 for( size_t c = r; c < hess[0].size(); ++c )
00364                     if( hess[0].get_block( r, c ) )
00365                         CPPUNIT_ASSERT_MATRICES_EQUAL( -*hess[0].get_block( r, c ), *hess[1].get_block( r, c ), 1e-6 );
00366         case DIAG:
00367             // NOTE: When case HESS: falls through to here, diag[0] and diag[1]
00368             // will be empty, making this a no-op.
00369             CPPUNIT_ASSERT_EQUAL( diag[0].size(), diag[1].size() );
00370             for( size_t j = 0; j < diag[0].size(); ++j )
00371                 CPPUNIT_ASSERT_MATRICES_EQUAL( -diag[0][j], diag[1][j], 1e-6 );
00372         case GRAD:
00373             CPPUNIT_ASSERT_EQUAL( grad[0].size(), grad[1].size() );
00374             for( size_t j = 0; j < grad[0].size(); ++j )
00375                 CPPUNIT_ASSERT_VECTORS_EQUAL( -grad[0][j], grad[1][j], 1e-6 );
00376         default:
00377             CPPUNIT_ASSERT_DOUBLES_EQUAL( -value[0], value[1], 1e-6 );
00378     }
00379 }
00380 
00381 void ObjectiveFunctionTests::compare_numerical_gradient( ObjectiveFunction* of )
00382 {
00383     MsqPrintError err( std::cout );
00384     PatchData pd;
00385     create_twelve_hex_patch( pd, err );
00386     ASSERT_NO_ERROR( err );
00387 
00388     std::vector< Vector3D > num_grad, ana_grad;
00389     double num_val, ana_val;
00390     bool valid;
00391 
00392     valid = of->evaluate_with_gradient( ObjectiveFunction::CALCULATE, pd, ana_val, ana_grad, err );
00393     ASSERT_NO_ERROR( err );
00394     CPPUNIT_ASSERT( valid );
00395     CPPUNIT_ASSERT_EQUAL( pd.num_free_vertices(), ana_grad.size() );
00396 
00397     valid = of->ObjectiveFunction::evaluate_with_gradient( ObjectiveFunction::CALCULATE, pd, num_val, num_grad, err );
00398     ASSERT_NO_ERROR( err );
00399     CPPUNIT_ASSERT( valid );
00400     CPPUNIT_ASSERT_EQUAL( pd.num_free_vertices(), num_grad.size() );
00401 
00402     CPPUNIT_ASSERT_DOUBLES_EQUAL( ana_val, num_val, 1e-6 );
00403     for( size_t i = 0; i < pd.num_free_vertices(); ++i )
00404     {
00405         CPPUNIT_ASSERT_VECTORS_EQUAL( num_grad[i], ana_grad[i], 1e-3 );
00406     }
00407 }
00408 
00409 void ObjectiveFunctionTests::compare_hessian_gradient( ObjectiveFunction* of )
00410 {
00411     MsqPrintError err( std::cout );
00412     PatchData pd;
00413     create_twelve_hex_patch( pd, err );
00414     ASSERT_NO_ERROR( err );
00415 
00416     std::vector< Vector3D > grad, hess_grad;
00417     MsqHessian hess;
00418     double grad_val, hess_val;
00419     bool valid;
00420 
00421     valid = of->evaluate_with_gradient( ObjectiveFunction::CALCULATE, pd, grad_val, grad, err );
00422     ASSERT_NO_ERROR( err );
00423     CPPUNIT_ASSERT( valid );
00424     CPPUNIT_ASSERT_EQUAL( pd.num_free_vertices(), grad.size() );
00425 
00426     hess.initialize( pd, err );
00427     ASSERT_NO_ERROR( err );
00428     valid = of->evaluate_with_Hessian( ObjectiveFunction::CALCULATE, pd, hess_val, hess_grad, hess, err );
00429     ASSERT_NO_ERROR( err );
00430     CPPUNIT_ASSERT( valid );
00431     CPPUNIT_ASSERT_EQUAL( pd.num_free_vertices(), hess_grad.size() );
00432 
00433     CPPUNIT_ASSERT_DOUBLES_EQUAL( grad_val, hess_val, 1e-6 );
00434     for( size_t i = 0; i < pd.num_free_vertices(); ++i )
00435     {
00436         CPPUNIT_ASSERT_VECTORS_EQUAL( grad[i], hess_grad[i], 1e-6 );
00437     }
00438 }
00439 
00440 void ObjectiveFunctionTests::compare_diagonal_gradient( ObjectiveFunction* of )
00441 {
00442     MsqPrintError err( std::cout );
00443     PatchData pd;
00444     create_twelve_hex_patch( pd, err );
00445     ASSERT_NO_ERROR( err );
00446 
00447     std::vector< Vector3D > grad, hess_grad;
00448     std::vector< SymMatrix3D > hess;
00449     double grad_val, hess_val;
00450     bool valid;
00451 
00452     valid = of->evaluate_with_gradient( ObjectiveFunction::CALCULATE, pd, grad_val, grad, err );
00453     ASSERT_NO_ERROR( err );
00454     CPPUNIT_ASSERT( valid );
00455     CPPUNIT_ASSERT_EQUAL( pd.num_free_vertices(), grad.size() );
00456 
00457     valid = of->evaluate_with_Hessian_diagonal( ObjectiveFunction::CALCULATE, pd, hess_val, hess_grad, hess, err );
00458     ASSERT_NO_ERROR( err );
00459     CPPUNIT_ASSERT( valid );
00460     CPPUNIT_ASSERT_EQUAL( pd.num_free_vertices(), hess_grad.size() );
00461 
00462     CPPUNIT_ASSERT_DOUBLES_EQUAL( grad_val, hess_val, 1e-6 );
00463     for( size_t i = 0; i < pd.num_free_vertices(); ++i )
00464     {
00465         CPPUNIT_ASSERT_VECTORS_EQUAL( grad[i], hess_grad[i], 1e-6 );
00466     }
00467 }
00468 
00469 void ObjectiveFunctionTests::compare_hessian_diagonal( ObjectiveFunction* of )
00470 {
00471     MsqPrintError err( std::cout );
00472     PatchData pd;
00473     create_twelve_hex_patch( pd, err );
00474     ASSERT_NO_ERROR( err );
00475 
00476     std::vector< Vector3D > diag_grad, hess_grad;
00477     std::vector< SymMatrix3D > diag;
00478     MsqHessian hess;
00479     double diag_val, hess_val;
00480     bool valid;
00481 
00482     valid = of->evaluate_with_Hessian_diagonal( ObjectiveFunction::CALCULATE, pd, diag_val, diag_grad, diag, err );
00483     ASSERT_NO_ERROR( err );
00484     CPPUNIT_ASSERT( valid );
00485     CPPUNIT_ASSERT_EQUAL( pd.num_free_vertices(), diag_grad.size() );
00486     CPPUNIT_ASSERT_EQUAL( pd.num_free_vertices(), diag.size() );
00487 
00488     hess.initialize( pd, err );
00489     ASSERT_NO_ERROR( err );
00490     valid = of->evaluate_with_Hessian( ObjectiveFunction::CALCULATE, pd, hess_val, hess_grad, hess, err );
00491     ASSERT_NO_ERROR( err );
00492     CPPUNIT_ASSERT( valid );
00493     CPPUNIT_ASSERT_EQUAL( pd.num_free_vertices(), hess_grad.size() );
00494     CPPUNIT_ASSERT_EQUAL( pd.num_free_vertices(), hess.size() );
00495 
00496     CPPUNIT_ASSERT_DOUBLES_EQUAL( hess_val, diag_val, 1e-6 );
00497     for( size_t i = 0; i < pd.num_free_vertices(); ++i )
00498     {
00499         CPPUNIT_ASSERT_VECTORS_EQUAL( hess_grad[i], diag_grad[i], 1e-6 );
00500         CPPUNIT_ASSERT_MATRICES_EQUAL( *hess.get_block( i, i ), diag[i], 1e-6 );
00501     }
00502 }
00503 
00504 double MAT_EPS( const Matrix3D& A )
00505 {
00506     return std::max( 0.001 * sqrt( Frobenius_2( A ) ), 0.001 );
00507 }
00508 
00509 #define CHECK_EQUAL_MATRICES( A, B ) CPPUNIT_ASSERT_MATRICES_EQUAL( ( A ), ( B ), MAT_EPS( A ) )
00510 
00511 void ObjectiveFunctionTests::compare_numerical_hessian( ObjectiveFunction* of, bool diagonal_only )
00512 {
00513     const double delta = 0.0001;
00514 
00515     MsqPrintError err( std::cout );
00516     PatchData pd;
00517     create_qm_two_tet_patch( pd, err );
00518     ASSERT_NO_ERROR( err );
00519     CPPUNIT_ASSERT( pd.num_free_vertices() != 0 );
00520 
00521     // get analytical Hessian from objective function
00522     std::vector< Vector3D > grad;
00523     std::vector< SymMatrix3D > diag;
00524     MsqHessian hess;
00525     hess.initialize( pd, err );
00526     ASSERT_NO_ERROR( err );
00527     double value;
00528     bool valid;
00529     if( diagonal_only )
00530         valid = of->evaluate_with_Hessian_diagonal( ObjectiveFunction::CALCULATE, pd, value, grad, diag, err );
00531     else
00532         valid = of->evaluate_with_Hessian( ObjectiveFunction::CALCULATE, pd, value, grad, hess, err );
00533     ASSERT_NO_ERROR( err );
00534     CPPUNIT_ASSERT( valid );
00535 
00536     // do numerical approximation of each block and compare to analytical value
00537     for( size_t i = 0; i < pd.num_free_vertices(); ++i )
00538     {
00539         const size_t j_end = diagonal_only ? i + 1 : pd.num_free_vertices();
00540         for( size_t j = i; j < j_end; ++j )
00541         {
00542             // do numerical approximation for block corresponding to
00543             // coorindates for ith and jth vertices.
00544             Matrix3D block;
00545             for( int k = 0; k < 3; ++k )
00546             {
00547                 for( int m = 0; m < 3; ++m )
00548                 {
00549                     double dk, dm, dkm;
00550                     Vector3D ik = pd.vertex_by_index( i );
00551                     Vector3D im = pd.vertex_by_index( j );
00552 
00553                     Vector3D delta_k( 0.0 );
00554                     delta_k[k] = delta;
00555                     pd.move_vertex( delta_k, i, err );
00556                     ASSERT_NO_ERROR( err );
00557                     valid = of->evaluate( ObjectiveFunction::CALCULATE, pd, dk, true, err );
00558                     ASSERT_NO_ERROR( err );
00559                     CPPUNIT_ASSERT( valid );
00560 
00561                     Vector3D delta_m( 0.0 );
00562                     delta_m[m] = delta;
00563                     pd.move_vertex( delta_m, j, err );
00564                     ASSERT_NO_ERROR( err );
00565                     valid = of->evaluate( ObjectiveFunction::CALCULATE, pd, dkm, true, err );
00566                     ASSERT_NO_ERROR( err );
00567                     CPPUNIT_ASSERT( valid );
00568 
00569                     // be careful here that we do the right thing if i==j
00570                     pd.set_vertex_coordinates( ik, i, err );
00571                     ASSERT_NO_ERROR( err );
00572                     pd.set_vertex_coordinates( im, j, err );
00573                     ASSERT_NO_ERROR( err );
00574                     pd.move_vertex( delta_m, j, err );
00575                     ASSERT_NO_ERROR( err );
00576                     valid = of->evaluate( ObjectiveFunction::CALCULATE, pd, dm, true, err );
00577                     ASSERT_NO_ERROR( err );
00578                     CPPUNIT_ASSERT( valid );
00579 
00580                     pd.set_vertex_coordinates( ik, i, err );
00581                     ASSERT_NO_ERROR( err );
00582                     pd.set_vertex_coordinates( im, j, err );
00583                     ASSERT_NO_ERROR( err );
00584 
00585                     block[k][m] = ( dkm - dk - dm + value ) / ( delta * delta );
00586                 }
00587             }
00588             // compare to analytical value
00589             if( diagonal_only )
00590             {
00591                 CPPUNIT_ASSERT( i == j );  // see j_end above
00592                 CPPUNIT_ASSERT( i < diag.size() );
00593                 CHECK_EQUAL_MATRICES( block, Matrix3D( diag[i] ) );
00594             }
00595             else
00596             {
00597                 Matrix3D* m  = hess.get_block( i, j );
00598                 Matrix3D* mt = hess.get_block( j, i );
00599                 if( NULL != m ) { CHECK_EQUAL_MATRICES( block, *m ); }
00600                 if( NULL != mt ) { CHECK_EQUAL_MATRICES( transpose( block ), *m ); }
00601                 if( NULL == mt && NULL == m ) { CHECK_EQUAL_MATRICES( Matrix3D( 0.0 ), block ); }
00602             }
00603         }
00604     }
00605 }
00606 
00607 const size_t HANDLE_VAL = 0xDEADBEEF;
00608 class OFTestBadQM : public QualityMetric
00609 {
00610   public:
00611     OFTestBadQM( bool return_error ) : mError( return_error ) {}
00612 
00613     virtual MetricType get_metric_type() const
00614     {
00615         return ELEMENT_BASED;
00616     }
00617 
00618     virtual string get_name() const
00619     {
00620         return "ObjectiveFunctionTests";
00621     }
00622 
00623     virtual int get_negate_flag() const
00624     {
00625         return 1;
00626     }
00627 
00628     virtual void get_evaluations( PatchData&, vector< size_t >& h, bool, MsqError& )
00629     {
00630         h.clear();
00631         h.push_back( HANDLE_VAL );
00632     }
00633 
00634     virtual bool evaluate( PatchData&, size_t h, double&, MsqError& err )
00635     {
00636         CPPUNIT_ASSERT_EQUAL( HANDLE_VAL, h );
00637         if( mError ) MSQ_SETERR( err )( MsqError::INVALID_STATE );
00638         return false;
00639     }
00640 
00641     virtual bool evaluate_with_indices( PatchData&, size_t h, double&, vector< size_t >&, MsqError& err )
00642     {
00643         CPPUNIT_ASSERT_EQUAL( HANDLE_VAL, h );
00644         if( mError ) MSQ_SETERR( err )( MsqError::INVALID_STATE );
00645         return false;
00646     }
00647 
00648     virtual bool evaluate_with_gradient( PatchData&, size_t h, double&, vector< size_t >&, vector< Vector3D >&,
00649                                          MsqError& err )
00650     {
00651         CPPUNIT_ASSERT_EQUAL( HANDLE_VAL, h );
00652         if( mError ) MSQ_SETERR( err )( MsqError::INVALID_STATE );
00653         return false;
00654     }
00655 
00656     virtual bool evaluate_with_Hessian( PatchData&, size_t h, double&, vector< size_t >&, vector< Vector3D >&,
00657                                         vector< Matrix3D >&, MsqError& err )
00658     {
00659         CPPUNIT_ASSERT_EQUAL( HANDLE_VAL, h );
00660         if( mError ) MSQ_SETERR( err )( MsqError::INVALID_STATE );
00661         return false;
00662     }
00663 
00664   private:
00665     bool mError;
00666 };
00667 
00668 void ObjectiveFunctionTests::test_handles_invalid_qm( OFTestMode test_mode, ObjectiveFunctionTemplate* of )
00669 {
00670     OFTestBadQM metric( false );
00671     of->set_quality_metric( &metric );
00672 
00673     MsqPrintError err( cout );
00674     vector< Vector3D > grad;
00675     vector< SymMatrix3D > diag;
00676     MsqHessian hess;
00677     double result;
00678     bool valid = false;
00679 
00680     switch( test_mode )
00681     {
00682         case EVAL:
00683             valid = of->evaluate( ObjectiveFunction::CALCULATE, patch(), result, OF_FREE_EVALS_ONLY, err );
00684             break;
00685         case GRAD:
00686             valid = of->evaluate_with_gradient( ObjectiveFunction::CALCULATE, patch(), result, grad, err );
00687             break;
00688         case DIAG:
00689             valid =
00690                 of->evaluate_with_Hessian_diagonal( ObjectiveFunction::CALCULATE, patch(), result, grad, diag, err );
00691             break;
00692         case HESS:
00693             hess.initialize( patch(), err );
00694             ASSERT_NO_ERROR( err );
00695             valid = of->evaluate_with_Hessian( ObjectiveFunction::CALCULATE, patch(), result, grad, hess, err );
00696             break;
00697         default:
00698             CPPUNIT_ASSERT( false );
00699     }
00700 
00701     ASSERT_NO_ERROR( err );
00702     CPPUNIT_ASSERT( !valid );
00703 }
00704 
00705 void ObjectiveFunctionTests::test_handles_qm_error( OFTestMode test_mode, ObjectiveFunctionTemplate* of )
00706 
00707 {
00708     OFTestBadQM metric( true );
00709     of->set_quality_metric( &metric );
00710 
00711     MsqError err;
00712     vector< Vector3D > grad;
00713     vector< SymMatrix3D > diag;
00714     MsqHessian hess;
00715     double result;
00716     bool valid;
00717 
00718     switch( test_mode )
00719     {
00720         case EVAL:
00721             valid = of->evaluate( ObjectiveFunction::CALCULATE, patch(), result, OF_FREE_EVALS_ONLY, err );
00722             break;
00723         case GRAD:
00724             valid = of->evaluate_with_gradient( ObjectiveFunction::CALCULATE, patch(), result, grad, err );
00725             break;
00726         case DIAG:
00727             valid =
00728                 of->evaluate_with_Hessian_diagonal( ObjectiveFunction::CALCULATE, patch(), result, grad, diag, err );
00729             break;
00730         case HESS:
00731             hess.initialize( patch(), err );
00732             ASSERT_NO_ERROR( err );
00733             valid = of->evaluate_with_Hessian( ObjectiveFunction::CALCULATE, patch(), result, grad, hess, err );
00734             break;
00735         default:
00736             CPPUNIT_ASSERT( false );
00737     }
00738 
00739     CPPUNIT_ASSERT( err );
00740 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines