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 Regres regres Regression
12 */
13
14 #include <cdi.h>
15
16 #include "process_int.h"
17 #include "cdo_vlist.h"
18 #include "cdo_options.h"
19 #include "datetime.h"
20 #include "pmlist.h"
21 #include "param_conversion.h"
22
23 // Same code as Trend!
24
25 static void
regresGetParameter(bool & tstepIsEqual)26 regresGetParameter(bool &tstepIsEqual)
27 {
28 const auto pargc = cdo_operator_argc();
29 if (pargc)
30 {
31 auto pargv = cdo_get_oper_argv();
32
33 KVList kvlist;
34 kvlist.name = "TREND";
35 if (kvlist.parse_arguments(pargc, pargv) != 0) cdo_abort("Parse error!");
36 if (Options::cdoVerbose) kvlist.print();
37
38 for (const auto &kv : kvlist)
39 {
40 const auto &key = kv.key;
41 if (kv.nvalues > 1) cdo_abort("Too many values for parameter key >%s<!", key.c_str());
42 if (kv.nvalues < 1) cdo_abort("Missing value for parameter key >%s<!", key.c_str());
43 const auto &value = kv.values[0];
44
45 // clang-format off
46 if (key == "equal") tstepIsEqual = parameter_to_bool(value);
47 else cdo_abort("Invalid parameter key >%s<!", key.c_str());
48 // clang-format on
49 }
50 }
51 }
52
53 void *
Regres(void * process)54 Regres(void *process)
55 {
56 cdo_initialize(process);
57
58 auto tstepIsEqual = true;
59 regresGetParameter(tstepIsEqual);
60
61 const auto streamID1 = cdo_open_read(0);
62
63 const auto vlistID1 = cdo_stream_inq_vlist(streamID1);
64 const auto vlistID2 = vlistDuplicate(vlistID1);
65
66 vlistDefNtsteps(vlistID2, 1);
67
68 const auto taxisID1 = vlistInqTaxis(vlistID1);
69 const auto taxisID2 = taxisDuplicate(taxisID1);
70 vlistDefTaxis(vlistID2, taxisID2);
71
72 // const auto streamID2 = cdo_open_write(1);
73 // cdo_def_vlist(streamID2, vlistID2);
74 const auto streamID3 = cdo_open_write(1);
75 cdo_def_vlist(streamID3, vlistID2);
76
77 VarList varList;
78 varListInit(varList, vlistID1);
79
80 const auto maxrecs = vlistNrecs(vlistID1);
81 std::vector<RecordInfo> recList(maxrecs);
82
83 const auto gridsizemax = vlistGridsizeMax(vlistID1);
84
85 Field field1, field2;
86 field1.resize(gridsizemax);
87 field2.resize(gridsizemax);
88
89 constexpr size_t nwork = 5;
90 FieldVector2D work[nwork];
91 for (auto &w : work) fields_from_vlist(vlistID1, w, FIELD_VEC, 0);
92
93 const auto calendar = taxisInqCalendar(taxisID1);
94
95 CheckTimeInc checkTimeInc;
96 JulianDate juldate0;
97 double deltat1 = 0;
98 int64_t vdate = 0;
99 int vtime = 0;
100 int tsID = 0;
101 while (true)
102 {
103 const auto nrecs = cdo_stream_inq_timestep(streamID1, tsID);
104 if (nrecs == 0) break;
105
106 vdate = taxisInqVdate(taxisID1);
107 vtime = taxisInqVtime(taxisID1);
108
109 if (tstepIsEqual) check_time_increment(tsID, calendar, vdate, vtime, checkTimeInc);
110 const auto zj = tstepIsEqual ? (double) tsID : delta_time_step_0(tsID, calendar, vdate, vtime, juldate0, deltat1);
111
112 for (int recID = 0; recID < nrecs; recID++)
113 {
114 int varID, levelID;
115 cdo_inq_record(streamID1, &varID, &levelID);
116
117 if (tsID == 0)
118 {
119 recList[recID].varID = varID;
120 recList[recID].levelID = levelID;
121 recList[recID].lconst = (varList[varID].timetype == TIME_CONSTANT);
122 }
123
124 size_t nmiss;
125 cdo_read_record(streamID1, field1.vec_d.data(), &nmiss);
126
127 const auto gridsize = varList[varID].gridsize;
128 const auto missval = varList[varID].missval;
129
130 auto &sumj = work[0][varID][levelID].vec_d;
131 auto &sumjj = work[1][varID][levelID].vec_d;
132 auto &sumjx = work[2][varID][levelID].vec_d;
133 auto &sumx = work[3][varID][levelID].vec_d;
134 auto &zn = work[4][varID][levelID].vec_d;
135
136 for (size_t i = 0; i < gridsize; i++)
137 if (!DBL_IS_EQUAL(field1.vec_d[i], missval))
138 {
139 sumj[i] += zj;
140 sumjj[i] += zj * zj;
141 sumjx[i] += zj * field1.vec_d[i];
142 sumx[i] += field1.vec_d[i];
143 zn[i]++;
144 }
145 }
146
147 tsID++;
148 }
149
150 taxisDefVdate(taxisID2, vdate);
151 taxisDefVtime(taxisID2, vtime);
152 // cdo_def_timestep(streamID2, 0);
153 cdo_def_timestep(streamID3, 0);
154
155 for (int recID = 0; recID < maxrecs; recID++)
156 {
157 const auto varID = recList[recID].varID;
158 const auto levelID = recList[recID].levelID;
159 const auto gridsize = varList[varID].gridsize;
160 const auto missval = varList[varID].missval;
161 const auto missval1 = missval;
162 const auto missval2 = missval;
163 field1.size = gridsize;
164 field1.missval = missval;
165 field2.size = gridsize;
166 field2.missval = missval;
167
168 const auto &sumj = work[0][varID][levelID].vec_d;
169 const auto &sumjj = work[1][varID][levelID].vec_d;
170 const auto &sumjx = work[2][varID][levelID].vec_d;
171 const auto &sumx = work[3][varID][levelID].vec_d;
172 const auto &zn = work[4][varID][levelID].vec_d;
173
174 for (size_t i = 0; i < gridsize; i++)
175 {
176 const auto temp1 = SUBMN(sumjx[i], DIVMN(MULMN(sumj[i], sumx[i]), zn[i]));
177 const auto temp2 = SUBMN(sumjj[i], DIVMN(MULMN(sumj[i], sumj[i]), zn[i]));
178
179 field2.vec_d[i] = DIVMN(temp1, temp2);
180 field1.vec_d[i] = SUBMN(DIVMN(sumx[i], zn[i]), MULMN(DIVMN(sumj[i], zn[i]), field2.vec_d[i]));
181 }
182
183 cdo_def_record(streamID3, varID, levelID);
184 cdo_write_record(streamID3, field2.vec_d.data(), field_num_miss(field2));
185 }
186
187 cdo_stream_close(streamID3);
188 cdo_stream_close(streamID1);
189
190 cdo_finish();
191
192 return nullptr;
193 }
194