MOAB: Mesh Oriented datABase  (version 5.4.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,
00410                                 size_t vtx_idx,
00411                                 double& value,
00412                                 vector< size_t >& indices,
00413                                 MsqError& err );
00414     // bool evaluate_with_gradient( PatchData& pd, size_t vtx_idx, double &value, vector<size_t>&
00415     // indices, vector<Vector3D>& grad, MsqError& err );
00416     bool falseEval;
00417     bool failEval;
00418 };
00419 
00420 bool DistTestMetric::evaluate( PatchData& pd, size_t vtx_idx, double& value, MsqError& err )
00421 {
00422     if( failEval )
00423     {
00424         MSQ_SETERR( err )( MsqError::INVALID_STATE );
00425         return true;
00426     }
00427 
00428     const MsqVertex& vtx = pd.vertex_by_index( vtx_idx );
00429     value                = vtx.length_squared();
00430     return !falseEval;
00431 }
00432 
00433 bool DistTestMetric::evaluate_with_indices( PatchData& pd,
00434                                             size_t vtx_idx,
00435                                             double& value,
00436                                             vector< size_t >& indices,
00437                                             MsqError& err )
00438 {
00439     if( failEval )
00440     {
00441         MSQ_SETERR( err )( MsqError::INVALID_STATE );
00442         return true;
00443     }
00444 
00445     indices.clear();
00446     if( vtx_idx < pd.num_free_vertices() ) indices.push_back( vtx_idx );
00447 
00448     const MsqVertex& vtx = pd.vertex_by_index( vtx_idx );
00449     value                = vtx.length_squared();
00450     return !falseEval;
00451 }
00452 
00453 void PMeanPTemplateTest::test_eval_type( OFTestMode eval_func, ObjectiveFunction::EvalType type )
00454 {
00455     PMeanPTemplate func( 1, NULL );
00456     ObjectiveFunctionTests::test_eval_type( type, eval_func, &func );
00457 }
00458 
00459 void PMeanPTemplateTest::test_evaluate( double power )
00460 {
00461     MsqPrintError err( cout );
00462     double value;
00463     bool rval;
00464 
00465     DistTestMetric metric;
00466     PMeanPTemplate func( power, &metric );
00467     rval = func.evaluate( ObjectiveFunction::CALCULATE, mPatch, value, OF_FREE_EVALS_ONLY, err );
00468     CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
00469     CPPUNIT_ASSERT( rval );
00470 
00471     check_result( mPatch, power, value );
00472 }
00473 
00474 void PMeanPTemplateTest::test_gradient( double power )
00475 {
00476     MsqPrintError err( cout );
00477     double value;
00478     bool rval;
00479     vector< Vector3D > grad;
00480 
00481     DistTestMetric metric;
00482     PMeanPTemplate func( power, &metric );
00483     rval = func.evaluate_with_gradient( ObjectiveFunction::CALCULATE, mPatch, value, grad, err );
00484     CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
00485     CPPUNIT_ASSERT( rval );
00486     CPPUNIT_ASSERT_EQUAL( mPatch.num_free_vertices(), grad.size() );
00487 
00488     if( !grad.empty() ) check_result( mPatch, power, value, arrptr( grad ) );
00489 }
00490 
00491 void PMeanPTemplateTest::test_diagonal( double power )
00492 {
00493     MsqPrintError err( cout );
00494     double value;
00495     bool rval;
00496     vector< Vector3D > grad;
00497     vector< SymMatrix3D > Hess;
00498     CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
00499 
00500     DistTestMetric metric;
00501     PMeanPTemplate func( power, &metric );
00502     rval = func.evaluate_with_Hessian_diagonal( ObjectiveFunction::CALCULATE, mPatch, value, grad, Hess, err );
00503     CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
00504     CPPUNIT_ASSERT( rval );
00505     size_t n = mPatch.num_free_vertices();
00506     CPPUNIT_ASSERT_EQUAL( n, grad.size() );
00507     CPPUNIT_ASSERT_EQUAL( n, Hess.size() );
00508 
00509     vector< Matrix3D > Hessians( n );
00510     for( size_t r = 0; r < n; ++r )
00511         Hessians[r] = Hess[r];
00512 
00513     if( !grad.empty() ) check_result( mPatch, power, value, arrptr( grad ), arrptr( Hessians ) );
00514 }
00515 
00516 void PMeanPTemplateTest::test_Hessian( double power )
00517 {
00518     MsqPrintError err( cout );
00519     double value;
00520     bool rval;
00521     vector< Vector3D > grad;
00522     MsqHessian Hess;
00523     Hess.initialize( mPatch, err );
00524     CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
00525 
00526     DistTestMetric metric;
00527     PMeanPTemplate func( power, &metric );
00528     rval = func.evaluate_with_Hessian( ObjectiveFunction::CALCULATE, mPatch, value, grad, Hess, err );
00529     CPPUNIT_ASSERT( !MSQ_CHKERR( err ) );
00530     CPPUNIT_ASSERT( rval );
00531     size_t n = mPatch.num_free_vertices();
00532     CPPUNIT_ASSERT_EQUAL( n, grad.size() );
00533     CPPUNIT_ASSERT_EQUAL( n, Hess.size() );
00534 
00535     Matrix3D zero( 0, 0, 0, 0, 0, 0, 0, 0, 0 );
00536     vector< Matrix3D > Hessians( n );
00537     for( size_t r = 0; r < n; ++r )
00538     {
00539         Matrix3D* mat = Hess.get_block( r, r );
00540         CPPUNIT_ASSERT( mat != 0 );
00541         Hessians[r] = *mat;
00542 
00543         for( size_t c = r + 1; c < n; ++c )
00544         {
00545             mat = Hess.get_block( r, c );
00546             if( mat ) CPPUNIT_ASSERT_MATRICES_EQUAL( zero, *mat, EPSILON );
00547         }
00548     }
00549 
00550     if( !grad.empty() ) check_result( mPatch, power, value, arrptr( grad ), arrptr( Hessians ) );
00551 }
00552 
00553 void PMeanPTemplateTest::check_result( PatchData& pd,
00554                                        double power,
00555                                        double value,
00556                                        Vector3D* gradient,
00557                                        Matrix3D* Hessian )
00558 {
00559     MsqPrintError err( cout );
00560     double mvalue, sum = 0;
00561     bool rval;
00562     vector< Vector3D > grads;
00563     vector< Matrix3D > Hess;
00564     vector< size_t > indices;
00565 
00566     DistTestMetric metric;
00567     vector< size_t > handles;
00568     metric.get_evaluations( pd, handles, OF_FREE_EVALS_ONLY, err );
00569     ASSERT_NO_ERROR( err );
00570 
00571     for( size_t i = 0; i < handles.size(); ++i )
00572     {
00573         rval = metric.evaluate_with_Hessian( pd, handles[i], mvalue, indices, grads, Hess, err );
00574         CPPUNIT_ASSERT( !MSQ_CHKERR( err ) && rval );
00575         sum += pow( mvalue, power );
00576 
00577         // if (!OF_FREE_EVALS_ONLY && indices.empty())
00578         //  continue;
00579 
00580         CPPUNIT_ASSERT_EQUAL( (size_t)1, indices.size() );
00581         CPPUNIT_ASSERT_EQUAL( handles[i], indices[0] );
00582         CPPUNIT_ASSERT_EQUAL( (size_t)1, grads.size() );
00583         CPPUNIT_ASSERT_EQUAL( (size_t)1, Hess.size() );
00584 
00585         if( gradient )
00586         {
00587             double f = power * pow( mvalue, power - 1 ) / handles.size();
00588             CPPUNIT_ASSERT_VECTORS_EQUAL( f * grads[0], gradient[i], EPSILON );
00589         }
00590         if( Hessian )
00591         {
00592             double f  = power / handles.size();
00593             double p2 = ( power - 1 ) * pow( mvalue, power - 2 );
00594             double p1 = pow( mvalue, power - 1 );
00595             Matrix3D m;
00596             m.outer_product( grads[0], grads[0] );
00597             m *= p2;
00598             m += p1 * Hess[0];
00599             m *= f;
00600             CPPUNIT_ASSERT_MATRICES_EQUAL( m, Hessian[i], EPSILON );
00601         }
00602     }
00603 
00604     CPPUNIT_ASSERT_DOUBLES_EQUAL( sum / handles.size(), value, EPSILON );
00605 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines