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"
#include "quickstats/quickstats.hpp"
 
namespace scran_blocks {
 
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) {
        quickstats::SingleQuantileVariableNumber<Output_, I<decltype(nblocks)> > 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.data());
        }
 
    } else {
        quickstats::SingleQuantileFixedNumber<Output_, I<decltype(nblocks)> > 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.data());
        }
    }
}
 
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