From 24fcb9c42944bd50d18001e0edd4681ad3dc2286 Mon Sep 17 00:00:00 2001
From: Thibault Hallouin <thibault.hallouin@inrae.fr>
Date: Thu, 20 Oct 2022 15:15:45 +0200
Subject: [PATCH] relocate std deviation in maths header

---
 src/determinist/evaluator.hpp | 21 ++++-----------------
 src/maths.hpp                 | 12 ++++++++++++
 2 files changed, 16 insertions(+), 17 deletions(-)

diff --git a/src/determinist/evaluator.hpp b/src/determinist/evaluator.hpp
index bec074d..0fdab32 100644
--- a/src/determinist/evaluator.hpp
+++ b/src/determinist/evaluator.hpp
@@ -6,23 +6,10 @@
 #include <xtensor/xexpression.hpp>
 #include <xtensor/xtensor.hpp>
 
+#include "../maths.hpp"
+
 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
     {
         class Evaluator
@@ -326,8 +313,8 @@ namespace evalhyd
                     auto prd = xt::view(prd_masked, xt::all(), b_exp[e]);
                     auto obs = xt::view(obs_masked, xt::all(), b_exp[e]);
                     xt::view(alpha, m, e) =
-                            detail::std_dev(prd, xt::view(mean_prd, m, e))
-                            / detail::std_dev(obs, xt::view(mean_obs, m, e));
+                            evalhyd::maths::nanstd(prd, xt::view(mean_prd, m, e))
+                            / evalhyd::maths::nanstd(obs, xt::view(mean_obs, m, e));
                 }
             }
         }
diff --git a/src/maths.hpp b/src/maths.hpp
index f0613fb..7fccf31 100644
--- a/src/maths.hpp
+++ b/src/maths.hpp
@@ -13,6 +13,18 @@ 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)
+        // 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)
         {
             std::size_t n = t.size();
-- 
GitLab