median.hpp

Go to the documentation of this file.

#ifndef TATAMI_STATS_MEDIAN_HPP
#define TATAMI_STATS_MEDIAN_HPP
 
#include "utils.hpp"
 
#include <cmath>
#include <vector>
#include <algorithm>
#include <limits>
 
#include "tatami/tatami.hpp"
#include "sanisizer/sanisizer.hpp"
#include "quickstats/quickstats.hpp"
 
namespace tatami_stats {
 
struct MedianOptions {
    bool skip_nan = false;
 
    int num_threads = 1;
};
 
template<typename Output_ = double, typename Value_, typename Index_>
Output_ median_direct(Value_* ptr, Index_ num, bool skip_nan) {
    nanable_ifelse<Value_>(
        skip_nan,
        [&]() -> void {
            num = shift_nans(ptr, num);
        },
        []() -> void {}
    );
 
    return quickstats::median<Output_>(num, ptr);
}
 
template<typename Output_ = double, typename Value_, typename Index_>
Output_ median_direct(Value_* value, Index_ num_nonzero, Index_ num_all, bool skip_nan) {
    nanable_ifelse<Value_>(
        skip_nan,
        [&]() -> void {
            auto new_nonzero = shift_nans(value, num_nonzero);
            num_all -= num_nonzero - new_nonzero;
            num_nonzero = new_nonzero;
        },
        []() -> void {}
    );
 
    return quickstats::median<Output_>(num_all, num_nonzero, value);
}
template<typename Value_, typename Index_, typename Output_>
void median(const bool row, const tatami::Matrix<Value_, Index_>& mat, Output_* const output, const MedianOptions& opt) {
    const auto dim = (row ? mat.nrow() : mat.ncol());
    const auto otherdim = (row ? mat.ncol() : mat.nrow());
 
    if (mat.sparse()) {
        tatami::Options topt;
        topt.sparse_extract_index = false;
        topt.sparse_ordered_index = false; // we'll be sorting by value anyway.
 
        tatami::parallelize([&](int, Index_ s, Index_ l) -> void {
            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);
                output[x + s] = median_direct<Output_>(vbuffer, range.number, otherdim, opt.skip_nan);
            }
        }, dim, opt.num_threads);
 
    } else {
        tatami::parallelize([&](int, Index_ s, Index_ l) -> void {
            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 ptr = ext->fetch(buffer.data());
                tatami::copy_n(ptr, otherdim, buffer.data());
                output[x + s] = median_direct<Output_>(buffer.data(), otherdim, opt.skip_nan);
            }
        }, dim, opt.num_threads);
    }
}
 
template<typename Output_ = double, typename Value_, typename Index_>
std::vector<Output_> median(const bool row, const tatami::Matrix<Value_, Index_>& mat, const MedianOptions& opt) {
    const auto dim = (row ? mat.nrow() : mat.ncol());
    auto output = sanisizer::create<std::vector<Output_> >(dim
#ifdef TATAMI_STATS_TEST_DIRTY
        , -1
#endif
    );
    median(row, mat, output.data(), opt);
    return output;
}
 
}
 
#endif