18namespace Calibration {
19Vector
calibration(
const Vector& pc,
const Vector& pa,
const Vector& ph, Numeric tc, Numeric th)
noexcept
21 const auto N = Index(pc.size());
23 for (Index i=0; i<N; i++) calib[i] =
calibration(pc[i], pa[i], ph[i], tc, th);
27Vector
systemtemp(
const Vector& pc,
const Vector& ph, Numeric tc, Numeric th)
noexcept
29 const auto N = Index(pc.size());
31 for (Index i=0; i<N; i++) systemp[i] =
systemtemp(Numeric(pc[i]), Numeric(ph[i]), tc, th);
35ArrayOfVector
caha(
const ArrayOfVector& data,
const Vector& tcvec,
const Vector& thvec,
const Index first_c_index)
37 assert(first_c_index >= 0);
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);
56 const Vector * pc =
nullptr;
57 const Vector * pa =
nullptr;
58 const Vector * ph =
nullptr;
65 for (Index i=start; i<data.nelem(); i++) {
67 const Index pos = i % 4;
73 }
else if (pos == hpos) {
84 if ((pos == cpos or pos == hpos) and pc and pa and ph) {
85 out.emplace_back(
calibration(*pc, *pa, *ph, tc, th));
94VectorView
avg(VectorView y,
const ArrayOfVector& ys,
const Index start,
const Index end_tmp)
97 const Index end = end_tmp >= 0 ? end_tmp : 1 + ys.nelem() + end_tmp;
100 const Numeric scale = 1 / Numeric(end - start);
106 for (Index k=start; k<end; k++)
107 for (Index i=0; i<y.nelem(); i++)
108 y[i] += ys[k][i] * scale;
112VectorView
nanavg(VectorView y,
const ArrayOfVector& ys,
const Index start,
const Index end_tmp)
115 const Index end = end_tmp >= 0 ? end_tmp : 1 + ys.nelem() + end_tmp;
121 for (Index i=0; i<y.nelem(); i++) {
123 for (Index k=start; k<end; k++) {
132 y[i] /= Numeric(numnormal);
134 y[i] = std::numeric_limits<Numeric>::quiet_NaN();
141VectorView
var(VectorView
var,
const Vector& y,
const ArrayOfVector& ys,
const Index start,
const Index end_tmp)
144 const Index end = end_tmp >= 0 ? end_tmp : 1 + ys.nelem() + end_tmp;
147 const Numeric scale = 1 / Numeric(end - start);
153 for (Index k=start; k<end; k++)
154 for (Index i=0; i<y.nelem(); i++)
159VectorView
nanvar(VectorView
var,
const Vector& y,
const ArrayOfVector& ys,
const Index start,
const Index end_tmp)
162 const Index end = end_tmp >= 0 ? end_tmp : 1 + ys.nelem() + end_tmp;
168 for (Index i=0; i<y.nelem(); i++) {
170 for (Index k=start; k<end; k++) {
179 var[i] /= Numeric(numnormal);
181 var[i] = std::numeric_limits<Numeric>::quiet_NaN();
187VectorView
std(VectorView
std,
const Vector& y,
const ArrayOfVector& ys,
const Index start,
const Index end_tmp)
189 var(
std, y, ys, start, end_tmp);
190 std::transform(
std.begin(),
std.end(),
std.begin(), [](
const auto& x){return std::sqrt(x);});
194VectorView
nanstd(VectorView
std,
const Vector& y,
const ArrayOfVector& ys,
const Index start,
const Index end_tmp)
197 std::transform(
std.begin(),
std.end(),
std.begin(), [](
const auto& x){return std::sqrt(x);});
201MatrixView
cov(MatrixView
cov,
const Vector& y,
const ArrayOfVector& ys,
const Index start,
const Index end_tmp)
204 const Index end = end_tmp >= 0 ? end_tmp : 1 + ys.nelem() + end_tmp;
207 const Numeric scale = 1 / Numeric(end - start - 1);
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;
223 const Index n = pos.
nelem() ? pos.
nelem() :
v.nelem();
227 for (Index i=0; i<n; i++)
234 std::sort(calc.begin(), calc.end());
239 return (calc[(n-1)/2] + calc[n/2]) / 2;
245 const Index n = pos.
nelem() ? pos.
nelem() :
v.nelem();
249 for (Index i=0; i<n; i++) {
258 std::sort(calc.begin(), calc.end());
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);
267 return calc[numnormal/2];
268 return (calc[(numnormal-1)/2] + calc[numnormal/2]) / 2;
271 return std::numeric_limits<Numeric>::quiet_NaN();
281 auto scaling = [](Numeric X, Numeric X0, Numeric DX) {
282 const Index p = 1 + Index(std::abs(X - X0) / DX);
283 Index this_scale = 0;
284 while ((p >> this_scale) > 1) this_scale++;
285 return 1 << this_scale;
288 for (Index i=0; i<x.size(); i++) scale[i] = scaling(x[i], x0, dx);
294 Index first=x.
nelem() - 1;
296 for (Index i=0; i<x.
nelem(); i++) {
299 first = std::min(i, first);
303 assert(first <=
last);
304 return {first,
last};
309 const Index N = x.size();
310 Index first_i, last_i;
316 std::vector<Numeric> out;
317 for (Index i=firstone; i<=lastone; i++) out.emplace_back(x[i]);
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)]));
329 std::vector<Numeric> fx(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)]));
343 const Index N = x.size();
344 Index first_i, last_i;
350 std::vector<Numeric> out;
351 for (Index i=firstone; i<=lastone; i++) {
353 out.emplace_back(x[i]);
355 out.emplace_back(std::numeric_limits<Numeric>::quiet_NaN());
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)]));
369 std::vector<Numeric> fx(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)]));
383std::vector<bool>
out_of_bounds(
const Vector& x,
const Numeric xmin,
const Numeric xmax)
385 std::vector<bool>
mask(x.size(),
false);
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]));
390VectorView
mask(VectorView x,
const std::vector<bool>& masking)
392 const Index N = x.nelem();
393 for (Index i=0; i<N; i++)
395 x[i] = std::numeric_limits<Numeric>::quiet_NaN();
400namespace Correction {
406 const Numeric t = std::exp(-tau);
407 if (std::isnormal(t)) {
409 bt += trop_bt * std::expm1(-tau) / t;
TimeStep mean(const ArrayOfTimeStep &dt)
Returns the mean time step.
Index nelem() const ARTS_NOEXCEPT
void Reduce(Numeric &o, const Vector &i, const Verbosity &)
WORKSPACE METHOD: Reduce.
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.
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.
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.
VectorView var(VectorView var, const Vector &y, const ArrayOfVector &ys, const Index start, const Index end_tmp)
Compute the variance of the ranged ys.
Numeric nanmedian(const ConstVectorView &v, const ArrayOfIndex &pos)
Get the median of the vector in the range ignoring all non-normal values.
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.
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.
Numeric median(const ConstVectorView &v, const ArrayOfIndex &pos)
Get the median of the vector in the range.
ArrayOfVector caha(const ArrayOfVector &data, const Vector &tcvec, const Vector &thvec, const Index first_c_index)
Calibrate the data by CAHA.
Vector systemtemp(const Vector &pc, const Vector &ph, Numeric tc, Numeric th) noexcept
Computes the linear receiver temperature formula elementwise.
Vector calibration(const Vector &pc, const Vector &pa, const Vector &ph, Numeric tc, Numeric th) noexcept
Computes the linear calibration formula elementwise.
VectorView naive_tropospheric(VectorView bt, const Numeric tau, const Numeric trop_bt)
Apply tropospheric correction.
Numeric naive_tropospheric_singletau_median(const ConstVectorView &bt, const Numeric trop_bt, const Numeric target_bt)
Naive tropospheric correction parameterization.
std::vector< bool > out_of_bounds(const Vector &x, const Numeric xmin, const Numeric xmax)
Masks values that are out of bounds in x.
VectorView mask(VectorView x, const std::vector< bool > &masking)
Masks all true entries of masking in x by turning them into NaN.
Vector nanfocus(const Vector &x, const ArrayOfIndex &scaling)
Returns the re-averaged values of x according to scaling ignoring all non-normal values.
std::pair< Index, Index > find_first_and_last_1(const ArrayOfIndex &x)
ArrayOfIndex focus_doublescaling(const Vector &x, const Numeric x0, const Numeric dx)
Returns the relative position scale for each value in x.
Vector focus(const Vector &x, const ArrayOfIndex &scaling)
Returns the re-averaged values of x according to scaling.
constexpr bool isnormal_or_zero(T x) noexcept
Stuff related to generating y-data from raw data.