17 #ifndef MOAB_HAVE_EIGEN3
18 #error compareMaps tool requires eigen3 configuration
27 #include <Eigen/Sparse>
30 printf( "Error: %s\n", nc_strerror( e ) ); \
39 #define GET_DIM1( ncdim, name, val ) \
41 int gdfail = nc_inq_dimid( ncFile1, name, &( ncdim ) ); \
42 if( NC_NOERR == gdfail ) \
45 gdfail = nc_inq_dimlen( ncFile1, ncdim, &tmp_val ); \
46 if( NC_NOERR != gdfail ) \
57 #define GET_VAR1( name, id, dims ) \
60 int gvfail = nc_inq_varid( ncFile1, name, &( id ) ); \
61 if( NC_NOERR == gvfail ) \
64 gvfail = nc_inq_varndims( ncFile1, id, &ndims ); \
65 if( NC_NOERR == gvfail ) \
67 ( dims ).resize( ndims ); \
68 gvfail = nc_inq_vardimid( ncFile1, id, &( dims )[0] ); \
69 if( NC_NOERR != gvfail ) \
77 #define GET_1D_INT_VAR1( name, id, vals ) \
79 GET_VAR1( name, id, vals ); \
83 int ivfail = nc_inq_dimlen( ncFile1, ( vals )[0], &ntmp ); \
84 if( NC_NOERR != ivfail ) \
88 ( vals ).resize( ntmp ); \
90 ivfail = nc_get_vara_int( ncFile1, id, &ntmp1, &ntmp, &( vals )[0] ); \
91 if( NC_NOERR != ivfail ) \
98 #define GET_1D_DBL_VAR1( name, id, vals ) \
100 std::vector< int > dum_dims; \
101 GET_VAR1( name, id, dum_dims ); \
105 int dvfail = nc_inq_dimlen( ncFile1, dum_dims[0], &ntmp ); \
106 if( NC_NOERR != dvfail ) \
110 ( vals ).resize( ntmp ); \
112 dvfail = nc_get_vara_double( ncFile1, id, &ntmp1, &ntmp, &( vals )[0] ); \
113 if( NC_NOERR != dvfail ) \
122 #define GET_DIM2( ncdim, name, val ) \
124 int gdfail = nc_inq_dimid( ncFile2, name, &( ncdim ) ); \
125 if( NC_NOERR == gdfail ) \
128 gdfail = nc_inq_dimlen( ncFile2, ncdim, &tmp_val ); \
129 if( NC_NOERR != gdfail ) \
140 #define GET_VAR2( name, id, dims ) \
143 int gvfail = nc_inq_varid( ncFile2, name, &( id ) ); \
144 if( NC_NOERR == gvfail ) \
147 gvfail = nc_inq_varndims( ncFile2, id, &ndims ); \
148 if( NC_NOERR == gvfail ) \
150 ( dims ).resize( ndims ); \
151 gvfail = nc_inq_vardimid( ncFile2, id, &( dims )[0] ); \
152 if( NC_NOERR != gvfail ) \
160 #define GET_1D_INT_VAR2( name, id, vals ) \
162 GET_VAR2( name, id, vals ); \
166 int ivfail = nc_inq_dimlen( ncFile2, ( vals )[0], &ntmp ); \
167 if( NC_NOERR != ivfail ) \
171 ( vals ).resize( ntmp ); \
173 ivfail = nc_get_vara_int( ncFile2, id, &ntmp1, &ntmp, &( vals )[0] ); \
174 if( NC_NOERR != ivfail ) \
181 #define GET_1D_DBL_VAR2( name, id, vals ) \
183 std::vector< int > dum_dims; \
184 GET_VAR2( name, id, dum_dims ); \
188 int dvfail = nc_inq_dimlen( ncFile2, dum_dims[0], &ntmp ); \
189 if( NC_NOERR != dvfail ) \
193 ( vals ).resize( ntmp ); \
195 dvfail = nc_get_vara_double( ncFile2, id, &ntmp1, &ntmp, &( vals )[0] ); \
196 if( NC_NOERR != dvfail ) \
203 #define GET_2D_DBL_VAR1( name, id, vals ) \
205 std::vector< int > dum_dims; \
206 GET_VAR1( name, id, dum_dims ); \
210 int dvfail = nc_inq_dimlen( ncFile1, dum_dims[0], &ntmp[0] ); \
211 if( NC_NOERR != dvfail ) \
215 dvfail = nc_inq_dimlen( ncFile1, dum_dims[1], &ntmp[1] ); \
216 if( NC_NOERR != dvfail ) \
220 ( vals ).resize( ntmp[0] * ntmp[1] ); \
221 size_t ntmp1[2] = { 0, 0 }; \
222 dvfail = nc_get_vara_double( ncFile1, id, ntmp1, ntmp, &( vals )[0] ); \
223 if( NC_NOERR != dvfail ) \
230 #define GET_2D_DBL_VAR2( name, id, vals ) \
232 std::vector< int > dum_dims; \
233 GET_VAR2( name, id, dum_dims ); \
237 int dvfail = nc_inq_dimlen( ncFile2, dum_dims[0], &ntmp[0] ); \
238 if( NC_NOERR != dvfail ) \
242 dvfail = nc_inq_dimlen( ncFile2, dum_dims[1], &ntmp[1] ); \
243 if( NC_NOERR != dvfail ) \
247 ( vals ).resize( ntmp[0] * ntmp[1] ); \
248 size_t ntmp1[2] = { 0, 0 }; \
249 dvfail = nc_get_vara_double( ncFile2, id, ntmp1, ntmp, &( vals )[0] ); \
250 if( NC_NOERR != dvfail ) \
257 typedef Eigen::Map< Eigen::VectorXd >
EigenV;
262 std::vector< double > fraca1( n ), fraca2( n );
265 EigenV fa1( fraca1.data(), n );
267 EigenV fa2( fraca2.data(), n );
269 EigenV diff( fraca2.data(), n );
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";
283 std::vector< double > fraca1( n ), fraca2( n );
286 EigenV fa1( fraca1.data(), n );
288 EigenV fa2( fraca2.data(), n );
289 std::cout << var_name <<
" diff norm: " << ( fa1 - fa2 ).norm() <<
"\n";
297 return fabs( a.value() ) > fabs( b.value() ) ;
301 int main(
int argc,
char* argv[] )
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 );
310 opts.
addOpt<
int >(
"print_differences,p",
"print differences ", &print_diff );
315 int fail = nc_open( inputfile1.c_str(), 0, &
ncFile1 );
316 if( NC_NOWRITE !=
fail )
321 std::cout <<
" opened " << inputfile1 <<
" for map 1 \n";
324 if( NC_NOWRITE !=
fail )
329 std::cout <<
" opened " << inputfile2 <<
" for map 2 \n";
331 int na1, nb1, ns1, na2, nb2, ns2, nv_a, nv_b, nv_a2, nv_b2;
339 GET_DIM2( temp_dim,
"nv_a", nv_a2 );
341 GET_DIM2( temp_dim,
"nv_b", nv_b2 );
342 if( nv_a != nv_a2 || nv_b != nv_b2 )
344 std::cout <<
" different nv dimensions:" << nv_a <<
" == " << nv_a2 <<
" or " << nv_b <<
" == " << nv_b2
348 std::cout <<
" n_a, n_b, n_s : " << na1 <<
", " << nb1 <<
", " << ns1 <<
" for map 1 \n";
350 std::cout <<
" n_a, n_b, n_s : " << na2 <<
", " << nb2 <<
", " << ns2 <<
" for map 2 \n";
351 if( na1 != na2 || nb1 != nb2 )
353 std::cout <<
" different dimensions bail out \n";
356 std::vector< int > col1( ns1 ), row1( ns1 );
357 std::vector< int > col2( ns2 ), row2( ns2 );
359 std::vector< double > val1( ns1 ), val2( ns2 );
361 int idrow1, idcol1, idrow2, idcol2, ids1, ids2;
371 std::vector< Triplet > tripletList;
372 tripletList.reserve( ns1 );
373 for(
int iv = 0; iv < ns1; iv++ )
375 tripletList.push_back(
Triplet( row1[iv] - 1, col1[iv] - 1, val1[iv] ) );
377 Eigen::SparseMatrix< double > weight1( nb1, na1 );
379 weight1.setFromTriplets( tripletList.begin(), tripletList.end() );
380 weight1.makeCompressed();
382 if( ns1 != ns2 ) tripletList.resize( ns2 );
383 for(
int iv = 0; iv < ns2; iv++ )
385 tripletList[iv] =
Triplet( row2[iv] - 1, col2[iv] - 1, val2[iv] );
387 Eigen::SparseMatrix< double > weight2( nb2, na2 );
389 weight2.setFromTriplets( tripletList.begin(), tripletList.end() );
390 weight2.makeCompressed();
393 Eigen::SparseMatrix< double > diff = weight1 - weight2;
394 diff.makeCompressed();
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";
402 std::priority_queue< Triplet, std::vector<Triplet>,
CompareTriplets > largestDiffs;
404 for(
int k = 0; ( k < diff.outerSize() ); ++k )
406 for( Eigen::SparseMatrix< double >::InnerIterator it( diff, k ); ( it ) ; ++it )
408 double val = it.value();
409 Triplet tp(it.row(), it.col(), val);
410 largestDiffs.push(tp);
411 if (largestDiffs.size() > print_diff)
415 std::cout << std::setprecision( 16 );
417 while (!largestDiffs.empty()) {
418 Triplet tp = largestDiffs.top();
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 )