diff --git a/CMakeLists.txt b/CMakeLists.txt
index 5ab913f99a660bfbd6cc26b07c37b595b9e7e98e..27df2280a43567e4bbbe2893bf6630b313f4fd81 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -18,7 +18,7 @@ project(
find_package(xtl 0.7.5 REQUIRED)
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")
# ------------------------------------------------------------------------------
diff --git a/include/evalhyd/detail/determinist/efficiencies.hpp b/include/evalhyd/detail/determinist/efficiencies.hpp
index af83745fa530c1e8f349348a6742c1ed7c4c1a3c..6ae18ff67ee925025ec18d90aa5252252def1757 100644
--- a/include/evalhyd/detail/determinist/efficiencies.hpp
+++ b/include/evalhyd/detail/determinist/efficiencies.hpp
@@ -9,6 +9,9 @@
#include <xtensor/xview.hpp>
#include <xtensor/xoperation.hpp>
+#include "../maths.hpp"
+
+
namespace evalhyd
{
namespace determinist
diff --git a/include/evalhyd/detail/maths.hpp b/include/evalhyd/detail/maths.hpp
index f8b44d55fc1ee539d3c33c3c8e5f8e58dce874df..49786ad67145a847f5b5faaf5b94b5d94ac80630 100644
--- a/include/evalhyd/detail/maths.hpp
+++ b/include/evalhyd/detail/maths.hpp
@@ -17,8 +17,8 @@ namespace evalhyd
{
namespace maths
{
- // TODO: substitute with `xt::stddev` when fixed for `xt::rtensor`
- // (see https://gitlab.irstea.fr/HYCAR-Hydro/evalhyd/evalhyd-r/-/issues/1)
+ // TODO: substitute with `xt::stddev` when performance fixed
+ // (see https://github.com/xtensor-stack/xtensor/pull/2656)
// function to calculate standard deviation on last axis of n-dim expressions
template <class A1, class A2>
inline auto nanstd(A1&& arr, A2&& mean_arr)
@@ -27,56 +27,6 @@ namespace evalhyd
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)();
- }
- }
}
}
diff --git a/include/evalhyd/detail/probabilist/intervals.hpp b/include/evalhyd/detail/probabilist/intervals.hpp
index ceca4a2531c961f57f18bcbc884e6529c125c8d2..0a2f4c2a2e3df0fe2fbc802be8c0fdcd4d0af515 100644
--- a/include/evalhyd/detail/probabilist/intervals.hpp
+++ b/include/evalhyd/detail/probabilist/intervals.hpp
@@ -11,9 +11,6 @@
#include <xtensor/xview.hpp>
#include <xtensor/xindex_view.hpp>
#include <xtensor/xsort.hpp>
-#include <xtensor/xio.hpp>
-
-#include "../maths.hpp"
namespace evalhyd
@@ -64,34 +61,16 @@ namespace evalhyd
xt::col(quantiles, 1) = 0.5 + c_lvl / 200.;
// compute predictive interval bounds from quantiles
-
- // 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 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;
}
@@ -217,8 +196,6 @@ namespace evalhyd
xt::view(q_obs_masked, xt::all(), xt::all(), b_exp[e]);
// compute "climatology" interval
-
- // 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++)
@@ -234,15 +211,16 @@ 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)}
+ )();
}
else
{
diff --git a/include/evalhyd/detail/uncertainty.hpp b/include/evalhyd/detail/uncertainty.hpp
index 36b459ea100830b62f75a74cfa82f00c3725d450..25b13af0957e7552f52fe648564f2abc67be0054 100644
--- a/include/evalhyd/detail/uncertainty.hpp
+++ b/include/evalhyd/detail/uncertainty.hpp
@@ -16,8 +16,8 @@
#include <xtensor/xtensor.hpp>
#include <xtensor/xadapt.hpp>
#include <xtensor/xrandom.hpp>
+#include <xtensor/xsort.hpp>
-#include "maths.hpp"
typedef std::chrono::time_point<
std::chrono::system_clock, std::chrono::minutes
diff --git a/include/evalhyd/evald.hpp b/include/evalhyd/evald.hpp
index 834583319a7852ecb503214fd2bcebfd4e07f226..f56db11eed76781139529369e982a5744a82e07c 100644
--- a/include/evalhyd/evald.hpp
+++ b/include/evalhyd/evald.hpp
@@ -15,7 +15,6 @@
#include "detail/utils.hpp"
#include "detail/masks.hpp"
-#include "detail/maths.hpp"
#include "detail/uncertainty.hpp"
#include "detail/determinist/evaluator.hpp"
diff --git a/include/evalhyd/evalp.hpp b/include/evalhyd/evalp.hpp
index fa5e7e752ac59d0fd1f58e156a54e291003c164e..5dfe36984077f92cff08fe908e7755721d5eab49 100644
--- a/include/evalhyd/evalp.hpp
+++ b/include/evalhyd/evalp.hpp
@@ -16,7 +16,6 @@
#include "detail/utils.hpp"
#include "detail/masks.hpp"
-#include "detail/maths.hpp"
#include "detail/uncertainty.hpp"
#include "detail/probabilist/evaluator.hpp"