maths.hpp

// Copyright (c) 2023, INRAE.
// Distributed under the terms of the GPL-3 Licence.
// The full licence is in the file LICENCE, distributed with this software.

#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
    {
        // 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();

            // 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