Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
Classes | Macros | Typedefs | Functions | Variables
compareMaps.cpp File Reference
#include "moab/MOABConfig.h"
#include "moab/ProgOptions.hpp"
#include "netcdf.h"
#include <cmath>
#include <iomanip>
#include <queue>
#include <Eigen/Sparse>
+ Include dependency graph for compareMaps.cpp:

Go to the source code of this file.

Classes

struct  CompareTriplets
 

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
 
typedef Eigen::Triplet< double > Triplet
 

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
 

Macro Definition Documentation

◆ ERR_NC

#define ERR_NC (   e)
Value:
{ \
printf( "Error: %s\n", nc_strerror( e ) ); \
exit( 2 ); \
}

Definition at line 28 of file compareMaps.cpp.

◆ GET_1D_DBL_VAR1

#define GET_1D_DBL_VAR1 (   name,
  id,
  vals 
)
Value:
{ \
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 98 of file compareMaps.cpp.

◆ GET_1D_DBL_VAR2

#define GET_1D_DBL_VAR2 (   name,
  id,
  vals 
)
Value:
{ \
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 181 of file compareMaps.cpp.

◆ GET_1D_INT_VAR1

#define GET_1D_INT_VAR1 (   name,
  id,
  vals 
)
Value:
{ \
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 77 of file compareMaps.cpp.

◆ GET_1D_INT_VAR2

#define GET_1D_INT_VAR2 (   name,
  id,
  vals 
)
Value:
{ \
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 160 of file compareMaps.cpp.

◆ GET_2D_DBL_VAR1

#define GET_2D_DBL_VAR1 (   name,
  id,
  vals 
)
Value:
{ \
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 203 of file compareMaps.cpp.

◆ GET_2D_DBL_VAR2

#define GET_2D_DBL_VAR2 (   name,
  id,
  vals 
)
Value:
{ \
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 230 of file compareMaps.cpp.

◆ GET_DIM1

#define GET_DIM1 (   ncdim,
  name,
  val 
)
Value:
{ \
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 39 of file compareMaps.cpp.

◆ GET_DIM2

#define GET_DIM2 (   ncdim,
  name,
  val 
)
Value:
{ \
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 122 of file compareMaps.cpp.

◆ GET_VAR1

#define GET_VAR1 (   name,
  id,
  dims 
)
Value:
{ \
( 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 57 of file compareMaps.cpp.

◆ GET_VAR2

#define GET_VAR2 (   name,
  id,
  dims 
)
Value:
{ \
( 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 140 of file compareMaps.cpp.

Typedef Documentation

◆ EigenV

typedef Eigen::Map< Eigen::VectorXd > EigenV

Definition at line 257 of file compareMaps.cpp.

◆ Triplet

typedef Eigen::Triplet< double > Triplet

Definition at line 293 of file compareMaps.cpp.

Function Documentation

◆ diff_2d_vect()

void diff_2d_vect ( const char *  var_name,
int  n 
)

Definition at line 280 of file compareMaps.cpp.

281 {
282  // compare frac_a between maps
283  std::vector< double > fraca1( n ), fraca2( n );
284  int idfa1, idfa2;
285  GET_2D_DBL_VAR1( var_name, idfa1, fraca1 );
286  EigenV fa1( fraca1.data(), n );
287  GET_2D_DBL_VAR2( var_name, idfa2, fraca2 );
288  EigenV fa2( fraca2.data(), n );
289  std::cout << var_name << " diff norm: " << ( fa1 - fa2 ).norm() << "\n";
290  return;
291 }

References GET_2D_DBL_VAR1, and GET_2D_DBL_VAR2.

Referenced by main().

◆ diff_vect()

void diff_vect ( const char *  var_name,
int  n 
)

Definition at line 259 of file compareMaps.cpp.

260 {
261  // compare frac_a between maps
262  std::vector< double > fraca1( n ), fraca2( n );
263  int idfa1, idfa2;
264  GET_1D_DBL_VAR1( var_name, idfa1, fraca1 );
265  EigenV fa1( fraca1.data(), n );
266  GET_1D_DBL_VAR2( var_name, idfa2, fraca2 );
267  EigenV fa2( fraca2.data(), n );
268 
269  EigenV diff( fraca2.data(), n );
270  diff = fa1 - fa2;
271 
272  int imin, imax;
273  double minV = diff.minCoeff( &imin );
274  double maxV = diff.maxCoeff( &imax );
275  std::cout << var_name << " diff norm: " << diff.norm() << " min at " << imin << " : " << minV << " max at " << imax
276  << " : " << maxV << "\n";
277  return;
278 }

References GET_1D_DBL_VAR1, and GET_1D_DBL_VAR2.

Referenced by main().

◆ main()

int main ( int  argc,
char *  argv[] 
)

Definition at line 301 of file compareMaps.cpp.

302 {
303 
304  ProgOptions opts;
305 
306  std::string inputfile1, inputfile2;
307  opts.addOpt< std::string >( "firstMap,i", "input filename 1", &inputfile1 );
308  opts.addOpt< std::string >( "secondMap,j", "input second map", &inputfile2 );
309  int print_diff = 10;
310  opts.addOpt< int >( "print_differences,p", "print differences ", &print_diff );
311 
312  opts.parseCommandLine( argc, argv );
313 
314  // Open netcdf/exodus file
315  int fail = nc_open( inputfile1.c_str(), 0, &ncFile1 );
316  if( NC_NOWRITE != fail )
317  {
318  ERR_NC( fail )
319  }
320 
321  std::cout << " opened " << inputfile1 << " for map 1 \n";
322 
323  fail = nc_open( inputfile2.c_str(), 0, &ncFile2 );
324  if( NC_NOWRITE != fail )
325  {
326  ERR_NC( fail )
327  }
328 
329  std::cout << " opened " << inputfile2 << " for map 2 \n";
330  int temp_dim;
331  int na1, nb1, ns1, na2, nb2, ns2, nv_a, nv_b, nv_a2, nv_b2;
332  GET_DIM1( temp_dim, "n_a", na1 );
333  GET_DIM2( temp_dim, "n_a", na2 );
334  GET_DIM1( temp_dim, "n_b", nb1 );
335  GET_DIM2( temp_dim, "n_b", nb2 );
336  GET_DIM1( temp_dim, "n_s", ns1 );
337  GET_DIM2( temp_dim, "n_s", ns2 );
338  GET_DIM1( temp_dim, "nv_a", nv_a );
339  GET_DIM2( temp_dim, "nv_a", nv_a2 );
340  GET_DIM1( temp_dim, "nv_b", nv_b );
341  GET_DIM2( temp_dim, "nv_b", nv_b2 );
342  if( nv_a != nv_a2 || nv_b != nv_b2 )
343  {
344  std::cout << " different nv dimensions:" << nv_a << " == " << nv_a2 << " or " << nv_b << " == " << nv_b2
345  << " bail out \n";
346  return 1;
347  }
348  std::cout << " n_a, n_b, n_s : " << na1 << ", " << nb1 << ", " << ns1 << " for map 1 \n";
349 
350  std::cout << " n_a, n_b, n_s : " << na2 << ", " << nb2 << ", " << ns2 << " for map 2 \n";
351  if( na1 != na2 || nb1 != nb2 )
352  {
353  std::cout << " different dimensions bail out \n";
354  return 1;
355  }
356  std::vector< int > col1( ns1 ), row1( ns1 );
357  std::vector< int > col2( ns2 ), row2( ns2 );
358 
359  std::vector< double > val1( ns1 ), val2( ns2 );
360 
361  int idrow1, idcol1, idrow2, idcol2, ids1, ids2;
362  GET_1D_INT_VAR1( "row", idrow1, row1 );
363  GET_1D_INT_VAR1( "col", idcol1, col1 );
364  GET_1D_DBL_VAR1( "S", ids1, val1 );
365  GET_1D_INT_VAR2( "row", idrow2, row2 );
366  GET_1D_INT_VAR2( "col", idcol2, col2 );
367  GET_1D_DBL_VAR2( "S", ids2, val2 );
368 
369  // first matrix
370 
371  std::vector< Triplet > tripletList;
372  tripletList.reserve( ns1 );
373  for( int iv = 0; iv < ns1; iv++ )
374  {
375  tripletList.push_back( Triplet( row1[iv] - 1, col1[iv] - 1, val1[iv] ) );
376  }
377  Eigen::SparseMatrix< double > weight1( nb1, na1 );
378 
379  weight1.setFromTriplets( tripletList.begin(), tripletList.end() );
380  weight1.makeCompressed();
381 
382  if( ns1 != ns2 ) tripletList.resize( ns2 );
383  for( int iv = 0; iv < ns2; iv++ )
384  {
385  tripletList[iv] = Triplet( row2[iv] - 1, col2[iv] - 1, val2[iv] );
386  }
387  Eigen::SparseMatrix< double > weight2( nb2, na2 );
388 
389  weight2.setFromTriplets( tripletList.begin(), tripletList.end() );
390  weight2.makeCompressed();
391 
392  // default storage type is column major
393  Eigen::SparseMatrix< double > diff = weight1 - weight2;
394  diff.makeCompressed(); // is it needed or not ?
395  auto coeffs = diff.coeffs();
396  double maxv = coeffs.maxCoeff();
397  double minv = coeffs.minCoeff();
398  std::cout << " euclidian norm for difference: " << diff.norm()
399  << " \n squared norm for difference: " << diff.squaredNorm() << "\n"
400  << " minv: " << minv << " maxv: " << maxv << "\n";
401  // print out the largest 20 absolute values and position in diff sparse matrix
402  std::priority_queue< Triplet, std::vector<Triplet>, CompareTriplets > largestDiffs;
403 
404  for( int k = 0; ( k < diff.outerSize() ); ++k ) // this is by column
405  {
406  for( Eigen::SparseMatrix< double >::InnerIterator it( diff, k ); ( it ) ; ++it )
407  {
408  double val = it.value();
409  Triplet tp(it.row(), it.col(), val);
410  largestDiffs.push(tp);
411  if (largestDiffs.size() > print_diff)
412  largestDiffs.pop();
413  }
414  }
415  std::cout << std::setprecision( 16 );
416  int counter=0;
417  while (!largestDiffs.empty()) {
418  Triplet tp = largestDiffs.top();
419  largestDiffs.pop();
420 
421  int row = tp.row();
422  int col = tp.col();
423  std::cout << " counter:" << counter << "\t col: " << col + 1 << "\t row: " << row + 1
424  << "\t diff: " << tp.value();
425  std::cout << "\t map1: " << weight1.coeffRef( row, col ) << "\t map2: " << weight2.coeffRef( row, col )
426  << "\n"; // row index
427  counter ++;
428  }
429 
430  // compare frac_a between maps
431  diff_vect( "frac_a", na1 );
432  diff_vect( "frac_b", nb1 );
433  diff_vect( "area_a", na1 );
434  diff_vect( "area_b", nb1 );
435  diff_vect( "yc_a", na1 );
436  diff_vect( "yc_b", nb1 );
437  diff_vect( "xc_a", na1 );
438  diff_vect( "xc_b", nb1 );
439  diff_2d_vect( "xv_a", na1 * nv_a );
440  diff_2d_vect( "yv_a", na1 * nv_a );
441  diff_2d_vect( "xv_b", nb1 * nv_b );
442  diff_2d_vect( "yv_b", nb1 * nv_b );
443 
444  return 0;
445 }

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().

Variable Documentation

◆ ncFile1

int ncFile1

Definition at line 37 of file compareMaps.cpp.

Referenced by main().

◆ ncFile2

int ncFile2

Definition at line 120 of file compareMaps.cpp.

Referenced by main().