Commit 24fcb9c4 authored by Thibault Hallouin's avatar Thibault Hallouin
Browse files

relocate std deviation in maths header

Showing with 16 additions and 17 deletions
+16 -17
...@@ -6,23 +6,10 @@ ...@@ -6,23 +6,10 @@
#include <xtensor/xexpression.hpp> #include <xtensor/xexpression.hpp>
#include <xtensor/xtensor.hpp> #include <xtensor/xtensor.hpp>
#include "../maths.hpp"
namespace evalhyd namespace evalhyd
{ {
namespace detail {
// function to calculate standard deviation on last axis
//
// TODO: substitute with `xt::stddev` when fixed for `xt::rtensor`
// (see https://gitlab.irstea.fr/HYCAR-Hydro/evalhyd/evalhyd-r/-/issues/1)
template <class A1, class A2>
inline auto std_dev(A1&& arr, A2&& mean_arr)
{
return xt::sqrt(
xt::nanmean(xt::square(arr - mean_arr), -1)
);
}
}
namespace determinist namespace determinist
{ {
class Evaluator class Evaluator
...@@ -326,8 +313,8 @@ namespace evalhyd ...@@ -326,8 +313,8 @@ namespace evalhyd
auto prd = xt::view(prd_masked, xt::all(), b_exp[e]); auto prd = xt::view(prd_masked, xt::all(), b_exp[e]);
auto obs = xt::view(obs_masked, xt::all(), b_exp[e]); auto obs = xt::view(obs_masked, xt::all(), b_exp[e]);
xt::view(alpha, m, e) = xt::view(alpha, m, e) =
detail::std_dev(prd, xt::view(mean_prd, m, e)) evalhyd::maths::nanstd(prd, xt::view(mean_prd, m, e))
/ detail::std_dev(obs, xt::view(mean_obs, m, e)); / evalhyd::maths::nanstd(obs, xt::view(mean_obs, m, e));
} }
} }
} }
......
...@@ -13,6 +13,18 @@ namespace evalhyd ...@@ -13,6 +13,18 @@ namespace evalhyd
{ {
namespace maths namespace maths
{ {
// TODO: substitute with `xt::stddev` when fixed for `xt::rtensor`
// (see https://gitlab.irstea.fr/HYCAR-Hydro/evalhyd/evalhyd-r/-/issues/1)
// 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)
{
return xt::sqrt(
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) inline double quantile(const xt::xtensor<double, 1>& t, double p)
{ {
std::size_t n = t.size(); std::size_t n = t.size();
......
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