MOAB: Mesh Oriented datABase  (version 5.3.1)
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 map1.nc 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 <sstream>
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 ) { ERR_NC( gdfail ) }       \
00048             else                                                \
00049                 ( val ) = tmp_val;                              \
00050         }                                                       \
00051         else                                                    \
00052             ( val ) = 0;                                        \
00053     }
00054 
00055 #define GET_DIMB1( ncdim, name, varname, id, val ) \
00056     INS_ID( name, varname, id );                   \
00057     GET_DIM1( ncdim, name, val );
00058 
00059 #define GET_VAR1( name, id, dims )                                     \
00060     {                                                                  \
00061         ( id )     = -1;                                               \
00062         int gvfail = nc_inq_varid( ncFile1, name, &( id ) );           \
00063         if( NC_NOERR == gvfail )                                       \
00064         {                                                              \
00065             int ndims;                                                 \
00066             gvfail = nc_inq_varndims( ncFile1, id, &ndims );           \
00067             if( NC_NOERR == gvfail )                                   \
00068             {                                                          \
00069                 ( dims ).resize( ndims );                              \
00070                 gvfail = nc_inq_vardimid( ncFile1, id, &( dims )[0] ); \
00071                 if( NC_NOERR != gvfail ) { ERR_NC( gvfail ) }          \
00072             }                                                          \
00073         }                                                              \
00074     }
00075 
00076 #define GET_1D_INT_VAR1( name, id, vals )                                               \
00077     {                                                                                   \
00078         GET_VAR1( name, id, vals );                                                     \
00079         if( -1 != ( id ) )                                                              \
00080         {                                                                               \
00081             size_t ntmp;                                                                \
00082             int ivfail = nc_inq_dimlen( ncFile1, ( vals )[0], &ntmp );                  \
00083             if( NC_NOERR != ivfail ) { ERR_NC( ivfail ) }                               \
00084             ( vals ).resize( ntmp );                                                    \
00085             size_t ntmp1 = 0;                                                           \
00086             ivfail       = nc_get_vara_int( ncFile1, id, &ntmp1, &ntmp, &( vals )[0] ); \
00087             if( NC_NOERR != ivfail ) { ERR_NC( ivfail ) }                               \
00088         }                                                                               \
00089     }
00090 
00091 #define GET_1D_DBL_VAR1( name, id, vals )                                                  \
00092     {                                                                                      \
00093         std::vector< int > dum_dims;                                                       \
00094         GET_VAR1( name, id, dum_dims );                                                    \
00095         if( -1 != ( id ) )                                                                 \
00096         {                                                                                  \
00097             size_t ntmp;                                                                   \
00098             int dvfail = nc_inq_dimlen( ncFile1, dum_dims[0], &ntmp );                     \
00099             if( NC_NOERR != dvfail ) { ERR_NC( dvfail ) }                                  \
00100             ( vals ).resize( ntmp );                                                       \
00101             size_t ntmp1 = 0;                                                              \
00102             dvfail       = nc_get_vara_double( ncFile1, id, &ntmp1, &ntmp, &( vals )[0] ); \
00103             if( NC_NOERR != dvfail ) { ERR_NC( dvfail ) }                                  \
00104         }                                                                                  \
00105     }
00106 
00107 int ncFile2;
00108 
00109 #define GET_DIM2( ncdim, name, val )                            \
00110     {                                                           \
00111         int gdfail = nc_inq_dimid( ncFile2, name, &( ncdim ) ); \
00112         if( NC_NOERR == gdfail )                                \
00113         {                                                       \
00114             size_t tmp_val;                                     \
00115             gdfail = nc_inq_dimlen( ncFile2, ncdim, &tmp_val ); \
00116             if( NC_NOERR != gdfail ) { ERR_NC( gdfail ) }       \
00117             else                                                \
00118                 ( val ) = tmp_val;                              \
00119         }                                                       \
00120         else                                                    \
00121             ( val ) = 0;                                        \
00122     }
00123 
00124 #define GET_DIMB2( ncdim, name, varname, id, val ) \
00125     INS_ID( name, varname, id );                   \
00126     GET_DIM2( ncdim, name, val );
00127 
00128 #define GET_VAR2( name, id, dims )                                     \
00129     {                                                                  \
00130         ( id )     = -1;                                               \
00131         int gvfail = nc_inq_varid( ncFile2, name, &( id ) );           \
00132         if( NC_NOERR == gvfail )                                       \
00133         {                                                              \
00134             int ndims;                                                 \
00135             gvfail = nc_inq_varndims( ncFile2, id, &ndims );           \
00136             if( NC_NOERR == gvfail )                                   \
00137             {                                                          \
00138                 ( dims ).resize( ndims );                              \
00139                 gvfail = nc_inq_vardimid( ncFile2, id, &( dims )[0] ); \
00140                 if( NC_NOERR != gvfail ) { ERR_NC( gvfail ) }          \
00141             }                                                          \
00142         }                                                              \
00143     }
00144 
00145 #define GET_1D_INT_VAR2( name, id, vals )                                               \
00146     {                                                                                   \
00147         GET_VAR2( name, id, vals );                                                     \
00148         if( -1 != ( id ) )                                                              \
00149         {                                                                               \
00150             size_t ntmp;                                                                \
00151             int ivfail = nc_inq_dimlen( ncFile2, ( vals )[0], &ntmp );                  \
00152             if( NC_NOERR != ivfail ) { ERR_NC( ivfail ) }                               \
00153             ( vals ).resize( ntmp );                                                    \
00154             size_t ntmp1 = 0;                                                           \
00155             ivfail       = nc_get_vara_int( ncFile2, id, &ntmp1, &ntmp, &( vals )[0] ); \
00156             if( NC_NOERR != ivfail ) { ERR_NC( ivfail ) }                               \
00157         }                                                                               \
00158     }
00159 
00160 #define GET_1D_DBL_VAR2( name, id, vals )                                                  \
00161     {                                                                                      \
00162         std::vector< int > dum_dims;                                                       \
00163         GET_VAR2( name, id, dum_dims );                                                    \
00164         if( -1 != ( id ) )                                                                 \
00165         {                                                                                  \
00166             size_t ntmp;                                                                   \
00167             int dvfail = nc_inq_dimlen( ncFile2, dum_dims[0], &ntmp );                     \
00168             if( NC_NOERR != dvfail ) { ERR_NC( dvfail ) }                                  \
00169             ( vals ).resize( ntmp );                                                       \
00170             size_t ntmp1 = 0;                                                              \
00171             dvfail       = nc_get_vara_double( ncFile2, id, &ntmp1, &ntmp, &( vals )[0] ); \
00172             if( NC_NOERR != dvfail ) { ERR_NC( dvfail ) }                                  \
00173         }                                                                                  \
00174     }
00175 
00176 #define GET_2D_DBL_VAR1( name, id, vals )                                                   \
00177     {                                                                                       \
00178         std::vector< int > dum_dims;                                                        \
00179         GET_VAR1( name, id, dum_dims );                                                     \
00180         if( -1 != ( id ) )                                                                  \
00181         {                                                                                   \
00182             size_t ntmp[2];                                                                 \
00183             int dvfail = nc_inq_dimlen( ncFile1, dum_dims[0], &ntmp[0] );                   \
00184             if( NC_NOERR != dvfail ) { ERR_NC( dvfail ) }                                   \
00185             dvfail = nc_inq_dimlen( ncFile1, dum_dims[1], &ntmp[1] );                       \
00186             if( NC_NOERR != dvfail ) { ERR_NC( dvfail ) }                                   \
00187             ( vals ).resize( ntmp[0] * ntmp[1] );                                           \
00188             size_t ntmp1[2] = { 0, 0 };                                                     \
00189             dvfail          = nc_get_vara_double( ncFile1, id, ntmp1, ntmp, &( vals )[0] ); \
00190             if( NC_NOERR != dvfail ) { ERR_NC( dvfail ) }                                   \
00191         }                                                                                   \
00192     }
00193 
00194 #define GET_2D_DBL_VAR2( name, id, vals )                                                   \
00195     {                                                                                       \
00196         std::vector< int > dum_dims;                                                        \
00197         GET_VAR2( name, id, dum_dims );                                                     \
00198         if( -1 != ( id ) )                                                                  \
00199         {                                                                                   \
00200             size_t ntmp[2];                                                                 \
00201             int dvfail = nc_inq_dimlen( ncFile2, dum_dims[0], &ntmp[0] );                   \
00202             if( NC_NOERR != dvfail ) { ERR_NC( dvfail ) }                                   \
00203             dvfail = nc_inq_dimlen( ncFile2, dum_dims[1], &ntmp[1] );                       \
00204             if( NC_NOERR != dvfail ) { ERR_NC( dvfail ) }                                   \
00205             ( vals ).resize( ntmp[0] * ntmp[1] );                                           \
00206             size_t ntmp1[2] = { 0, 0 };                                                     \
00207             dvfail          = nc_get_vara_double( ncFile2, id, ntmp1, ntmp, &( vals )[0] ); \
00208             if( NC_NOERR != dvfail ) { ERR_NC( dvfail ) }                                   \
00209         }                                                                                   \
00210     }
00211 
00212 typedef Eigen::Map< Eigen::VectorXd > EigenV;
00213 
00214 void diff_vect( const char* var_name, int n )
00215 {
00216     // compare frac_a between maps
00217     std::vector< double > fraca1( n ), fraca2( n );
00218     int idfa1, idfa2;
00219     GET_1D_DBL_VAR1( var_name, idfa1, fraca1 );
00220     EigenV fa1( fraca1.data(), n );
00221     GET_1D_DBL_VAR2( var_name, idfa2, fraca2 );
00222     EigenV fa2( fraca2.data(), n );
00223 
00224     EigenV diff( fraca2.data(), n );
00225     diff = fa1 - fa2;
00226 
00227     int imin, imax;
00228     double minV = diff.minCoeff( &imin );
00229     double maxV = diff.maxCoeff( &imax );
00230     std::cout << var_name << " diff norm: " << diff.norm() << " min at " << imin << " : " << minV << " max at " << imax
00231               << " : " << maxV << "\n";
00232     return;
00233 }
00234 
00235 void diff_2d_vect( const char* var_name, int n )
00236 {
00237     // compare frac_a between maps
00238     std::vector< double > fraca1( n ), fraca2( n );
00239     int idfa1, idfa2;
00240     GET_2D_DBL_VAR1( var_name, idfa1, fraca1 );
00241     EigenV fa1( fraca1.data(), n );
00242     GET_2D_DBL_VAR2( var_name, idfa2, fraca2 );
00243     EigenV fa2( fraca2.data(), n );
00244     std::cout << var_name << " diff norm: " << ( fa1 - fa2 ).norm() << "\n";
00245     return;
00246 }
00247 int main( int argc, char* argv[] )
00248 {
00249 
00250     ProgOptions opts;
00251 
00252     std::string inputfile1, inputfile2;
00253     opts.addOpt< std::string >( "firstMap,i", "input filename 1", &inputfile1 );
00254     opts.addOpt< std::string >( "secondMap,j", "input second map", &inputfile2 );
00255 
00256     opts.parseCommandLine( argc, argv );
00257 
00258     // Open netcdf/exodus file
00259     int fail = nc_open( inputfile1.c_str(), 0, &ncFile1 );
00260     if( NC_NOWRITE != fail ) { ERR_NC( fail ) }
00261 
00262     std::cout << " opened " << inputfile1 << " for map 1 \n";
00263 
00264     fail = nc_open( inputfile2.c_str(), 0, &ncFile2 );
00265     if( NC_NOWRITE != fail ) { ERR_NC( fail ) }
00266 
00267     std::cout << " opened " << inputfile2 << " for map 2 \n";
00268     int temp_dim;
00269     int na1, nb1, ns1, na2, nb2, ns2, nv_a, nv_b, nv_a2, nv_b2;
00270     GET_DIM1( temp_dim, "n_a", na1 );
00271     GET_DIM2( temp_dim, "n_a", na2 );
00272     GET_DIM1( temp_dim, "n_b", nb1 );
00273     GET_DIM2( temp_dim, "n_b", nb2 );
00274     GET_DIM1( temp_dim, "n_s", ns1 );
00275     GET_DIM2( temp_dim, "n_s", ns2 );
00276     GET_DIM1( temp_dim, "nv_a", nv_a );
00277     GET_DIM2( temp_dim, "nv_a", nv_a2 );
00278     GET_DIM1( temp_dim, "nv_b", nv_b );
00279     GET_DIM2( temp_dim, "nv_b", nv_b2 );
00280     if( nv_a != nv_a2 || nv_b != nv_b2 )
00281     {
00282         std::cout << " different nv dimensions:" << nv_a << " == " << nv_a2 << "  or " << nv_b << " == " << nv_b2
00283                   << "  bail out \n";
00284         return 1;
00285     }
00286     std::cout << " n_a, n_b, n_s : " << na1 << ", " << nb1 << ", " << ns1 << " for map 1 \n";
00287 
00288     std::cout << " n_a, n_b, n_s : " << na2 << ", " << nb2 << ", " << ns2 << " for map 2 \n";
00289     if( na1 != na2 || nb1 != nb2 )
00290     {
00291         std::cout << " different dimensions bail out \n";
00292         return 1;
00293     }
00294     std::vector< int > col1( ns1 ), row1( ns1 );
00295     std::vector< int > col2( ns2 ), row2( ns2 );
00296 
00297     std::vector< double > val1( ns1 ), val2( ns2 );
00298 
00299     int idrow1, idcol1, idrow2, idcol2, ids1, ids2;
00300     GET_1D_INT_VAR1( "row", idrow1, row1 );
00301     GET_1D_INT_VAR1( "col", idcol1, col1 );
00302     GET_1D_DBL_VAR1( "S", ids1, val1 );
00303     GET_1D_INT_VAR2( "row", idrow2, row2 );
00304     GET_1D_INT_VAR2( "col", idcol2, col2 );
00305     GET_1D_DBL_VAR2( "S", ids2, val2 );
00306 
00307     // first matrix
00308     typedef Eigen::Triplet< double > Triplet;
00309     std::vector< Triplet > tripletList;
00310     tripletList.reserve( ns1 );
00311     for( int iv = 0; iv < ns1; iv++ )
00312     {
00313         tripletList.push_back( Triplet( row1[iv] - 1, col1[iv] - 1, val1[iv] ) );
00314     }
00315     Eigen::SparseMatrix< double > weight1( nb1, na1 );
00316 
00317     weight1.setFromTriplets( tripletList.begin(), tripletList.end() );
00318     weight1.makeCompressed();
00319 
00320     if( ns1 != ns2 ) tripletList.resize( ns2 );
00321     for( int iv = 0; iv < ns2; iv++ )
00322     {
00323         tripletList[iv] = Triplet( row2[iv] - 1, col2[iv] - 1, val2[iv] );
00324     }
00325     Eigen::SparseMatrix< double > weight2( nb2, na2 );
00326 
00327     weight2.setFromTriplets( tripletList.begin(), tripletList.end() );
00328     weight2.makeCompressed();
00329 
00330     Eigen::SparseMatrix< double > diff = weight1 - weight2;
00331     double maxv                        = diff.coeffs().maxCoeff();
00332     double minv                        = diff.coeffs().minCoeff();
00333     std::cout << " euclidian norm for difference: " << diff.norm()
00334               << " \n squared norm for difference: " << diff.squaredNorm() << "\n"
00335               << " minv: " << minv << " maxv: " << maxv << "\n";
00336 
00337     // compare frac_a between maps
00338     diff_vect( "frac_a", na1 );
00339     diff_vect( "frac_b", nb1 );
00340     diff_vect( "area_a", na1 );
00341     diff_vect( "area_b", nb1 );
00342     diff_vect( "yc_a", na1 );
00343     diff_vect( "yc_b", nb1 );
00344     diff_vect( "xc_a", na1 );
00345     diff_vect( "xc_b", nb1 );
00346     diff_2d_vect( "xv_a", na1 * nv_a );
00347     diff_2d_vect( "yv_a", na1 * nv_a );
00348     diff_2d_vect( "xv_b", nb1 * nv_b );
00349     diff_2d_vect( "yv_b", nb1 * nv_b );
00350 
00351     return 0;
00352 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines