Commit 01e05adb authored by Thibault Hallouin's avatar Thibault Hallouin
Browse files

replace custom-made quantile function with new `xt::quantile`

1 merge request!3release v0.1.0
Pipeline #43888 failed with stage
in 3 minutes and 54 seconds
Showing with 17 additions and 85 deletions
+17 -85
......@@ -19,8 +19,6 @@
#include <xtensor/xsort.hpp>
#include <xtensor/xindex_view.hpp>
#include "maths.hpp"
typedef std::map<std::string, std::vector<std::vector<std::string>>> msk_tree;
......@@ -234,6 +232,8 @@ namespace evalhyd
auto q = get_q();
// define lambda function to precompute mean/median/quantile
auto get_val =
[&](const std::string& str, const std::string& num)
{
......@@ -251,9 +251,8 @@ namespace evalhyd
}
else // (str == "quantile")
{
return evalhyd::maths::quantile(q, std::stod(num));
return xt::quantile(q, {std::stod(num)})();
}
};
// preprocess conditions to identify special cases
......
......@@ -27,54 +27,6 @@ namespace evalhyd
xt::nanmean(xt::square(xt::abs(arr - mean_arr)), -1)
);
}
// function to calculate quantile on 1-dim expressions
inline double quantile(const xt::xtensor<double, 1>& t, double p)
{
std::size_t n = t.size();
// compute virtual index
auto virtual_idx = (double(n) - 1) * p;
if (std::fmod(virtual_idx, 1) == 0)
// virtual index is an integer
{
std::size_t i = {static_cast<unsigned long long>(virtual_idx)};
std::array<std::size_t, 1> kth = {i};
auto values = xt::partition(t, kth);
return xt::mean(xt::view(values, xt::range(i, i + 1)))();
}
else
// virtual index is not an integer
{
// determine indices before and after virtual index
std::size_t prv_idx = std::floor(virtual_idx);
std::size_t nxt_idx = prv_idx + 1;
// deal with edge cases
if (virtual_idx > double(n) - 1)
{
prv_idx = n - 1;
nxt_idx = n - 1;
}
if (virtual_idx < 0)
{
prv_idx = 0;
nxt_idx = 0;
}
std::array<std::size_t, 2> kth = {prv_idx, nxt_idx};
auto values = xt::partition(t, kth);
auto vw = xt::view(values, xt::range(prv_idx, nxt_idx + 1));
// perform linear interpolation to determine quantile
auto gamma = virtual_idx - double(prv_idx);
auto prv = xt::amin(vw);
auto nxt = xt::amax(vw);
return (prv + (nxt - prv) * gamma)();
}
}
}
}
......
......@@ -10,8 +10,6 @@
#include <xtensor/xindex_view.hpp>
#include <xtensor/xsort.hpp>
#include "../maths.hpp"
namespace evalhyd
{
......@@ -61,31 +59,14 @@ namespace evalhyd
xt::col(quantiles, 1) = 0.5 + c_lvl / 200.;
// compute predictive interval bounds from quantiles
// TODO: replace with `xt::quantile` when available
for (std::size_t s = 0; s < n_sit; s++)
for (std::size_t i = 0; i < n_itv; i++)
{
for (std::size_t l = 0; l < n_ldt; l++)
{
for (std::size_t i = 0; i < n_itv; i++)
{
for (std::size_t t = 0; t < n_tim; t++)
{
// lower bound
xt::view(itv_bnds, s, l, i, 0, t) =
evalhyd::maths::quantile(
xt::view(q_prd, s, l, xt::all(), t),
quantiles(i, 0)
);
// upper bound
xt::view(itv_bnds, s, l, i, 1, t) =
evalhyd::maths::quantile(
xt::view(q_prd, s, l, xt::all(), t),
quantiles(i, 1)
);
}
}
}
auto q = xt::quantile(q_prd, xt::view(quantiles, i), 2);
xt::view(itv_bnds, xt::all(), xt::all(), i, 0, xt::all()) =
xt::view(q, 0);
xt::view(itv_bnds, xt::all(), xt::all(), i, 1, xt::all()) =
xt::view(q, 1);
}
return itv_bnds;
......@@ -214,7 +195,7 @@ namespace evalhyd
// compute "climatology" interval
// TODO: replace with `xt::quantile` when available
// TODO: replace with `xt::nanquantile` when available
for (std::size_t s = 0; s < n_sit; s++)
{
for (std::size_t l = 0; l < n_ldt; l++)
......@@ -228,14 +209,14 @@ namespace evalhyd
// lower bound
xt::view(clim_bnds, s, l, m, e, i, 0) =
evalhyd::maths::quantile(
obs_filtered, quantiles(i, 0)
);
xt::quantile(
obs_filtered, {quantiles(i, 0)}
)();
// upper bound
xt::view(clim_bnds, s, l, m, e, i, 1) =
evalhyd::maths::quantile(
obs_filtered, quantiles(i, 1)
);
xt::quantile(
obs_filtered, {quantiles(i, 1)}
)();
}
}
}
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment