Mesh Oriented datABase  (version 5.4.1)
Array-based unstructured mesh datastructure
compareMaps.cpp
Go to the documentation of this file.
00001 /*
00002  * compareMaps.cpp
00003  * this tool will take 2 existing map files in nc format, and compare their sparse matrices
00004  * we can compare xc, yc, areas, fractions with ncdiff from nco
00005  * maybe there is another utility in nco, need to ask Charlie Zender
00006  *
00007  * example of usage:
00008  * ./mbcmpmaps -i map1.nc -j map2.nc
00009  * will look for row, col, S entries, and use eigen3 sparse matrix constructor
00010  *
00011  * can be built only if netcdf and eigen3 are available
00012  *
00013  *
00014  */
00015 #include "moab/MOABConfig.h"
00016 
00017 #ifndef MOAB_HAVE_EIGEN3
00018 #error compareMaps tool requires eigen3 configuration
00019 #endif
00020 
00021 #include "moab/ProgOptions.hpp"
00022 
00023 #include "netcdf.h"
00024 #include <cmath>
00025 #include <iomanip>
00026 #include <Eigen/Sparse>
00027 #define ERR_NC( e )                                \
00028     {                                              \
00029         printf( "Error: %s\n", nc_strerror( e ) ); \
00030         exit( 2 );                                 \
00031     }
00032 
00033 // copy from ReadNCDF.cpp some useful macros for reading from a netcdf file (exodus?)
00034 // ncFile is an integer initialized when opening the nc file in read mode
00035 
00036 int ncFile1;
00037 
00038 #define INS_ID( stringvar, prefix, id ) sprintf( stringvar, prefix, id )
00039 
00040 #define GET_DIM1( ncdim, name, val )                            \
00041     {                                                           \
00042         int gdfail = nc_inq_dimid( ncFile1, name, &( ncdim ) ); \
00043         if( NC_NOERR == gdfail )                                \
00044         {                                                       \
00045             size_t tmp_val;                                     \
00046             gdfail = nc_inq_dimlen( ncFile1, ncdim, &tmp_val ); \
00047             if( NC_NOERR != gdfail )                            \
00048             {                                                   \
00049                 ERR_NC( gdfail )                                \
00050             }                                                   \
00051             else                                                \
00052                 ( val ) = tmp_val;                              \
00053         }                                                       \
00054         else                                                    \
00055             ( val ) = 0;                                        \
00056     }
00057 
00058 #define GET_DIMB1( ncdim, name, varname, id, val ) \
00059     INS_ID( name, varname, id );                   \
00060     GET_DIM1( ncdim, name, val );
00061 
00062 #define GET_VAR1( name, id, dims )                                     \
00063     {                                                                  \
00064         ( id )     = -1;                                               \
00065         int gvfail = nc_inq_varid( ncFile1, name, &( id ) );           \
00066         if( NC_NOERR == gvfail )                                       \
00067         {                                                              \
00068             int ndims;                                                 \
00069             gvfail = nc_inq_varndims( ncFile1, id, &ndims );           \
00070             if( NC_NOERR == gvfail )                                   \
00071             {                                                          \
00072                 ( dims ).resize( ndims );                              \
00073                 gvfail = nc_inq_vardimid( ncFile1, id, &( dims )[0] ); \
00074                 if( NC_NOERR != gvfail )                               \
00075                 {                                                      \
00076                     ERR_NC( gvfail )                                   \
00077                 }                                                      \
00078             }                                                          \
00079         }                                                              \
00080     }
00081 
00082 #define GET_1D_INT_VAR1( name, id, vals )                                               \
00083     {                                                                                   \
00084         GET_VAR1( name, id, vals );                                                     \
00085         if( -1 != ( id ) )                                                              \
00086         {                                                                               \
00087             size_t ntmp;                                                                \
00088             int ivfail = nc_inq_dimlen( ncFile1, ( vals )[0], &ntmp );                  \
00089             if( NC_NOERR != ivfail )                                                    \
00090             {                                                                           \
00091                 ERR_NC( ivfail )                                                        \
00092             }                                                                           \
00093             ( vals ).resize( ntmp );                                                    \
00094             size_t ntmp1 = 0;                                                           \
00095             ivfail       = nc_get_vara_int( ncFile1, id, &ntmp1, &ntmp, &( vals )[0] ); \
00096             if( NC_NOERR != ivfail )                                                    \
00097             {                                                                           \
00098                 ERR_NC( ivfail )                                                        \
00099             }                                                                           \
00100         }                                                                               \
00101     }
00102 
00103 #define GET_1D_DBL_VAR1( name, id, vals )                                                  \
00104     {                                                                                      \
00105         std::vector< int > dum_dims;                                                       \
00106         GET_VAR1( name, id, dum_dims );                                                    \
00107         if( -1 != ( id ) )                                                                 \
00108         {                                                                                  \
00109             size_t ntmp;                                                                   \
00110             int dvfail = nc_inq_dimlen( ncFile1, dum_dims[0], &ntmp );                     \
00111             if( NC_NOERR != dvfail )                                                       \
00112             {                                                                              \
00113                 ERR_NC( dvfail )                                                           \
00114             }                                                                              \
00115             ( vals ).resize( ntmp );                                                       \
00116             size_t ntmp1 = 0;                                                              \
00117             dvfail       = nc_get_vara_double( ncFile1, id, &ntmp1, &ntmp, &( vals )[0] ); \
00118             if( NC_NOERR != dvfail )                                                       \
00119             {                                                                              \
00120                 ERR_NC( dvfail )                                                           \
00121             }                                                                              \
00122         }                                                                                  \
00123     }
00124 
00125 int ncFile2;
00126 
00127 #define GET_DIM2( ncdim, name, val )                            \
00128     {                                                           \
00129         int gdfail = nc_inq_dimid( ncFile2, name, &( ncdim ) ); \
00130         if( NC_NOERR == gdfail )                                \
00131         {                                                       \
00132             size_t tmp_val;                                     \
00133             gdfail = nc_inq_dimlen( ncFile2, ncdim, &tmp_val ); \
00134             if( NC_NOERR != gdfail )                            \
00135             {                                                   \
00136                 ERR_NC( gdfail )                                \
00137             }                                                   \
00138             else                                                \
00139                 ( val ) = tmp_val;                              \
00140         }                                                       \
00141         else                                                    \
00142             ( val ) = 0;                                        \
00143     }
00144 
00145 #define GET_DIMB2( ncdim, name, varname, id, val ) \
00146     INS_ID( name, varname, id );                   \
00147     GET_DIM2( ncdim, name, val );
00148 
00149 #define GET_VAR2( name, id, dims )                                     \
00150     {                                                                  \
00151         ( id )     = -1;                                               \
00152         int gvfail = nc_inq_varid( ncFile2, name, &( id ) );           \
00153         if( NC_NOERR == gvfail )                                       \
00154         {                                                              \
00155             int ndims;                                                 \
00156             gvfail = nc_inq_varndims( ncFile2, id, &ndims );           \
00157             if( NC_NOERR == gvfail )                                   \
00158             {                                                          \
00159                 ( dims ).resize( ndims );                              \
00160                 gvfail = nc_inq_vardimid( ncFile2, id, &( dims )[0] ); \
00161                 if( NC_NOERR != gvfail )                               \
00162                 {                                                      \
00163                     ERR_NC( gvfail )                                   \
00164                 }                                                      \
00165             }                                                          \
00166         }                                                              \
00167     }
00168 
00169 #define GET_1D_INT_VAR2( name, id, vals )                                               \
00170     {                                                                                   \
00171         GET_VAR2( name, id, vals );                                                     \
00172         if( -1 != ( id ) )                                                              \
00173         {                                                                               \
00174             size_t ntmp;                                                                \
00175             int ivfail = nc_inq_dimlen( ncFile2, ( vals )[0], &ntmp );                  \
00176             if( NC_NOERR != ivfail )                                                    \
00177             {                                                                           \
00178                 ERR_NC( ivfail )                                                        \
00179             }                                                                           \
00180             ( vals ).resize( ntmp );                                                    \
00181             size_t ntmp1 = 0;                                                           \
00182             ivfail       = nc_get_vara_int( ncFile2, id, &ntmp1, &ntmp, &( vals )[0] ); \
00183             if( NC_NOERR != ivfail )                                                    \
00184             {                                                                           \
00185                 ERR_NC( ivfail )                                                        \
00186             }                                                                           \
00187         }                                                                               \
00188     }
00189 
00190 #define GET_1D_DBL_VAR2( name, id, vals )                                                  \
00191     {                                                                                      \
00192         std::vector< int > dum_dims;                                                       \
00193         GET_VAR2( name, id, dum_dims );                                                    \
00194         if( -1 != ( id ) )                                                                 \
00195         {                                                                                  \
00196             size_t ntmp;                                                                   \
00197             int dvfail = nc_inq_dimlen( ncFile2, dum_dims[0], &ntmp );                     \
00198             if( NC_NOERR != dvfail )                                                       \
00199             {                                                                              \
00200                 ERR_NC( dvfail )                                                           \
00201             }                                                                              \
00202             ( vals ).resize( ntmp );                                                       \
00203             size_t ntmp1 = 0;                                                              \
00204             dvfail       = nc_get_vara_double( ncFile2, id, &ntmp1, &ntmp, &( vals )[0] ); \
00205             if( NC_NOERR != dvfail )                                                       \
00206             {                                                                              \
00207                 ERR_NC( dvfail )                                                           \
00208             }                                                                              \
00209         }                                                                                  \
00210     }
00211 
00212 #define GET_2D_DBL_VAR1( name, id, vals )                                                   \
00213     {                                                                                       \
00214         std::vector< int > dum_dims;                                                        \
00215         GET_VAR1( name, id, dum_dims );                                                     \
00216         if( -1 != ( id ) )                                                                  \
00217         {                                                                                   \
00218             size_t ntmp[2];                                                                 \
00219             int dvfail = nc_inq_dimlen( ncFile1, dum_dims[0], &ntmp[0] );                   \
00220             if( NC_NOERR != dvfail )                                                        \
00221             {                                                                               \
00222                 ERR_NC( dvfail )                                                            \
00223             }                                                                               \
00224             dvfail = nc_inq_dimlen( ncFile1, dum_dims[1], &ntmp[1] );                       \
00225             if( NC_NOERR != dvfail )                                                        \
00226             {                                                                               \
00227                 ERR_NC( dvfail )                                                            \
00228             }                                                                               \
00229             ( vals ).resize( ntmp[0] * ntmp[1] );                                           \
00230             size_t ntmp1[2] = { 0, 0 };                                                     \
00231             dvfail          = nc_get_vara_double( ncFile1, id, ntmp1, ntmp, &( vals )[0] ); \
00232             if( NC_NOERR != dvfail )                                                        \
00233             {                                                                               \
00234                 ERR_NC( dvfail )                                                            \
00235             }                                                                               \
00236         }                                                                                   \
00237     }
00238 
00239 #define GET_2D_DBL_VAR2( name, id, vals )                                                   \
00240     {                                                                                       \
00241         std::vector< int > dum_dims;                                                        \
00242         GET_VAR2( name, id, dum_dims );                                                     \
00243         if( -1 != ( id ) )                                                                  \
00244         {                                                                                   \
00245             size_t ntmp[2];                                                                 \
00246             int dvfail = nc_inq_dimlen( ncFile2, dum_dims[0], &ntmp[0] );                   \
00247             if( NC_NOERR != dvfail )                                                        \
00248             {                                                                               \
00249                 ERR_NC( dvfail )                                                            \
00250             }                                                                               \
00251             dvfail = nc_inq_dimlen( ncFile2, dum_dims[1], &ntmp[1] );                       \
00252             if( NC_NOERR != dvfail )                                                        \
00253             {                                                                               \
00254                 ERR_NC( dvfail )                                                            \
00255             }                                                                               \
00256             ( vals ).resize( ntmp[0] * ntmp[1] );                                           \
00257             size_t ntmp1[2] = { 0, 0 };                                                     \
00258             dvfail          = nc_get_vara_double( ncFile2, id, ntmp1, ntmp, &( vals )[0] ); \
00259             if( NC_NOERR != dvfail )                                                        \
00260             {                                                                               \
00261                 ERR_NC( dvfail )                                                            \
00262             }                                                                               \
00263         }                                                                                   \
00264     }
00265 
00266 typedef Eigen::Map< Eigen::VectorXd > EigenV;
00267 
00268 void diff_vect( const char* var_name, int n )
00269 {
00270     // compare frac_a between maps
00271     std::vector< double > fraca1( n ), fraca2( n );
00272     int idfa1, idfa2;
00273     GET_1D_DBL_VAR1( var_name, idfa1, fraca1 );
00274     EigenV fa1( fraca1.data(), n );
00275     GET_1D_DBL_VAR2( var_name, idfa2, fraca2 );
00276     EigenV fa2( fraca2.data(), n );
00277 
00278     EigenV diff( fraca2.data(), n );
00279     diff = fa1 - fa2;
00280 
00281     int imin, imax;
00282     double minV = diff.minCoeff( &imin );
00283     double maxV = diff.maxCoeff( &imax );
00284     std::cout << var_name << " diff norm: " << diff.norm() << " min at " << imin << " : " << minV << " max at " << imax
00285               << " : " << maxV << "\n";
00286     return;
00287 }
00288 
00289 void diff_2d_vect( const char* var_name, int n )
00290 {
00291     // compare frac_a between maps
00292     std::vector< double > fraca1( n ), fraca2( n );
00293     int idfa1, idfa2;
00294     GET_2D_DBL_VAR1( var_name, idfa1, fraca1 );
00295     EigenV fa1( fraca1.data(), n );
00296     GET_2D_DBL_VAR2( var_name, idfa2, fraca2 );
00297     EigenV fa2( fraca2.data(), n );
00298     std::cout << var_name << " diff norm: " << ( fa1 - fa2 ).norm() << "\n";
00299     return;
00300 }
00301 int main( int argc, char* argv[] )
00302 {
00303 
00304     ProgOptions opts;
00305 
00306     std::string inputfile1, inputfile2;
00307     opts.addOpt< std::string >( "firstMap,i", "input filename 1", &inputfile1 );
00308     opts.addOpt< std::string >( "secondMap,j", "input second map", &inputfile2 );
00309     int print_diff = 20;
00310     opts.addOpt< int >( "print_differences,p", "print differences ", &print_diff );
00311     double fraction = 0.9;
00312     opts.addOpt< double >( "fraction_diff,f", "fraction threshold", &fraction );
00313 
00314     opts.parseCommandLine( argc, argv );
00315 
00316     // Open netcdf/exodus file
00317     int fail = nc_open( inputfile1.c_str(), 0, &ncFile1 );
00318     if( NC_NOWRITE != fail )
00319     {
00320         ERR_NC( fail )
00321     }
00322 
00323     std::cout << " opened " << inputfile1 << " for map 1 \n";
00324 
00325     fail = nc_open( inputfile2.c_str(), 0, &ncFile2 );
00326     if( NC_NOWRITE != fail )
00327     {
00328         ERR_NC( fail )
00329     }
00330 
00331     std::cout << " opened " << inputfile2 << " for map 2 \n";
00332     int temp_dim;
00333     int na1, nb1, ns1, na2, nb2, ns2, nv_a, nv_b, nv_a2, nv_b2;
00334     GET_DIM1( temp_dim, "n_a", na1 );
00335     GET_DIM2( temp_dim, "n_a", na2 );
00336     GET_DIM1( temp_dim, "n_b", nb1 );
00337     GET_DIM2( temp_dim, "n_b", nb2 );
00338     GET_DIM1( temp_dim, "n_s", ns1 );
00339     GET_DIM2( temp_dim, "n_s", ns2 );
00340     GET_DIM1( temp_dim, "nv_a", nv_a );
00341     GET_DIM2( temp_dim, "nv_a", nv_a2 );
00342     GET_DIM1( temp_dim, "nv_b", nv_b );
00343     GET_DIM2( temp_dim, "nv_b", nv_b2 );
00344     if( nv_a != nv_a2 || nv_b != nv_b2 )
00345     {
00346         std::cout << " different nv dimensions:" << nv_a << " == " << nv_a2 << "  or " << nv_b << " == " << nv_b2
00347                   << "  bail out \n";
00348         return 1;
00349     }
00350     std::cout << " n_a, n_b, n_s : " << na1 << ", " << nb1 << ", " << ns1 << " for map 1 \n";
00351 
00352     std::cout << " n_a, n_b, n_s : " << na2 << ", " << nb2 << ", " << ns2 << " for map 2 \n";
00353     if( na1 != na2 || nb1 != nb2 )
00354     {
00355         std::cout << " different dimensions bail out \n";
00356         return 1;
00357     }
00358     std::vector< int > col1( ns1 ), row1( ns1 );
00359     std::vector< int > col2( ns2 ), row2( ns2 );
00360 
00361     std::vector< double > val1( ns1 ), val2( ns2 );
00362 
00363     int idrow1, idcol1, idrow2, idcol2, ids1, ids2;
00364     GET_1D_INT_VAR1( "row", idrow1, row1 );
00365     GET_1D_INT_VAR1( "col", idcol1, col1 );
00366     GET_1D_DBL_VAR1( "S", ids1, val1 );
00367     GET_1D_INT_VAR2( "row", idrow2, row2 );
00368     GET_1D_INT_VAR2( "col", idcol2, col2 );
00369     GET_1D_DBL_VAR2( "S", ids2, val2 );
00370 
00371     // first matrix
00372     typedef Eigen::Triplet< double > Triplet;
00373     std::vector< Triplet > tripletList;
00374     tripletList.reserve( ns1 );
00375     for( int iv = 0; iv < ns1; iv++ )
00376     {
00377         tripletList.push_back( Triplet( row1[iv] - 1, col1[iv] - 1, val1[iv] ) );
00378     }
00379     Eigen::SparseMatrix< double > weight1( nb1, na1 );
00380 
00381     weight1.setFromTriplets( tripletList.begin(), tripletList.end() );
00382     weight1.makeCompressed();
00383 
00384     if( ns1 != ns2 ) tripletList.resize( ns2 );
00385     for( int iv = 0; iv < ns2; iv++ )
00386     {
00387         tripletList[iv] = Triplet( row2[iv] - 1, col2[iv] - 1, val2[iv] );
00388     }
00389     Eigen::SparseMatrix< double > weight2( nb2, na2 );
00390 
00391     weight2.setFromTriplets( tripletList.begin(), tripletList.end() );
00392     weight2.makeCompressed();
00393 
00394     // default storage type is column major
00395     Eigen::SparseMatrix< double > diff = weight1 - weight2;
00396     diff.makeCompressed();  // is it needed or not ?
00397     auto coeffs = diff.coeffs();
00398     double maxv = coeffs.maxCoeff();
00399     double minv = coeffs.minCoeff();
00400     std::cout << " euclidian norm for difference: " << diff.norm()
00401               << " \n squared norm for difference: " << diff.squaredNorm() << "\n"
00402               << " minv: " << minv << " maxv: " << maxv << "\n";
00403     // print out the first 10 positions for which the value is outside 90% of min/max values
00404     double min_threshold = fraction * minv;
00405     double max_threshold = fraction * maxv;
00406     int counter          = 0;
00407     std::cout << std::setprecision( 9 );
00408     for( int k = 0; ( k < diff.outerSize() ); ++k )  // this is by column
00409     {
00410         for( Eigen::SparseMatrix< double >::InnerIterator it( diff, k ); ( it ) && ( counter < print_diff ); ++it )
00411         {
00412             double val = it.value();
00413             if( val <= min_threshold || val >= max_threshold )
00414             {
00415                 int row = it.row();
00416                 int col = it.col();
00417                 std::cout << " counter:" << counter << "\t col: " << col + 1 << "\t row: " << row + 1
00418                           << "\t diff: " << val;
00419                 std::cout << "\t map1: " << weight1.coeffRef( row, col ) << "\t map2: " << weight2.coeffRef( row, col )
00420                           << "\n";  // row index
00421                 counter++;
00422             }
00423         }
00424     }
00425 
00426     // compare frac_a between maps
00427     diff_vect( "frac_a", na1 );
00428     diff_vect( "frac_b", nb1 );
00429     diff_vect( "area_a", na1 );
00430     diff_vect( "area_b", nb1 );
00431     diff_vect( "yc_a", na1 );
00432     diff_vect( "yc_b", nb1 );
00433     diff_vect( "xc_a", na1 );
00434     diff_vect( "xc_b", nb1 );
00435     diff_2d_vect( "xv_a", na1 * nv_a );
00436     diff_2d_vect( "yv_a", na1 * nv_a );
00437     diff_2d_vect( "xv_b", nb1 * nv_b );
00438     diff_2d_vect( "yv_b", nb1 * nv_b );
00439 
00440     return 0;
00441 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines