 |
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
|
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
#include "moab/MOABConfig.h"#include "moab/ProgOptions.hpp"#include "netcdf.h"#include <cmath>#include <iomanip>#include <Eigen/Sparse>
Include dependency graph for compareMaps.cpp:Go to the source code of this file.
Macros | |
| #define | ERR_NC(e) |
| #define | GET_DIM1(ncdim, name, val) |
| #define | GET_VAR1(name, id, dims) |
| #define | GET_1D_INT_VAR1(name, id, vals) |
| #define | GET_1D_DBL_VAR1(name, id, vals) |
| #define | GET_DIM2(ncdim, name, val) |
| #define | GET_VAR2(name, id, dims) |
| #define | GET_1D_INT_VAR2(name, id, vals) |
| #define | GET_1D_DBL_VAR2(name, id, vals) |
| #define | GET_2D_DBL_VAR1(name, id, vals) |
| #define | GET_2D_DBL_VAR2(name, id, vals) |
Typedefs | |
| typedef Eigen::Map< Eigen::VectorXd > | EigenV |
Functions | |
| void | diff_vect (const char *var_name, int n) |
| void | diff_2d_vect (const char *var_name, int n) |
| int | main (int argc, char *argv[]) |
Variables | |
| int | ncFile1 |
| int | ncFile2 |
| #define ERR_NC | ( | e | ) |
{ \
printf( "Error: %s\n", nc_strerror( e ) ); \
exit( 2 ); \
}
Definition at line 27 of file compareMaps.cpp.
| #define GET_1D_DBL_VAR1 | ( | name, | |
| id, | |||
| vals | |||
| ) |
{ \
std::vector< int > dum_dims; \
GET_VAR1( name, id, dum_dims ); \
if( -1 != ( id ) ) \
{ \
size_t ntmp; \
int dvfail = nc_inq_dimlen( ncFile1, dum_dims[0], &ntmp ); \
if( NC_NOERR != dvfail ) \
{ \
ERR_NC( dvfail ) \
} \
( vals ).resize( ntmp ); \
size_t ntmp1 = 0; \
dvfail = nc_get_vara_double( ncFile1, id, &ntmp1, &ntmp, &( vals )[0] ); \
if( NC_NOERR != dvfail ) \
{ \
ERR_NC( dvfail ) \
} \
} \
}
Definition at line 97 of file compareMaps.cpp.
| #define GET_1D_DBL_VAR2 | ( | name, | |
| id, | |||
| vals | |||
| ) |
{ \
std::vector< int > dum_dims; \
GET_VAR2( name, id, dum_dims ); \
if( -1 != ( id ) ) \
{ \
size_t ntmp; \
int dvfail = nc_inq_dimlen( ncFile2, dum_dims[0], &ntmp ); \
if( NC_NOERR != dvfail ) \
{ \
ERR_NC( dvfail ) \
} \
( vals ).resize( ntmp ); \
size_t ntmp1 = 0; \
dvfail = nc_get_vara_double( ncFile2, id, &ntmp1, &ntmp, &( vals )[0] ); \
if( NC_NOERR != dvfail ) \
{ \
ERR_NC( dvfail ) \
} \
} \
}
Definition at line 180 of file compareMaps.cpp.
| #define GET_1D_INT_VAR1 | ( | name, | |
| id, | |||
| vals | |||
| ) |
{ \
GET_VAR1( name, id, vals ); \
if( -1 != ( id ) ) \
{ \
size_t ntmp; \
int ivfail = nc_inq_dimlen( ncFile1, ( vals )[0], &ntmp ); \
if( NC_NOERR != ivfail ) \
{ \
ERR_NC( ivfail ) \
} \
( vals ).resize( ntmp ); \
size_t ntmp1 = 0; \
ivfail = nc_get_vara_int( ncFile1, id, &ntmp1, &ntmp, &( vals )[0] ); \
if( NC_NOERR != ivfail ) \
{ \
ERR_NC( ivfail ) \
} \
} \
}
Definition at line 76 of file compareMaps.cpp.
| #define GET_1D_INT_VAR2 | ( | name, | |
| id, | |||
| vals | |||
| ) |
{ \
GET_VAR2( name, id, vals ); \
if( -1 != ( id ) ) \
{ \
size_t ntmp; \
int ivfail = nc_inq_dimlen( ncFile2, ( vals )[0], &ntmp ); \
if( NC_NOERR != ivfail ) \
{ \
ERR_NC( ivfail ) \
} \
( vals ).resize( ntmp ); \
size_t ntmp1 = 0; \
ivfail = nc_get_vara_int( ncFile2, id, &ntmp1, &ntmp, &( vals )[0] ); \
if( NC_NOERR != ivfail ) \
{ \
ERR_NC( ivfail ) \
} \
} \
}
Definition at line 159 of file compareMaps.cpp.
| #define GET_2D_DBL_VAR1 | ( | name, | |
| id, | |||
| vals | |||
| ) |
{ \
std::vector< int > dum_dims; \
GET_VAR1( name, id, dum_dims ); \
if( -1 != ( id ) ) \
{ \
size_t ntmp[2]; \
int dvfail = nc_inq_dimlen( ncFile1, dum_dims[0], &ntmp[0] ); \
if( NC_NOERR != dvfail ) \
{ \
ERR_NC( dvfail ) \
} \
dvfail = nc_inq_dimlen( ncFile1, dum_dims[1], &ntmp[1] ); \
if( NC_NOERR != dvfail ) \
{ \
ERR_NC( dvfail ) \
} \
( vals ).resize( ntmp[0] * ntmp[1] ); \
size_t ntmp1[2] = { 0, 0 }; \
dvfail = nc_get_vara_double( ncFile1, id, ntmp1, ntmp, &( vals )[0] ); \
if( NC_NOERR != dvfail ) \
{ \
ERR_NC( dvfail ) \
} \
} \
}
Definition at line 202 of file compareMaps.cpp.
| #define GET_2D_DBL_VAR2 | ( | name, | |
| id, | |||
| vals | |||
| ) |
{ \
std::vector< int > dum_dims; \
GET_VAR2( name, id, dum_dims ); \
if( -1 != ( id ) ) \
{ \
size_t ntmp[2]; \
int dvfail = nc_inq_dimlen( ncFile2, dum_dims[0], &ntmp[0] ); \
if( NC_NOERR != dvfail ) \
{ \
ERR_NC( dvfail ) \
} \
dvfail = nc_inq_dimlen( ncFile2, dum_dims[1], &ntmp[1] ); \
if( NC_NOERR != dvfail ) \
{ \
ERR_NC( dvfail ) \
} \
( vals ).resize( ntmp[0] * ntmp[1] ); \
size_t ntmp1[2] = { 0, 0 }; \
dvfail = nc_get_vara_double( ncFile2, id, ntmp1, ntmp, &( vals )[0] ); \
if( NC_NOERR != dvfail ) \
{ \
ERR_NC( dvfail ) \
} \
} \
}
Definition at line 229 of file compareMaps.cpp.
| #define GET_DIM1 | ( | ncdim, | |
| name, | |||
| val | |||
| ) |
{ \
int gdfail = nc_inq_dimid( ncFile1, name, &( ncdim ) ); \
if( NC_NOERR == gdfail ) \
{ \
size_t tmp_val; \
gdfail = nc_inq_dimlen( ncFile1, ncdim, &tmp_val ); \
if( NC_NOERR != gdfail ) \
{ \
ERR_NC( gdfail ) \
} \
else \
( val ) = tmp_val; \
} \
else \
( val ) = 0; \
}
Definition at line 38 of file compareMaps.cpp.
| #define GET_DIM2 | ( | ncdim, | |
| name, | |||
| val | |||
| ) |
{ \
int gdfail = nc_inq_dimid( ncFile2, name, &( ncdim ) ); \
if( NC_NOERR == gdfail ) \
{ \
size_t tmp_val; \
gdfail = nc_inq_dimlen( ncFile2, ncdim, &tmp_val ); \
if( NC_NOERR != gdfail ) \
{ \
ERR_NC( gdfail ) \
} \
else \
( val ) = tmp_val; \
} \
else \
( val ) = 0; \
}
Definition at line 121 of file compareMaps.cpp.
| #define GET_VAR1 | ( | name, | |
| id, | |||
| dims | |||
| ) |
{ \
( id ) = -1; \
int gvfail = nc_inq_varid( ncFile1, name, &( id ) ); \
if( NC_NOERR == gvfail ) \
{ \
int ndims; \
gvfail = nc_inq_varndims( ncFile1, id, &ndims ); \
if( NC_NOERR == gvfail ) \
{ \
( dims ).resize( ndims ); \
gvfail = nc_inq_vardimid( ncFile1, id, &( dims )[0] ); \
if( NC_NOERR != gvfail ) \
{ \
ERR_NC( gvfail ) \
} \
} \
} \
}
Definition at line 56 of file compareMaps.cpp.
| #define GET_VAR2 | ( | name, | |
| id, | |||
| dims | |||
| ) |
{ \
( id ) = -1; \
int gvfail = nc_inq_varid( ncFile2, name, &( id ) ); \
if( NC_NOERR == gvfail ) \
{ \
int ndims; \
gvfail = nc_inq_varndims( ncFile2, id, &ndims ); \
if( NC_NOERR == gvfail ) \
{ \
( dims ).resize( ndims ); \
gvfail = nc_inq_vardimid( ncFile2, id, &( dims )[0] ); \
if( NC_NOERR != gvfail ) \
{ \
ERR_NC( gvfail ) \
} \
} \
} \
}
Definition at line 139 of file compareMaps.cpp.
| typedef Eigen::Map< Eigen::VectorXd > EigenV |
Definition at line 256 of file compareMaps.cpp.
| void diff_2d_vect | ( | const char * | var_name, |
| int | n | ||
| ) |
Definition at line 279 of file compareMaps.cpp.
280 {
281 // compare frac_a between maps
282 std::vector< double > fraca1( n ), fraca2( n );
283 int idfa1, idfa2;
284 GET_2D_DBL_VAR1( var_name, idfa1, fraca1 );
285 EigenV fa1( fraca1.data(), n );
286 GET_2D_DBL_VAR2( var_name, idfa2, fraca2 );
287 EigenV fa2( fraca2.data(), n );
288 std::cout << var_name << " diff norm: " << ( fa1 - fa2 ).norm() << "\n";
289 return;
290 }
References GET_2D_DBL_VAR1, and GET_2D_DBL_VAR2.
Referenced by main().
| void diff_vect | ( | const char * | var_name, |
| int | n | ||
| ) |
Definition at line 258 of file compareMaps.cpp.
259 {
260 // compare frac_a between maps
261 std::vector< double > fraca1( n ), fraca2( n );
262 int idfa1, idfa2;
263 GET_1D_DBL_VAR1( var_name, idfa1, fraca1 );
264 EigenV fa1( fraca1.data(), n );
265 GET_1D_DBL_VAR2( var_name, idfa2, fraca2 );
266 EigenV fa2( fraca2.data(), n );
267
268 EigenV diff( fraca2.data(), n );
269 diff = fa1 - fa2;
270
271 int imin, imax;
272 double minV = diff.minCoeff( &imin );
273 double maxV = diff.maxCoeff( &imax );
274 std::cout << var_name << " diff norm: " << diff.norm() << " min at " << imin << " : " << minV << " max at " << imax
275 << " : " << maxV << "\n";
276 return;
277 }
References GET_1D_DBL_VAR1, and GET_1D_DBL_VAR2.
Referenced by main().
| int main | ( | int | argc, |
| char * | argv[] | ||
| ) |
Definition at line 291 of file compareMaps.cpp.
292 {
293
294 ProgOptions opts;
295
296 std::string inputfile1, inputfile2;
297 opts.addOpt< std::string >( "firstMap,i", "input filename 1", &inputfile1 );
298 opts.addOpt< std::string >( "secondMap,j", "input second map", &inputfile2 );
299 int print_diff = 20;
300 opts.addOpt< int >( "print_differences,p", "print differences ", &print_diff );
301 double fraction = 0.9;
302 opts.addOpt< double >( "fraction_diff,f", "fraction threshold", &fraction );
303
304 opts.parseCommandLine( argc, argv );
305
306 // Open netcdf/exodus file
307 int fail = nc_open( inputfile1.c_str(), 0, &ncFile1 );
308 if( NC_NOWRITE != fail )
309 {
310 ERR_NC( fail )
311 }
312
313 std::cout << " opened " << inputfile1 << " for map 1 \n";
314
315 fail = nc_open( inputfile2.c_str(), 0, &ncFile2 );
316 if( NC_NOWRITE != fail )
317 {
318 ERR_NC( fail )
319 }
320
321 std::cout << " opened " << inputfile2 << " for map 2 \n";
322 int temp_dim;
323 int na1, nb1, ns1, na2, nb2, ns2, nv_a, nv_b, nv_a2, nv_b2;
324 GET_DIM1( temp_dim, "n_a", na1 );
325 GET_DIM2( temp_dim, "n_a", na2 );
326 GET_DIM1( temp_dim, "n_b", nb1 );
327 GET_DIM2( temp_dim, "n_b", nb2 );
328 GET_DIM1( temp_dim, "n_s", ns1 );
329 GET_DIM2( temp_dim, "n_s", ns2 );
330 GET_DIM1( temp_dim, "nv_a", nv_a );
331 GET_DIM2( temp_dim, "nv_a", nv_a2 );
332 GET_DIM1( temp_dim, "nv_b", nv_b );
333 GET_DIM2( temp_dim, "nv_b", nv_b2 );
334 if( nv_a != nv_a2 || nv_b != nv_b2 )
335 {
336 std::cout << " different nv dimensions:" << nv_a << " == " << nv_a2 << " or " << nv_b << " == " << nv_b2
337 << " bail out \n";
338 return 1;
339 }
340 std::cout << " n_a, n_b, n_s : " << na1 << ", " << nb1 << ", " << ns1 << " for map 1 \n";
341
342 std::cout << " n_a, n_b, n_s : " << na2 << ", " << nb2 << ", " << ns2 << " for map 2 \n";
343 if( na1 != na2 || nb1 != nb2 )
344 {
345 std::cout << " different dimensions bail out \n";
346 return 1;
347 }
348 std::vector< int > col1( ns1 ), row1( ns1 );
349 std::vector< int > col2( ns2 ), row2( ns2 );
350
351 std::vector< double > val1( ns1 ), val2( ns2 );
352
353 int idrow1, idcol1, idrow2, idcol2, ids1, ids2;
354 GET_1D_INT_VAR1( "row", idrow1, row1 );
355 GET_1D_INT_VAR1( "col", idcol1, col1 );
356 GET_1D_DBL_VAR1( "S", ids1, val1 );
357 GET_1D_INT_VAR2( "row", idrow2, row2 );
358 GET_1D_INT_VAR2( "col", idcol2, col2 );
359 GET_1D_DBL_VAR2( "S", ids2, val2 );
360
361 // first matrix
362 typedef Eigen::Triplet< double > Triplet;
363 std::vector< Triplet > tripletList;
364 tripletList.reserve( ns1 );
365 for( int iv = 0; iv < ns1; iv++ )
366 {
367 tripletList.push_back( Triplet( row1[iv] - 1, col1[iv] - 1, val1[iv] ) );
368 }
369 Eigen::SparseMatrix< double > weight1( nb1, na1 );
370
371 weight1.setFromTriplets( tripletList.begin(), tripletList.end() );
372 weight1.makeCompressed();
373
374 if( ns1 != ns2 ) tripletList.resize( ns2 );
375 for( int iv = 0; iv < ns2; iv++ )
376 {
377 tripletList[iv] = Triplet( row2[iv] - 1, col2[iv] - 1, val2[iv] );
378 }
379 Eigen::SparseMatrix< double > weight2( nb2, na2 );
380
381 weight2.setFromTriplets( tripletList.begin(), tripletList.end() );
382 weight2.makeCompressed();
383
384 // default storage type is column major
385 Eigen::SparseMatrix< double > diff = weight1 - weight2;
386 diff.makeCompressed(); // is it needed or not ?
387 auto coeffs = diff.coeffs();
388 double maxv = coeffs.maxCoeff();
389 double minv = coeffs.minCoeff();
390 std::cout << " euclidian norm for difference: " << diff.norm()
391 << " \n squared norm for difference: " << diff.squaredNorm() << "\n"
392 << " minv: " << minv << " maxv: " << maxv << "\n";
393 // print out the first 10 positions for which the value is outside 90% of min/max values
394 double min_threshold = fraction * minv;
395 double max_threshold = fraction * maxv;
396 int counter = 0;
397 std::cout << std::setprecision( 9 );
398 for( int k = 0; ( k < diff.outerSize() ); ++k ) // this is by column
399 {
400 for( Eigen::SparseMatrix< double >::InnerIterator it( diff, k ); ( it ) && ( counter < print_diff ); ++it )
401 {
402 double val = it.value();
403 if( val <= min_threshold || val >= max_threshold )
404 {
405 int row = it.row();
406 int col = it.col();
407 std::cout << " counter:" << counter << "\t col: " << col + 1 << "\t row: " << row + 1
408 << "\t diff: " << val;
409 std::cout << "\t map1: " << weight1.coeffRef( row, col ) << "\t map2: " << weight2.coeffRef( row, col )
410 << "\n"; // row index
411 counter++;
412 }
413 }
414 }
415
416 // compare frac_a between maps
417 diff_vect( "frac_a", na1 );
418 diff_vect( "frac_b", nb1 );
419 diff_vect( "area_a", na1 );
420 diff_vect( "area_b", nb1 );
421 diff_vect( "yc_a", na1 );
422 diff_vect( "yc_b", nb1 );
423 diff_vect( "xc_a", na1 );
424 diff_vect( "xc_b", nb1 );
425 diff_2d_vect( "xv_a", na1 * nv_a );
426 diff_2d_vect( "yv_a", na1 * nv_a );
427 diff_2d_vect( "xv_b", nb1 * nv_b );
428 diff_2d_vect( "yv_b", nb1 * nv_b );
429
430 return 0;
431 }
References ProgOptions::addOpt(), diff_2d_vect(), diff_vect(), ERR_NC, fail(), GET_1D_DBL_VAR1, GET_1D_DBL_VAR2, GET_1D_INT_VAR1, GET_1D_INT_VAR2, GET_DIM1, GET_DIM2, ncFile1, ncFile2, and ProgOptions::parseCommandLine().
| int ncFile1 |
Definition at line 36 of file compareMaps.cpp.
Referenced by main().
| int ncFile2 |
Definition at line 119 of file compareMaps.cpp.
Referenced by main().
1.9.1.