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