MOAB: Mesh Oriented datABase  (version 5.4.1)
CompareQM.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) [email protected]
00024 
00025   ***************************************************************** */
00026 
00027 /** \file CompareQM.cpp
00028  *  \brief
00029  *  \author Jason Kraftcheck
00030  */
00031 
00032 #include "CompareQM.hpp"
00033 #include "MsqError.hpp"
00034 #include <iostream>
00035 #include <cstdio>
00036 #include <cstdlib>
00037 
00038 namespace MBMesquite
00039 {
00040 
00041 CompareQM::CompareQM( QualityMetric* primary, QualityMetric* other, const char* primary_name, const char* other_name )
00042     : primaryMetric( primary ), otherMetric( other ), abortOnMismatch( false ), toleranceFactor( 1e-6 )
00043 {
00044     if( primary_name ) primaryName = primary_name;
00045     if( other_name ) otherName = other_name;
00046 }
00047 
00048 CompareQM::~CompareQM() {}
00049 
00050 void CompareQM::abort_on_mismatch( double tolerance )
00051 {
00052     abortOnMismatch = true;
00053     toleranceFactor = tolerance;
00054 }
00055 
00056 void CompareQM::do_not_abort()
00057 {
00058     abortOnMismatch = false;
00059 }
00060 
00061 bool CompareQM::will_abort_on_mismatch() const
00062 {
00063     return abortOnMismatch;
00064 }
00065 
00066 QualityMetric::MetricType CompareQM::get_metric_type() const
00067 {
00068     MetricType type1, type2;
00069     type1 = primaryMetric->get_metric_type();
00070     type2 = otherMetric->get_metric_type();
00071     if( type1 != type2 )
00072     {
00073         std::cerr << "Incompatible metric types in CompareQM" << std::endl << __FILE__ << ':' << __LINE__ << std::endl;
00074         if( abortOnMismatch ) abort();
00075         // otherwise just return some type because this function can't
00076         // flag an error.  The mismatch should cause get_evaluations
00077         // to fail anyway.
00078     }
00079     return type1;
00080 }
00081 
00082 std::string CompareQM::get_name() const
00083 {
00084     return primaryMetric->get_name();
00085 }
00086 
00087 int CompareQM::get_negate_flag() const
00088 {
00089     return primaryMetric->get_negate_flag();
00090 }
00091 
00092 void CompareQM::get_evaluations( PatchData& pd, std::vector< size_t >& handles, bool free_only, MsqError& err )
00093 {
00094     primaryMetric->get_evaluations( pd, handles, free_only, err );MSQ_ERRRTN( err );
00095 
00096     std::vector< size_t > handles2;
00097     otherMetric->get_evaluations( pd, handles2, free_only, err );MSQ_ERRRTN( err );
00098 
00099     std::vector< size_t > handles1( handles );
00100     std::sort( handles1.begin(), handles1.end() );
00101     std::sort( handles2.begin(), handles2.end() );
00102     if( handles1 != handles2 )
00103     {
00104         MSQ_SETERR( err )( "Incompatible metrics cannot be compared", MsqError::INVALID_STATE );
00105     }
00106 }
00107 
00108 bool CompareQM::check_valid( size_t handle, bool rval1, bool rval2 )
00109 {
00110     if( rval1 != rval2 )
00111     {
00112         std::cerr << "Metrics returned conflicting validity at location " << std::hex << handle << std::dec << std::endl
00113                   << __FILE__ << ":" << __LINE__ << std::endl;
00114         if( abortOnMismatch )
00115             abort();
00116         else
00117             return false;
00118     }
00119 
00120     return rval1;
00121 }
00122 
00123 double CompareQM::epsilon( double a, double b )
00124 {
00125     return toleranceFactor * std::max( 1.0, std::max( a, b ) );
00126 }
00127 
00128 void CompareQM::check_value( size_t handle, double value, double value2 )
00129 {
00130     valPrimary.add_value( value );
00131     valOther.add_value( value2 );
00132     valDiff.add_value( fabs( value - value2 ) );
00133 
00134     if( abortOnMismatch )
00135     {
00136         if( fabs( value - value2 ) > epsilon( value, value2 ) )
00137         {
00138             std::cerr << "Metric values to not match at location " << std::hex << handle << std::dec << std::endl
00139                       << "Primary: " << value << std::endl
00140                       << "Other  : " << value2 << std::endl
00141                       << __FILE__ << ":" << __LINE__ << std::endl;
00142             abort();
00143         }
00144     }
00145 }
00146 
00147 bool CompareQM::evaluate( PatchData& pd, size_t handle, double& value, MsqError& err )
00148 {
00149     double value2;
00150     bool rval1, rval2;
00151     rval1 = primaryMetric->evaluate( pd, handle, value, err );
00152     MSQ_ERRZERO( err );
00153     rval2 = otherMetric->evaluate( pd, handle, value2, err );
00154     MSQ_ERRZERO( err );
00155     if( !check_valid( handle, rval1, rval2 ) ) return false;
00156 
00157     check_value( handle, value, value2 );
00158     return true;
00159 }
00160 
00161 void CompareQM::index_mismatch( size_t handle,
00162                                 const std::vector< size_t >& idx1,
00163                                 const std::vector< size_t >& idx2,
00164                                 MsqError& err )
00165 {
00166     std::vector< size_t >::const_iterator i;
00167 
00168     std::cerr << "Metrics cannot be compared at location " << std::hex << handle << std::dec
00169               << " because they are incompatible." << std::endl
00170               << "Primary metric vertices: ";
00171     if( idx1.empty() )
00172         std::cerr << "(empty)";
00173     else
00174     {
00175         i = idx1.begin();
00176         std::cerr << *i;
00177         for( ++i; i != idx1.end(); ++i )
00178             std::cerr << ',' << *i;
00179     }
00180     std::cerr << std::endl << "Other metric vertices: ";
00181     if( idx2.empty() )
00182         std::cerr << "(empty)";
00183     else
00184     {
00185         i = idx2.begin();
00186         std::cerr << *i;
00187         for( ++i; i != idx2.end(); ++i )
00188             std::cerr << ',' << *i;
00189     }
00190     std::cerr << std::endl;
00191 
00192     if( abortOnMismatch ) abort();
00193     MSQ_SETERR( err )( "Cannot compare incompatible metrics", MsqError::INVALID_STATE );
00194 }
00195 
00196 void CompareQM::check_indices( size_t handle,
00197                                const std::vector< size_t >& idx1,
00198                                const std::vector< size_t >& idx2,
00199                                std::vector< size_t >& map_out,
00200                                MsqError& err )
00201 {
00202     if( idx1.size() != idx2.size() )
00203     {
00204         index_mismatch( handle, idx1, idx2, err );MSQ_ERRRTN( err );
00205     }
00206 
00207     std::vector< size_t >::const_iterator i, j;
00208     map_out.clear();
00209     for( i = idx1.begin(); i != idx1.end(); ++i )
00210     {
00211         j = std::find( idx2.begin(), idx2.end(), *i );
00212         if( j == idx2.end() )
00213         {
00214             index_mismatch( handle, idx1, idx2, err );MSQ_ERRRTN( err );
00215         }
00216         map_out.push_back( j - idx2.begin() );
00217     }
00218 }
00219 
00220 bool CompareQM::evaluate_with_indices( PatchData& pd,
00221                                        size_t handle,
00222                                        double& value,
00223                                        std::vector< size_t >& indices,
00224                                        MsqError& err )
00225 {
00226     bool valid1, valid2;
00227     double value2;
00228     std::vector< size_t > indices2, junk;
00229     valid1 = primaryMetric->evaluate_with_indices( pd, handle, value, indices, err );
00230     MSQ_ERRZERO( err );
00231     valid2 = otherMetric->evaluate_with_indices( pd, handle, value2, indices2, err );
00232     MSQ_ERRZERO( err );
00233     if( !check_valid( handle, valid1, valid2 ) ) return false;
00234 
00235     check_value( handle, value, value2 );
00236     check_indices( handle, indices, indices2, junk, err );
00237     MSQ_ERRZERO( err );
00238     return true;
00239 }
00240 
00241 void CompareQM::check_grad( size_t handle,
00242                             const std::vector< size_t >& indices,
00243                             const std::vector< size_t >& index_map,
00244                             const std::vector< Vector3D >& grad1,
00245                             const std::vector< Vector3D >& grad2 )
00246 {
00247     assert( index_map.size() == indices.size() );
00248     assert( index_map.size() == grad1.size() );
00249     assert( index_map.size() == grad2.size() );
00250     for( size_t i = 0; i < index_map.size(); ++i )
00251     {
00252         gradPrimary.add( grad1[i] );
00253         gradOther.add( grad2[i] );
00254         gradDiff.add_diff( grad1[i], grad2[index_map[i]] );
00255 
00256         if( abortOnMismatch )
00257         {
00258             if( ( grad1[i] - grad2[index_map[i]] ).length() >
00259                 epsilon( grad1[i].length(), grad2[index_map[i]].length() ) )
00260             {
00261                 std::cerr << "Gradients differ for metric evaluation at " << std::hex << handle << std::dec << std::endl
00262                           << "Primary metric derivs with respect to vertex " << indices[i] << ": " << grad1[i]
00263                           << std::endl
00264                           << "Other metric derivs with presect to vertex " << indices[i] << ": " << grad2[index_map[i]]
00265                           << std::endl
00266                           << __FILE__ << ":" << __LINE__ << std::endl;
00267                 abort();
00268             }
00269         }
00270     }
00271 }
00272 
00273 bool CompareQM::evaluate_with_gradient( PatchData& pd,
00274                                         size_t handle,
00275                                         double& value,
00276                                         std::vector< size_t >& indices,
00277                                         std::vector< Vector3D >& grad,
00278                                         MsqError& err )
00279 {
00280     bool valid1, valid2;
00281     double value2;
00282     std::vector< size_t > indices2, map;
00283     std::vector< Vector3D > grad2;
00284     valid1 = primaryMetric->evaluate_with_gradient( pd, handle, value, indices, grad, err );
00285     MSQ_ERRZERO( err );
00286     valid2 = otherMetric->evaluate_with_gradient( pd, handle, value2, indices2, grad2, err );
00287     MSQ_ERRZERO( err );
00288     if( !check_valid( handle, valid1, valid2 ) ) return false;
00289 
00290     check_value( handle, value, value2 );
00291     check_indices( handle, indices, indices2, map, err );
00292     MSQ_ERRZERO( err );
00293     check_grad( handle, indices, map, grad, grad2 );
00294     return true;
00295 }
00296 
00297 void CompareQM::check_hess_diag( size_t handle,
00298                                  const std::vector< size_t >& indices,
00299                                  const std::vector< size_t >& index_map,
00300                                  const std::vector< SymMatrix3D >& hess1,
00301                                  const std::vector< SymMatrix3D >& hess2 )
00302 {
00303     assert( index_map.size() == indices.size() );
00304     assert( index_map.size() == hess1.size() );
00305     assert( index_map.size() == hess2.size() );
00306     for( size_t i = 0; i < index_map.size(); ++i )
00307     {
00308         hessPrimary.add_diag( hess1[i] );
00309         hessOther.add_diag( hess2[i] );
00310         hessDiff.add_diag_diff( hess1[i], hess2[index_map[i]] );
00311 
00312         if( abortOnMismatch )
00313         {
00314             double eps = epsilon( Frobenius( hess1[i] ), Frobenius( hess2[index_map[i]] ) );
00315             if( Frobenius( hess1[i] - hess2[index_map[i]] ) > eps )
00316             {
00317                 std::cerr << "Hessian diagonal blocks differ for metric evaluation at " << std::hex << handle
00318                           << std::dec << std::endl
00319                           << "For second derivatives with repsect to vertex " << indices[i] << "twice" << std::endl
00320                           << "Primary metric derivs: " << hess1[i] << std::endl
00321                           << "Other metric derivs: " << hess2[index_map[i]] << std::endl
00322                           << __FILE__ << ":" << __LINE__ << std::endl;
00323                 abort();
00324             }
00325         }
00326     }
00327 }
00328 
00329 bool CompareQM::evaluate_with_Hessian_diagonal( PatchData& pd,
00330                                                 size_t handle,
00331                                                 double& value,
00332                                                 std::vector< size_t >& indices,
00333                                                 std::vector< Vector3D >& grad,
00334                                                 std::vector< SymMatrix3D >& hess,
00335                                                 MsqError& err )
00336 {
00337     bool valid1, valid2;
00338     double value2;
00339     std::vector< size_t > indices2, map;
00340     std::vector< Vector3D > grad2;
00341     std::vector< SymMatrix3D > hess2;
00342     valid1 = primaryMetric->evaluate_with_Hessian_diagonal( pd, handle, value, indices, grad, hess, err );
00343     MSQ_ERRZERO( err );
00344     valid2 = otherMetric->evaluate_with_Hessian_diagonal( pd, handle, value2, indices2, grad2, hess2, err );
00345     MSQ_ERRZERO( err );
00346     if( !check_valid( handle, valid1, valid2 ) ) return false;
00347 
00348     check_value( handle, value, value2 );
00349     check_indices( handle, indices, indices2, map, err );
00350     MSQ_ERRZERO( err );
00351     check_grad( handle, indices, map, grad, grad2 );
00352     check_hess_diag( handle, indices, map, hess, hess2 );
00353     return true;
00354 }
00355 
00356 void CompareQM::check_hess( size_t handle,
00357                             const std::vector< size_t >& indices,
00358                             const std::vector< size_t >& index_map,
00359                             const std::vector< Matrix3D >& hess1,
00360                             const std::vector< Matrix3D >& hess2 )
00361 {
00362     const size_t n = index_map.size();
00363     const size_t N = ( n + 1 ) * n / 2;
00364     assert( n == indices.size() );
00365     assert( N == hess1.size() );
00366     assert( N == hess2.size() );
00367 
00368     for( size_t r = 0; r < n; ++r )
00369     {
00370         const size_t r2 = index_map[r];
00371         for( size_t c = r; c < n; ++c )
00372         {
00373             const size_t c2 = index_map[c];
00374             size_t idx1     = n * r - r * ( r + 1 ) / 2 + c;
00375             Matrix3D h2;
00376             if( r2 <= c2 )
00377             {
00378                 size_t idx2 = n * r2 - r2 * ( r2 + 1 ) / 2 + c2;
00379                 h2          = hess2[idx2];
00380             }
00381             else
00382             {
00383                 size_t idx2 = n * c2 - c2 * ( c2 + 1 ) / 2 + r2;
00384                 h2          = transpose( hess2[idx2] );
00385             }
00386 
00387             if( r == c )
00388             {
00389                 hessPrimary.add_diag( hess1[idx1] );
00390                 hessOther.add_diag( h2 );
00391                 hessDiff.add_diag_diff( hess1[idx1], h2 );
00392             }
00393             else
00394             {
00395                 hessPrimary.add_nondiag( hess1[idx1] );
00396                 hessOther.add_nondiag( h2 );
00397                 hessDiff.add_nondiag_diff( hess1[idx1], h2 );
00398             }
00399 
00400             if( abortOnMismatch )
00401             {
00402                 double eps = epsilon( sqrt( Frobenius_2( hess1[idx1] ) ), sqrt( Frobenius_2( h2 ) ) );
00403                 if( sqrt( Frobenius_2( hess1[idx1] - h2 ) ) > eps )
00404                 {
00405                     std::cerr << "Hessian blocks differ for metric evaluation at " << std::hex << handle << std::dec
00406                               << std::endl
00407                               << "For second derivatives with repsect to vertices " << indices[r] << " and "
00408                               << indices[c] << std::endl
00409                               << "Primary metric derivs: " << hess1[idx1] << std::endl
00410                               << "Other metric derivs: " << h2 << std::endl
00411                               << __FILE__ << ":" << __LINE__ << std::endl;
00412                     abort();
00413                 }
00414             }
00415         }
00416     }
00417 }
00418 
00419 bool CompareQM::evaluate_with_Hessian( PatchData& pd,
00420                                        size_t handle,
00421                                        double& value,
00422                                        std::vector< size_t >& indices,
00423                                        std::vector< Vector3D >& grad,
00424                                        std::vector< Matrix3D >& hess,
00425                                        MsqError& err )
00426 {
00427     bool valid1, valid2;
00428     double value2;
00429     std::vector< size_t > indices2, map;
00430     std::vector< Vector3D > grad2;
00431     std::vector< Matrix3D > hess2;
00432     valid1 = primaryMetric->evaluate_with_Hessian( pd, handle, value, indices, grad, hess, err );
00433     MSQ_ERRZERO( err );
00434     valid2 = otherMetric->evaluate_with_Hessian( pd, handle, value2, indices2, grad2, hess2, err );
00435     MSQ_ERRZERO( err );
00436     if( !check_valid( handle, valid1, valid2 ) ) return false;
00437 
00438     check_value( handle, value, value2 );
00439     check_indices( handle, indices, indices2, map, err );
00440     MSQ_ERRZERO( err );
00441     check_grad( handle, indices, map, grad, grad2 );
00442     check_hess( handle, indices, map, hess, hess2 );
00443     return true;
00444 }
00445 
00446 void CompareQM::GradStat::add( Vector3D grad )
00447 {
00448     x.add_value( grad[0] );
00449     y.add_value( grad[1] );
00450     z.add_value( grad[2] );
00451 }
00452 
00453 void CompareQM::GradStat::add_diff( Vector3D grad1, Vector3D grad2 )
00454 {
00455     x.add_value( fabs( grad1[0] - grad2[0] ) );
00456     y.add_value( fabs( grad1[1] - grad2[1] ) );
00457     z.add_value( fabs( grad1[2] - grad2[2] ) );
00458 }
00459 
00460 void CompareQM::HessStat::add_diag( Matrix3D hess )
00461 {
00462     xx.add_value( hess[0][0] );
00463     xy.add_value( ( hess[0][1] + hess[1][0] ) / 2 );
00464     xz.add_value( ( hess[0][2] + hess[2][0] ) / 2 );
00465     yy.add_value( hess[1][1] );
00466     yz.add_value( ( hess[1][2] + hess[2][1] ) / 2 );
00467     zz.add_value( hess[2][2] );
00468 }
00469 
00470 void CompareQM::HessStat::add_diag( SymMatrix3D hess )
00471 {
00472     xx.add_value( hess[0] );
00473     xy.add_value( hess[1] );
00474     xz.add_value( hess[2] );
00475     yy.add_value( hess[3] );
00476     yz.add_value( hess[4] );
00477     zz.add_value( hess[5] );
00478 }
00479 
00480 void CompareQM::HessStat::add_diag_diff( Matrix3D hess1, Matrix3D hess2 )
00481 {
00482     Matrix3D d = hess1 - hess2;
00483     Matrix3D hess( fabs( d[0][0] ), fabs( d[0][1] ), fabs( d[0][2] ), fabs( d[1][0] ), fabs( d[1][1] ), fabs( d[1][2] ),
00484                    fabs( d[2][0] ), fabs( d[2][1] ), fabs( d[2][2] ) );
00485     add_diag( hess );
00486 }
00487 
00488 void CompareQM::HessStat::add_diag_diff( SymMatrix3D hess1, SymMatrix3D hess2 )
00489 {
00490     SymMatrix3D d = hess1 - hess2;
00491     SymMatrix3D hess( fabs( d[0] ), fabs( d[1] ), fabs( d[2] ), fabs( d[3] ), fabs( d[4] ), fabs( d[5] ) );
00492     add_diag( hess );
00493 }
00494 
00495 void CompareQM::HessStat::add_nondiag( Matrix3D hess )
00496 {
00497     xx.add_value( hess[0][0] );
00498     xy.add_value( hess[0][1] );
00499     xy.add_value( hess[1][0] );
00500     xz.add_value( hess[0][2] );
00501     xz.add_value( hess[2][0] );
00502     yy.add_value( hess[1][1] );
00503     yz.add_value( hess[1][2] );
00504     yz.add_value( hess[2][1] );
00505     zz.add_value( hess[2][2] );
00506 }
00507 
00508 void CompareQM::HessStat::add_nondiag_diff( Matrix3D hess1, Matrix3D hess2 )
00509 {
00510     Matrix3D d = hess1 - hess2;
00511     Matrix3D hess( fabs( d[0][0] ), fabs( d[0][1] ), fabs( d[0][2] ), fabs( d[1][0] ), fabs( d[1][1] ), fabs( d[1][2] ),
00512                    fabs( d[2][0] ), fabs( d[2][1] ), fabs( d[2][2] ) );
00513     add_nondiag( hess );
00514 }
00515 
00516 static void print( const char* title,
00517                    const char* name1,
00518                    const char* name2,
00519                    const SimpleStats& s1,
00520                    const SimpleStats& s2,
00521                    const SimpleStats& sd )
00522 {
00523     const char named[] = "difference";
00524     int len = std::max( std::max( strlen( named ), strlen( title ) ), std::max( strlen( name1 ), strlen( name2 ) ) );
00525     std::vector< char > dashes( len, '-' );
00526     dashes.push_back( '\0' );
00527     printf( "%*s %12s %12s %12s %12s %12s\n", len, title, "minimum", "average", "rms", "maximum", "std.dev." );
00528     printf( "%s ------------ ------------ ------------ ------------ ------------\n", &dashes[0] );
00529     printf( "%*s % 12g % 12g % 12g % 12g % 12g\n", len, name1, s1.minimum(), s1.average(), s1.rms(), s1.maximum(),
00530             s1.standard_deviation() );
00531     printf( "%*s % 12g % 12g % 12g % 12g % 12g\n", len, name2, s2.minimum(), s2.average(), s2.rms(), s2.maximum(),
00532             s2.standard_deviation() );
00533     printf( "%*s % 12g % 12g % 12g % 12g % 12g\n", len, named, sd.minimum(), sd.average(), sd.rms(), sd.maximum(),
00534             sd.standard_deviation() );
00535     printf( "\n" );
00536 }
00537 
00538 void CompareQM::print_stats() const
00539 {
00540     std::string name1 = primaryName.empty() ? primaryMetric->get_name() : primaryName;
00541     std::string name2 = otherName.empty() ? otherMetric->get_name() : otherName;
00542     print( "Values", name1.c_str(), name2.c_str(), valPrimary, valOther, valDiff );
00543     print( "Gradient X", name1.c_str(), name2.c_str(), gradPrimary.x, gradOther.x, gradDiff.x );
00544     print( "Gradient Y", name1.c_str(), name2.c_str(), gradPrimary.y, gradOther.y, gradDiff.y );
00545     print( "Gradient Z", name1.c_str(), name2.c_str(), gradPrimary.z, gradOther.z, gradDiff.z );
00546     print( "Hessian XX", name1.c_str(), name2.c_str(), hessPrimary.xx, hessOther.xx, hessDiff.xx );
00547     print( "Hessian XY", name1.c_str(), name2.c_str(), hessPrimary.xy, hessOther.xy, hessDiff.xy );
00548     print( "Hessian XZ", name1.c_str(), name2.c_str(), hessPrimary.xz, hessOther.xz, hessDiff.xz );
00549     print( "Hessian YY", name1.c_str(), name2.c_str(), hessPrimary.yy, hessOther.yy, hessDiff.yy );
00550     print( "Hessian YZ", name1.c_str(), name2.c_str(), hessPrimary.yz, hessOther.yz, hessDiff.yz );
00551     print( "Hessian ZZ", name1.c_str(), name2.c_str(), hessPrimary.zz, hessOther.zz, hessDiff.zz );
00552 }
00553 
00554 }  // namespace MBMesquite
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines