 |
Mesh Oriented datABase
(version 5.4.1)
Array-based unstructured mesh datastructure
|
|
Mesh Oriented datABase
(version 5.4.1)
Array-based unstructured mesh datastructure
|
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
00025 #include
00026 #include
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 }
1.7.6.1.
Dark theme by Tilen Majerle. All rights reserved.