Diff.cc 12 KB
Newer Older
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1
2
3
4
/*
  This file is part of CDO. CDO is a collection of Operators to
  manipulate and analyse Climate model Data.

Uwe Schulzweida's avatar
Uwe Schulzweida committed
5
  Copyright (C) 2003-2020 Uwe Schulzweida, <uwe.schulzweida AT mpimet.mpg.de>
Uwe Schulzweida's avatar
Uwe Schulzweida committed
6
7
8
9
10
11
12
13
14
15
16
17
18
  See COPYING file for copying and redistribution conditions.

  This program is free software; you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
  the Free Software Foundation; version 2 of the License.

  This program is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU General Public License for more details.
*/

/*
19
   This module contains the following operators:
Uwe Schulzweida's avatar
Uwe Schulzweida committed
20

21
22
      Diff       diff            Compare two datasets
*/
Uwe Schulzweida's avatar
Uwe Schulzweida committed
23

Uwe Schulzweida's avatar
Uwe Schulzweida committed
24
#include <map>
25
#include <algorithm>
26

27
28
#include <cdi.h>

29
#include "functs.h"
30
#include "process_int.h"
31
#include "cdo_vlist.h"
32
#include "param_conversion.h"
33
#include "mpmo_color.h"
34
#include "cdo_options.h"
Oliver Heidmann's avatar
Oliver Heidmann committed
35
#include "printinfo.h"
Oliver Heidmann's avatar
Oliver Heidmann committed
36
#include "pmlist.h"
37
#include "cdo_zaxis.h"
Uwe Schulzweida's avatar
Uwe Schulzweida committed
38

39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111

static void
varListSynchronizeMemtype(VarList &varList1, VarList &varList2, const std::map<int, int> &mapOfVarIDs)
{
  int nvars = varList1.size();
  for (int varID1 = 0; varID1 < nvars; varID1++)
    {
      auto it = mapOfVarIDs.find(varID1);
      if (it == mapOfVarIDs.end()) continue;

      auto varID2 = it->second;

      if (varList1[varID1].memType == MemType::Float && varList2[varID2].memType == MemType::Double)
        varList1[varID1].memType = MemType::Double;
      else if (varList1[varID1].memType == MemType::Double && varList2[varID2].memType == MemType::Float)
        varList2[varID2].memType = MemType::Double;
    }
}

template <typename T>
static void
diff_kernel(bool lmiss, const T v1, const T v2, const T missval1, const T missval2,
            size_t &ndiff, bool &dsgn, bool &zero, double &absm, double &relm)
{
  const auto v1isnan = std::isnan(v1);
  const auto v2isnan = std::isnan(v2);
  const auto v1ismissval = lmiss ? DBL_IS_EQUAL(v1, missval1) : false;
  const auto v2ismissval = lmiss ? DBL_IS_EQUAL(v2, missval2) : false;
  if ((v1isnan && !v2isnan) || (!v1isnan && v2isnan))
    {
      ndiff++;
      relm = 1.0;
    }
  else if (lmiss == false || (!v1ismissval && !v2ismissval))
    {
      const double absdiff = std::fabs(v1 - v2);
      if (absdiff > 0.) ndiff++;

      absm = std::max(absm, absdiff);

      const auto vv = v1 * v2;
      if (vv < 0.)
        dsgn = true;
      else if (IS_EQUAL(vv, 0.))
        zero = true;
      else
        relm = std::max(relm, absdiff / std::max(std::fabs(v1), std::fabs(v2)));
    }
  else if ((v1ismissval && !v2ismissval) || (!v1ismissval && v2ismissval))
    {
      ndiff++;
      relm = 1.0;
    }
}

static void
diff_kernel(size_t i, bool lmiss, const Field &field1, const Field &field2, size_t &ndiff, bool &dsgn, bool &zero, double &absm, double &relm)
{
  if (field1.memType == MemType::Float)
    diff_kernel(lmiss, field1.vec_f[i], field2.vec_f[i], (float)field1.missval, (float)field2.missval, ndiff, dsgn, zero, absm, relm);
  else
    diff_kernel(lmiss, field1.vec_d[i], field2.vec_d[i], field1.missval, field2.missval, ndiff, dsgn, zero, absm, relm);
}

static void
use_real_part(size_t gridsize, Field &field)
{
  if (field.memType == MemType::Float)
    for (size_t i = 0; i < gridsize; ++i) field.vec_f[i] = field.vec_f[i * 2];
  else
    for (size_t i = 0; i < gridsize; ++i) field.vec_d[i] = field.vec_d[i * 2];
}

112
static void
113
diffGetParameter(double &abslim, double &abslim2, double &rellim, int &mapflag, int &maxcount)
114
{
115
  const auto pargc = operatorArgc();
116
117
  if (pargc)
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
118
      auto pargv = cdoGetOperArgv();
119

Uwe Schulzweida's avatar
Uwe Schulzweida committed
120
121
122
123
      KVList kvlist;
      kvlist.name = "DIFF";
      if (kvlist.parseArguments(pargc, pargv) != 0) cdoAbort("Parse error!");
      if (Options::cdoVerbose) kvlist.print();
124

Uwe Schulzweida's avatar
Uwe Schulzweida committed
125
      for (const auto &kv : kvlist)
126
        {
127
128
129
130
          const auto &key = kv.key;
          if (kv.nvalues > 1) cdoAbort("Too many values for parameter key >%s<!", key.c_str());
          if (kv.nvalues < 1) cdoAbort("Missing value for parameter key >%s<!", key.c_str());
          const auto &value = kv.values[0];
131
132

          // clang-format off
133
134
135
136
          if      (key == "abslim")   abslim = parameter2double(value);
          else if (key == "abslim2")  abslim2 = parameter2double(value);
          else if (key == "rellim")   rellim = parameter2double(value);
          else if (key == "maxcount") maxcount = parameter2int(value);
137
          else if (key == "names")
138
            {
139
140
141
142
              if      (value == "left")      mapflag = 1;
              else if (value == "right")     mapflag = 2;
              else if (value == "intersect") mapflag = 3;
              else cdoAbort("Invalid value for key >%s< (names=<left/right/intersect>)", key.c_str(), value.c_str());
143
            }
144
          else cdoAbort("Invalid parameter key >%s<!", key.c_str());
145
146
147
148
149
          // clang-format on
        }
    }
}

150
151
void *
Diff(void *process)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
152
{
153
  bool lhead = true;
154
  int nrecs, nrecs2;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
155
  int varID1, varID2 = -1;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
156
  int levelID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
157
  int ndrec = 0, nd2rec = 0, ngrec = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
158
  char paramstr[32];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
159

160
  cdoInitialize(process);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
161

Uwe Schulzweida's avatar
Uwe Schulzweida committed
162
  // clang-format off
163
164
165
166
  const auto DIFF  = cdoOperatorAdd("diff",  0, 0, nullptr);
  const auto DIFFP = cdoOperatorAdd("diffp", 0, 0, nullptr);
  const auto DIFFN = cdoOperatorAdd("diffn", 0, 0, nullptr);
  const auto DIFFC = cdoOperatorAdd("diffc", 0, 0, nullptr);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
167
  // clang-format on
Uwe Schulzweida's avatar
Uwe Schulzweida committed
168

169
  const auto operatorID = cdoOperatorID();
Uwe Schulzweida's avatar
Uwe Schulzweida committed
170

171
172
173
  int mapflag = 0, maxcount = 0;
  double abslim = 0., abslim2 = 1.e-3, rellim = 1.0;
  diffGetParameter(abslim, abslim2, rellim, mapflag, maxcount);
174

175
  if (rellim < -1.e33 || rellim > 1.e+33) cdoAbort("Rel. limit out of range!");
176
177
  if (abslim < -1.e33 || abslim > 1.e+33) cdoAbort("Abs. limit out of range!");
  if (abslim2 < -1.e33 || abslim2 > 1.e+33) cdoAbort("Abs2. limit out of range!");
178

179
180
  const auto streamID1 = cdoOpenRead(0);
  const auto streamID2 = cdoOpenRead(1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
181

182
183
  const auto vlistID1 = cdoStreamInqVlist(streamID1);
  const auto vlistID2 = cdoStreamInqVlist(streamID2);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
184

185
  const auto nvars = vlistNvars(vlistID1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
186
187
  std::map<int, int> mapOfVarIDs;

Uwe Schulzweida's avatar
Uwe Schulzweida committed
188
  if (mapflag == 0)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
189
190
191
192
    {
      vlistCompare(vlistID1, vlistID2, CMP_ALL);
      for (int varID = 0; varID < nvars; ++varID) mapOfVarIDs[varID] = varID;
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
193
194
195
196
  else
    {
      vlistMap(vlistID1, vlistID2, CMP_ALL, mapflag, mapOfVarIDs);
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
197

198
199
200
201
  VarList varList1, varList2;
  varListInit(varList1, vlistID1);
  varListInit(varList2, vlistID2);
  varListSynchronizeMemtype(varList1, varList2, mapOfVarIDs);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
202

203
  Field field1, field2;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
204

205
  const auto taxisID = vlistInqTaxis(vlistID1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
206

207
208
  int indg = 0;
  int tsID = 0;
209
  while (true)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
210
    {
211
212
      bool lstop = false;

213
      nrecs = cdoStreamInqTimestep(streamID1, tsID);
214
215
      const auto vdateString = dateToString(taxisInqVdate(taxisID));
      const auto vtimeString = timeToString(taxisInqVtime(taxisID));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
216

217
      nrecs2 = cdoStreamInqTimestep(streamID2, tsID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
218

219
      if (nrecs == 0 || nrecs2 == 0) break;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
220

Uwe Schulzweida's avatar
Uwe Schulzweida committed
221
222
      int recID2next = 0;

223
224
      for (int recID = 0; recID < nrecs; recID++)
        {
225
          cdoInqRecord(streamID1, &varID1, &levelID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
226

Uwe Schulzweida's avatar
Uwe Schulzweida committed
227
228
229
          auto it = mapOfVarIDs.find(varID1);
          if (it == mapOfVarIDs.end())
            {
230
              if (mapflag == 2 || mapflag == 3) continue;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
231
232
233
234
235
              cdoAbort("Internal problem: varID1=%d not found!", varID1);
            }

          for (; recID2next < nrecs2; ++recID2next)
            {
236
              cdoInqRecord(streamID2, &varID2, &levelID);
237
238
239
240
241
              if (it->second == varID2)
                {
                  ++recID2next;
                  break;
                }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
242
            }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
243

244
245
          if (it->second != varID2 && recID2next == nrecs2)
            cdoAbort("Internal problem: varID2=%d not found in second stream!", it->second);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
246

247
          indg += 1;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
248

249
          const auto gridsize = varList1[varID1].gridsize;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
250

251
          // checkrel = gridInqType(gridID) != GRID_SPECTRAL;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
252
          const bool checkrel = true;
253

254
          cdiParamToString(varList1[varID1].param, paramstr, sizeof(paramstr));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
255

256
257
258
          field1.init(varList1[varID1]);
          cdoReadRecord(streamID1, field1);
          if (varList1[varID1].nwpv == CDI_COMP) use_real_part(gridsize, field1);
259

260
261
262
          field2.init(varList2[varID2]);
          cdoReadRecord(streamID2, field2);
          if (varList2[varID2].nwpv == CDI_COMP) use_real_part(gridsize, field2);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
263

264
265
          bool lmiss = field1.nmiss || field2.nmiss;
          size_t ndiff = 0;
266
267
          bool dsgn = false, zero = false;
          double absm = 0.0, relm = 0.0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
268

269
270
          for (size_t i = 0; i < gridsize; i++)
            {
271
              diff_kernel(i, lmiss, field1, field2, ndiff, dsgn, zero, absm, relm);
272
273
            }

274
          if (!Options::silentMode || Options::cdoVerbose)
275
            {
276
              if (absm > abslim || (checkrel && relm >= rellim) || Options::cdoVerbose)
277
278
279
280
281
                {
                  if (lhead)
                    {
                      lhead = false;

282
                      fprintf(stdout, "               Date     Time   Level Gridsize    Miss ");
283
                      fprintf(stdout, "   Diff ");
284
                      fprintf(stdout, ": S Z  Max_Absdiff Max_Reldiff : ");
285
286

                      if (operatorID == DIFFN)
287
                        fprintf(stdout, "Parameter name");
288
                      else if (operatorID == DIFF || operatorID == DIFFP)
289
                        fprintf(stdout, "Parameter ID");
290
                      else if (operatorID == DIFFC)
291
                        fprintf(stdout, "Code number");
292
293
294
295
296
297
298

                      fprintf(stdout, "\n");
                    }

                  fprintf(stdout, "%6d ", indg);
                  fprintf(stdout, ":");

299
                  set_text_color(stdout, MAGENTA);
300
                  fprintf(stdout, "%s %s ", vdateString.c_str(), vtimeString.c_str());
301
                  reset_text_color(stdout);
302
                  set_text_color(stdout, GREEN);
303
304
305
                  fprintf(stdout, "%7g ", cdoZaxisInqLevel(varList1[varID1].zaxisID, levelID));
                  fprintf(stdout, "%8zu %7zu ", gridsize, std::max(field1.nmiss, field2.nmiss));
                  fprintf(stdout, "%7zu ", ndiff);
306
307
308
                  reset_text_color(stdout);

                  fprintf(stdout, ":");
Uwe Schulzweida's avatar
Uwe Schulzweida committed
309
                  fprintf(stdout, " %c %c ", dsgn ? 'T' : 'F', zero ? 'T' : 'F');
310
                  set_text_color(stdout, BLUE);
311
312
                  fprintf(stdout, "%#12.5g%#12.5g", absm, relm);
                  reset_text_color(stdout);
313
                  fprintf(stdout, " : ");
314
315
316

                  set_text_color(stdout, BRIGHT, GREEN);
                  if (operatorID == DIFFN)
317
                    fprintf(stdout, "%-11s", varList1[varID1].name);
318
319
320
                  else if (operatorID == DIFF || operatorID == DIFFP)
                    fprintf(stdout, "%-11s", paramstr);
                  else if (operatorID == DIFFC)
321
                    fprintf(stdout, "%4d", varList1[varID1].code);
322
323
324
325
326
327
328
329
330
                  reset_text_color(stdout);

                  fprintf(stdout, "\n");
                }
            }

          ngrec++;
          if (absm > abslim || (checkrel && relm >= rellim)) ndrec++;
          if (absm > abslim2 || (checkrel && relm >= rellim)) nd2rec++;
331
332
333
334
335
336

          if (maxcount > 0 && ndrec >= maxcount)
            {
              lstop = true;
              break;
            }
337
        }
338
339
340

      if (lstop) break;

Uwe Schulzweida's avatar
Uwe Schulzweida committed
341
342
343
      tsID++;
    }

344
  if (ndrec > 0)
345
    {
346
347
      Options::cdoExitStatus = 1;

Uwe Schulzweida's avatar
Uwe Schulzweida committed
348
349
350
351
352
      set_text_color(stdout, BRIGHT, RED);
      fprintf(stdout, "  %d of %d records differ", ndrec, ngrec);
      reset_text_color(stdout);
      fprintf(stdout, "\n");

353
      if (ndrec != nd2rec && abslim < abslim2) fprintf(stdout, "  %d of %d records differ more than %g\n", nd2rec, ngrec, abslim2);
354
      //  fprintf(stdout, "  %d of %d records differ more then one thousandth\n", nprec, ngrec);
355
356
    }

357
358
  if (nrecs == 0 && nrecs2 > 0) cdoWarning("stream2 has more time steps than stream1!");
  if (nrecs > 0 && nrecs2 == 0) cdoWarning("stream1 has more time steps than stream2!");
Uwe Schulzweida's avatar
Uwe Schulzweida committed
359

360
361
  cdoStreamClose(streamID1);
  cdoStreamClose(streamID2);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
362
363
364

  cdoFinish();

Uwe Schulzweida's avatar
Uwe Schulzweida committed
365
  return nullptr;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
366
}