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