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