blocked.hpp

#ifndef MEDIAN_MARKERS_BLOCKED_HPP
#define MEDIAN_MARKERS_BLOCKED_HPP
 
#include <cstddef>
#include <optional>
#include <vector>
#include <limits>
#include <cmath>
 
#include "sanisizer/sanisizer.hpp"
#include "tatami/tatami.hpp"
#include "tatami_stats/tatami_stats.hpp"
 
#include "queue.hpp"
#include "scan.hpp"
#include "number.hpp"
 
namespace singler_classic_markers {
 
struct ChooseBlockedOptions {
    std::optional<std::size_t> number;
 
    bool keep_ties = false;
 
    bool use_minimum = false;
 
    int num_threads = 1;
};
 
template<bool include_stat_, typename Stat_, typename Value_, typename Index_, typename Label_, typename Block_>
Markers<include_stat_, Index_, Stat_> choose_blocked_raw(
    const tatami::Matrix<Value_, Index_>& matrix, 
    const Label_* label,
    const Block_* block,
    const ChooseBlockedOptions& options
) {
    const auto NC = matrix.ncol();
    const std::size_t ngroups = tatami_stats::total_groups/*<std::size_t>*/(label, NC);
    const std::size_t nblocks = tatami_stats::total_groups/*<std::size_t>*/(block, NC);
 
    const auto num_keep = get_num_keep<Index_>(ngroups, options.number);
    auto pqueues = sanisizer::create<std::vector<std::optional<PairwiseTopQueues<Stat_, Index_> > > >(options.num_threads);
 
    // Creating the combinations between block and not.
    const auto ncombos = sanisizer::product<std::size_t>(ngroups, nblocks); // check that all producs below are safe.
    auto combinations = sanisizer::create<std::vector<std::size_t> >(NC);
    for (I<decltype(NC)> c = 0; c < NC; ++c) {
        combinations[c] = sanisizer::nd_offset<std::size_t>(label[c], ngroups, block[c]); // group is the faster changing dimension.
    }
    auto combo_sizes = tatami_stats::tabulate_groups(combinations.data(), NC);
 
    const auto num_used = scan_matrix<Stat_>(
        matrix,
        sanisizer::cast<std::size_t>(ncombos),
        combinations.data(),
        combo_sizes,
 
        /* setup = */ [&]() -> PairwiseTopQueues<Stat_, Index_> {
            PairwiseTopQueues<Stat_, Index_> output;
            allocate_pairwise_queues(output, num_keep, ngroups, options.keep_ties, /* check_nan = */ false); // we'll check it ourselves.
            return output;
        },
 
        /* fun = */ [&](const Index_ r, const std::vector<Stat_>& medians, PairwiseTopQueues<Stat_, Index_>& curqueues) -> void {
            for (I<decltype(ngroups)> g1 = 1; g1 < ngroups; ++g1) {
                for (I<decltype(ngroups)> g2 = 0; g2 < g1; ++g2) {
 
                    if (options.use_minimum) {
                        Stat_ xval = std::numeric_limits<Stat_>::infinity();
                        Stat_ yval = std::numeric_limits<Stat_>::infinity();
                        for (I<decltype(nblocks)> b = 0; b < nblocks; ++b) {
                            const auto delta = medians[sanisizer::nd_offset<std::size_t>(g1, ngroups, b)] - medians[sanisizer::nd_offset<std::size_t>(g2, ngroups, b)];
                            if (!std::isnan(delta)) {
                                xval = std::min(xval, delta);
                                yval = std::min(yval, -delta);
                            }
                        }
 
                        if (std::isfinite(xval)) {
                            curqueues[g1][g2].emplace(xval, r); 
                            curqueues[g2][g1].emplace(yval, r); 
                        }
 
                    } else {
                        Stat_ val = 0;
                        std::size_t denom = 0;
                        for (I<decltype(nblocks)> b = 0; b < nblocks; ++b) {
                            const auto delta = medians[sanisizer::nd_offset<std::size_t>(g1, ngroups, b)] - medians[sanisizer::nd_offset<std::size_t>(g2, ngroups, b)];
                            if (!std::isnan(delta)) {
                                ++denom;
                                val += delta; 
                            }
                        }
 
                        if (denom) {
                            val /= denom;
                            curqueues[g1][g2].emplace(val, r); 
                            curqueues[g2][g1].emplace(-val, r); 
                        }
                    }
 
                }
            }
        },
 
        /* finalize = */ [&](const int t, PairwiseTopQueues<Stat_, Index_>& curqueues) -> void {
            pqueues[t] = std::move(curqueues);
        },
 
        options.num_threads
    );
 
    pqueues.resize(num_used); 
    Markers<include_stat_, Index_, Stat_> output;
    report_best_top_queues<include_stat_>(pqueues, ngroups, output);
    return output;
}
template<typename Stat_ = double, typename Value_, typename Index_, typename Label_, typename Block_>
std::vector<std::vector<std::vector<std::pair<Index_, Stat_> > > > choose_blocked(
    const tatami::Matrix<Value_, Index_>& matrix, 
    const Label_* label,
    const Block_* block,
    const ChooseBlockedOptions& options
) {
    return choose_blocked_raw<true, Stat_>(matrix, label, block, options);
}
 
template<typename Stat_ = double, typename Value_, typename Index_, typename Label_, typename Block_>
std::vector<std::vector<std::vector<Index_> > > choose_blocked_index(
    const tatami::Matrix<Value_, Index_>& matrix, 
    const Label_* label,
    const Block_* block,
    const ChooseBlockedOptions& options
) {
    return choose_blocked_raw<false, Stat_>(matrix, label, block, options);
}
 
}
 
#endif