1 /*
2 This file is part of CDO. CDO is a collection of Operators to manipulate and analyse Climate model Data.
3
4 Author: Uwe Schulzweida
5
6 */
7
8 /*
9 This module contains the following operators:
10
11 Diff diff Compare two datasets
12 */
13
14 #include <map>
15 #include <algorithm>
16
17 #include <cdi.h>
18
19 #include "process_int.h"
20 #include "cdo_vlist.h"
21 #include "param_conversion.h"
22 #include "mpmo_color.h"
23 #include "cdo_options.h"
24 #include "printinfo.h"
25 #include "pmlist.h"
26 #include "cdo_zaxis.h"
27
28 static void
varListSynchronizeMemtype(VarList & varList1,VarList & varList2,const std::map<int,int> & mapOfVarIDs)29 varListSynchronizeMemtype(VarList &varList1, VarList &varList2, const std::map<int, int> &mapOfVarIDs)
30 {
31 const int nvars = varList1.size();
32 for (int varID1 = 0; varID1 < nvars; varID1++)
33 {
34 const auto it = mapOfVarIDs.find(varID1);
35 if (it == mapOfVarIDs.end()) continue;
36
37 const auto varID2 = it->second;
38
39 if (varList1[varID1].memType == MemType::Float && varList2[varID2].memType == MemType::Double)
40 varList1[varID1].memType = MemType::Double;
41 else if (varList1[varID1].memType == MemType::Double && varList2[varID2].memType == MemType::Float)
42 varList2[varID2].memType = MemType::Double;
43 }
44 }
45
46 template <typename T>
47 static void
diff_kernel(bool hasMissvals,const T v1,const T v2,const T missval1,const T missval2,size_t & ndiff,bool & dsgn,bool & zero,double & absm,double & relm)48 diff_kernel(bool hasMissvals, const T v1, const T v2, const T missval1, const T missval2, size_t &ndiff, bool &dsgn, bool &zero,
49 double &absm, double &relm)
50 {
51 const auto v1isnan = std::isnan(v1);
52 const auto v2isnan = std::isnan(v2);
53 const auto v1ismissval = hasMissvals ? DBL_IS_EQUAL(v1, missval1) : false;
54 const auto v2ismissval = hasMissvals ? DBL_IS_EQUAL(v2, missval2) : false;
55 if ((v1isnan && !v2isnan) || (!v1isnan && v2isnan))
56 {
57 ndiff++;
58 relm = 1.0;
59 }
60 else if (hasMissvals == false || (!v1ismissval && !v2ismissval))
61 {
62 const double absdiff = std::fabs(v1 - v2);
63 if (absdiff > 0.0) ndiff++;
64
65 absm = std::max(absm, absdiff);
66
67 const auto vv = v1 * v2;
68 if (vv < 0.0)
69 dsgn = true;
70 else if (IS_EQUAL(vv, 0.0))
71 zero = true;
72 else
73 relm = std::max(relm, absdiff / std::max(std::fabs(v1), std::fabs(v2)));
74 }
75 else if ((v1ismissval && !v2ismissval) || (!v1ismissval && v2ismissval))
76 {
77 ndiff++;
78 relm = 1.0;
79 }
80 }
81
82 static void
diff_kernel(size_t i,bool hasMissvals,const Field & field1,const Field & field2,size_t & ndiff,bool & dsgn,bool & zero,double & absm,double & relm)83 diff_kernel(size_t i, bool hasMissvals, const Field &field1, const Field &field2, size_t &ndiff, bool &dsgn, bool &zero,
84 double &absm, double &relm)
85 {
86 if (field1.memType == MemType::Float)
87 diff_kernel(hasMissvals, field1.vec_f[i], field2.vec_f[i], (float) field1.missval, (float) field2.missval, ndiff, dsgn, zero,
88 absm, relm);
89 else
90 diff_kernel(hasMissvals, field1.vec_d[i], field2.vec_d[i], field1.missval, field2.missval, ndiff, dsgn, zero, absm, relm);
91 }
92
93 static void
use_real_part(const size_t gridsize,Field & field)94 use_real_part(const size_t gridsize, Field &field)
95 {
96 if (field.memType == MemType::Float)
97 for (size_t i = 0; i < gridsize; ++i) field.vec_f[i] = field.vec_f[i * 2];
98 else
99 for (size_t i = 0; i < gridsize; ++i) field.vec_d[i] = field.vec_d[i * 2];
100 }
101
102 static void
diffGetParameter(double & abslim,double & abslim2,double & rellim,int & mapflag,int & maxcount)103 diffGetParameter(double &abslim, double &abslim2, double &rellim, int &mapflag, int &maxcount)
104 {
105 const auto pargc = cdo_operator_argc();
106 if (pargc)
107 {
108 const auto pargv = cdo_get_oper_argv();
109
110 KVList kvlist;
111 kvlist.name = "DIFF";
112 if (kvlist.parse_arguments(pargc, pargv) != 0) cdo_abort("Parse error!");
113 if (Options::cdoVerbose) kvlist.print();
114
115 for (const auto &kv : kvlist)
116 {
117 const auto &key = kv.key;
118 if (kv.nvalues > 1) cdo_abort("Too many values for parameter key >%s<!", key.c_str());
119 if (kv.nvalues < 1) cdo_abort("Missing value for parameter key >%s<!", key.c_str());
120 const auto &value = kv.values[0];
121
122 // clang-format off
123 if (key == "abslim") abslim = parameter_to_double(value);
124 else if (key == "abslim2") abslim2 = parameter_to_double(value);
125 else if (key == "rellim") rellim = parameter_to_double(value);
126 else if (key == "maxcount") maxcount = parameter_to_int(value);
127 else if (key == "names")
128 {
129 if (value == "left") mapflag = 1;
130 else if (value == "right") mapflag = 2;
131 else if (value == "intersect") mapflag = 3;
132 else cdo_abort("Invalid value for key >%s< (names=<left/right/intersect>)", key.c_str(), value.c_str());
133 }
134 else cdo_abort("Invalid parameter key >%s<!", key.c_str());
135 // clang-format on
136 }
137 }
138 }
139
140 void *
Diff(void * process)141 Diff(void *process)
142 {
143 auto printHeader = true;
144 int varID1, varID2 = -1;
145 int levelID;
146 int ndrec = 0, nd2rec = 0, ngrec = 0;
147 char paramstr[32];
148
149 cdo_initialize(process);
150
151 // clang-format off
152 const auto DIFF = cdo_operator_add("diff", 0, 0, nullptr);
153 const auto DIFFP = cdo_operator_add("diffp", 0, 0, nullptr);
154 const auto DIFFN = cdo_operator_add("diffn", 0, 0, nullptr);
155 const auto DIFFC = cdo_operator_add("diffc", 0, 0, nullptr);
156 // clang-format on
157
158 const auto operatorID = cdo_operator_id();
159
160 int mapflag = 0, maxcount = 0;
161 double abslim = 0.0, abslim2 = 1.e-3, rellim = 1.0;
162 diffGetParameter(abslim, abslim2, rellim, mapflag, maxcount);
163
164 if (rellim < -1.e33 || rellim > 1.e+33) cdo_abort("Rel. limit out of range!");
165 if (abslim < -1.e33 || abslim > 1.e+33) cdo_abort("Abs. limit out of range!");
166 if (abslim2 < -1.e33 || abslim2 > 1.e+33) cdo_abort("Abs2. limit out of range!");
167
168 const auto streamID1 = cdo_open_read(0);
169 const auto streamID2 = cdo_open_read(1);
170
171 const auto vlistID1 = cdo_stream_inq_vlist(streamID1);
172 const auto vlistID2 = cdo_stream_inq_vlist(streamID2);
173
174 const auto nvars = vlistNvars(vlistID1);
175 std::map<int, int> mapOfVarIDs;
176
177 if (mapflag == 0)
178 {
179 vlist_compare(vlistID1, vlistID2, CMP_ALL);
180 for (int varID = 0; varID < nvars; ++varID) mapOfVarIDs[varID] = varID;
181 }
182 else
183 {
184 vlist_map(vlistID1, vlistID2, CMP_ALL, mapflag, mapOfVarIDs);
185 }
186
187 VarList varList1, varList2;
188 varListInit(varList1, vlistID1);
189 varListInit(varList2, vlistID2);
190 varListSynchronizeMemtype(varList1, varList2, mapOfVarIDs);
191
192 Field field1, field2;
193
194 const auto taxisID = vlistInqTaxis(vlistID1);
195
196 int nrecs, nrecs2;
197 int indg = 0;
198 int tsID = 0;
199 while (true)
200 {
201 auto stopLoop = false;
202
203 nrecs = cdo_stream_inq_timestep(streamID1, tsID);
204 const auto vdateString = date_to_string(taxisInqVdate(taxisID));
205 const auto vtimeString = time_to_string(taxisInqVtime(taxisID));
206
207 nrecs2 = cdo_stream_inq_timestep(streamID2, tsID);
208
209 if (nrecs == 0 || nrecs2 == 0) break;
210
211 int recID2next = 0;
212
213 for (int recID = 0; recID < nrecs; recID++)
214 {
215 cdo_inq_record(streamID1, &varID1, &levelID);
216
217 auto it = mapOfVarIDs.find(varID1);
218 if (it == mapOfVarIDs.end())
219 {
220 if (mapflag == 2 || mapflag == 3) continue;
221 cdo_abort("Internal problem: varID1=%d not found!", varID1);
222 }
223
224 for (; recID2next < nrecs2; ++recID2next)
225 {
226 cdo_inq_record(streamID2, &varID2, &levelID);
227 if (it->second == varID2)
228 {
229 ++recID2next;
230 break;
231 }
232 }
233
234 if (it->second != varID2 && recID2next == nrecs2)
235 cdo_abort("Internal problem: varID2=%d not found in second stream!", it->second);
236
237 indg += 1;
238
239 const auto gridsize = varList1[varID1].gridsize;
240
241 // checkrel = gridInqType(gridID) != GRID_SPECTRAL;
242 const auto checkrel = true;
243
244 cdiParamToString(varList1[varID1].param, paramstr, sizeof(paramstr));
245
246 field1.init(varList1[varID1]);
247 cdo_read_record(streamID1, field1);
248 if (varList1[varID1].nwpv == CDI_COMP) use_real_part(gridsize, field1);
249
250 field2.init(varList2[varID2]);
251 cdo_read_record(streamID2, field2);
252 if (varList2[varID2].nwpv == CDI_COMP) use_real_part(gridsize, field2);
253
254 const auto hasMissvals = (field1.nmiss || field2.nmiss);
255 size_t ndiff = 0;
256 auto dsgn = false, zero = false;
257 double absm = 0.0, relm = 0.0;
258
259 for (size_t i = 0; i < gridsize; i++)
260 {
261 diff_kernel(i, hasMissvals, field1, field2, ndiff, dsgn, zero, absm, relm);
262 }
263
264 if (!Options::silentMode || Options::cdoVerbose)
265 {
266 if (absm > abslim || (checkrel && relm >= rellim) || Options::cdoVerbose)
267 {
268 if (printHeader)
269 {
270 printHeader = false;
271
272 fprintf(stdout, " Date Time Level Gridsize Miss ");
273 fprintf(stdout, " Diff ");
274 fprintf(stdout, ": S Z Max_Absdiff Max_Reldiff : ");
275
276 if (operatorID == DIFFN)
277 fprintf(stdout, "Parameter name");
278 else if (operatorID == DIFF || operatorID == DIFFP)
279 fprintf(stdout, "Parameter ID");
280 else if (operatorID == DIFFC)
281 fprintf(stdout, "Code number");
282
283 fprintf(stdout, "\n");
284 }
285
286 fprintf(stdout, "%6d ", indg);
287 fprintf(stdout, ":");
288
289 set_text_color(stdout, MAGENTA);
290 fprintf(stdout, "%s %s ", vdateString.c_str(), vtimeString.c_str());
291 reset_text_color(stdout);
292 set_text_color(stdout, GREEN);
293 fprintf(stdout, "%7g ", cdo_zaxis_inq_level(varList1[varID1].zaxisID, levelID));
294 fprintf(stdout, "%8zu %7zu ", gridsize, std::max(field1.nmiss, field2.nmiss));
295 fprintf(stdout, "%7zu ", ndiff);
296 reset_text_color(stdout);
297
298 fprintf(stdout, ":");
299 fprintf(stdout, " %c %c ", dsgn ? 'T' : 'F', zero ? 'T' : 'F');
300 set_text_color(stdout, BLUE);
301 fprintf(stdout, "%#12.5g%#12.5g", absm, relm);
302 reset_text_color(stdout);
303 fprintf(stdout, " : ");
304
305 set_text_color(stdout, BRIGHT, GREEN);
306 if (operatorID == DIFFN)
307 fprintf(stdout, "%-11s", varList1[varID1].name);
308 else if (operatorID == DIFF || operatorID == DIFFP)
309 fprintf(stdout, "%-11s", paramstr);
310 else if (operatorID == DIFFC)
311 fprintf(stdout, "%4d", varList1[varID1].code);
312 reset_text_color(stdout);
313
314 fprintf(stdout, "\n");
315 }
316 }
317
318 ngrec++;
319 if (absm > abslim || (checkrel && relm >= rellim)) ndrec++;
320 if (absm > abslim2 || (checkrel && relm >= rellim)) nd2rec++;
321
322 if (maxcount > 0 && ndrec >= maxcount)
323 {
324 stopLoop = true;
325 break;
326 }
327 }
328
329 if (stopLoop) break;
330
331 tsID++;
332 }
333
334 if (ndrec > 0)
335 {
336 Options::cdoExitStatus = 1;
337
338 set_text_color(stdout, BRIGHT, RED);
339 fprintf(stdout, " %d of %d records differ", ndrec, ngrec);
340 reset_text_color(stdout);
341 fprintf(stdout, "\n");
342
343 if (ndrec != nd2rec && abslim < abslim2) fprintf(stdout, " %d of %d records differ more than %g\n", nd2rec, ngrec, abslim2);
344 // fprintf(stdout, " %d of %d records differ more then one thousandth\n", nprec, ngrec);
345 }
346
347 if (nrecs == 0 && nrecs2 > 0) cdo_warning("stream2 has more time steps than stream1!");
348 if (nrecs > 0 && nrecs2 == 0) cdo_warning("stream1 has more time steps than stream2!");
349
350 cdo_stream_close(streamID1);
351 cdo_stream_close(streamID2);
352
353 cdo_finish();
354
355 return nullptr;
356 }
357