Diff.cc 7.35 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-2017 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
24


Ralf Mueller's avatar
Ralf Mueller committed
25
#include <cdi.h>
26
27
28
#include "cdo.h"
#include "cdo_int.h"
#include "pstream.h"
Uwe Schulzweida's avatar
Uwe Schulzweida committed
29
30
31
32


void *Diff(void *argument)
{
33
  bool lhead = true;
34
  int nrecs, nrecs2;
35
  int varID1, varID2;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
36
37
  int levelID;
  int nmiss1, nmiss2;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
38
  int ndrec = 0, nd2rec = 0, ngrec = 0;
39
  char varname[CDI_MAX_NAME];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
40
  char paramstr[32];
41
  char vdatestr[32], vtimestr[32];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
42
  double abslim = 0., abslim2 = 1.e-3, rellim = 1.0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
43
44
45

  cdoInitialize(argument);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
46
  // clang-format off
47
48
49
50
  int DIFF  = cdoOperatorAdd("diff",  0, 0, NULL);
  int DIFFP = cdoOperatorAdd("diffp", 0, 0, NULL);
  int DIFFN = cdoOperatorAdd("diffn", 0, 0, NULL);
  int DIFFC = cdoOperatorAdd("diffc", 0, 0, NULL);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
51
  // clang-format on
Uwe Schulzweida's avatar
Uwe Schulzweida committed
52

53
  int operatorID = cdoOperatorID();
Uwe Schulzweida's avatar
Uwe Schulzweida committed
54

Uwe Schulzweida's avatar
Uwe Schulzweida committed
55
  if ( operatorArgc() >= 1 ) abslim = parameter2double(operatorArgv()[0]);
56
  if ( abslim < -1.e33 || abslim > 1.e+33 ) cdoAbort("Abs. limit out of range!");
Uwe Schulzweida's avatar
Uwe Schulzweida committed
57
  if ( operatorArgc() == 2 ) rellim = parameter2double(operatorArgv()[1]);
58
  if ( rellim < -1.e33 || rellim > 1.e+33 ) cdoAbort("Rel. limit out of range!");
59

60
61
  int streamID1 = pstreamOpenRead(cdoStreamName(0));
  int streamID2 = pstreamOpenRead(cdoStreamName(1));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
62

63
64
  int vlistID1 = pstreamInqVlist(streamID1);
  int vlistID2 = pstreamInqVlist(streamID2);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
65

66
  vlistCompare(vlistID1, vlistID2, CMP_ALL);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
67

68
  int gridsize = vlistGridsizeMax(vlistID1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
69

70
71
  double *array1 = (double*) Malloc(gridsize*sizeof(double));
  double *array2 = (double*) Malloc(gridsize*sizeof(double));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
72

73
74
75
  int indg = 0;
  int tsID = 0;
  int taxisID = vlistInqTaxis(vlistID1);
76
  while ( TRUE )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
77
    {
78
      nrecs = pstreamInqTimestep(streamID1, tsID);
79
80
      if ( nrecs > 0 )
	{
81
82
	  date2str(taxisInqVdate(taxisID), vdatestr, sizeof(vdatestr));
	  time2str(taxisInqVtime(taxisID), vtimestr, sizeof(vtimestr));
83
	}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
84

85
      nrecs2 = pstreamInqTimestep(streamID2, tsID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
86

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

89
      for ( int recID = 0; recID < nrecs; recID++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
90
	{
91
92
	  pstreamInqRecord(streamID1, &varID1, &levelID);
	  pstreamInqRecord(streamID2, &varID2, &levelID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
93
94
95

	  indg += 1;

96
97
98
99
100
101
102
	  int param    = vlistInqVarParam(vlistID1, varID1);
	  int code     = vlistInqVarCode(vlistID1, varID1);
	  int gridID   = vlistInqVarGrid(vlistID1, varID1);
	  int zaxisID  = vlistInqVarZaxis(vlistID1, varID1);
	  int gridsize = gridInqSize(gridID);
	  double missval1 = vlistInqVarMissval(vlistID1, varID1);
	  double missval2 = vlistInqVarMissval(vlistID2, varID2);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
103

104
	  //checkrel = gridInqType(gridID) != GRID_SPECTRAL;
105
          bool checkrel = true;
106

Uwe Schulzweida's avatar
Uwe Schulzweida committed
107
	  cdiParamToString(param, paramstr, sizeof(paramstr));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
108

109
110
	  pstreamReadRecord(streamID1, array1, &nmiss1);
	  pstreamReadRecord(streamID2, array2, &nmiss2);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
111

112
113
114
115
116
117
	  int ndiff = 0;
	  bool dsgn = false;
          bool zero = false;
          double absm = 0.0;
	  double relm = 0.0;
          double absdiff;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
118

119
	  for ( int i = 0; i < gridsize; i++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
120
	    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
121
122
123
124
125
126
127
	      if ( (DBL_IS_NAN(array1[i]) && !DBL_IS_NAN(array2[i])) ||
		  (!DBL_IS_NAN(array1[i]) &&  DBL_IS_NAN(array2[i])) )
		{
		  ndiff++;
		  relm = 1.0;
		}
	      else if ( !DBL_IS_EQUAL(array1[i], missval1) && !DBL_IS_EQUAL(array2[i], missval2) )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
128
		{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
129
130
		  absdiff = fabs(array1[i] - array2[i]);
		  if ( absdiff > 0. ) ndiff++;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
131

Uwe Schulzweida's avatar
Uwe Schulzweida committed
132
133
134
		  absm = MAX(absm, absdiff);

		  if ( array1[i]*array2[i] < 0. )
135
		    dsgn = true;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
136
		  else if ( IS_EQUAL(array1[i]*array2[i], 0.) )
137
		    zero = true;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
138
		  else
Uwe Schulzweida's avatar
Uwe Schulzweida committed
139
		    relm = MAX(relm, absdiff / MAX(fabs(array1[i]), fabs(array2[i])));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
140
141
		}
	      else if ( (DBL_IS_EQUAL(array1[i], missval1) && !DBL_IS_EQUAL(array2[i], missval2)) || 
Uwe Schulzweida's avatar
Uwe Schulzweida committed
142
		       (!DBL_IS_EQUAL(array1[i], missval1) &&  DBL_IS_EQUAL(array2[i], missval2)) )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
143
		{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
144
		  ndiff++;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
145
146
147
148
		  relm = 1.0;
		}
	    }

149
	  if ( ! cdoSilentMode || cdoVerbose )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
150
	    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
151
	      if ( absm > abslim || (checkrel && relm >= rellim) || cdoVerbose )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
152
		{
153
154
		  if ( lhead )
		    {
155
		      lhead = false;
156

Uwe Schulzweida's avatar
Uwe Schulzweida committed
157
		      set_text_color(stdout, BRIGHT, BLACK);
158
		      fprintf(stdout, "               Date     Time   Level Gridsize    Miss ");
159
		      fprintf(stdout, "   Diff ");
160
161
		      fprintf(stdout, ": S Z  Max_Absdiff Max_Reldiff");

162
		      if ( operatorID == DIFFN )
163
			fprintf(stdout, " : Parameter name");
Uwe Schulzweida's avatar
Uwe Schulzweida committed
164
		      else if ( operatorID == DIFF || operatorID == DIFFP )
165
166
167
			fprintf(stdout, " : Parameter ID");
		      else if ( operatorID == DIFFC )
			fprintf(stdout, " : Code number");
Uwe Schulzweida's avatar
Uwe Schulzweida committed
168
		      reset_text_color(stdout);
169
170
171
172
173

		      fprintf(stdout, "\n");
		    }

		  if ( operatorID == DIFFN ) vlistInqVarName(vlistID1, varID1, varname);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
174
175
176
177
178
179
180
		  
		  set_text_color(stdout, BRIGHT, BLACK);
		  fprintf(stdout, "%6d ", indg);
		  reset_text_color(stdout);
		  set_text_color(stdout, RESET, BLACK);
		  fprintf(stdout, ":");
		  reset_text_color(stdout);
181
		
Uwe Schulzweida's avatar
Uwe Schulzweida committed
182
		  set_text_color(stdout, RESET, MAGENTA);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
183
		  fprintf(stdout, "%s %s ", vdatestr, vtimestr);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
184
185
		  reset_text_color(stdout);
		  set_text_color(stdout, RESET, GREEN);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
186
                  double level = cdoZaxisInqLevel(zaxisID, levelID);
187
		  fprintf(stdout, "%7g ", level);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
188
		  fprintf(stdout, "%8d %7d ", gridsize, MAX(nmiss1, nmiss2));
189
		  fprintf(stdout, "%7d ", ndiff);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
190
		  reset_text_color(stdout);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
191
		
Uwe Schulzweida's avatar
Uwe Schulzweida committed
192
193
194
195
		  set_text_color(stdout, RESET, BLACK);
		  fprintf(stdout, ":");
		  reset_text_color(stdout);
		  fprintf(stdout, " %c %c ", dsgn ? 'T' : 'F', zero ? 'T' : 'F');
Uwe Schulzweida's avatar
Uwe Schulzweida committed
196
		  set_text_color(stdout, RESET, BLUE);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
197
		  fprintf(stdout, "%#12.5g%#12.5g", absm, relm);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
198
199
200
		  set_text_color(stdout, RESET, BLACK);
		  fprintf(stdout, " : ");
		  reset_text_color(stdout);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
201

Uwe Schulzweida's avatar
Uwe Schulzweida committed
202
		  set_text_color(stdout, BRIGHT, GREEN);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
203
		  if ( operatorID == DIFFN )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
204
		    fprintf(stdout, "%-11s", varname);
205
		  else if ( operatorID == DIFF || operatorID == DIFFP )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
206
		    fprintf(stdout, "%-11s", paramstr);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
207
		  else if ( operatorID == DIFFC )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
208
209
		    fprintf(stdout, "%4d", code);
		  reset_text_color(stdout);
210
		      
Uwe Schulzweida's avatar
Uwe Schulzweida committed
211
212
213
		  fprintf(stdout, "\n");
		}
	    }
214

Uwe Schulzweida's avatar
Uwe Schulzweida committed
215
	  ngrec++;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
216
217
	  if ( absm > abslim  || (checkrel && relm >= rellim) ) ndrec++;
	  if ( absm > abslim2 || (checkrel && relm >= rellim) ) nd2rec++;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
218
219
220
221
	}
      tsID++;
    }

222
223
  if ( ndrec > 0 )
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
224
225
226
227
228
      set_text_color(stdout, BRIGHT, RED);
      fprintf(stdout, "  %d of %d records differ", ndrec, ngrec);
      reset_text_color(stdout);
      fprintf(stdout, "\n");

Uwe Schulzweida's avatar
Uwe Schulzweida committed
229
      if ( ndrec != nd2rec && abslim < abslim2 )
230
231
232
233
	fprintf(stdout, "  %d of %d records differ more than 0.001\n", nd2rec, ngrec);
      /*  fprintf(stdout, "  %d of %d records differ more then one thousandth\n", nprec, ngrec); */
    }

234
235
236
237
  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
238

239
240
  pstreamClose(streamID1);
  pstreamClose(streamID2);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
241

242
243
  if ( array1 ) Free(array1);
  if ( array2 ) Free(array2);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
244
245
246

  cdoFinish();

247
  return 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
248
}