Skip to content
GitLab
Commit 5fe62314 authored by Thibault Hallouin's avatar Thibault Hallouin
Browse files

use newly fixed xtensor quantile function 

available in xtensor 0.24., see https://github.com/xtensor-stack/xtensor/issues/2651
parent cbe1588d
main 11-check-the-computation-of-the-ranks-for-ties dev v0.1.2 v0.1.1 v0.1.0
1 merge request!3release v0.1.0
Pipeline #45464 failed with stage
in 3 minutes and 44 seconds
Hide whitespace changes
Inline Side-by-side
Showing
with 23 additions and 94 deletions
+23 -94
CMakeLists.txt
+ 1
− 1
View file @ 5fe62314
...@@ -18,7 +18,7 @@ project( ...@@ -18,7 +18,7 @@ project(
find_package(xtl 0.7.5 REQUIRED) find_package(xtl 0.7.5 REQUIRED)
message(STATUS "Found xtl: ${xtl_INCLUDE_DIRS}/xtl") message(STATUS "Found xtl: ${xtl_INCLUDE_DIRS}/xtl")
find_package(xtensor 0.24.4 REQUIRED) find_package(xtensor 0.24.6 REQUIRED)
message(STATUS "Found xtensor: ${xtensor_INCLUDE_DIRS}/xtensor") message(STATUS "Found xtensor: ${xtensor_INCLUDE_DIRS}/xtensor")
# ------------------------------------------------------------------------------ # ------------------------------------------------------------------------------
......
include/evalhyd/detail/determinist/efficiencies.hpp
+ 3
− 0
View file @ 5fe62314
...@@ -9,6 +9,9 @@ ...@@ -9,6 +9,9 @@
#include <xtensor/xview.hpp> #include <xtensor/xview.hpp>
#include <xtensor/xoperation.hpp> #include <xtensor/xoperation.hpp>
#include "../maths.hpp"
namespace evalhyd namespace evalhyd
{ {
namespace determinist namespace determinist
......
include/evalhyd/detail/maths.hpp
+ 2
− 52
View file @ 5fe62314
...@@ -17,8 +17,8 @@ namespace evalhyd ...@@ -17,8 +17,8 @@ namespace evalhyd
{ {
namespace maths namespace maths
{ {
// TODO: substitute with `xt::stddev` when fixed for `xt::rtensor` // TODO: substitute with `xt::stddev` when performance fixed
// (see https://gitlab.irstea.fr/HYCAR-Hydro/evalhyd/evalhyd-r/-/issues/1) // (see https://github.com/xtensor-stack/xtensor/pull/2656)
// function to calculate standard deviation on last axis of n-dim expressions // function to calculate standard deviation on last axis of n-dim expressions
template <class A1, class A2> template <class A1, class A2>
inline auto nanstd(A1&& arr, A2&& mean_arr) inline auto nanstd(A1&& arr, A2&& mean_arr)
...@@ -27,56 +27,6 @@ namespace evalhyd ...@@ -27,56 +27,6 @@ namespace evalhyd
xt::nanmean(xt::square(xt::abs(arr - mean_arr)), -1) xt::nanmean(xt::square(xt::abs(arr - mean_arr)), -1)
); );
} }
// TODO: drop in favour of xt::quantile when xtensor issue fixed
// https://github.com/xtensor-stack/xtensor/issues/2651
// 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)();
}
}
} }
} }
......
include/evalhyd/detail/probabilist/intervals.hpp
+ 16
− 38
View file @ 5fe62314
...@@ -11,9 +11,6 @@ ...@@ -11,9 +11,6 @@
#include <xtensor/xview.hpp> #include <xtensor/xview.hpp>
#include <xtensor/xindex_view.hpp> #include <xtensor/xindex_view.hpp>
#include <xtensor/xsort.hpp> #include <xtensor/xsort.hpp>
#include <xtensor/xio.hpp>
#include "../maths.hpp"
namespace evalhyd namespace evalhyd
...@@ -64,34 +61,16 @@ namespace evalhyd ...@@ -64,34 +61,16 @@ namespace evalhyd
xt::col(quantiles, 1) = 0.5 + c_lvl / 200.; xt::col(quantiles, 1) = 0.5 + c_lvl / 200.;
// compute predictive interval bounds from quantiles // compute predictive interval bounds from quantiles
for (std::size_t i = 0; i < n_itv; i++)
// TODO: replace with xt::quantile when xtensor issue fixed
// https://github.com/xtensor-stack/xtensor/issues/2651
for (std::size_t s = 0; s < n_sit; s++)
{ {
for (std::size_t l = 0; l < n_ldt; l++) auto q = xt::quantile(q_prd, xt::view(quantiles, i), 2);
{
for (std::size_t i = 0; i < n_itv; i++) xt::view(itv_bnds, xt::all(), xt::all(), i, 0, xt::all()) =
{ xt::view(q, 0);
for (std::size_t t = 0; t < n_tim; t++) xt::view(itv_bnds, xt::all(), xt::all(), i, 1, xt::all()) =
{ xt::view(q, 1);
// 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)
);
}
}
}
} }
return itv_bnds; return itv_bnds;
} }
...@@ -217,8 +196,6 @@ namespace evalhyd ...@@ -217,8 +196,6 @@ namespace evalhyd
xt::view(q_obs_masked, xt::all(), xt::all(), b_exp[e]); xt::view(q_obs_masked, xt::all(), xt::all(), b_exp[e]);
// compute "climatology" interval // compute "climatology" interval
// TODO: replace with `xt::nanquantile` when available
for (std::size_t s = 0; s < n_sit; s++) for (std::size_t s = 0; s < n_sit; s++)
{ {
for (std::size_t l = 0; l < n_ldt; l++) for (std::size_t l = 0; l < n_ldt; l++)
...@@ -234,15 +211,16 @@ namespace evalhyd ...@@ -234,15 +211,16 @@ namespace evalhyd
{ {
// lower bound // lower bound
xt::view(clim_bnds, s, l, m, e, i, 0) = xt::view(clim_bnds, s, l, m, e, i, 0) =
evalhyd::maths::quantile( xt::quantile(
obs_filtered, quantiles(i, 0) obs_filtered,
); {quantiles(i, 0)}
)();
// upper bound // upper bound
xt::view(clim_bnds, s, l, m, e, i, 1) = xt::view(clim_bnds, s, l, m, e, i, 1) =
evalhyd::maths::quantile( xt::quantile(
obs_filtered, quantiles(i, 1) obs_filtered,
); {quantiles(i, 1)}
)();
} }
else else
{ {
......
include/evalhyd/detail/uncertainty.hpp
+ 1
− 1
View file @ 5fe62314
...@@ -16,8 +16,8 @@ ...@@ -16,8 +16,8 @@
#include <xtensor/xtensor.hpp> #include <xtensor/xtensor.hpp>
#include <xtensor/xadapt.hpp> #include <xtensor/xadapt.hpp>
#include <xtensor/xrandom.hpp> #include <xtensor/xrandom.hpp>
#include <xtensor/xsort.hpp>
#include "maths.hpp"
typedef std::chrono::time_point< typedef std::chrono::time_point<
std::chrono::system_clock, std::chrono::minutes std::chrono::system_clock, std::chrono::minutes
......
include/evalhyd/evald.hpp
+ 0
− 1
View file @ 5fe62314
...@@ -15,7 +15,6 @@ ...@@ -15,7 +15,6 @@
#include "detail/utils.hpp" #include "detail/utils.hpp"
#include "detail/masks.hpp" #include "detail/masks.hpp"
#include "detail/maths.hpp"
#include "detail/uncertainty.hpp" #include "detail/uncertainty.hpp"
#include "detail/determinist/evaluator.hpp" #include "detail/determinist/evaluator.hpp"
......
include/evalhyd/evalp.hpp
+ 0
− 1
View file @ 5fe62314
...@@ -16,7 +16,6 @@ ...@@ -16,7 +16,6 @@
#include "detail/utils.hpp" #include "detail/utils.hpp"
#include "detail/masks.hpp" #include "detail/masks.hpp"
#include "detail/maths.hpp"
#include "detail/uncertainty.hpp" #include "detail/uncertainty.hpp"
#include "detail/probabilist/evaluator.hpp" #include "detail/probabilist/evaluator.hpp"
......
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