Mesh Oriented datABase  (version 5.5.1)
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 <Eigen/Sparse>
27 #define ERR_NC( e ) \
28  { \
29  printf( "Error: %s\n", nc_strerror( e ) ); \
30  exit( 2 ); \
31  }
32 
33 // copy from ReadNCDF.cpp some useful macros for reading from a netcdf file (exodus?)
34 // ncFile is an integer initialized when opening the nc file in read mode
35 
36 int ncFile1;
37 
38 #define GET_DIM1( ncdim, name, val ) \
39  { \
40  int gdfail = nc_inq_dimid( ncFile1, name, &( ncdim ) ); \
41  if( NC_NOERR == gdfail ) \
42  { \
43  size_t tmp_val; \
44  gdfail = nc_inq_dimlen( ncFile1, ncdim, &tmp_val ); \
45  if( NC_NOERR != gdfail ) \
46  { \
47  ERR_NC( gdfail ) \
48  } \
49  else \
50  ( val ) = tmp_val; \
51  } \
52  else \
53  ( val ) = 0; \
54  }
55 
56 #define GET_VAR1( name, id, dims ) \
57  { \
58  ( id ) = -1; \
59  int gvfail = nc_inq_varid( ncFile1, name, &( id ) ); \
60  if( NC_NOERR == gvfail ) \
61  { \
62  int ndims; \
63  gvfail = nc_inq_varndims( ncFile1, id, &ndims ); \
64  if( NC_NOERR == gvfail ) \
65  { \
66  ( dims ).resize( ndims ); \
67  gvfail = nc_inq_vardimid( ncFile1, id, &( dims )[0] ); \
68  if( NC_NOERR != gvfail ) \
69  { \
70  ERR_NC( gvfail ) \
71  } \
72  } \
73  } \
74  }
75 
76 #define GET_1D_INT_VAR1( name, id, vals ) \
77  { \
78  GET_VAR1( name, id, vals ); \
79  if( -1 != ( id ) ) \
80  { \
81  size_t ntmp; \
82  int ivfail = nc_inq_dimlen( ncFile1, ( vals )[0], &ntmp ); \
83  if( NC_NOERR != ivfail ) \
84  { \
85  ERR_NC( ivfail ) \
86  } \
87  ( vals ).resize( ntmp ); \
88  size_t ntmp1 = 0; \
89  ivfail = nc_get_vara_int( ncFile1, id, &ntmp1, &ntmp, &( vals )[0] ); \
90  if( NC_NOERR != ivfail ) \
91  { \
92  ERR_NC( ivfail ) \
93  } \
94  } \
95  }
96 
97 #define GET_1D_DBL_VAR1( name, id, vals ) \
98  { \
99  std::vector< int > dum_dims; \
100  GET_VAR1( name, id, dum_dims ); \
101  if( -1 != ( id ) ) \
102  { \
103  size_t ntmp; \
104  int dvfail = nc_inq_dimlen( ncFile1, dum_dims[0], &ntmp ); \
105  if( NC_NOERR != dvfail ) \
106  { \
107  ERR_NC( dvfail ) \
108  } \
109  ( vals ).resize( ntmp ); \
110  size_t ntmp1 = 0; \
111  dvfail = nc_get_vara_double( ncFile1, id, &ntmp1, &ntmp, &( vals )[0] ); \
112  if( NC_NOERR != dvfail ) \
113  { \
114  ERR_NC( dvfail ) \
115  } \
116  } \
117  }
118 
119 int ncFile2;
120 
121 #define GET_DIM2( ncdim, name, val ) \
122  { \
123  int gdfail = nc_inq_dimid( ncFile2, name, &( ncdim ) ); \
124  if( NC_NOERR == gdfail ) \
125  { \
126  size_t tmp_val; \
127  gdfail = nc_inq_dimlen( ncFile2, ncdim, &tmp_val ); \
128  if( NC_NOERR != gdfail ) \
129  { \
130  ERR_NC( gdfail ) \
131  } \
132  else \
133  ( val ) = tmp_val; \
134  } \
135  else \
136  ( val ) = 0; \
137  }
138 
139 #define GET_VAR2( name, id, dims ) \
140  { \
141  ( id ) = -1; \
142  int gvfail = nc_inq_varid( ncFile2, name, &( id ) ); \
143  if( NC_NOERR == gvfail ) \
144  { \
145  int ndims; \
146  gvfail = nc_inq_varndims( ncFile2, id, &ndims ); \
147  if( NC_NOERR == gvfail ) \
148  { \
149  ( dims ).resize( ndims ); \
150  gvfail = nc_inq_vardimid( ncFile2, id, &( dims )[0] ); \
151  if( NC_NOERR != gvfail ) \
152  { \
153  ERR_NC( gvfail ) \
154  } \
155  } \
156  } \
157  }
158 
159 #define GET_1D_INT_VAR2( name, id, vals ) \
160  { \
161  GET_VAR2( name, id, vals ); \
162  if( -1 != ( id ) ) \
163  { \
164  size_t ntmp; \
165  int ivfail = nc_inq_dimlen( ncFile2, ( vals )[0], &ntmp ); \
166  if( NC_NOERR != ivfail ) \
167  { \
168  ERR_NC( ivfail ) \
169  } \
170  ( vals ).resize( ntmp ); \
171  size_t ntmp1 = 0; \
172  ivfail = nc_get_vara_int( ncFile2, id, &ntmp1, &ntmp, &( vals )[0] ); \
173  if( NC_NOERR != ivfail ) \
174  { \
175  ERR_NC( ivfail ) \
176  } \
177  } \
178  }
179 
180 #define GET_1D_DBL_VAR2( name, id, vals ) \
181  { \
182  std::vector< int > dum_dims; \
183  GET_VAR2( name, id, dum_dims ); \
184  if( -1 != ( id ) ) \
185  { \
186  size_t ntmp; \
187  int dvfail = nc_inq_dimlen( ncFile2, dum_dims[0], &ntmp ); \
188  if( NC_NOERR != dvfail ) \
189  { \
190  ERR_NC( dvfail ) \
191  } \
192  ( vals ).resize( ntmp ); \
193  size_t ntmp1 = 0; \
194  dvfail = nc_get_vara_double( ncFile2, id, &ntmp1, &ntmp, &( vals )[0] ); \
195  if( NC_NOERR != dvfail ) \
196  { \
197  ERR_NC( dvfail ) \
198  } \
199  } \
200  }
201 
202 #define GET_2D_DBL_VAR1( name, id, vals ) \
203  { \
204  std::vector< int > dum_dims; \
205  GET_VAR1( name, id, dum_dims ); \
206  if( -1 != ( id ) ) \
207  { \
208  size_t ntmp[2]; \
209  int dvfail = nc_inq_dimlen( ncFile1, dum_dims[0], &ntmp[0] ); \
210  if( NC_NOERR != dvfail ) \
211  { \
212  ERR_NC( dvfail ) \
213  } \
214  dvfail = nc_inq_dimlen( ncFile1, dum_dims[1], &ntmp[1] ); \
215  if( NC_NOERR != dvfail ) \
216  { \
217  ERR_NC( dvfail ) \
218  } \
219  ( vals ).resize( ntmp[0] * ntmp[1] ); \
220  size_t ntmp1[2] = { 0, 0 }; \
221  dvfail = nc_get_vara_double( ncFile1, id, ntmp1, ntmp, &( vals )[0] ); \
222  if( NC_NOERR != dvfail ) \
223  { \
224  ERR_NC( dvfail ) \
225  } \
226  } \
227  }
228 
229 #define GET_2D_DBL_VAR2( name, id, vals ) \
230  { \
231  std::vector< int > dum_dims; \
232  GET_VAR2( name, id, dum_dims ); \
233  if( -1 != ( id ) ) \
234  { \
235  size_t ntmp[2]; \
236  int dvfail = nc_inq_dimlen( ncFile2, dum_dims[0], &ntmp[0] ); \
237  if( NC_NOERR != dvfail ) \
238  { \
239  ERR_NC( dvfail ) \
240  } \
241  dvfail = nc_inq_dimlen( ncFile2, dum_dims[1], &ntmp[1] ); \
242  if( NC_NOERR != dvfail ) \
243  { \
244  ERR_NC( dvfail ) \
245  } \
246  ( vals ).resize( ntmp[0] * ntmp[1] ); \
247  size_t ntmp1[2] = { 0, 0 }; \
248  dvfail = nc_get_vara_double( ncFile2, id, ntmp1, ntmp, &( vals )[0] ); \
249  if( NC_NOERR != dvfail ) \
250  { \
251  ERR_NC( dvfail ) \
252  } \
253  } \
254  }
255 
256 typedef Eigen::Map< Eigen::VectorXd > EigenV;
257 
258 void diff_vect( const char* var_name, int n )
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 }
278 
279 void diff_2d_vect( const char* var_name, int n )
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 }
291 int main( int argc, char* argv[] )
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 }