MOAB: Mesh Oriented datABase  (version 5.3.0)
PMeanPTemplateTest.cpp
Go to the documentation of this file.
00001 /* *****************************************************************
00002     MESQUITE -- The Mesh Quality Improvement Toolkit
00003 
00004     Copyright 2006 Lawrence Livermore National Laboratory.  Under
00005     the terms of Contract B545069 with the University of Wisconsin --
00006     Madison, Lawrence Livermore National Laboratory retains certain
00007     rights in 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     (2006) kraftche@cae.wisc.edu
00024 
00025   ***************************************************************** */
00026 
00027 /** \file PMeanPTemplateTest.cpp
00028  *  \brief previous name: PowerMeanPTest.cpp
00029  *  \author Jason Kraftcheck
00030  */
00031 
00032 #include "Mesquite.hpp"
00033 #include "PMeanPTemplate.hpp"
00034 #include "MsqError.hpp"
00035 #include "PatchData.hpp"
00036 #include "UnitUtil.hpp"
00037 #include "VertexQM.hpp"
00038 #include "MsqHessian.hpp"
00039 
00040 #include "ObjectiveFunctionTests.hpp"
00041 
00042 using namespace MBMesquite;
00043 using namespace std;
00044 
00045 const double EPSILON = 1e-4;
00046 
00047 class PMeanPTemplateTest : public CppUnit::TestFixture, public ObjectiveFunctionTests
00048 {
00049   private:
00050     void test_eval_type( OFTestMode eval_func, ObjectiveFunction::EvalType type );
00051     void test_evaluate( double power );
00052     void test_gradient( double power );
00053     void test_diagonal( double power );
00054     void test_Hessian( double power );
00055 
00056     void check_result( PatchData& pd, double power, double value, Vector3D* gradient = 0, Matrix3D* Hessian = 0 );
00057 
00058     CPPUNIT_TEST_SUITE( PMeanPTemplateTest );
00059 
00060     CPPUNIT_TEST( test_eval_calc );
00061     CPPUNIT_TEST( test_eval_accum );
00062     CPPUNIT_TEST( test_eval_save );
00063     CPPUNIT_TEST( test_eval_update );
00064     CPPUNIT_TEST( test_eval_temp );
00065 
00066     CPPUNIT_TEST( test_grad_calc );
00067     CPPUNIT_TEST( test_grad_save );
00068     CPPUNIT_TEST( test_grad_update );
00069     CPPUNIT_TEST( test_grad_temp );
00070 
00071     CPPUNIT_TEST( test_diag_calc );
00072     CPPUNIT_TEST( test_diag_save );
00073     CPPUNIT_TEST( test_diag_update );
00074     CPPUNIT_TEST( test_diag_temp );
00075 
00076     CPPUNIT_TEST( test_Hess_calc );
00077     CPPUNIT_TEST( test_Hess_save );
00078     CPPUNIT_TEST( test_Hess_update );
00079     CPPUNIT_TEST( test_Hess_temp );
00080 
00081     CPPUNIT_TEST( test_clone );
00082 
00083     CPPUNIT_TEST( test_failed_metric_in_eval );
00084     CPPUNIT_TEST( test_failed_metric_in_grad );
00085     CPPUNIT_TEST( test_failed_metric_in_diag );
00086     CPPUNIT_TEST( test_failed_metric_in_Hess );
00087 
00088     CPPUNIT_TEST( test_false_metric_in_eval );
00089     CPPUNIT_TEST( test_false_metric_in_grad );
00090     CPPUNIT_TEST( test_false_metric_in_diag );
00091     CPPUNIT_TEST( test_false_metric_in_Hess );
00092 
00093     CPPUNIT_TEST( test_evaluate_arithmatic );
00094     CPPUNIT_TEST( test_evaluate_rms );
00095 
00096     CPPUNIT_TEST( test_gradient_arithmatic );
00097     CPPUNIT_TEST( test_gradient_rms );
00098 
00099     CPPUNIT_TEST( test_Hessian_arithmatic );
00100     CPPUNIT_TEST( test_Hessian_rms );
00101 
00102     CPPUNIT_TEST( compare_gradient_arithmatic );
00103     CPPUNIT_TEST( compare_gradient_rms );
00104 
00105     CPPUNIT_TEST( compare_diagonal_gradient_arithmatic );
00106     CPPUNIT_TEST( compare_diagonal_gradient_rms );
00107 
00108     CPPUNIT_TEST( compare_hessian_gradient_arithmatic );
00109     CPPUNIT_TEST( compare_hessian_gradient_rms );
00110 
00111     CPPUNIT_TEST( compare_hessian_diagonal_arithmatic );
00112     CPPUNIT_TEST( compare_hessian_diagonal_rms );
00113 
00114     CPPUNIT_TEST( compare_hessian_arithmatic );
00115     CPPUNIT_TEST( compare_hessian_rms );
00116 
00117     CPPUNIT_TEST( compare_hessian_diag_arithmatic );
00118     CPPUNIT_TEST( compare_hessian_diag_rms );
00119 
00120     CPPUNIT_TEST( test_negate_eval );
00121     CPPUNIT_TEST( test_negate_grad );
00122     CPPUNIT_TEST( test_negate_diag );
00123     CPPUNIT_TEST( test_negate_hess );
00124 
00125     CPPUNIT_TEST_SUITE_END();
00126     PatchData mPatch;
00127 
00128   public:
00129     void setUp();
00130 
00131     void test_eval_calc()
00132     {
00133         test_eval_type( EVAL, ObjectiveFunction::CALCULATE );
00134     }
00135     void test_eval_accum()
00136     {
00137         test_eval_type( EVAL, ObjectiveFunction::ACCUMULATE );
00138     }
00139     void test_eval_save()
00140     {
00141         test_eval_type( EVAL, ObjectiveFunction::SAVE );
00142     }
00143     void test_eval_update()
00144     {
00145         test_eval_type( EVAL, ObjectiveFunction::UPDATE );
00146     }
00147     void test_eval_temp()
00148     {
00149         test_eval_type( EVAL, ObjectiveFunction::TEMPORARY );
00150     }
00151 
00152     void test_grad_calc()
00153     {
00154         test_eval_type( GRAD, ObjectiveFunction::CALCULATE );
00155     }
00156     void test_grad_save()
00157     {
00158         test_eval_type( GRAD, ObjectiveFunction::SAVE );
00159     }
00160     void test_grad_update()
00161     {
00162         test_eval_type( GRAD, ObjectiveFunction::UPDATE );
00163     }
00164     void test_grad_temp()
00165     {
00166         test_eval_type( GRAD, ObjectiveFunction::TEMPORARY );
00167     }
00168 
00169     void test_diag_calc()
00170     {
00171         test_eval_type( DIAG, ObjectiveFunction::CALCULATE );
00172     }
00173     void test_diag_save()
00174     {
00175         test_eval_type( DIAG, ObjectiveFunction::SAVE );
00176     }
00177     void test_diag_update()
00178     {
00179         test_eval_type( DIAG, ObjectiveFunction::UPDATE );
00180     }
00181     void test_diag_temp()
00182     {
00183         test_eval_type( DIAG, ObjectiveFunction::TEMPORARY );
00184     }
00185 
00186     void test_Hess_calc()
00187     {
00188         test_eval_type( HESS, ObjectiveFunction::CALCULATE );
00189     }
00190     void test_Hess_save()
00191     {
00192         test_eval_type( HESS, ObjectiveFunction::SAVE );
00193     }
00194     void test_Hess_update()
00195     {
00196         test_eval_type( HESS, ObjectiveFunction::UPDATE );
00197     }
00198     void test_Hess_temp()
00199     {
00200         test_eval_type( HESS, ObjectiveFunction::TEMPORARY );
00201     }
00202 
00203     void test_clone()
00204     {
00205         PMeanPTemplate of( 1, NULL );
00206         ObjectiveFunctionTests::test_clone( &of );
00207     }
00208 
00209     void test_failed_metric_in_eval()
00210     {
00211         PMeanPTemplate of( 1, NULL );
00212         test_handles_qm_error( EVAL, &of );
00213     }
00214     void test_failed_metric_in_grad()
00215     {
00216         PMeanPTemplate of( 1, NULL );
00217         test_handles_qm_error( GRAD, &of );
00218     }
00219     void test_failed_metric_in_diag()
00220     {
00221         PMeanPTemplate of( 1, NULL );
00222         test_handles_qm_error( DIAG, &of );
00223     }
00224     void test_failed_metric_in_Hess()
00225     {
00226         PMeanPTemplate of( 1, NULL );
00227         test_handles_qm_error( HESS, &of );
00228     }
00229 
00230     void test_false_metric_in_eval()
00231     {
00232         PMeanPTemplate of( 1, NULL );
00233         test_handles_invalid_qm( EVAL, &of );
00234     }
00235     void test_false_metric_in_grad()
00236     {
00237         PMeanPTemplate of( 1, NULL );
00238         test_handles_invalid_qm( GRAD, &of );
00239     }
00240     void test_false_metric_in_diag()
00241     {
00242         PMeanPTemplate of( 1, NULL );
00243         test_handles_invalid_qm( DIAG, &of );
00244     }
00245     void test_false_metric_in_Hess()
00246     {
00247         PMeanPTemplate of( 1, NULL );
00248         test_handles_invalid_qm( HESS, &of );
00249     }
00250 
00251     void test_evaluate_arithmatic()
00252     {
00253         test_evaluate( 1 );
00254     }
00255     void test_evaluate_rms()
00256     {
00257         test_evaluate( 2 );
00258     }
00259 
00260     void test_gradient_arithmatic()
00261     {
00262         test_gradient( 1 );
00263     }
00264     void test_gradient_rms()
00265     {
00266         test_gradient( 2 );
00267     }
00268 
00269     void test_diagonal_arithmatic()
00270     {
00271         test_diagonal( 1 );
00272     }
00273     void test_diagonal_rms()
00274     {
00275         test_diagonal( 2 );
00276     }
00277 
00278     void test_Hessian_arithmatic()
00279     {
00280         test_Hessian( 1 );
00281     }
00282     void test_Hessian_rms()
00283     {
00284         test_Hessian( 2 );
00285     }
00286 
00287     void compare_gradient_arithmatic()
00288     {
00289         PMeanPTemplate of( 1, NULL );
00290         compare_numerical_gradient( &of );
00291     }
00292     void compare_gradient_rms()
00293     {
00294         PMeanPTemplate of( 2, NULL );
00295         compare_numerical_gradient( &of );
00296     }
00297 
00298     void compare_diagonal_gradient_arithmatic()
00299     {
00300         PMeanPTemplate of( 1, NULL );
00301         compare_diagonal_gradient( &of );
00302     }
00303     void compare_diagonal_gradient_rms()
00304     {
00305         PMeanPTemplate of( 2, NULL );
00306         compare_diagonal_gradient( &of );
00307     }
00308 
00309     void compare_hessian_gradient_arithmatic()
00310     {
00311         PMeanPTemplate of( 1, NULL );
00312         compare_hessian_gradient( &of );
00313     }
00314     void compare_hessian_gradient_rms()
00315     {
00316         PMeanPTemplate of( 2, NULL );
00317         compare_hessian_gradient( &of );
00318     }
00319 
00320     void compare_hessian_diagonal_arithmatic()
00321     {
00322         PMeanPTemplate of( 1, NULL );
00323         compare_hessian_diagonal( &of );
00324     }
00325     void compare_hessian_diagonal_rms()
00326     {
00327         PMeanPTemplate of( 2, NULL );
00328         compare_hessian_diagonal( &of );
00329     }
00330 
00331     void compare_hessian_arithmatic()
00332     {
00333         PMeanPTemplate of( 1, NULL );
00334         compare_numerical_hessian( &of );
00335     }
00336     void compare_hessian_rms()
00337     {
00338         PMeanPTemplate of( 2, NULL );
00339         compare_numerical_hessian( &of );
00340     }
00341 
00342     void compare_hessian_diag_arithmatic()
00343     {
00344         PMeanPTemplate of( 1, NULL );
00345         compare_numerical_hessian_diagonal( &of );
00346     }
00347     void compare_hessian_diag_rms()
00348     {
00349         PMeanPTemplate of( 2, NULL );
00350         compare_numerical_hessian_diagonal( &of );
00351     }
00352 
00353     void test_negate_eval()
00354     {
00355         PMeanPTemplate of( 2, NULL );
00356         test_negate_flag( EVAL, &of );
00357     }
00358     void test_negate_grad()
00359     {
00360         PMeanPTemplate of( 2, NULL );
00361         test_negate_flag( GRAD, &of );
00362     }
00363     void test_negate_diag()
00364     {
00365         PMeanPTemplate of( 2, NULL );
00366         test_negate_flag( DIAG, &of );
00367     }
00368     void test_negate_hess()
00369     {
00370         PMeanPTemplate of( 2, NULL );
00371         test_negate_flag( HESS, &of );
00372     }
00373 };
00374 
00375 CPPUNIT_TEST_SUITE_NAMED_REGISTRATION( PMeanPTemplateTest, "PMeanPTemplateTest" );
00376 CPPUNIT_TEST_SUITE_NAMED_REGISTRATION( PMeanPTemplateTest, "Unit" );
00377 
00378 void PMeanPTemplateTest::setUp()
00379 {
00380     MsqPrintError err( std::cout );
00381 
00382     // Create a triangle mesh with three free vertices
00383     const double coords[]   = { 1, 1, 0, 2, 1, 0, 3, 1, 0, 2, 2, 0, 1, 3, 0, 1, 2, 0 };
00384     const bool fixed_vtx[]  = { true, false, true, false, true, false };
00385     const size_t tri_conn[] = { 0, 1, 5, 1, 2, 3, 3, 4, 5, 1, 3, 5 };
00386     mPatch.fill( 6, coords, 4, TRIANGLE, tri_conn, fixed_vtx, err );
00387     CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
00388 }
00389 
00390 /** Define a fake quality metric for testing the objective function
00391  *
00392  * Returns, for each vertex in the patch, the distance of that
00393  * vertex from the origin as the quality.  Each evaluation depends
00394  * only on a single vertex.
00395  */
00396 class DistTestMetric : public VertexQM
00397 {
00398   public:
00399     DistTestMetric() : falseEval( false ), failEval( false ) {}
00400     string get_name() const
00401     {
00402         return "Fake metric for testing objective function";
00403     }
00404     int get_negate_flag() const
00405     {
00406         return 1;
00407     }
00408     bool evaluate( PatchData& pd, size_t vtx_idx, double& value, MsqError& err );
00409     bool evaluate_with_indices( PatchData& pd, size_t vtx_idx, double& value, vector< size_t >& indices,
00410                                 MsqError& err );
00411     // bool evaluate_with_gradient( PatchData& pd, size_t vtx_idx, double &value, vector<size_t>&
00412     // indices, vector<Vector3D>& grad, MsqError& err );
00413     bool falseEval;
00414     bool failEval;
00415 };
00416 
00417 bool DistTestMetric::evaluate( PatchData& pd, size_t vtx_idx, double& value, MsqError& err )
00418 {
00419     if( failEval )
00420     {
00421         MSQ_SETERR( err )( MsqError::INVALID_STATE );
00422         return true;
00423     }
00424 
00425     const MsqVertex& vtx = pd.vertex_by_index( vtx_idx );
00426     value                = vtx.length_squared();
00427     return !falseEval;
00428 }
00429 
00430 bool DistTestMetric::evaluate_with_indices( PatchData& pd, size_t vtx_idx, double& value, vector< size_t >& indices,
00431                                             MsqError& err )
00432 {
00433     if( failEval )
00434     {
00435         MSQ_SETERR( err )( MsqError::INVALID_STATE );
00436         return true;
00437     }
00438 
00439     indices.clear();
00440     if( vtx_idx < pd.num_free_vertices() ) indices.push_back( vtx_idx );
00441 
00442     const MsqVertex& vtx = pd.vertex_by_index( vtx_idx );
00443     value                = vtx.length_squared();
00444     return !falseEval;
00445 }
00446 
00447 void PMeanPTemplateTest::test_eval_type( OFTestMode eval_func, ObjectiveFunction::EvalType type )
00448 {
00449     PMeanPTemplate func( 1, NULL );
00450     ObjectiveFunctionTests::test_eval_type( type, eval_func, &func );
00451 }
00452 
00453 void PMeanPTemplateTest::test_evaluate( double power )
00454 {
00455     MsqPrintError err( cout );
00456     double value;
00457     bool rval;
00458 
00459     DistTestMetric metric;
00460     PMeanPTemplate func( power, &metric );
00461     rval = func.evaluate( ObjectiveFunction::CALCULATE, mPatch, value, OF_FREE_EVALS_ONLY, err );
00462     CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
00463     CPPUNIT_ASSERT( rval );
00464 
00465     check_result( mPatch, power, value );
00466 }
00467 
00468 void PMeanPTemplateTest::test_gradient( double power )
00469 {
00470     MsqPrintError err( cout );
00471     double value;
00472     bool rval;
00473     vector< Vector3D > grad;
00474 
00475     DistTestMetric metric;
00476     PMeanPTemplate func( power, &metric );
00477     rval = func.evaluate_with_gradient( ObjectiveFunction::CALCULATE, mPatch, value, grad, err );
00478     CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
00479     CPPUNIT_ASSERT( rval );
00480     CPPUNIT_ASSERT_EQUAL( mPatch.num_free_vertices(), grad.size() );
00481 
00482     if( !grad.empty() ) check_result( mPatch, power, value, arrptr( grad ) );
00483 }
00484 
00485 void PMeanPTemplateTest::test_diagonal( double power )
00486 {
00487     MsqPrintError err( cout );
00488     double value;
00489     bool rval;
00490     vector< Vector3D > grad;
00491     vector< SymMatrix3D > Hess;
00492     CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
00493 
00494     DistTestMetric metric;
00495     PMeanPTemplate func( power, &metric );
00496     rval = func.evaluate_with_Hessian_diagonal( ObjectiveFunction::CALCULATE, mPatch, value, grad, Hess, err );
00497     CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
00498     CPPUNIT_ASSERT( rval );
00499     size_t n = mPatch.num_free_vertices();
00500     CPPUNIT_ASSERT_EQUAL( n, grad.size() );
00501     CPPUNIT_ASSERT_EQUAL( n, Hess.size() );
00502 
00503     vector< Matrix3D > Hessians( n );
00504     for( size_t r = 0; r < n; ++r )
00505         Hessians[r] = Hess[r];
00506 
00507     if( !grad.empty() ) check_result( mPatch, power, value, arrptr( grad ), arrptr( Hessians ) );
00508 }
00509 
00510 void PMeanPTemplateTest::test_Hessian( double power )
00511 {
00512     MsqPrintError err( cout );
00513     double value;
00514     bool rval;
00515     vector< Vector3D > grad;
00516     MsqHessian Hess;
00517     Hess.initialize( mPatch, err );
00518     CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
00519 
00520     DistTestMetric metric;
00521     PMeanPTemplate func( power, &metric );
00522     rval = func.evaluate_with_Hessian( ObjectiveFunction::CALCULATE, mPatch, value, grad, Hess, err );
00523     CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
00524     CPPUNIT_ASSERT( rval );
00525     size_t n = mPatch.num_free_vertices();
00526     CPPUNIT_ASSERT_EQUAL( n, grad.size() );
00527     CPPUNIT_ASSERT_EQUAL( n, Hess.size() );
00528 
00529     Matrix3D zero( 0, 0, 0, 0, 0, 0, 0, 0, 0 );
00530     vector< Matrix3D > Hessians( n );
00531     for( size_t r = 0; r < n; ++r )
00532     {
00533         Matrix3D* mat = Hess.get_block( r, r );
00534         CPPUNIT_ASSERT( mat != 0 );
00535         Hessians[r] = *mat;
00536 
00537         for( size_t c = r + 1; c < n; ++c )
00538         {
00539             mat = Hess.get_block( r, c );
00540             if( mat ) CPPUNIT_ASSERT_MATRICES_EQUAL( zero, *mat, EPSILON );
00541         }
00542     }
00543 
00544     if( !grad.empty() ) check_result( mPatch, power, value, arrptr( grad ), arrptr( Hessians ) );
00545 }
00546 
00547 void PMeanPTemplateTest::check_result( PatchData& pd, double power, double value, Vector3D* gradient,
00548                                        Matrix3D* Hessian )
00549 {
00550     MsqPrintError err( cout );
00551     double mvalue, sum = 0;
00552     bool rval;
00553     vector< Vector3D > grads;
00554     vector< Matrix3D > Hess;
00555     vector< size_t > indices;
00556 
00557     DistTestMetric metric;
00558     vector< size_t > handles;
00559     metric.get_evaluations( pd, handles, OF_FREE_EVALS_ONLY, err );
00560     ASSERT_NO_ERROR( err );
00561 
00562     for( size_t i = 0; i < handles.size(); ++i )
00563     {
00564         rval = metric.evaluate_with_Hessian( pd, handles[i], mvalue, indices, grads, Hess, err );
00565         CPPUNIT_ASSERT( !MSQ_CHKERR( err ) && rval );
00566         sum += pow( mvalue, power );
00567 
00568         // if (!OF_FREE_EVALS_ONLY && indices.empty())
00569         //  continue;
00570 
00571         CPPUNIT_ASSERT_EQUAL( (size_t)1, indices.size() );
00572         CPPUNIT_ASSERT_EQUAL( handles[i], indices[0] );
00573         CPPUNIT_ASSERT_EQUAL( (size_t)1, grads.size() );
00574         CPPUNIT_ASSERT_EQUAL( (size_t)1, Hess.size() );
00575 
00576         if( gradient )
00577         {
00578             double f = power * pow( mvalue, power - 1 ) / handles.size();
00579             CPPUNIT_ASSERT_VECTORS_EQUAL( f * grads[0], gradient[i], EPSILON );
00580         }
00581         if( Hessian )
00582         {
00583             double f  = power / handles.size();
00584             double p2 = ( power - 1 ) * pow( mvalue, power - 2 );
00585             double p1 = pow( mvalue, power - 1 );
00586             Matrix3D m;
00587             m.outer_product( grads[0], grads[0] );
00588             m *= p2;
00589             m += p1 * Hess[0];
00590             m *= f;
00591             CPPUNIT_ASSERT_MATRICES_EQUAL( m, Hessian[i], EPSILON );
00592         }
00593     }
00594 
00595     CPPUNIT_ASSERT_DOUBLES_EQUAL( sum / handles.size(), value, EPSILON );
00596 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines