allow mean/median/quantile# in masking conditions

src/masks.hpp

@@ -15,6 +15,10 @@
 #include <xtensor/xsort.hpp>
 #include <xtensor/xindex_view.hpp>
 
+#include "maths.hpp"
+
+namespace eh = evalhyd;
+
 typedef std::map<std::string, std::vector<std::vector<std::string>>> msk_tree;
 
 namespace evalhyd
@@ -30,7 +34,7 @@ namespace evalhyd
             // observed or predicted (median or mean for probabilist) streamflow
             // e.g. q{>9.} q{<9} q{>=99.0} q{<=99} q{>9,<99} q{==9} q{!=9}
             std::regex exp_q (
-                    R"((q_obs|q_prd_median|q_prd_mean)\{((([><!=]?=?[0-9]+\.?[0-9]*),*)+)\})"
+                    R"((q_obs|q_prd_median|q_prd_mean)\{((([><!=]?=?(mean|median|quantile[0-9]+\.?[0-9]*|[0-9]+\.?[0-9]*)),*)+)\})"
             );
 
             for (std::sregex_iterator i =
@@ -40,33 +44,41 @@ namespace evalhyd
                 const std::smatch & mtc = *i;
 
                 std::string var = mtc[1];
-                std::string s = mtc[2];
+                std::string str = mtc[2];
 
                 // process masking conditions on streamflow
                 std::vector<std::vector<std::string>> conditions;
 
                 // pattern supported to specify masking conditions based on streamflow
-                std::regex e (R"(([><!=]?=?)([0-9]+\.?[0-9]*))");
+                std::regex ex (R"(([><!=]?=?)(mean|median|quantile|[0-9]+\.?[0-9]*)([0-9]+\.?[0-9]*)?)");
 
                 for (std::sregex_iterator j =
-                        std::sregex_iterator(s.begin(), s.end(), e);
+                        std::sregex_iterator(str.begin(), str.end(), ex);
                      j != std::sregex_iterator(); j++)
                 {
-                    const std::smatch & m = *j;
+                    const std::smatch & mt = *j;
 
                     // check that operator is provided and is supported
                     std::set<std::string> supported_op =
                             {"<", ">", "<=", ">=", "!=", "=="};
-                    if (supported_op.find(m[1]) != supported_op.end())
-                        conditions.push_back({m[1].str(), m[2].str()});
-                    else if (m[1].str().empty())
+                    if (mt[1].str().empty())
                         throw std::runtime_error(
                                 "missing operator for streamflow masking condition"
                         );
+                    else if (supported_op.find(mt[1]) != supported_op.end())
+                    {
+                        if ((mt[2].str() == "median")
+                            or (mt[2].str() == "mean")
+                            or (mt[2].str() == "quantile"))
+                            conditions.push_back({mt[1].str(), mt[2].str(), mt[3].str()});
+                        else
+                            // it is a simple numerical value, swap last two
+                            conditions.push_back({mt[1].str(), mt[3].str(), mt[2].str()});
+                    }
                     else
                         throw std::runtime_error(
                                 "invalid operator for streamflow masking "
-                                "condition: " + m[1].str()
+                                "condition: " + mt[1].str()
                         );
                 }
 
@@ -193,6 +205,19 @@ 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) {
+                        if (str.empty())
+                            // it is a simple numerical value
+                            return std::stod(num);
+                        else if (str == "median")
+                            return xt::median(q);
+                        else if (str == "mean")
+                            return xt::mean(q)();
+                        else  // (str == "quantile")
+                            return eh::maths::quantile(q, std::stod(num));
+                    };
+
                     // preprocess conditions to identify special cases
                     // within/without
                     bool within = false;
@@ -204,9 +229,9 @@ namespace evalhyd
                     if (cond.size() == 2)
                     {
                         opr1 = cond[0][0];
-                        val1 = std::stod(cond[0][1]);
+                        val1= get_val(cond[0][1], cond[0][2]);
                         opr2 = cond[1][0];
-                        val2 = std::stod(cond[1][1]);
+                        val2 = get_val(cond[1][1], cond[1][2]);
 
                         if ((opr1 == "<") or (opr1 == "<="))
                         {
@@ -291,9 +316,7 @@ namespace evalhyd
                         for (const auto & opr_val : cond)
                         {
                             auto opr = opr_val[0];
-
-                            // convert comparison value from string to double
-                            double val = std::stod(opr_val[1]);
+                            double val = get_val(opr_val[1], opr_val[2]);
 
                             // apply masking condition to given subset
                             if (opr == "<")

src/maths.hpp

+#ifndef EVALHYD_MATHS_HPP
+#define EVALHYD_MATHS_HPP
+
+#include <xtensor/xtensor.hpp>
+#include <xtensor/xview.hpp>
+#include <xtensor/xsort.hpp>
+#include <xtensor/xbuilder.hpp>
+#include <xtensor/xutils.hpp>
+
+#include <cmath>
+
+namespace evalhyd
+{
+    namespace maths
+    {
+        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)();
+            }
+        }
+    }
+}
+
+#endif //EVALHYD_MATHS_HPP