quantile.hpp

Go to the documentation of this file.

#ifndef TATAMI_STATS_QUANTILE_HPP
#define TATAMI_STATS_QUANTILE_HPP
 
#include "utils.hpp"
 
#include <cmath>
#include <vector>
#include <algorithm>
#include <limits>
#include <type_traits>
 
#include "tatami/tatami.hpp"
#include "sanisizer/sanisizer.hpp"
#include "quickstats/quickstats.hpp"
 
namespace tatami_stats {
 
struct QuantileOptions {
    bool skip_nan = false;
 
    int num_threads = 1;
};
 
template<typename Value_, typename Index_, typename Output_>
void quantile(
    const bool row,
    const tatami::Matrix<Value_, Index_>& mat,
    const double prob,
    Output_* const output,
    const QuantileOptions& opt
) {
    const auto dim = (row ? mat.nrow() : mat.ncol());
    const auto otherdim = (row ? mat.ncol() : mat.nrow());
    if (otherdim == 0) {
        // Prevent initalize() from constructing with a fixed instance with otherdim == 0.
        std::fill_n(output, dim, std::numeric_limits<Output_>::quiet_NaN());
        return;
    }
 
    tatami::parallelize([&](int, Index_ s, Index_ l) -> void {
        std::optional<quickstats::SingleQuantileFixedNumber<Output_> > qcalcs_fixed;
        std::optional<quickstats::SingleQuantileVariableNumber<Output_> > qcalcs_var;
        // Index_ is safe to cast to std::size_t as that's part of the tatami contract.
        nanable_ifelse<Value_>(
            opt.skip_nan,
            [&]() -> void {
                qcalcs_var.emplace(otherdim, prob);
            },
            [&]() -> void {
                qcalcs_fixed.emplace(otherdim, prob);
            }
        );
 
        if (mat.sparse()) {
            tatami::Options topt;
            topt.sparse_extract_index = false;
            topt.sparse_ordered_index = false; // we'll be sorting by value anyway.
            auto ext = tatami::consecutive_extractor<true>(mat, row, s, l, topt);
            auto buffer = tatami::create_container_of_Index_size<std::vector<Value_> >(otherdim);
            auto vbuffer = buffer.data();
 
            for (Index_ x = 0; x < l; ++x) {
                auto range = ext->fetch(vbuffer, NULL);
                tatami::copy_n(range.value, range.number, vbuffer);
 
                nanable_ifelse<Value_>(
                    opt.skip_nan,
                    [&]() -> void {
                        const auto new_non_zeros = shift_nans(vbuffer, range.number);
                        output[x + s] = (*qcalcs_var)(otherdim - (range.number - new_non_zeros), new_non_zeros, vbuffer);
                    },
                    [&]() -> void {
                        output[x + s] = (*qcalcs_fixed)(range.number, vbuffer);
                    }
                );
            }
 
        } else {
            auto buffer = tatami::create_container_of_Index_size<std::vector<Value_> >(otherdim);
            auto ext = tatami::consecutive_extractor<false>(mat, row, s, l);
 
            for (Index_ x = 0; x < l; ++x) {
                auto bufptr = buffer.data();
                auto raw = ext->fetch(bufptr);
                tatami::copy_n(raw, otherdim, bufptr);
 
                nanable_ifelse<Value_>(
                    opt.skip_nan,
                    [&]() -> void {
                        const auto new_total = shift_nans(bufptr, otherdim);
                        output[x + s] = (*qcalcs_var)(new_total, bufptr);
                    },
                    [&]() -> void {
                        output[x + s] = (*qcalcs_fixed)(bufptr);
                    }
                );
            }
        }
    }, dim, opt.num_threads);
}
 
template<typename Output_ = double, typename Value_, typename Index_>
std::vector<Output_> quantile(
    const bool row,
    const tatami::Matrix<Value_, Index_>& mat,
    const double prob,
    const QuantileOptions& opt
) {
    const auto dim = (row ? mat.nrow() : mat.ncol());
    auto output = sanisizer::create<std::vector<Output_> >(dim
#ifdef TATAMI_STATS_TEST_DIRTY
        , -1
#endif
    );
    quantile(row, mat, prob, output.data(), opt);
    return output;
}
 
}
 
#endif