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         Timstat3        varquot2test
12         Timstat3        meandiff2test
13 */
14 
15 #include <cdi.h>
16 
17 #include "process_int.h"
18 #include "cdo_vlist.h"
19 #include "param_conversion.h"
20 #include "statistic.h"
21 
22 #define NIN 2
23 #define NOUT 1
24 #define NFWORK 4
25 #define NIWORK 2
26 
27 void *
Timstat3(void * process)28 Timstat3(void *process)
29 {
30   CdoStreamID streamID[NIN];
31   int vlistID[NIN], vlistID2 = -1;
32   int64_t vdate = 0;
33   int vtime = 0;
34   int is;
35   Varray3D<int> iwork[NIWORK];
36   FieldVector2D fwork[NFWORK];
37   int reached_eof[NIN];
38   constexpr int n_in = NIN;
39 
40   cdo_initialize(process);
41 
42   // clang-format off
43   const auto VARQUOT2TEST  = cdo_operator_add("varquot2test",  0, 0, nullptr);
44   const auto MEANDIFF2TEST = cdo_operator_add("meandiff2test", 0, 0, nullptr);
45   // clang-format on
46 
47   const auto operatorID = cdo_operator_id();
48 
49   operator_input_arg("constant and risk (e.g. 0.05)");
50   operator_check_argc(2);
51   const auto rconst = parameter_to_double(cdo_operator_argv(0));
52   const auto risk = parameter_to_double(cdo_operator_argv(1));
53 
54   if (operatorID == VARQUOT2TEST)
55     {
56       if (rconst <= 0) cdo_abort("Constant must be positive!");
57       if (risk <= 0 || risk >= 1) cdo_abort("Risk must be greater than 0 and lower than 1!");
58     }
59 
60   for (is = 0; is < NIN; ++is)
61     {
62       streamID[is] = cdo_open_read(is);
63       vlistID[is] = cdo_stream_inq_vlist(streamID[is]);
64       if (is > 0)
65         {
66           vlistID2 = cdo_stream_inq_vlist(streamID[is]);
67           vlist_compare(vlistID[0], vlistID2, CMP_ALL);
68         }
69     }
70 
71   const auto vlistID3 = vlistDuplicate(vlistID[0]);
72 
73   const auto gridsizemax = vlistGridsizeMax(vlistID[0]);
74   const auto nvars = vlistNvars(vlistID[0]);
75 
76   const auto maxrecs = vlistNrecs(vlistID[0]);
77   std::vector<RecordInfo> recList(maxrecs);
78 
79   const auto taxisID1 = vlistInqTaxis(vlistID[0]);
80   const auto taxisID3 = taxisDuplicate(taxisID1);
81 
82   vlistDefTaxis(vlistID3, taxisID3);
83   const auto streamID3 = cdo_open_write(2);
84   cdo_def_vlist(streamID3, vlistID3);
85 
86   for (int i = 0; i < NIN; ++i) reached_eof[i] = 0;
87 
88   Field in[NIN], out[NOUT];
89   for (int i = 0; i < NIN; ++i) in[i].resize(gridsizemax);
90   for (int i = 0; i < NOUT; ++i) out[i].resize(gridsizemax);
91 
92   for (int iw = 0; iw < NFWORK; ++iw) fwork[iw].resize(nvars);
93   for (int iw = 0; iw < NIWORK; ++iw) iwork[iw].resize(nvars);
94 
95   for (int varID = 0; varID < nvars; ++varID)
96     {
97       const auto gridID = vlistInqVarGrid(vlistID[0], varID);
98       const auto gridsize = vlistGridsizeMax(vlistID[0]);
99       const auto nlevels = zaxisInqSize(vlistInqVarZaxis(vlistID[0], varID));
100       const auto missval = vlistInqVarMissval(vlistID[0], varID);
101 
102       for (int iw = 0; iw < NFWORK; ++iw) fwork[iw][varID].resize(nlevels);
103       for (int iw = 0; iw < NIWORK; ++iw) iwork[iw][varID].resize(nlevels);
104 
105       for (int levelID = 0; levelID < nlevels; ++levelID)
106         {
107           for (int iw = 0; iw < NFWORK; ++iw)
108             {
109               fwork[iw][varID][levelID].grid = gridID;
110               fwork[iw][varID][levelID].missval = missval;
111               fwork[iw][varID][levelID].resize(gridsize);
112             }
113 
114           for (int iw = 0; iw < NIWORK; ++iw) iwork[iw][varID][levelID].resize(gridsize, 0);
115         }
116     }
117 
118   int tsID = 0;
119   while (true)
120     {
121       for (is = 0; is < NIN; ++is)
122         {
123           if (reached_eof[is]) continue;
124 
125           const auto nrecs = cdo_stream_inq_timestep(streamID[is], tsID);
126           if (nrecs == 0)
127             {
128               reached_eof[is] = 1;
129               continue;
130             }
131 
132           vdate = taxisInqVdate(taxisID1);
133           vtime = taxisInqVtime(taxisID1);
134 
135           for (int recID = 0; recID < nrecs; recID++)
136             {
137               int varID, levelID;
138               cdo_inq_record(streamID[is], &varID, &levelID);
139 
140               const auto gridsize = gridInqSize(vlistInqVarGrid(vlistID[is], varID));
141 
142               in[is].missval = vlistInqVarMissval(vlistID[is], varID);
143 
144               if (tsID == 0 && is == 0)
145                 {
146                   recList[recID].varID = varID;
147                   recList[recID].levelID = levelID;
148                   recList[recID].lconst = (vlistInqVarTimetype(vlistID[0], varID) == TIME_CONSTANT);
149                 }
150 
151               cdo_read_record(streamID[is], in[is].vec_d.data(), &in[is].nmiss);
152 
153               for (size_t i = 0; i < gridsize; ++i)
154                 {
155                   // if ( ( ! DBL_IS_EQUAL(array1[i], missval1) ) && ( ! DBL_IS_EQUAL(array2[i], missval2) ) )
156                   {
157                     fwork[NIN * is + 0][varID][levelID].vec_d[i] += in[is].vec_d[i];
158                     fwork[NIN * is + 1][varID][levelID].vec_d[i] += in[is].vec_d[i] * in[is].vec_d[i];
159                     iwork[is][varID][levelID][i]++;
160                   }
161                 }
162             }
163         }
164 
165       for (is = 0; is < NIN; ++is)
166         if (!reached_eof[is]) break;
167 
168       if (is == NIN) break;
169 
170       tsID++;
171     }
172 
173   taxisDefVdate(taxisID3, vdate);
174   taxisDefVtime(taxisID3, vtime);
175   cdo_def_timestep(streamID3, 0);
176 
177   for (int recID = 0; recID < maxrecs; recID++)
178     {
179       const auto varID = recList[recID].varID;
180       const auto levelID = recList[recID].levelID;
181 
182       const auto missval1 = fwork[0][varID][levelID].missval;
183       const auto missval2 = missval1;
184 
185       if (operatorID == VARQUOT2TEST)
186         {
187           const auto fwork0 = fwork[0][varID][levelID];
188           const auto fwork1 = fwork[1][varID][levelID];
189           const auto fwork2 = fwork[2][varID][levelID];
190           const auto fwork3 = fwork[3][varID][levelID];
191           for (size_t i = 0; i < gridsizemax; ++i)
192             {
193               const auto fnvals0 = iwork[0][varID][levelID][i];
194               const auto fnvals1 = iwork[1][varID][levelID][i];
195 
196               const auto temp0 = DIVMN(MULMN(fwork0.vec_d[i], fwork0.vec_d[i]), fnvals0);
197               const auto temp1 = DIVMN(MULMN(fwork2.vec_d[i], fwork2.vec_d[i]), fnvals1);
198               const auto temp2 = SUBMN(fwork1.vec_d[i], temp0);
199               const auto temp3 = SUBMN(fwork3.vec_d[i], temp1);
200               const auto statistic = DIVMN(temp2, ADDMN(temp2, MULMN(rconst, temp3)));
201 
202               double fractil_1 = missval1, fractil_2 = missval1;
203               if (fnvals0 > 1 && fnvals1 > 1)
204                 cdo::beta_distr_constants((fnvals0 - 1) / 2, (fnvals1 - 1) / 2, 1 - risk, &fractil_1, &fractil_2, __func__);
205 
206               out[0].vec_d[i]
207                   = DBL_IS_EQUAL(statistic, missval1) ? missval1 : (statistic <= fractil_1 || statistic >= fractil_2) ? 1.0 : 0.0;
208             }
209         }
210       else if (operatorID == MEANDIFF2TEST)
211         {
212           double mean_factor[NIN], var_factor[NIN];
213 
214           mean_factor[0] = 1.0;
215           mean_factor[1] = -1.0;
216           var_factor[0] = var_factor[1] = 1.0;
217 
218           for (size_t i = 0; i < gridsizemax; ++i)
219             {
220               double temp0 = 0.0;
221               double deg_of_freedom = -n_in;
222               for (int j = 0; j < n_in; j++)
223                 {
224                   const auto fnvals = iwork[j][varID][levelID][i];
225                   auto tmp = DIVMN(MULMN(fwork[2 * j][varID][levelID].vec_d[i], fwork[2 * j][varID][levelID].vec_d[i]), fnvals);
226                   temp0 = ADDMN(temp0, DIVMN(SUBMN(fwork[2 * j + 1][varID][levelID].vec_d[i], tmp), var_factor[j]));
227                   deg_of_freedom = ADDMN(deg_of_freedom, fnvals);
228                 }
229 
230               if (!DBL_IS_EQUAL(temp0, missval1) && temp0 < 0) temp0 = 0;  // This is possible because of rounding errors
231 
232               const auto stddev_estimator = SQRTMN(DIVMN(temp0, deg_of_freedom));
233               auto mean_estimator = -rconst;
234               for (int j = 0; j < n_in; j++)
235                 {
236                   const auto fnvals = iwork[j][varID][levelID][i];
237                   mean_estimator
238                       = ADDMN(mean_estimator, MULMN(mean_factor[j], DIVMN(fwork[2 * j][varID][levelID].vec_d[i], fnvals)));
239                 }
240 
241               double temp1 = 0.0;
242               for (int j = 0; j < n_in; j++)
243                 {
244                   const auto fnvals = iwork[j][varID][levelID][i];
245                   temp1 = ADDMN(temp1, DIVMN(MUL(MUL(mean_factor[j], mean_factor[j]), var_factor[j]), fnvals));
246                 }
247 
248               const auto norm = SQRTMN(temp1);
249 
250               const auto temp2 = DIVMN(DIVMN(mean_estimator, norm), stddev_estimator);
251               const auto fractil = (deg_of_freedom < 1) ? missval1 : cdo::student_t_inv(deg_of_freedom, 1 - risk / 2, __func__);
252 
253               out[0].vec_d[i]
254                   = DBL_IS_EQUAL(temp2, missval1) || DBL_IS_EQUAL(fractil, missval1) ? missval1 : std::fabs(temp2) >= fractil;
255             }
256         }
257 
258       cdo_def_record(streamID3, varID, levelID);
259       cdo_write_record(streamID3, out[0].vec_d.data(), field_num_miss(out[0]));
260     }
261 
262   cdo_stream_close(streamID3);
263   for (is = 0; is < NIN; ++is) cdo_stream_close(streamID[is]);
264 
265   cdo_finish();
266 
267   return nullptr;
268 }
269