ARTS 2.5.0 (git: 9ee3ac6c)
raw.cc
Go to the documentation of this file.
1/* Copyright (C) 2020
2 * Richard Larsson <ric.larsson@gmail.com>
3 *
4 * This program is free software; you can redistribute it and/or modify it
5 * under the terms of the GNU General Public License as published by the
6 * Free Software Foundation; either version 2, or (at your option) any
7 * later version.
8 *
9 * This program is distributed in the hope that it will be useful,
10 * but WITHOUT ANY WARRANTY; without even the implied warranty of
11 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 * GNU General Public License for more details.
13 *
14 * You should have received a copy of the GNU General Public License
15 * along with this program; if not, write to the Free Software
16 * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
17 * USA. */
18
27#include <algorithm>
28
29#include "constants.h"
30#include "raw.h"
31
32template <typename T>
33constexpr bool isnormal_or_zero(T x) noexcept {return std::isnormal(x) or x == 0;}
34
35namespace Raw {
36namespace Calibration {
37Vector calibration(const Vector& pc, const Vector& pa, const Vector& ph, Numeric tc, Numeric th) noexcept
38{
39 const auto N = Index(pc.size());
40 Vector calib(N);
41 for (Index i=0; i<N; i++) calib[i] = calibration(pc[i], pa[i], ph[i], tc, th);
42 return calib;
43}
44
45Vector systemtemp(const Vector& pc, const Vector& ph, Numeric tc, Numeric th) noexcept
46{
47 const auto N = Index(pc.size());
48 Vector systemp(N);
49 for (Index i=0; i<N; i++) systemp[i] = systemtemp(Numeric(pc[i]), Numeric(ph[i]), tc, th);
50 return systemp;
51}
52
53ArrayOfVector caha(const ArrayOfVector& data, const Vector& tcvec, const Vector& thvec, const Index first_c_index)
54{
55 assert(first_c_index >= 0);
56
57 /* Must compute the positions due to c-Index offset
58 *
59 * The start position should be on C or H
60 * if first_c_index > 1 then start at H two steps early
61 *
62 * The cold position is always a multiple of 4 from
63 * first_c_index. In case first_c_index is > 3, we
64 * must take the modulus of the position
65 *
66 * The hot position is either two indices before or
67 * two indices behind the cold position, whichever
68 * is lowest and positive and below 4 */
69 const Index start = first_c_index - ((first_c_index > 1) ? 2 : 0);
70 const Index cpos = first_c_index % 4;
71 const Index hpos = cpos + ((cpos < 2) ? 2 : -2);
72
73 // Initialize pointers and values of interest
74 const Vector * pc = nullptr;
75 const Vector * pa = nullptr;
76 const Vector * ph = nullptr;
77 Numeric tc=0, th=0;
78
79 // Output vector
80 ArrayOfVector out;
81
82 // For all data entries
83 for (Index i=start; i<data.nelem(); i++) {
84 // Cycle is CAHA → 4 steps
85 const Index pos = i % 4;
86
87 // Set the values when we are at the correct position
88 if (pos == cpos) {
89 pc = &data[i];
90 tc = tcvec[i];
91 } else if (pos == hpos) {
92 ph = &data[i];
93 th = thvec[i];
94 } else {
95 pa = &data[i];
96 }
97
98 /* Add new measurement if we are at C or H since we have
99 * then completed a new calibration cycle
100 *
101 * Note that nullptr is false */
102 if ((pos == cpos or pos == hpos) and pc and pa and ph) {
103 out.emplace_back(calibration(*pc, *pa, *ph, tc, th));
104 }
105 }
106
107 return out;
108}
109} // namespace Calibration
110
111namespace Average {
112VectorView avg(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 // Scale
118 const Numeric scale = 1 / Numeric(end - start);
119
120 // Reset
121 y = 0;
122
123 // Compute
124 for (Index k=start; k<end; k++)
125 for (Index i=0; i<y.nelem(); i++)
126 y[i] += ys[k][i] * scale;
127 return y;
128}
129
130VectorView nanavg(VectorView y, const ArrayOfVector& ys, const Index start, const Index end_tmp)
131{
132 // True end
133 const Index end = end_tmp >= 0 ? end_tmp : 1 + ys.nelem() + end_tmp;
134
135 // Reset y
136 y = 0;
137
138 // Compute the averages ignoring all NaN
139 for (Index i=0; i<y.nelem(); i++) {
140 Index numnormal=0;
141 for (Index k=start; k<end; k++) {
142 if (isnormal_or_zero(ys[k][i])) {
143 y[i] += ys[k][i];
144 numnormal++;
145 }
146 }
147
148 // Compute the average or set to NaN
149 if (numnormal) {
150 y[i] /= Numeric(numnormal);
151 } else {
152 y[i] = std::numeric_limits<Numeric>::quiet_NaN();
153 }
154 }
155
156 return y;
157}
158
159VectorView var(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 // Scale
165 const Numeric scale = 1 / Numeric(end - start);
166
167 // Reset
168 var = 0;
169
170 // Compute
171 for (Index k=start; k<end; k++)
172 for (Index i=0; i<y.nelem(); i++)
173 var[i] += Constant::pow2(ys[k][i] - y[i]) * scale;
174 return var;
175}
176
177VectorView nanvar(VectorView var, const Vector& y, const ArrayOfVector& ys, const Index start, const Index end_tmp)
178{
179 // True end
180 const Index end = end_tmp >= 0 ? end_tmp : 1 + ys.nelem() + end_tmp;
181
182 // Reset
183 var = 0;
184
185 // Compute
186 for (Index i=0; i<y.nelem(); i++) {
187 Index numnormal=0;
188 for (Index k=start; k<end; k++) {
189 if (isnormal_or_zero(ys[k][i])) {
190 var[i] += Constant::pow2(ys[k][i] - y[i]);
191 numnormal++;
192 }
193 }
194
195 // Compute the average or set to NaN
196 if (numnormal > 0) {
197 var[i] /= Numeric(numnormal);
198 } else {
199 var[i] = std::numeric_limits<Numeric>::quiet_NaN();
200 }
201 }
202 return var;
203}
204
205VectorView std(VectorView std, const Vector& y, const ArrayOfVector& ys, const Index start, const Index end_tmp)
206{
207 var(std, y, ys, start, end_tmp);
208 std::transform(std.begin(), std.end(), std.begin(), [](const auto& x){return std::sqrt(x);});
209 return std;
210}
211
212VectorView nanstd(VectorView std, const Vector& y, const ArrayOfVector& ys, const Index start, const Index end_tmp)
213{
214 nanvar(std, y, ys, start, end_tmp);
215 std::transform(std.begin(), std.end(), std.begin(), [](const auto& x){return std::sqrt(x);});
216 return std;
217}
218
219MatrixView cov(MatrixView cov, const Vector& y, const ArrayOfVector& ys, const Index start, const Index end_tmp)
220{
221 // True end
222 const Index end = end_tmp >= 0 ? end_tmp : 1 + ys.nelem() + end_tmp;
223
224 // Scale
225 const Numeric scale = 1 / Numeric(end - start - 1);
226
227 // Reset
228 cov = 0;
229
230 // Compute
231 for (Index k=start; k<end; k++)
232 for (Index i=0; i<y.nelem(); i++)
233 for (Index j=0; j<y.nelem(); j++)
234 cov(i, j) += (ys[k][i] - y[i]) * (ys[k][j] - y[j]) * scale;
235 return cov;
236}
237
239{
240 // Size of problem
241 const Index n = pos.nelem() ? pos.nelem() : v.nelem();
242
243 // Create a vector to sort
244 ArrayOfNumeric calc(n);
245 for (Index i=0; i<n; i++)
246 if (pos.nelem())
247 calc[i] = v[pos[i]];
248 else
249 calc[i] = v[i];
250
251 // Sort the vector
252 std::sort(calc.begin(), calc.end());
253
254 // Return the median
255 if (n % 2)
256 return calc[n/2];
257 return (calc[(n-1)/2] + calc[n/2]) / 2;
258}
259
261{
262 // Size of problem
263 const Index n = pos.nelem() ? pos.nelem() : v.nelem();
264
265 // Create a vector to sort
266 ArrayOfNumeric calc(n);
267 for (Index i=0; i<n; i++) {
268 if (pos.nelem()) {
269 calc[i] = v[pos[i]];
270 } else {
271 calc[i] = v[i];
272 }
273 }
274
275 // Sort the vector
276 std::sort(calc.begin(), calc.end());
277
278 // Remove non-normal
279 auto newend = std::remove_if(calc.begin(), calc.end(), [](auto x){return not isnormal_or_zero(x);});
280 const Index numnormal = n - Index(calc.end() - newend);
281
282 // Return the median
283 if (numnormal) {
284 if (numnormal % 2)
285 return calc[numnormal/2];
286 return (calc[(numnormal-1)/2] + calc[numnormal/2]) / 2;
287 }
288
289 return std::numeric_limits<Numeric>::quiet_NaN();
290}
291} // namespace Average
292
293namespace Reduce {
295{
296 ArrayOfIndex scale(x.size());
297
298 // Scaling function, rounds down an Index
299 auto scaling = [](Numeric X, Numeric X0, Numeric DX) {
300 const Index p = 1 + Index(std::abs(X - X0) / DX); // floor + 1 is close enough
301 Index this_scale = 0;
302 while ((p >> this_scale) > 1) this_scale++;
303 return 1 << this_scale;
304 };
305
306 for (Index i=0; i<x.size(); i++) scale[i] = scaling(x[i], x0, dx);
307 return scale;
308}
309
310std::pair<Index, Index> find_first_and_last_1(const ArrayOfIndex& x)
311{
312 Index first=x.nelem() - 1;
313 Index last=0;
314 for (Index i=0; i<x.nelem(); i++) {
315 if (x[i] == 1) {
316 last = std::max(i, last);
317 first = std::min(i, first);
318 }
319 }
320
321 assert(first <= last);
322 return {first, last};
323}
324
325Vector focus(const Vector& x, const ArrayOfIndex& scaling)
326{
327 const Index N = x.size();
328 Index first_i, last_i;
329
330 // Find the first and last one in scaling
331 const auto [firstone, lastone] = find_first_and_last_1(scaling);
332
333 // All the central numbers
334 std::vector<Numeric> out;
335 for (Index i=firstone; i<=lastone; i++) out.emplace_back(x[i]);
336
337 // All the last numbers appended
338 first_i = lastone + 1;
339 while (first_i not_eq N) {
340 last_i = first_i + scaling[first_i];
341 if (last_i > N) last_i = N;
342 out.emplace_back(mean(x[Range(first_i, last_i-first_i)]));
343 first_i = last_i;
344 }
345
346 // All the first numbers prepended
347 std::vector<Numeric> fx(0);
348 last_i = firstone;
349 while (last_i > 0) {
350 first_i = last_i - scaling[last_i];
351 if (first_i < 0) first_i = 0;
352 out.insert(out.begin(), mean(x[Range(first_i, last_i-first_i)]));
353 last_i = first_i;
354 }
355
356 return out;
357}
358
359Vector nanfocus(const Vector& x, const ArrayOfIndex& scaling)
360{
361 const Index N = x.size();
362 Index first_i, last_i;
363
364 // Find the first and last one in scaling
365 const auto [firstone, lastone] = find_first_and_last_1(scaling);
366
367 // All the central numbers
368 std::vector<Numeric> out;
369 for (Index i=firstone; i<=lastone; i++) {
370 if (isnormal_or_zero(x[i])) {
371 out.emplace_back(x[i]);
372 } else {
373 out.emplace_back(std::numeric_limits<Numeric>::quiet_NaN());
374 }
375 }
376
377 // All the last numbers appended
378 first_i = lastone + 1;
379 while (first_i not_eq N) {
380 last_i = first_i + scaling[first_i];
381 if (last_i > N) last_i = N;
382 out.emplace_back(nanmean(x[Range(first_i, last_i-first_i)]));
383 first_i = last_i;
384 }
385
386 // All the first numbers prepended
387 std::vector<Numeric> fx(0);
388 last_i = firstone;
389 while (last_i > 0) {
390 first_i = last_i - scaling[last_i];
391 if (first_i < 0) first_i = 0;
392 out.insert(out.begin(), nanmean(x[Range(first_i, last_i-first_i)]));
393 last_i = first_i;
394 }
395
396 return out;
397}
398} // namespace Reduce
399
400namespace Mask {
401std::vector<bool> out_of_bounds(const Vector& x, const Numeric xmin, const Numeric xmax)
402{
403 std::vector<bool> mask(x.size(), false);
404 for (Index i=0; i<x.size(); i++) mask[i] = (xmin > x[i]) or (xmax < x[i]) or (not isnormal_or_zero(x[i]));
405 return mask;
406}
407
408VectorView mask(VectorView x, const std::vector<bool>& masking)
409{
410 const Index N = x.nelem();
411 for (Index i=0; i<N; i++)
412 if (masking[i])
413 x[i] = std::numeric_limits<Numeric>::quiet_NaN();
414 return x;
415}
416} // namespace Mask
417
418namespace Correction {
419 Numeric naive_tropospheric_singletau_median(const ConstVectorView& bt, const Numeric trop_bt, const Numeric target_bt) {
420 return - std::log((trop_bt - Average::nanmedian(bt)) / (trop_bt - target_bt));
421}
422
424 const Numeric t = std::exp(-tau);
425 if (std::isnormal(t)) {
426 bt /= t;
427 bt += trop_bt * std::expm1(-tau) / t;
428 }
429 return bt;
430}
431} // namespace Correction
432} // namespace Raw
void * data
TimeStep mean(const ArrayOfTimeStep &dt)
Returns the mean time step.
Definition: artstime.cc:190
void Reduce(Numeric &o, const Vector &i, const Verbosity &verbosity)
WORKSPACE METHOD: Reduce.
Definition: m_reduce.h:132
Index nelem() const ARTS_NOEXCEPT
Number of elements.
Definition: array.h:195
A constant view of a Vector.
Definition: matpackI.h:489
Index size() const ARTS_NOEXCEPT
Definition: matpackI.cc:50
Index nelem() const ARTS_NOEXCEPT
Returns the number of elements.
Definition: matpackI.cc:48
The MatrixView class.
Definition: matpackI.h:1125
The range class.
Definition: matpackI.h:165
The VectorView class.
Definition: matpackI.h:626
Iterator1D begin() ARTS_NOEXCEPT
Return iterator to first element.
Definition: matpackI.cc:143
Iterator1D end() ARTS_NOEXCEPT
Return iterator behind last element.
Definition: matpackI.cc:147
The Vector class.
Definition: matpackI.h:876
Constants of physical expressions as constexpr.
#define abs(x)
#define min(a, b)
#define dx
#define max(a, b)
Numeric last(ConstVectorView x)
last
Definition: math_funcs.cc:159
Numeric nanmean(const ConstVectorView &x) ARTS_NOEXCEPT
Mean function, vector version ignoring nans and infs.
Definition: matpackI.cc:1600
void transform(VectorView y, double(&my_func)(double), ConstVectorView x)
A generic transform function for vectors, which can be used to implement mathematical functions opera...
Definition: matpackI.cc:1472
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:39
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:33
constexpr auto pow2(T x) -> decltype(x *x)
power of two
Definition: constants.h:65
VectorView avg(VectorView y, const ArrayOfVector &ys, const Index start, const Index end_tmp)
Compute the average of the ranged ys.
Definition: raw.cc:112
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:219
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:177
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:159
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:260
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:130
VectorView std(VectorView std, const Vector &y, const ArrayOfVector &ys, const Index start, const Index end_tmp)
Compute the standard deviation of the ranged ys.
Definition: raw.cc:205
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:212
Numeric median(const ConstVectorView &v, const ArrayOfIndex &pos)
Get the median of the vector in the range.
Definition: raw.cc:238
ArrayOfVector caha(const ArrayOfVector &data, const Vector &tcvec, const Vector &thvec, const Index first_c_index)
Calibrate the data by CAHA.
Definition: raw.cc:53
Vector systemtemp(const Vector &pc, const Vector &ph, Numeric tc, Numeric th) noexcept
Computes the linear receiver temperature formula elementwise.
Definition: raw.cc:45
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:37
VectorView naive_tropospheric(VectorView bt, const Numeric tau, const Numeric trop_bt)
Apply tropospheric correction.
Definition: raw.cc:423
Numeric naive_tropospheric_singletau_median(const ConstVectorView &bt, const Numeric trop_bt, const Numeric target_bt)
Naive tropospheric correction parameterization.
Definition: raw.cc:419
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:401
VectorView mask(VectorView x, const std::vector< bool > &masking)
Masks all true entries of masking in x by turning them into NaN.
Definition: raw.cc:408
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:359
std::pair< Index, Index > find_first_and_last_1(const ArrayOfIndex &x)
Definition: raw.cc:310
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:294
Vector focus(const Vector &x, const ArrayOfIndex &scaling)
Returns the re-averaged values of x according to scaling.
Definition: raw.cc:325
Definition: raw.cc:35
constexpr Rational start(Rational Ju, Rational Jl, Polarization type) noexcept
Gives the lowest M for a polarization type of this transition.
Definition: zeemandata.h:78
constexpr Rational end(Rational Ju, Rational Jl, Polarization type) noexcept
Gives the largest M for a polarization type of this transition.
Definition: zeemandata.h:109
#define v
constexpr bool isnormal_or_zero(T x) noexcept
Definition: raw.cc:33
Stuff related to generating y-data from raw data.
#define N
Definition: rng.cc:164