parallel_quantiles.hpp

Go to the documentation of this file.

#ifndef SCRAN_BLOCKS_PARALLEL_QUANTILES_HPP
#define SCRAN_BLOCKS_PARALLEL_QUANTILES_HPP
 
#include <vector>
#include <limits>
#include <algorithm>
#include <vector>
#include <cstddef>
#include <cmath>
#include <optional>
#include <cassert>
 
#include "utils.hpp"
 
#include "sanisizer/sanisizer.hpp"
 
namespace scran_blocks {
 
template<typename Output_, class Iterator_>
class SingleQuantile {
public:
    SingleQuantile(const std::size_t len, const double quantile) {
        assert(len > 0);
        sanisizer::can_ptrdiff<Iterator_>(len);
        const Output_ frac = static_cast<Output_>(len - 1) * static_cast<Output_>(quantile);
        const Output_ base = std::floor(frac);
        my_lower = base; // cast is known-safe if can_ptrdiff passes and 0 <= quantile <= 1.
        my_upper_fraction = frac - base;
        my_lower_fraction = static_cast<Output_>(1) - my_upper_fraction;
        my_skip_upper = my_upper_fraction == 0;
    }
 
private:
    std::size_t my_lower;
    Output_ my_upper_fraction;
    Output_ my_lower_fraction;
    bool my_skip_upper;
 
public:
    Output_ operator()(Iterator_ begin, Iterator_ end) const {
        assert(sanisizer::is_less_than(my_lower, end - begin));
        auto target = begin + my_lower;
        std::nth_element(begin, target, end);
 
        // Avoid looking at target+1 if don't need it - in particular, we'd get a out-of-bounds access if quantile = 1.
        if (my_skip_upper) {
            return *target;
        } else {
            const auto next = std::min_element(target + 1, end); 
            return static_cast<Output_>(*target) * my_lower_fraction + static_cast<Output_>(*next) * my_upper_fraction;
        }
    }
};
 
template<typename Output_, class Iterator_>
class SingleQuantileVariable {
public:
    SingleQuantileVariable(const std::size_t max_len, const double quantile) : my_quantile(quantile) {
        if (max_len >= 2) {
            sanisizer::resize(my_choices, max_len - 1);
        }
    }
 
private:
    std::vector<std::optional<SingleQuantile<Output_, Iterator_> > > my_choices;
    double my_quantile;
 
public:
    Output_ operator()(const std::size_t len, Iterator_ begin, Iterator_ end) {
        assert(sanisizer::is_equal(len, end - begin));
        if (len == 0) {
            return std::numeric_limits<Output_>::quiet_NaN();
        } else if (len == 1) {
            return *begin;
        } else {
            auto& ocalc = my_choices[len - 2];
            if (!ocalc.has_value()) { // Instantiating them on demand.
                ocalc.emplace(len, my_quantile);
            }
            return (*ocalc)(begin, end);
        }
    }
 
    Output_ operator()(Iterator_ begin, Iterator_ end) {
        return this->operator()(end - begin, begin, end);
    }
};
 
template<typename Stat_, typename Output_>
void parallel_quantiles(const std::size_t n, const std::vector<Stat_*>& in, const double quantile, Output_* const out, const bool skip_nan) {
    const auto nblocks = in.size();
    if (nblocks == 0) {
        std::fill_n(out, n, std::numeric_limits<Output_>::quiet_NaN());
        return;
    } else if (nblocks == 1) {
        std::copy_n(in[0], n, out);
        return;
    }
 
    std::vector<Stat_> tmp_buffer;
    tmp_buffer.reserve(nblocks);
 
    if (skip_nan) {
        SingleQuantileVariable<Output_, I<decltype(tmp_buffer.begin())> > calcs(nblocks, quantile);
        for (std::size_t g = 0; g < n; ++g) {
            tmp_buffer.clear();
            for (I<decltype(nblocks)> b = 0; b < nblocks; ++b) {
                const auto val = in[b][g];
                if (!std::isnan(val)) {
                    tmp_buffer.push_back(val);
                }
            }
            out[g] = calcs(tmp_buffer.size(), tmp_buffer.begin(), tmp_buffer.end());
        }
 
    } else {
 
        SingleQuantile<Output_, I<decltype(tmp_buffer.begin())> > calc(nblocks, quantile);
        for (std::size_t g = 0; g < n; ++g) {
            tmp_buffer.clear();
            for (I<decltype(nblocks)> b = 0; b < nblocks; ++b) {
                tmp_buffer.push_back(in[b][g]);
            }
            out[g] = calc(tmp_buffer.begin(), tmp_buffer.end());
        }
    }
}
 
template<typename Output_ = double, typename Stat_>
std::vector<Output_> parallel_quantiles(const std::size_t n, const std::vector<Stat_*>& in, const double quantile, const bool skip_nan) {
    auto out = sanisizer::create<std::vector<Output_> >(n);
    parallel_quantiles(n, in, quantile, out.data(), skip_nan);
    return out;
}
 
}
 
#endif