MOAB: Mesh Oriented datABase  (version 5.3.0)
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) kraftche@cae.wisc.edu
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, const std::vector< size_t >& idx1, const std::vector< size_t >& idx2,
00162                                 MsqError& err )
00163 {
00164     std::vector< size_t >::const_iterator i;
00165 
00166     std::cerr << "Metrics cannot be compared at location " << std::hex << handle << std::dec
00167               << " because they are incompatible." << std::endl
00168               << "Primary metric vertices: ";
00169     if( idx1.empty() )
00170         std::cerr << "(empty)";
00171     else
00172     {
00173         i = idx1.begin();
00174         std::cerr << *i;
00175         for( ++i; i != idx1.end(); ++i )
00176             std::cerr << ',' << *i;
00177     }
00178     std::cerr << std::endl << "Other metric vertices: ";
00179     if( idx2.empty() )
00180         std::cerr << "(empty)";
00181     else
00182     {
00183         i = idx2.begin();
00184         std::cerr << *i;
00185         for( ++i; i != idx2.end(); ++i )
00186             std::cerr << ',' << *i;
00187     }
00188     std::cerr << std::endl;
00189 
00190     if( abortOnMismatch ) abort();
00191     MSQ_SETERR( err )( "Cannot compare incompatible metrics", MsqError::INVALID_STATE );
00192 }
00193 
00194 void CompareQM::check_indices( size_t handle, const std::vector< size_t >& idx1, const std::vector< size_t >& idx2,
00195                                std::vector< size_t >& map_out, MsqError& err )
00196 {
00197     if( idx1.size() != idx2.size() )
00198     {
00199         index_mismatch( handle, idx1, idx2, err );MSQ_ERRRTN( err );
00200     }
00201 
00202     std::vector< size_t >::const_iterator i, j;
00203     map_out.clear();
00204     for( i = idx1.begin(); i != idx1.end(); ++i )
00205     {
00206         j = std::find( idx2.begin(), idx2.end(), *i );
00207         if( j == idx2.end() )
00208         {
00209             index_mismatch( handle, idx1, idx2, err );MSQ_ERRRTN( err );
00210         }
00211         map_out.push_back( j - idx2.begin() );
00212     }
00213 }
00214 
00215 bool CompareQM::evaluate_with_indices( PatchData& pd, size_t handle, double& value, std::vector< size_t >& indices,
00216                                        MsqError& err )
00217 {
00218     bool valid1, valid2;
00219     double value2;
00220     std::vector< size_t > indices2, junk;
00221     valid1 = primaryMetric->evaluate_with_indices( pd, handle, value, indices, err );
00222     MSQ_ERRZERO( err );
00223     valid2 = otherMetric->evaluate_with_indices( pd, handle, value2, indices2, err );
00224     MSQ_ERRZERO( err );
00225     if( !check_valid( handle, valid1, valid2 ) ) return false;
00226 
00227     check_value( handle, value, value2 );
00228     check_indices( handle, indices, indices2, junk, err );
00229     MSQ_ERRZERO( err );
00230     return true;
00231 }
00232 
00233 void CompareQM::check_grad( size_t handle, const std::vector< size_t >& indices, const std::vector< size_t >& index_map,
00234                             const std::vector< Vector3D >& grad1, const std::vector< Vector3D >& grad2 )
00235 {
00236     assert( index_map.size() == indices.size() );
00237     assert( index_map.size() == grad1.size() );
00238     assert( index_map.size() == grad2.size() );
00239     for( size_t i = 0; i < index_map.size(); ++i )
00240     {
00241         gradPrimary.add( grad1[i] );
00242         gradOther.add( grad2[i] );
00243         gradDiff.add_diff( grad1[i], grad2[index_map[i]] );
00244 
00245         if( abortOnMismatch )
00246         {
00247             if( ( grad1[i] - grad2[index_map[i]] ).length() >
00248                 epsilon( grad1[i].length(), grad2[index_map[i]].length() ) )
00249             {
00250                 std::cerr << "Gradients differ for metric evaluation at " << std::hex << handle << std::dec << std::endl
00251                           << "Primary metric derivs with respect to vertex " << indices[i] << ": " << grad1[i]
00252                           << std::endl
00253                           << "Other metric derivs with presect to vertex " << indices[i] << ": " << grad2[index_map[i]]
00254                           << std::endl
00255                           << __FILE__ << ":" << __LINE__ << std::endl;
00256                 abort();
00257             }
00258         }
00259     }
00260 }
00261 
00262 bool CompareQM::evaluate_with_gradient( PatchData& pd, size_t handle, double& value, std::vector< size_t >& indices,
00263                                         std::vector< Vector3D >& grad, MsqError& err )
00264 {
00265     bool valid1, valid2;
00266     double value2;
00267     std::vector< size_t > indices2, map;
00268     std::vector< Vector3D > grad2;
00269     valid1 = primaryMetric->evaluate_with_gradient( pd, handle, value, indices, grad, err );
00270     MSQ_ERRZERO( err );
00271     valid2 = otherMetric->evaluate_with_gradient( pd, handle, value2, indices2, grad2, err );
00272     MSQ_ERRZERO( err );
00273     if( !check_valid( handle, valid1, valid2 ) ) return false;
00274 
00275     check_value( handle, value, value2 );
00276     check_indices( handle, indices, indices2, map, err );
00277     MSQ_ERRZERO( err );
00278     check_grad( handle, indices, map, grad, grad2 );
00279     return true;
00280 }
00281 
00282 void CompareQM::check_hess_diag( size_t handle, const std::vector< size_t >& indices,
00283                                  const std::vector< size_t >& index_map, const std::vector< SymMatrix3D >& hess1,
00284                                  const std::vector< SymMatrix3D >& hess2 )
00285 {
00286     assert( index_map.size() == indices.size() );
00287     assert( index_map.size() == hess1.size() );
00288     assert( index_map.size() == hess2.size() );
00289     for( size_t i = 0; i < index_map.size(); ++i )
00290     {
00291         hessPrimary.add_diag( hess1[i] );
00292         hessOther.add_diag( hess2[i] );
00293         hessDiff.add_diag_diff( hess1[i], hess2[index_map[i]] );
00294 
00295         if( abortOnMismatch )
00296         {
00297             double eps = epsilon( Frobenius( hess1[i] ), Frobenius( hess2[index_map[i]] ) );
00298             if( Frobenius( hess1[i] - hess2[index_map[i]] ) > eps )
00299             {
00300                 std::cerr << "Hessian diagonal blocks differ for metric evaluation at " << std::hex << handle
00301                           << std::dec << std::endl
00302                           << "For second derivatives with repsect to vertex " << indices[i] << "twice" << std::endl
00303                           << "Primary metric derivs: " << hess1[i] << std::endl
00304                           << "Other metric derivs: " << hess2[index_map[i]] << std::endl
00305                           << __FILE__ << ":" << __LINE__ << std::endl;
00306                 abort();
00307             }
00308         }
00309     }
00310 }
00311 
00312 bool CompareQM::evaluate_with_Hessian_diagonal( PatchData& pd, size_t handle, double& value,
00313                                                 std::vector< size_t >& indices, std::vector< Vector3D >& grad,
00314                                                 std::vector< SymMatrix3D >& hess, MsqError& err )
00315 {
00316     bool valid1, valid2;
00317     double value2;
00318     std::vector< size_t > indices2, map;
00319     std::vector< Vector3D > grad2;
00320     std::vector< SymMatrix3D > hess2;
00321     valid1 = primaryMetric->evaluate_with_Hessian_diagonal( pd, handle, value, indices, grad, hess, err );
00322     MSQ_ERRZERO( err );
00323     valid2 = otherMetric->evaluate_with_Hessian_diagonal( pd, handle, value2, indices2, grad2, hess2, err );
00324     MSQ_ERRZERO( err );
00325     if( !check_valid( handle, valid1, valid2 ) ) return false;
00326 
00327     check_value( handle, value, value2 );
00328     check_indices( handle, indices, indices2, map, err );
00329     MSQ_ERRZERO( err );
00330     check_grad( handle, indices, map, grad, grad2 );
00331     check_hess_diag( handle, indices, map, hess, hess2 );
00332     return true;
00333 }
00334 
00335 void CompareQM::check_hess( size_t handle, const std::vector< size_t >& indices, const std::vector< size_t >& index_map,
00336                             const std::vector< Matrix3D >& hess1, const std::vector< Matrix3D >& hess2 )
00337 {
00338     const size_t n = index_map.size();
00339     const size_t N = ( n + 1 ) * n / 2;
00340     assert( n == indices.size() );
00341     assert( N == hess1.size() );
00342     assert( N == hess2.size() );
00343 
00344     for( size_t r = 0; r < n; ++r )
00345     {
00346         const size_t r2 = index_map[r];
00347         for( size_t c = r; c < n; ++c )
00348         {
00349             const size_t c2 = index_map[c];
00350             size_t idx1     = n * r - r * ( r + 1 ) / 2 + c;
00351             Matrix3D h2;
00352             if( r2 <= c2 )
00353             {
00354                 size_t idx2 = n * r2 - r2 * ( r2 + 1 ) / 2 + c2;
00355                 h2          = hess2[idx2];
00356             }
00357             else
00358             {
00359                 size_t idx2 = n * c2 - c2 * ( c2 + 1 ) / 2 + r2;
00360                 h2          = transpose( hess2[idx2] );
00361             }
00362 
00363             if( r == c )
00364             {
00365                 hessPrimary.add_diag( hess1[idx1] );
00366                 hessOther.add_diag( h2 );
00367                 hessDiff.add_diag_diff( hess1[idx1], h2 );
00368             }
00369             else
00370             {
00371                 hessPrimary.add_nondiag( hess1[idx1] );
00372                 hessOther.add_nondiag( h2 );
00373                 hessDiff.add_nondiag_diff( hess1[idx1], h2 );
00374             }
00375 
00376             if( abortOnMismatch )
00377             {
00378                 double eps = epsilon( sqrt( Frobenius_2( hess1[idx1] ) ), sqrt( Frobenius_2( h2 ) ) );
00379                 if( sqrt( Frobenius_2( hess1[idx1] - h2 ) ) > eps )
00380                 {
00381                     std::cerr << "Hessian blocks differ for metric evaluation at " << std::hex << handle << std::dec
00382                               << std::endl
00383                               << "For second derivatives with repsect to vertices " << indices[r] << " and "
00384                               << indices[c] << std::endl
00385                               << "Primary metric derivs: " << hess1[idx1] << std::endl
00386                               << "Other metric derivs: " << h2 << std::endl
00387                               << __FILE__ << ":" << __LINE__ << std::endl;
00388                     abort();
00389                 }
00390             }
00391         }
00392     }
00393 }
00394 
00395 bool CompareQM::evaluate_with_Hessian( PatchData& pd, size_t handle, double& value, std::vector< size_t >& indices,
00396                                        std::vector< Vector3D >& grad, std::vector< Matrix3D >& hess, MsqError& err )
00397 {
00398     bool valid1, valid2;
00399     double value2;
00400     std::vector< size_t > indices2, map;
00401     std::vector< Vector3D > grad2;
00402     std::vector< Matrix3D > hess2;
00403     valid1 = primaryMetric->evaluate_with_Hessian( pd, handle, value, indices, grad, hess, err );
00404     MSQ_ERRZERO( err );
00405     valid2 = otherMetric->evaluate_with_Hessian( pd, handle, value2, indices2, grad2, hess2, err );
00406     MSQ_ERRZERO( err );
00407     if( !check_valid( handle, valid1, valid2 ) ) return false;
00408 
00409     check_value( handle, value, value2 );
00410     check_indices( handle, indices, indices2, map, err );
00411     MSQ_ERRZERO( err );
00412     check_grad( handle, indices, map, grad, grad2 );
00413     check_hess( handle, indices, map, hess, hess2 );
00414     return true;
00415 }
00416 
00417 void CompareQM::GradStat::add( Vector3D grad )
00418 {
00419     x.add_value( grad[0] );
00420     y.add_value( grad[1] );
00421     z.add_value( grad[2] );
00422 }
00423 
00424 void CompareQM::GradStat::add_diff( Vector3D grad1, Vector3D grad2 )
00425 {
00426     x.add_value( fabs( grad1[0] - grad2[0] ) );
00427     y.add_value( fabs( grad1[1] - grad2[1] ) );
00428     z.add_value( fabs( grad1[2] - grad2[2] ) );
00429 }
00430 
00431 void CompareQM::HessStat::add_diag( Matrix3D hess )
00432 {
00433     xx.add_value( hess[0][0] );
00434     xy.add_value( ( hess[0][1] + hess[1][0] ) / 2 );
00435     xz.add_value( ( hess[0][2] + hess[2][0] ) / 2 );
00436     yy.add_value( hess[1][1] );
00437     yz.add_value( ( hess[1][2] + hess[2][1] ) / 2 );
00438     zz.add_value( hess[2][2] );
00439 }
00440 
00441 void CompareQM::HessStat::add_diag( SymMatrix3D hess )
00442 {
00443     xx.add_value( hess[0] );
00444     xy.add_value( hess[1] );
00445     xz.add_value( hess[2] );
00446     yy.add_value( hess[3] );
00447     yz.add_value( hess[4] );
00448     zz.add_value( hess[5] );
00449 }
00450 
00451 void CompareQM::HessStat::add_diag_diff( Matrix3D hess1, Matrix3D hess2 )
00452 {
00453     Matrix3D d = hess1 - hess2;
00454     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] ),
00455                    fabs( d[2][0] ), fabs( d[2][1] ), fabs( d[2][2] ) );
00456     add_diag( hess );
00457 }
00458 
00459 void CompareQM::HessStat::add_diag_diff( SymMatrix3D hess1, SymMatrix3D hess2 )
00460 {
00461     SymMatrix3D d = hess1 - hess2;
00462     SymMatrix3D hess( fabs( d[0] ), fabs( d[1] ), fabs( d[2] ), fabs( d[3] ), fabs( d[4] ), fabs( d[5] ) );
00463     add_diag( hess );
00464 }
00465 
00466 void CompareQM::HessStat::add_nondiag( Matrix3D hess )
00467 {
00468     xx.add_value( hess[0][0] );
00469     xy.add_value( hess[0][1] );
00470     xy.add_value( hess[1][0] );
00471     xz.add_value( hess[0][2] );
00472     xz.add_value( hess[2][0] );
00473     yy.add_value( hess[1][1] );
00474     yz.add_value( hess[1][2] );
00475     yz.add_value( hess[2][1] );
00476     zz.add_value( hess[2][2] );
00477 }
00478 
00479 void CompareQM::HessStat::add_nondiag_diff( Matrix3D hess1, Matrix3D hess2 )
00480 {
00481     Matrix3D d = hess1 - hess2;
00482     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] ),
00483                    fabs( d[2][0] ), fabs( d[2][1] ), fabs( d[2][2] ) );
00484     add_nondiag( hess );
00485 }
00486 
00487 static void print( const char* title, const char* name1, const char* name2, const SimpleStats& s1,
00488                    const SimpleStats& s2, const SimpleStats& sd )
00489 {
00490     const char named[] = "difference";
00491     int len = std::max( std::max( strlen( named ), strlen( title ) ), std::max( strlen( name1 ), strlen( name2 ) ) );
00492     std::vector< char > dashes( len, '-' );
00493     dashes.push_back( '\0' );
00494     printf( "%*s %12s %12s %12s %12s %12s\n", len, title, "minimum", "average", "rms", "maximum", "std.dev." );
00495     printf( "%s ------------ ------------ ------------ ------------ ------------\n", &dashes[0] );
00496     printf( "%*s % 12g % 12g % 12g % 12g % 12g\n", len, name1, s1.minimum(), s1.average(), s1.rms(), s1.maximum(),
00497             s1.standard_deviation() );
00498     printf( "%*s % 12g % 12g % 12g % 12g % 12g\n", len, name2, s2.minimum(), s2.average(), s2.rms(), s2.maximum(),
00499             s2.standard_deviation() );
00500     printf( "%*s % 12g % 12g % 12g % 12g % 12g\n", len, named, sd.minimum(), sd.average(), sd.rms(), sd.maximum(),
00501             sd.standard_deviation() );
00502     printf( "\n" );
00503 }
00504 
00505 void CompareQM::print_stats() const
00506 {
00507     std::string name1 = primaryName.empty() ? primaryMetric->get_name() : primaryName;
00508     std::string name2 = otherName.empty() ? otherMetric->get_name() : otherName;
00509     print( "Values", name1.c_str(), name2.c_str(), valPrimary, valOther, valDiff );
00510     print( "Gradient X", name1.c_str(), name2.c_str(), gradPrimary.x, gradOther.x, gradDiff.x );
00511     print( "Gradient Y", name1.c_str(), name2.c_str(), gradPrimary.y, gradOther.y, gradDiff.y );
00512     print( "Gradient Z", name1.c_str(), name2.c_str(), gradPrimary.z, gradOther.z, gradDiff.z );
00513     print( "Hessian XX", name1.c_str(), name2.c_str(), hessPrimary.xx, hessOther.xx, hessDiff.xx );
00514     print( "Hessian XY", name1.c_str(), name2.c_str(), hessPrimary.xy, hessOther.xy, hessDiff.xy );
00515     print( "Hessian XZ", name1.c_str(), name2.c_str(), hessPrimary.xz, hessOther.xz, hessDiff.xz );
00516     print( "Hessian YY", name1.c_str(), name2.c_str(), hessPrimary.yy, hessOther.yy, hessDiff.yy );
00517     print( "Hessian YZ", name1.c_str(), name2.c_str(), hessPrimary.yz, hessOther.yz, hessDiff.yz );
00518     print( "Hessian ZZ", name1.c_str(), name2.c_str(), hessPrimary.zz, hessOther.zz, hessDiff.zz );
00519 }
00520 
00521 }  // namespace MBMesquite
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines