|
MOAB: Mesh Oriented datABase
(version 5.4.1)
|
|
MOAB: Mesh Oriented datABase
(version 5.4.1)
|
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 }
1.7.6.1 
