![]() |
Mesh Oriented datABase
(version 5.4.1)
Array-based unstructured mesh datastructure
|
#include "moab/MOABConfig.h"
#include "moab/ProgOptions.hpp"
#include "netcdf.h"
#include <cmath>
#include <iomanip>
#include <Eigen/Sparse>
Go to the source code of this file.
Defines | |
#define | ERR_NC(e) |
#define | INS_ID(stringvar, prefix, id) sprintf( stringvar, prefix, id ) |
#define | GET_DIM1(ncdim, name, val) |
#define | GET_DIMB1(ncdim, name, varname, id, 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_DIMB2(ncdim, name, varname, id, 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.
Referenced by main().
#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 103 of file compareMaps.cpp.
Referenced by diff_vect(), and main().
#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 190 of file compareMaps.cpp.
Referenced by diff_vect(), and main().
#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 82 of file compareMaps.cpp.
Referenced by main().
#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 169 of file compareMaps.cpp.
Referenced by main().
#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 212 of file compareMaps.cpp.
Referenced by diff_2d_vect().
#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 239 of file compareMaps.cpp.
Referenced by diff_2d_vect().
#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 40 of file compareMaps.cpp.
Referenced by main().
#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 127 of file compareMaps.cpp.
Referenced by main().
#define GET_DIMB1 | ( | ncdim, | |
name, | |||
varname, | |||
id, | |||
val | |||
) |
INS_ID( name, varname, id ); \
GET_DIM1( ncdim, name, val );
Definition at line 58 of file compareMaps.cpp.
#define GET_DIMB2 | ( | ncdim, | |
name, | |||
varname, | |||
id, | |||
val | |||
) |
INS_ID( name, varname, id ); \
GET_DIM2( ncdim, name, val );
Definition at line 145 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 62 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 149 of file compareMaps.cpp.
#define INS_ID | ( | stringvar, | |
prefix, | |||
id | |||
) | sprintf( stringvar, prefix, id ) |
Definition at line 38 of file compareMaps.cpp.
typedef Eigen::Map< Eigen::VectorXd > EigenV |
Definition at line 266 of file compareMaps.cpp.
void diff_2d_vect | ( | const char * | var_name, |
int | n | ||
) |
Definition at line 289 of file compareMaps.cpp.
References GET_2D_DBL_VAR1, and GET_2D_DBL_VAR2.
Referenced by main().
{
// compare frac_a between maps
std::vector< double > fraca1( n ), fraca2( n );
int idfa1, idfa2;
GET_2D_DBL_VAR1( var_name, idfa1, fraca1 );
EigenV fa1( fraca1.data(), n );
GET_2D_DBL_VAR2( var_name, idfa2, fraca2 );
EigenV fa2( fraca2.data(), n );
std::cout << var_name << " diff norm: " << ( fa1 - fa2 ).norm() << "\n";
return;
}
void diff_vect | ( | const char * | var_name, |
int | n | ||
) |
Definition at line 268 of file compareMaps.cpp.
References GET_1D_DBL_VAR1, and GET_1D_DBL_VAR2.
Referenced by main().
{
// compare frac_a between maps
std::vector< double > fraca1( n ), fraca2( n );
int idfa1, idfa2;
GET_1D_DBL_VAR1( var_name, idfa1, fraca1 );
EigenV fa1( fraca1.data(), n );
GET_1D_DBL_VAR2( var_name, idfa2, fraca2 );
EigenV fa2( fraca2.data(), n );
EigenV diff( fraca2.data(), n );
diff = fa1 - fa2;
int imin, imax;
double minV = diff.minCoeff( &imin );
double maxV = diff.maxCoeff( &imax );
std::cout << var_name << " diff norm: " << diff.norm() << " min at " << imin << " : " << minV << " max at " << imax
<< " : " << maxV << "\n";
return;
}
int main | ( | int | argc, |
char * | argv[] | ||
) |
Definition at line 301 of file compareMaps.cpp.
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().
{
ProgOptions opts;
std::string inputfile1, inputfile2;
opts.addOpt< std::string >( "firstMap,i", "input filename 1", &inputfile1 );
opts.addOpt< std::string >( "secondMap,j", "input second map", &inputfile2 );
int print_diff = 20;
opts.addOpt< int >( "print_differences,p", "print differences ", &print_diff );
double fraction = 0.9;
opts.addOpt< double >( "fraction_diff,f", "fraction threshold", &fraction );
opts.parseCommandLine( argc, argv );
// Open netcdf/exodus file
int fail = nc_open( inputfile1.c_str(), 0, &ncFile1 );
if( NC_NOWRITE != fail )
{
ERR_NC( fail )
}
std::cout << " opened " << inputfile1 << " for map 1 \n";
fail = nc_open( inputfile2.c_str(), 0, &ncFile2 );
if( NC_NOWRITE != fail )
{
ERR_NC( fail )
}
std::cout << " opened " << inputfile2 << " for map 2 \n";
int temp_dim;
int na1, nb1, ns1, na2, nb2, ns2, nv_a, nv_b, nv_a2, nv_b2;
GET_DIM1( temp_dim, "n_a", na1 );
GET_DIM2( temp_dim, "n_a", na2 );
GET_DIM1( temp_dim, "n_b", nb1 );
GET_DIM2( temp_dim, "n_b", nb2 );
GET_DIM1( temp_dim, "n_s", ns1 );
GET_DIM2( temp_dim, "n_s", ns2 );
GET_DIM1( temp_dim, "nv_a", nv_a );
GET_DIM2( temp_dim, "nv_a", nv_a2 );
GET_DIM1( temp_dim, "nv_b", nv_b );
GET_DIM2( temp_dim, "nv_b", nv_b2 );
if( nv_a != nv_a2 || nv_b != nv_b2 )
{
std::cout << " different nv dimensions:" << nv_a << " == " << nv_a2 << " or " << nv_b << " == " << nv_b2
<< " bail out \n";
return 1;
}
std::cout << " n_a, n_b, n_s : " << na1 << ", " << nb1 << ", " << ns1 << " for map 1 \n";
std::cout << " n_a, n_b, n_s : " << na2 << ", " << nb2 << ", " << ns2 << " for map 2 \n";
if( na1 != na2 || nb1 != nb2 )
{
std::cout << " different dimensions bail out \n";
return 1;
}
std::vector< int > col1( ns1 ), row1( ns1 );
std::vector< int > col2( ns2 ), row2( ns2 );
std::vector< double > val1( ns1 ), val2( ns2 );
int idrow1, idcol1, idrow2, idcol2, ids1, ids2;
GET_1D_INT_VAR1( "row", idrow1, row1 );
GET_1D_INT_VAR1( "col", idcol1, col1 );
GET_1D_DBL_VAR1( "S", ids1, val1 );
GET_1D_INT_VAR2( "row", idrow2, row2 );
GET_1D_INT_VAR2( "col", idcol2, col2 );
GET_1D_DBL_VAR2( "S", ids2, val2 );
// first matrix
typedef Eigen::Triplet< double > Triplet;
std::vector< Triplet > tripletList;
tripletList.reserve( ns1 );
for( int iv = 0; iv < ns1; iv++ )
{
tripletList.push_back( Triplet( row1[iv] - 1, col1[iv] - 1, val1[iv] ) );
}
Eigen::SparseMatrix< double > weight1( nb1, na1 );
weight1.setFromTriplets( tripletList.begin(), tripletList.end() );
weight1.makeCompressed();
if( ns1 != ns2 ) tripletList.resize( ns2 );
for( int iv = 0; iv < ns2; iv++ )
{
tripletList[iv] = Triplet( row2[iv] - 1, col2[iv] - 1, val2[iv] );
}
Eigen::SparseMatrix< double > weight2( nb2, na2 );
weight2.setFromTriplets( tripletList.begin(), tripletList.end() );
weight2.makeCompressed();
// default storage type is column major
Eigen::SparseMatrix< double > diff = weight1 - weight2;
diff.makeCompressed(); // is it needed or not ?
auto coeffs = diff.coeffs();
double maxv = coeffs.maxCoeff();
double minv = coeffs.minCoeff();
std::cout << " euclidian norm for difference: " << diff.norm()
<< " \n squared norm for difference: " << diff.squaredNorm() << "\n"
<< " minv: " << minv << " maxv: " << maxv << "\n";
// print out the first 10 positions for which the value is outside 90% of min/max values
double min_threshold = fraction * minv;
double max_threshold = fraction * maxv;
int counter = 0;
std::cout << std::setprecision( 9 );
for( int k = 0; ( k < diff.outerSize() ); ++k ) // this is by column
{
for( Eigen::SparseMatrix< double >::InnerIterator it( diff, k ); ( it ) && ( counter < print_diff ); ++it )
{
double val = it.value();
if( val <= min_threshold || val >= max_threshold )
{
int row = it.row();
int col = it.col();
std::cout << " counter:" << counter << "\t col: " << col + 1 << "\t row: " << row + 1
<< "\t diff: " << val;
std::cout << "\t map1: " << weight1.coeffRef( row, col ) << "\t map2: " << weight2.coeffRef( row, col )
<< "\n"; // row index
counter++;
}
}
}
// compare frac_a between maps
diff_vect( "frac_a", na1 );
diff_vect( "frac_b", nb1 );
diff_vect( "area_a", na1 );
diff_vect( "area_b", nb1 );
diff_vect( "yc_a", na1 );
diff_vect( "yc_b", nb1 );
diff_vect( "xc_a", na1 );
diff_vect( "xc_b", nb1 );
diff_2d_vect( "xv_a", na1 * nv_a );
diff_2d_vect( "yv_a", na1 * nv_a );
diff_2d_vect( "xv_b", nb1 * nv_b );
diff_2d_vect( "yv_b", nb1 * nv_b );
return 0;
}
int ncFile1 |
Definition at line 36 of file compareMaps.cpp.
Referenced by main().
int ncFile2 |
Definition at line 125 of file compareMaps.cpp.
Referenced by main().