 |
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
|
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
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 }
1.9.1.