Mesh Oriented datABase  (version 5.6.0)
An array-based unstructured mesh library
compareMaps.cpp
Go to the documentation of this file.
1 /*
2  * compareMaps.cpp
3  * this tool will take 2 existing map files in nc format, and compare their sparse matrices
4  * we can compare xc, yc, areas, fractions with ncdiff from nco
5  * maybe there is another utility in nco, need to ask Charlie Zender
6  *
7  * example of usage:
8  * ./mbcmpmaps -i map1.nc -j map2.nc
9  * will look for row, col, S entries, and use eigen3 sparse matrix constructor
10  *
11  * can be built only if netcdf and eigen3 are available
12  *
13  *
14  */
15 #include "moab/MOABConfig.h"
16 
17 #ifndef MOAB_HAVE_EIGEN3
18 #error compareMaps tool requires eigen3 configuration
19 #endif
20 
21 #include "moab/ProgOptions.hpp"
22 
23 #include "netcdf.h"
24 #include <cmath>
25 #include <iomanip>
26 #include <queue>
27 #include <Eigen/Sparse>
28 #define ERR_NC( e ) \
29  { \
30  printf( "Error: %s\n", nc_strerror( e ) ); \
31  exit( 2 ); \
32  }
33 
34 // copy from ReadNCDF.cpp some useful macros for reading from a netcdf file (exodus?)
35 // ncFile is an integer initialized when opening the nc file in read mode
36 
37 int ncFile1;
38 
39 #define GET_DIM1( ncdim, name, val ) \
40  { \
41  int gdfail = nc_inq_dimid( ncFile1, name, &( ncdim ) ); \
42  if( NC_NOERR == gdfail ) \
43  { \
44  size_t tmp_val; \
45  gdfail = nc_inq_dimlen( ncFile1, ncdim, &tmp_val ); \
46  if( NC_NOERR != gdfail ) \
47  { \
48  ERR_NC( gdfail ) \
49  } \
50  else \
51  ( val ) = tmp_val; \
52  } \
53  else \
54  ( val ) = 0; \
55  }
56 
57 #define GET_VAR1( name, id, dims ) \
58  { \
59  ( id ) = -1; \
60  int gvfail = nc_inq_varid( ncFile1, name, &( id ) ); \
61  if( NC_NOERR == gvfail ) \
62  { \
63  int ndims; \
64  gvfail = nc_inq_varndims( ncFile1, id, &ndims ); \
65  if( NC_NOERR == gvfail ) \
66  { \
67  ( dims ).resize( ndims ); \
68  gvfail = nc_inq_vardimid( ncFile1, id, &( dims )[0] ); \
69  if( NC_NOERR != gvfail ) \
70  { \
71  ERR_NC( gvfail ) \
72  } \
73  } \
74  } \
75  }
76 
77 #define GET_1D_INT_VAR1( name, id, vals ) \
78  { \
79  GET_VAR1( name, id, vals ); \
80  if( -1 != ( id ) ) \
81  { \
82  size_t ntmp; \
83  int ivfail = nc_inq_dimlen( ncFile1, ( vals )[0], &ntmp ); \
84  if( NC_NOERR != ivfail ) \
85  { \
86  ERR_NC( ivfail ) \
87  } \
88  ( vals ).resize( ntmp ); \
89  size_t ntmp1 = 0; \
90  ivfail = nc_get_vara_int( ncFile1, id, &ntmp1, &ntmp, &( vals )[0] ); \
91  if( NC_NOERR != ivfail ) \
92  { \
93  ERR_NC( ivfail ) \
94  } \
95  } \
96  }
97 
98 #define GET_1D_DBL_VAR1( name, id, vals ) \
99  { \
100  std::vector< int > dum_dims; \
101  GET_VAR1( name, id, dum_dims ); \
102  if( -1 != ( id ) ) \
103  { \
104  size_t ntmp; \
105  int dvfail = nc_inq_dimlen( ncFile1, dum_dims[0], &ntmp ); \
106  if( NC_NOERR != dvfail ) \
107  { \
108  ERR_NC( dvfail ) \
109  } \
110  ( vals ).resize( ntmp ); \
111  size_t ntmp1 = 0; \
112  dvfail = nc_get_vara_double( ncFile1, id, &ntmp1, &ntmp, &( vals )[0] ); \
113  if( NC_NOERR != dvfail ) \
114  { \
115  ERR_NC( dvfail ) \
116  } \
117  } \
118  }
119 
120 int ncFile2;
121 
122 #define GET_DIM2( ncdim, name, val ) \
123  { \
124  int gdfail = nc_inq_dimid( ncFile2, name, &( ncdim ) ); \
125  if( NC_NOERR == gdfail ) \
126  { \
127  size_t tmp_val; \
128  gdfail = nc_inq_dimlen( ncFile2, ncdim, &tmp_val ); \
129  if( NC_NOERR != gdfail ) \
130  { \
131  ERR_NC( gdfail ) \
132  } \
133  else \
134  ( val ) = tmp_val; \
135  } \
136  else \
137  ( val ) = 0; \
138  }
139 
140 #define GET_VAR2( name, id, dims ) \
141  { \
142  ( id ) = -1; \
143  int gvfail = nc_inq_varid( ncFile2, name, &( id ) ); \
144  if( NC_NOERR == gvfail ) \
145  { \
146  int ndims; \
147  gvfail = nc_inq_varndims( ncFile2, id, &ndims ); \
148  if( NC_NOERR == gvfail ) \
149  { \
150  ( dims ).resize( ndims ); \
151  gvfail = nc_inq_vardimid( ncFile2, id, &( dims )[0] ); \
152  if( NC_NOERR != gvfail ) \
153  { \
154  ERR_NC( gvfail ) \
155  } \
156  } \
157  } \
158  }
159 
160 #define GET_1D_INT_VAR2( name, id, vals ) \
161  { \
162  GET_VAR2( name, id, vals ); \
163  if( -1 != ( id ) ) \
164  { \
165  size_t ntmp; \
166  int ivfail = nc_inq_dimlen( ncFile2, ( vals )[0], &ntmp ); \
167  if( NC_NOERR != ivfail ) \
168  { \
169  ERR_NC( ivfail ) \
170  } \
171  ( vals ).resize( ntmp ); \
172  size_t ntmp1 = 0; \
173  ivfail = nc_get_vara_int( ncFile2, id, &ntmp1, &ntmp, &( vals )[0] ); \
174  if( NC_NOERR != ivfail ) \
175  { \
176  ERR_NC( ivfail ) \
177  } \
178  } \
179  }
180 
181 #define GET_1D_DBL_VAR2( name, id, vals ) \
182  { \
183  std::vector< int > dum_dims; \
184  GET_VAR2( name, id, dum_dims ); \
185  if( -1 != ( id ) ) \
186  { \
187  size_t ntmp; \
188  int dvfail = nc_inq_dimlen( ncFile2, dum_dims[0], &ntmp ); \
189  if( NC_NOERR != dvfail ) \
190  { \
191  ERR_NC( dvfail ) \
192  } \
193  ( vals ).resize( ntmp ); \
194  size_t ntmp1 = 0; \
195  dvfail = nc_get_vara_double( ncFile2, id, &ntmp1, &ntmp, &( vals )[0] ); \
196  if( NC_NOERR != dvfail ) \
197  { \
198  ERR_NC( dvfail ) \
199  } \
200  } \
201  }
202 
203 #define GET_2D_DBL_VAR1( name, id, vals ) \
204  { \
205  std::vector< int > dum_dims; \
206  GET_VAR1( name, id, dum_dims ); \
207  if( -1 != ( id ) ) \
208  { \
209  size_t ntmp[2]; \
210  int dvfail = nc_inq_dimlen( ncFile1, dum_dims[0], &ntmp[0] ); \
211  if( NC_NOERR != dvfail ) \
212  { \
213  ERR_NC( dvfail ) \
214  } \
215  dvfail = nc_inq_dimlen( ncFile1, dum_dims[1], &ntmp[1] ); \
216  if( NC_NOERR != dvfail ) \
217  { \
218  ERR_NC( dvfail ) \
219  } \
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 ) \
224  { \
225  ERR_NC( dvfail ) \
226  } \
227  } \
228  }
229 
230 #define GET_2D_DBL_VAR2( name, id, vals ) \
231  { \
232  std::vector< int > dum_dims; \
233  GET_VAR2( name, id, dum_dims ); \
234  if( -1 != ( id ) ) \
235  { \
236  size_t ntmp[2]; \
237  int dvfail = nc_inq_dimlen( ncFile2, dum_dims[0], &ntmp[0] ); \
238  if( NC_NOERR != dvfail ) \
239  { \
240  ERR_NC( dvfail ) \
241  } \
242  dvfail = nc_inq_dimlen( ncFile2, dum_dims[1], &ntmp[1] ); \
243  if( NC_NOERR != dvfail ) \
244  { \
245  ERR_NC( dvfail ) \
246  } \
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 ) \
251  { \
252  ERR_NC( dvfail ) \
253  } \
254  } \
255  }
256 
257 typedef Eigen::Map< Eigen::VectorXd > EigenV;
258 
259 void diff_vect( const char* var_name, int n )
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 }
279 
280 void diff_2d_vect( const char* var_name, int n )
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 }
292 
293 typedef Eigen::Triplet< double > Triplet;
294 
295 struct CompareTriplets
296 {
297  bool operator()( Triplet& a, Triplet& b )
298  {
299  return fabs( a.value() ) > fabs( b.value() );
300  }
301 };
302 
303 int main( int argc, char* argv[] )
304 {
305 
306  ProgOptions opts;
307 
308  std::string inputfile1, inputfile2;
309  opts.addOpt< std::string >( "firstMap,i", "input filename 1", &inputfile1 );
310  opts.addOpt< std::string >( "secondMap,j", "input second map", &inputfile2 );
311  int print_diff = 10;
312  opts.addOpt< int >( "print_differences,p", "print differences ", &print_diff );
313 
314  opts.parseCommandLine( argc, argv );
315 
316  // Open netcdf/exodus file
317  int fail = nc_open( inputfile1.c_str(), 0, &ncFile1 );
318  if( NC_NOWRITE != fail )
319  {
320  ERR_NC( fail )
321  }
322 
323  std::cout << " opened " << inputfile1 << " for map 1 \n";
324 
325  fail = nc_open( inputfile2.c_str(), 0, &ncFile2 );
326  if( NC_NOWRITE != fail )
327  {
328  ERR_NC( fail )
329  }
330 
331  std::cout << " opened " << inputfile2 << " for map 2 \n";
332  int temp_dim;
333  int na1, nb1, ns1, na2, nb2, ns2, nv_a, nv_b, nv_a2, nv_b2;
334  GET_DIM1( temp_dim, "n_a", na1 );
335  GET_DIM2( temp_dim, "n_a", na2 );
336  GET_DIM1( temp_dim, "n_b", nb1 );
337  GET_DIM2( temp_dim, "n_b", nb2 );
338  GET_DIM1( temp_dim, "n_s", ns1 );
339  GET_DIM2( temp_dim, "n_s", ns2 );
340  GET_DIM1( temp_dim, "nv_a", nv_a );
341  GET_DIM2( temp_dim, "nv_a", nv_a2 );
342  GET_DIM1( temp_dim, "nv_b", nv_b );
343  GET_DIM2( temp_dim, "nv_b", nv_b2 );
344  if( nv_a != nv_a2 || nv_b != nv_b2 )
345  {
346  std::cout << " different nv dimensions:" << nv_a << " == " << nv_a2 << " or " << nv_b << " == " << nv_b2
347  << " bail out \n";
348  return 1;
349  }
350  std::cout << " n_a, n_b, n_s : " << na1 << ", " << nb1 << ", " << ns1 << " for map 1 \n";
351 
352  std::cout << " n_a, n_b, n_s : " << na2 << ", " << nb2 << ", " << ns2 << " for map 2 \n";
353  if( na1 != na2 || nb1 != nb2 )
354  {
355  std::cout << " different dimensions bail out \n";
356  return 1;
357  }
358  std::vector< int > col1( ns1 ), row1( ns1 );
359  std::vector< int > col2( ns2 ), row2( ns2 );
360 
361  std::vector< double > val1( ns1 ), val2( ns2 );
362 
363  int idrow1, idcol1, idrow2, idcol2, ids1, ids2;
364  GET_1D_INT_VAR1( "row", idrow1, row1 );
365  GET_1D_INT_VAR1( "col", idcol1, col1 );
366  GET_1D_DBL_VAR1( "S", ids1, val1 );
367  GET_1D_INT_VAR2( "row", idrow2, row2 );
368  GET_1D_INT_VAR2( "col", idcol2, col2 );
369  GET_1D_DBL_VAR2( "S", ids2, val2 );
370 
371  // first matrix
372 
373  std::vector< Triplet > tripletList;
374  tripletList.reserve( ns1 );
375  for( int iv = 0; iv < ns1; iv++ )
376  {
377  tripletList.push_back( Triplet( row1[iv] - 1, col1[iv] - 1, val1[iv] ) );
378  }
379  Eigen::SparseMatrix< double > weight1( nb1, na1 );
380 
381  weight1.setFromTriplets( tripletList.begin(), tripletList.end() );
382  weight1.makeCompressed();
383 
384  if( ns1 != ns2 ) tripletList.resize( ns2 );
385  for( int iv = 0; iv < ns2; iv++ )
386  {
387  tripletList[iv] = Triplet( row2[iv] - 1, col2[iv] - 1, val2[iv] );
388  }
389  Eigen::SparseMatrix< double > weight2( nb2, na2 );
390 
391  weight2.setFromTriplets( tripletList.begin(), tripletList.end() );
392  weight2.makeCompressed();
393 
394  // default storage type is column major
395  Eigen::SparseMatrix< double > diff = weight1 - weight2;
396  diff.makeCompressed(); // is it needed or not ?
397  auto coeffs = diff.coeffs();
398  double maxv = coeffs.maxCoeff();
399  double minv = coeffs.minCoeff();
400  std::cout << " euclidian norm for difference: " << diff.norm()
401  << " \n squared norm for difference: " << diff.squaredNorm() << "\n"
402  << " minv: " << minv << " maxv: " << maxv << "\n";
403 
404  // print out the largest 20 absolute values and position in diff sparse matrix
405  std::priority_queue< Triplet, std::vector< Triplet >, CompareTriplets > largestDiffs;
406 
407  for( int k = 0; ( k < diff.outerSize() ); ++k ) // this is by column
408  {
409  for( Eigen::SparseMatrix< double >::InnerIterator it( diff, k ); ( it ); ++it )
410  {
411  double val = it.value();
412  Triplet tp( it.row(), it.col(), val );
413  largestDiffs.push( tp );
414  if( largestDiffs.size() > print_diff ) largestDiffs.pop();
415  }
416  }
417  std::cout << std::setprecision( 16 );
418  int counter = 0;
419  while( !largestDiffs.empty() )
420  {
421  Triplet tp = largestDiffs.top();
422  largestDiffs.pop();
423 
424  int row = tp.row();
425  int col = tp.col();
426  std::cout << " counter:" << counter << "\t col: " << col + 1 << "\t row: " << row + 1
427  << "\t diff: " << tp.value();
428  std::cout << "\t map1: " << weight1.coeffRef( row, col ) << "\t map2: " << weight2.coeffRef( row, col )
429  << "\n"; // row index
430  counter++;
431  }
432 
433  // compare frac_a between maps
434  diff_vect( "frac_a", na1 );
435  diff_vect( "frac_b", nb1 );
436  diff_vect( "area_a", na1 );
437  diff_vect( "area_b", nb1 );
438  diff_vect( "yc_a", na1 );
439  diff_vect( "yc_b", nb1 );
440  diff_vect( "xc_a", na1 );
441  diff_vect( "xc_b", nb1 );
442  diff_2d_vect( "xv_a", na1 * nv_a );
443  diff_2d_vect( "yv_a", na1 * nv_a );
444  diff_2d_vect( "xv_b", nb1 * nv_b );
445  diff_2d_vect( "yv_b", nb1 * nv_b );
446 
447  return 0;
448 }