ARTS 2.5.11 (git: 725533f0)
raw.cc
Go to the documentation of this file.
1
9#include <algorithm>
10
11#include "arts_conversions.h"
12#include "raw.h"
13
14template <typename T>
15constexpr bool isnormal_or_zero(T x) noexcept {return std::isnormal(x) or x == 0;}
16
17namespace Raw {
18namespace Calibration {
19Vector calibration(const Vector& pc, const Vector& pa, const Vector& ph, Numeric tc, Numeric th) noexcept
20{
21 const auto N = Index(pc.size());
22 Vector calib(N);
23 for (Index i=0; i<N; i++) calib[i] = calibration(pc[i], pa[i], ph[i], tc, th);
24 return calib;
25}
26
27Vector systemtemp(const Vector& pc, const Vector& ph, Numeric tc, Numeric th) noexcept
28{
29 const auto N = Index(pc.size());
30 Vector systemp(N);
31 for (Index i=0; i<N; i++) systemp[i] = systemtemp(Numeric(pc[i]), Numeric(ph[i]), tc, th);
32 return systemp;
33}
34
35ArrayOfVector caha(const ArrayOfVector& data, const Vector& tcvec, const Vector& thvec, const Index first_c_index)
36{
37 assert(first_c_index >= 0);
38
39 /* Must compute the positions due to c-Index offset
40 *
41 * The start position should be on C or H
42 * if first_c_index > 1 then start at H two steps early
43 *
44 * The cold position is always a multiple of 4 from
45 * first_c_index. In case first_c_index is > 3, we
46 * must take the modulus of the position
47 *
48 * The hot position is either two indices before or
49 * two indices behind the cold position, whichever
50 * is lowest and positive and below 4 */
51 const Index start = first_c_index - ((first_c_index > 1) ? 2 : 0);
52 const Index cpos = first_c_index % 4;
53 const Index hpos = cpos + ((cpos < 2) ? 2 : -2);
54
55 // Initialize pointers and values of interest
56 const Vector * pc = nullptr;
57 const Vector * pa = nullptr;
58 const Vector * ph = nullptr;
59 Numeric tc=0, th=0;
60
61 // Output vector
62 ArrayOfVector out;
63
64 // For all data entries
65 for (Index i=start; i<data.nelem(); i++) {
66 // Cycle is CAHA → 4 steps
67 const Index pos = i % 4;
68
69 // Set the values when we are at the correct position
70 if (pos == cpos) {
71 pc = &data[i];
72 tc = tcvec[i];
73 } else if (pos == hpos) {
74 ph = &data[i];
75 th = thvec[i];
76 } else {
77 pa = &data[i];
78 }
79
80 /* Add new measurement if we are at C or H since we have
81 * then completed a new calibration cycle
82 *
83 * Note that nullptr is false */
84 if ((pos == cpos or pos == hpos) and pc and pa and ph) {
85 out.emplace_back(calibration(*pc, *pa, *ph, tc, th));
86 }
87 }
88
89 return out;
90}
91} // namespace Calibration
92
93namespace Average {
94VectorView avg(VectorView y, const ArrayOfVector& ys, const Index start, const Index end_tmp)
95{
96 // True end
97 const Index end = end_tmp >= 0 ? end_tmp : 1 + ys.nelem() + end_tmp;
98
99 // Scale
100 const Numeric scale = 1 / Numeric(end - start);
101
102 // Reset
103 y = 0;
104
105 // Compute
106 for (Index k=start; k<end; k++)
107 for (Index i=0; i<y.nelem(); i++)
108 y[i] += ys[k][i] * scale;
109 return y;
110}
111
112VectorView nanavg(VectorView y, const ArrayOfVector& ys, const Index start, const Index end_tmp)
113{
114 // True end
115 const Index end = end_tmp >= 0 ? end_tmp : 1 + ys.nelem() + end_tmp;
116
117 // Reset y
118 y = 0;
119
120 // Compute the averages ignoring all NaN
121 for (Index i=0; i<y.nelem(); i++) {
122 Index numnormal=0;
123 for (Index k=start; k<end; k++) {
124 if (isnormal_or_zero(ys[k][i])) {
125 y[i] += ys[k][i];
126 numnormal++;
127 }
128 }
129
130 // Compute the average or set to NaN
131 if (numnormal) {
132 y[i] /= Numeric(numnormal);
133 } else {
134 y[i] = std::numeric_limits<Numeric>::quiet_NaN();
135 }
136 }
137
138 return y;
139}
140
141VectorView var(VectorView var, const Vector& y, const ArrayOfVector& ys, const Index start, const Index end_tmp)
142{
143 // True end
144 const Index end = end_tmp >= 0 ? end_tmp : 1 + ys.nelem() + end_tmp;
145
146 // Scale
147 const Numeric scale = 1 / Numeric(end - start);
148
149 // Reset
150 var = 0;
151
152 // Compute
153 for (Index k=start; k<end; k++)
154 for (Index i=0; i<y.nelem(); i++)
155 var[i] += Math::pow2(ys[k][i] - y[i]) * scale;
156 return var;
157}
158
159VectorView nanvar(VectorView var, const Vector& y, const ArrayOfVector& ys, const Index start, const Index end_tmp)
160{
161 // True end
162 const Index end = end_tmp >= 0 ? end_tmp : 1 + ys.nelem() + end_tmp;
163
164 // Reset
165 var = 0;
166
167 // Compute
168 for (Index i=0; i<y.nelem(); i++) {
169 Index numnormal=0;
170 for (Index k=start; k<end; k++) {
171 if (isnormal_or_zero(ys[k][i])) {
172 var[i] += Math::pow2(ys[k][i] - y[i]);
173 numnormal++;
174 }
175 }
176
177 // Compute the average or set to NaN
178 if (numnormal > 0) {
179 var[i] /= Numeric(numnormal);
180 } else {
181 var[i] = std::numeric_limits<Numeric>::quiet_NaN();
182 }
183 }
184 return var;
185}
186
187VectorView std(VectorView std, const Vector& y, const ArrayOfVector& ys, const Index start, const Index end_tmp)
188{
189 var(std, y, ys, start, end_tmp);
190 std::transform(std.begin(), std.end(), std.begin(), [](const auto& x){return std::sqrt(x);});
191 return std;
192}
193
194VectorView nanstd(VectorView std, const Vector& y, const ArrayOfVector& ys, const Index start, const Index end_tmp)
195{
196 nanvar(std, y, ys, start, end_tmp);
197 std::transform(std.begin(), std.end(), std.begin(), [](const auto& x){return std::sqrt(x);});
198 return std;
199}
200
201MatrixView cov(MatrixView cov, const Vector& y, const ArrayOfVector& ys, const Index start, const Index end_tmp)
202{
203 // True end
204 const Index end = end_tmp >= 0 ? end_tmp : 1 + ys.nelem() + end_tmp;
205
206 // Scale
207 const Numeric scale = 1 / Numeric(end - start - 1);
208
209 // Reset
210 cov = 0;
211
212 // Compute
213 for (Index k=start; k<end; k++)
214 for (Index i=0; i<y.nelem(); i++)
215 for (Index j=0; j<y.nelem(); j++)
216 cov(i, j) += (ys[k][i] - y[i]) * (ys[k][j] - y[j]) * scale;
217 return cov;
218}
219
220Numeric median(const ConstVectorView& v, const ArrayOfIndex& pos)
221{
222 // Size of problem
223 const Index n = pos.nelem() ? pos.nelem() : v.nelem();
224
225 // Create a vector to sort
226 ArrayOfNumeric calc(n);
227 for (Index i=0; i<n; i++)
228 if (pos.nelem())
229 calc[i] = v[pos[i]];
230 else
231 calc[i] = v[i];
232
233 // Sort the vector
234 std::sort(calc.begin(), calc.end());
235
236 // Return the median
237 if (n % 2)
238 return calc[n/2];
239 return (calc[(n-1)/2] + calc[n/2]) / 2;
240}
241
242Numeric nanmedian(const ConstVectorView& v, const ArrayOfIndex& pos)
243{
244 // Size of problem
245 const Index n = pos.nelem() ? pos.nelem() : v.nelem();
246
247 // Create a vector to sort
248 ArrayOfNumeric calc(n);
249 for (Index i=0; i<n; i++) {
250 if (pos.nelem()) {
251 calc[i] = v[pos[i]];
252 } else {
253 calc[i] = v[i];
254 }
255 }
256
257 // Sort the vector
258 std::sort(calc.begin(), calc.end());
259
260 // Remove non-normal
261 auto newend = std::remove_if(calc.begin(), calc.end(), [](auto x){return not isnormal_or_zero(x);});
262 const Index numnormal = n - Index(calc.end() - newend);
263
264 // Return the median
265 if (numnormal) {
266 if (numnormal % 2)
267 return calc[numnormal/2];
268 return (calc[(numnormal-1)/2] + calc[numnormal/2]) / 2;
269 }
270
271 return std::numeric_limits<Numeric>::quiet_NaN();
272}
273} // namespace Average
274
275namespace Reduce {
276ArrayOfIndex focus_doublescaling(const Vector& x, const Numeric x0, const Numeric dx)
277{
278 ArrayOfIndex scale(x.size());
279
280 // Scaling function, rounds down an Index
281 auto scaling = [](Numeric X, Numeric X0, Numeric DX) {
282 const Index p = 1 + Index(std::abs(X - X0) / DX); // floor + 1 is close enough
283 Index this_scale = 0;
284 while ((p >> this_scale) > 1) this_scale++;
285 return 1 << this_scale;
286 };
287
288 for (Index i=0; i<x.size(); i++) scale[i] = scaling(x[i], x0, dx);
289 return scale;
290}
291
292std::pair<Index, Index> find_first_and_last_1(const ArrayOfIndex& x)
293{
294 Index first=x.nelem() - 1;
295 Index last=0;
296 for (Index i=0; i<x.nelem(); i++) {
297 if (x[i] == 1) {
298 last = std::max(i, last);
299 first = std::min(i, first);
300 }
301 }
302
303 assert(first <= last);
304 return {first, last};
305}
306
307Vector focus(const Vector& x, const ArrayOfIndex& scaling)
308{
309 const Index N = x.size();
310 Index first_i, last_i;
311
312 // Find the first and last one in scaling
313 const auto [firstone, lastone] = find_first_and_last_1(scaling);
314
315 // All the central numbers
316 std::vector<Numeric> out;
317 for (Index i=firstone; i<=lastone; i++) out.emplace_back(x[i]);
318
319 // All the last numbers appended
320 first_i = lastone + 1;
321 while (first_i not_eq N) {
322 last_i = first_i + scaling[first_i];
323 if (last_i > N) last_i = N;
324 out.emplace_back(mean(x[Range(first_i, last_i-first_i)]));
325 first_i = last_i;
326 }
327
328 // All the first numbers prepended
329 std::vector<Numeric> fx(0);
330 last_i = firstone;
331 while (last_i > 0) {
332 first_i = last_i - scaling[last_i];
333 if (first_i < 0) first_i = 0;
334 out.insert(out.begin(), mean(x[Range(first_i, last_i-first_i)]));
335 last_i = first_i;
336 }
337
338 return out;
339}
340
341Vector nanfocus(const Vector& x, const ArrayOfIndex& scaling)
342{
343 const Index N = x.size();
344 Index first_i, last_i;
345
346 // Find the first and last one in scaling
347 const auto [firstone, lastone] = find_first_and_last_1(scaling);
348
349 // All the central numbers
350 std::vector<Numeric> out;
351 for (Index i=firstone; i<=lastone; i++) {
352 if (isnormal_or_zero(x[i])) {
353 out.emplace_back(x[i]);
354 } else {
355 out.emplace_back(std::numeric_limits<Numeric>::quiet_NaN());
356 }
357 }
358
359 // All the last numbers appended
360 first_i = lastone + 1;
361 while (first_i not_eq N) {
362 last_i = first_i + scaling[first_i];
363 if (last_i > N) last_i = N;
364 out.emplace_back(nanmean(x[Range(first_i, last_i-first_i)]));
365 first_i = last_i;
366 }
367
368 // All the first numbers prepended
369 std::vector<Numeric> fx(0);
370 last_i = firstone;
371 while (last_i > 0) {
372 first_i = last_i - scaling[last_i];
373 if (first_i < 0) first_i = 0;
374 out.insert(out.begin(), nanmean(x[Range(first_i, last_i-first_i)]));
375 last_i = first_i;
376 }
377
378 return out;
379}
380} // namespace Reduce
381
383std::vector<bool> out_of_bounds(const Vector& x, const Numeric xmin, const Numeric xmax)
384{
386 for (Index i=0; i<x.size(); i++) mask[i] = (xmin > x[i]) or (xmax < x[i]) or (not isnormal_or_zero(x[i]));
388}
389
391{
392 const Index N = x.nelem();
393 for (Index i=0; i<N; i++)
395 x[i] = std::numeric_limits<Numeric>::quiet_NaN();
396 return x;
397}
399
400namespace Correction {
401 Numeric naive_tropospheric_singletau_median(const ConstVectorView& bt, const Numeric trop_bt, const Numeric target_bt) {
402 return - std::log((trop_bt - Average::nanmedian(bt)) / (trop_bt - target_bt));
403}
404
405VectorView naive_tropospheric(VectorView bt, const Numeric tau, const Numeric trop_bt) {
406 const Numeric t = std::exp(-tau);
407 if (std::isnormal(t)) {
408 bt /= t;
409 bt += trop_bt * std::expm1(-tau) / t;
410 }
411 return bt;
412}
413} // namespace Correction
414} // namespace Raw
Common ARTS conversions.
TimeStep mean(const ArrayOfTimeStep &dt)
Returns the mean time step.
Definition artstime.cc:172
Index nelem() const ARTS_NOEXCEPT
Definition array.h:75
void Reduce(Numeric &o, const Vector &i, const Verbosity &)
WORKSPACE METHOD: Reduce.
Definition m_reduce.h:115
Numeric last(ConstVectorView x)
last
constexpr auto pow2(auto x) noexcept
power of two
VectorView avg(VectorView y, const ArrayOfVector &ys, const Index start, const Index end_tmp)
Compute the average of the ranged ys.
Definition raw.cc:94
MatrixView cov(MatrixView cov, const Vector &y, const ArrayOfVector &ys, const Index start, const Index end_tmp)
Compute the covariance matrix of the ranged ys.
Definition raw.cc:201
VectorView nanvar(VectorView var, const Vector &y, const ArrayOfVector &ys, const Index start, const Index end_tmp)
Compute the variance of the ranged ys ignoring all non-normal values.
Definition raw.cc:159
VectorView var(VectorView var, const Vector &y, const ArrayOfVector &ys, const Index start, const Index end_tmp)
Compute the variance of the ranged ys.
Definition raw.cc:141
Numeric nanmedian(const ConstVectorView &v, const ArrayOfIndex &pos)
Get the median of the vector in the range ignoring all non-normal values.
Definition raw.cc:242
VectorView nanavg(VectorView y, const ArrayOfVector &ys, const Index start, const Index end_tmp)
Compute the average of the ranged ys ignoring all non-normal values.
Definition raw.cc:112
VectorView nanstd(VectorView std, const Vector &y, const ArrayOfVector &ys, const Index start, const Index end_tmp)
Compute the standard deviation of the ranged ys ignoring all non-normal values.
Definition raw.cc:194
Numeric median(const ConstVectorView &v, const ArrayOfIndex &pos)
Get the median of the vector in the range.
Definition raw.cc:220
ArrayOfVector caha(const ArrayOfVector &data, const Vector &tcvec, const Vector &thvec, const Index first_c_index)
Calibrate the data by CAHA.
Definition raw.cc:35
Vector systemtemp(const Vector &pc, const Vector &ph, Numeric tc, Numeric th) noexcept
Computes the linear receiver temperature formula elementwise.
Definition raw.cc:27
Vector calibration(const Vector &pc, const Vector &pa, const Vector &ph, Numeric tc, Numeric th) noexcept
Computes the linear calibration formula elementwise.
Definition raw.cc:19
VectorView naive_tropospheric(VectorView bt, const Numeric tau, const Numeric trop_bt)
Apply tropospheric correction.
Definition raw.cc:405
Numeric naive_tropospheric_singletau_median(const ConstVectorView &bt, const Numeric trop_bt, const Numeric target_bt)
Naive tropospheric correction parameterization.
Definition raw.cc:401
std::vector< bool > out_of_bounds(const Vector &x, const Numeric xmin, const Numeric xmax)
Masks values that are out of bounds in x.
Definition raw.cc:383
Masks all true entries of masking in x by turning them into NaN.
Definition raw.cc:390
Vector nanfocus(const Vector &x, const ArrayOfIndex &scaling)
Returns the re-averaged values of x according to scaling ignoring all non-normal values.
Definition raw.cc:341
std::pair< Index, Index > find_first_and_last_1(const ArrayOfIndex &x)
Definition raw.cc:292
ArrayOfIndex focus_doublescaling(const Vector &x, const Numeric x0, const Numeric dx)
Returns the relative position scale for each value in x.
Definition raw.cc:276
Vector focus(const Vector &x, const ArrayOfIndex &scaling)
Returns the re-averaged values of x according to scaling.
Definition raw.cc:307
Definition raw.cc:17
constexpr bool isnormal_or_zero(T x) noexcept
Definition raw.cc:15
Stuff related to generating y-data from raw data.
#define v