grouped_medians.hpp

Go to the documentation of this file.

#ifndef TATAMI_STATS_GROUPED_MEDIANS_HPP
#define TATAMI_STATS_GROUPED_MEDIANS_HPP
 
#include "utils.hpp"
#include "tatami/tatami.hpp"
#include "medians.hpp"
#include <vector>
#include <algorithm>
 
namespace tatami_stats {
 
namespace grouped_medians {
 
struct Options {
    bool skip_nan = false;
 
    int num_threads = 1;
};
 
template<typename Value_, typename Index_, typename Group_, class GroupSizes_, typename Output_>
void apply(bool row, const tatami::Matrix<Value_, Index_>* p, const Group_* group, const GroupSizes_& group_sizes, Output_** output, const Options& mopt) {
    Index_ dim = (row ? p->nrow() : p->ncol());
    Index_ otherdim = (row ? p->ncol() : p->nrow());
 
    tatami::parallelize([&](int, Index_ start, Index_ len) -> void {
        std::vector<Value_> xbuffer(otherdim);
 
        size_t ngroups = group_sizes.size();
        std::vector<std::vector<double> > workspace(ngroups);
        for (size_t g = 0; g < ngroups; ++g) {
            workspace[g].reserve(group_sizes[g]);
        }
 
        if (p->sparse()) {
            tatami::Options opt;
            opt.sparse_ordered_index = false;
 
            auto ext = tatami::consecutive_extractor<true>(p, row, start, len, opt);
            std::vector<Index_> ibuffer(otherdim);
            for (Index_ i = 0; i < len; ++i) {
                auto range = ext->fetch(xbuffer.data(), ibuffer.data());
                for (Index_ j = 0; j < range.number; ++j) {
                    workspace[group[range.index[j]]].push_back(range.value[j]);
                }
 
                for (size_t g = 0; g < ngroups; ++g) {
                    auto& w = workspace[g];
                    output[g][i + start] = medians::direct(w.data(), w.size(), static_cast<size_t>(group_sizes[g]), mopt.skip_nan);
                    w.clear();
                }
            }
 
        } else {
            auto ext = tatami::consecutive_extractor<false>(p, row, start, len);
            for (Index_ i = 0; i < len; ++i) {
                auto ptr = ext->fetch(xbuffer.data());
                for (Index_ j = 0; j < otherdim; ++j) {
                    workspace[group[j]].push_back(ptr[j]);
                }
 
                for (size_t g = 0; g < ngroups; ++g) {
                    auto& w = workspace[g];
                    output[g][i + start] = medians::direct(w.data(), w.size(), mopt.skip_nan);
                    w.clear();
                }
            }
        }
    }, dim, mopt.num_threads);
}
 
template<typename Output_ = double, typename Value_, typename Index_, typename Group_>
std::vector<std::vector<Output_> > by_row(const tatami::Matrix<Value_, Index_>* p, const Group_* group, const Options& mopt) {
    size_t mydim = p->nrow();
    auto group_sizes = tabulate_groups(group, p->ncol());
 
    std::vector<std::vector<Output_> > output(group_sizes.size());
    std::vector<Output_*> ptrs;
    ptrs.reserve(output.size());
    for (auto& o : output) {
        o.resize(mydim);
        ptrs.push_back(o.data());
    }
 
    apply(true, p, group, group_sizes, ptrs.data(), mopt);
    return output;
}
 
template<typename Output_ = double, typename Value_, typename Index_, typename Group_>
std::vector<std::vector<Output_> > by_row(const tatami::Matrix<Value_, Index_>* p, const Group_* group) {
    return by_row(p, group, Options());
}
 
template<typename Output_ = double, typename Value_, typename Index_, typename Group_>
std::vector<std::vector<Output_> > by_column(const tatami::Matrix<Value_, Index_>* p, const Group_* group, const Options& mopt) {
    size_t mydim = p->ncol();
    auto group_sizes = tabulate_groups(group, p->nrow());
 
    std::vector<std::vector<Output_> > output(group_sizes.size());
    std::vector<Output_*> ptrs;
    ptrs.reserve(output.size());
    for (auto& o : output) {
        o.resize(mydim);
        ptrs.push_back(o.data());
    }
 
    apply(false, p, group, group_sizes, ptrs.data(), mopt);
    return output;
}
 
template<typename Output_ = double, typename Value_, typename Index_, typename Group_>
std::vector<std::vector<Output_> > by_column(const tatami::Matrix<Value_, Index_>* p, const Group_* group) {
    return by_column(p, group, Options());
}
 
}
 
}
 
#endif