score_markers_pairwise.hpp

Go to the documentation of this file.

#ifndef SCRAN_MARKERS_SCORE_MARKERS_PAIRWISE_HPP
#define SCRAN_MARKERS_SCORE_MARKERS_PAIRWISE_HPP
 
#include <vector>
#include <cstddef>
#include <optional>
 
#include "scran_blocks/scran_blocks.hpp"
#include "tatami/tatami.hpp"
#include "tatami_stats/tatami_stats.hpp"
#include "sanisizer/sanisizer.hpp"
 
#include "scan_matrix.hpp"
#include "average_group_stats.hpp"
#include "block_averages.hpp"
#include "create_combinations.hpp"
#include "cohens_d.hpp"
#include "simple_diff.hpp"
#include "utils.hpp"
 
namespace scran_markers {
 
struct ScoreMarkersPairwiseOptions {
    double threshold = 0;
 
    int num_threads = 1;
 
    bool compute_group_mean = true;
 
    bool compute_group_detected = true;
 
    bool compute_cohens_d = true;
 
    bool compute_auc = true;
 
    bool compute_delta_mean = true;
 
    bool compute_delta_detected = true;
 
    BlockAveragePolicy block_average_policy = BlockAveragePolicy::MEAN;
 
    scran_blocks::WeightPolicy block_weight_policy = scran_blocks::WeightPolicy::VARIABLE;
 
    scran_blocks::VariableWeightParameters variable_block_weight_parameters;
 
    double block_quantile = 0.5;
};
 
template<typename Stat_>
struct ScoreMarkersPairwiseBuffers {
    std::vector<Stat_*> mean;
 
    std::vector<Stat_*> detected;
 
    Stat_* cohens_d = NULL;
 
    Stat_* auc = NULL;
 
    Stat_* delta_mean = NULL;
 
    Stat_* delta_detected = NULL;
};
 
template<typename Stat_>
struct ScoreMarkersPairwiseResults {
    std::vector<std::vector<Stat_> > mean;
 
    std::vector<std::vector<Stat_> > detected;
 
    std::vector<Stat_> cohens_d;
 
    std::vector<Stat_> auc;
 
    std::vector<Stat_> delta_mean;
 
    std::vector<Stat_> delta_detected;
};
 
namespace internal {
 
template<typename Index_, typename Stat_>
void process_simple_pairwise_effects(
    const Index_ ngenes,
    const std::size_t ngroups,
    const std::size_t nblocks,
    const std::size_t ncombos,
    const std::vector<Stat_>& combo_means,
    const std::vector<Stat_>& combo_vars,
    const std::vector<Stat_>& combo_detected,
    const double threshold,
    const BlockAverageInfo<Stat_>& average_info, 
    const ScoreMarkersPairwiseBuffers<Stat_>& output,
    const int num_threads
) {
    const Stat_* total_weights_ptr = NULL;
    std::optional<std::vector<Stat_> > total_weights_per_group;
    if (average_info.use_mean()) {
        if (!output.mean.empty() || !output.detected.empty()) {
            if (nblocks > 1) {
                total_weights_per_group = compute_total_weight_per_group(ngroups, nblocks, average_info.combo_weights().data());
                total_weights_ptr = total_weights_per_group->data();
            } else {
                total_weights_ptr = average_info.combo_weights().data();
            }
        }
    }
 
    std::optional<PrecomputedPairwiseWeights<Stat_> > preweights;
    if (average_info.use_mean()) {
        if (output.cohens_d != NULL || output.delta_mean != NULL || output.delta_detected != NULL) {
            preweights.emplace(ngroups, nblocks, average_info.combo_weights().data());
        }
    }
 
    tatami::parallelize([&](const int, const Index_ start, const Index_ length) -> void {
        std::optional<std::vector<Stat_> > qbuffer, qrevbuffer;
        std::optional<scran_blocks::SingleQuantileVariable<Stat_, typename std::vector<Stat_>::iterator> > qcalc;
        if (!average_info.use_mean()) {
            qbuffer.emplace();
            qrevbuffer.emplace();
            qcalc.emplace(nblocks, average_info.quantile());
        }
 
        for (Index_ gene = start, end = start + length; gene < end; ++gene) {
            auto in_offset = sanisizer::product_unsafe<std::size_t>(gene, ncombos);
 
            if (!output.mean.empty()) {
                const auto tmp_means = combo_means.data() + in_offset;
                if (average_info.use_mean()) {
                    average_group_stats_blockmean(gene, ngroups, nblocks, tmp_means, average_info.combo_weights().data(), total_weights_ptr, output.mean);
                } else {
                    average_group_stats_blockquantile(gene, ngroups, nblocks, tmp_means, *qbuffer, *qcalc, output.mean);
                }
            }
 
            if (!output.detected.empty()) {
                const auto tmp_detected = combo_detected.data() + in_offset;
                if (average_info.use_mean()) {
                    average_group_stats_blockmean(gene, ngroups, nblocks, tmp_detected, average_info.combo_weights().data(), total_weights_ptr, output.detected);
                } else {
                    average_group_stats_blockquantile(gene, ngroups, nblocks, tmp_detected, *qbuffer, *qcalc, output.detected);
                }
            }
 
            // Computing the effect sizes.
            const auto out_offset = sanisizer::product_unsafe<std::size_t>(gene, ngroups, ngroups);
 
            if (output.cohens_d != NULL) {
                const auto tmp_means = combo_means.data() + in_offset;
                const auto tmp_variances = combo_vars.data() + in_offset;
                const auto outptr = output.cohens_d + out_offset;
                if (average_info.use_mean()) {
                    compute_pairwise_cohens_d_blockmean(tmp_means, tmp_variances, ngroups, nblocks, threshold, *preweights, outptr);
                } else {
                    compute_pairwise_cohens_d_blockquantile(tmp_means, tmp_variances, ngroups, nblocks, threshold, *qbuffer, *qrevbuffer, *qcalc, outptr);
                }
            }
 
            if (output.delta_detected != NULL) {
                const auto tmp_detected = combo_detected.data() + in_offset;
                const auto outptr = output.delta_detected + out_offset;
                if (average_info.use_mean()) {
                    compute_pairwise_simple_diff_blockmean(tmp_detected, ngroups, nblocks, *preweights, outptr);
                } else {
                    compute_pairwise_simple_diff_blockquantile(tmp_detected, ngroups, nblocks, *qbuffer, *qcalc, outptr);
                }
            }
 
            if (output.delta_mean != NULL) {
                const auto tmp_means = combo_means.data() + in_offset;
                const auto outptr = output.delta_mean + out_offset;
                if (average_info.use_mean()) {
                    compute_pairwise_simple_diff_blockmean(tmp_means, ngroups, nblocks, *preweights, outptr);
                } else {
                    compute_pairwise_simple_diff_blockquantile(tmp_means, ngroups, nblocks, *qbuffer, *qcalc, outptr);
                }
            }
        }
    }, ngenes, num_threads);
}
 
template<typename Index_, typename Stat_>
ScoreMarkersPairwiseBuffers<Stat_> preallocate_pairwise_results(
    const Index_ ngenes,
    const std::size_t ngroups,
    ScoreMarkersPairwiseResults<Stat_>& store,
    const ScoreMarkersPairwiseOptions& opt
) {
    ScoreMarkersPairwiseBuffers<Stat_> output;
 
    if (opt.compute_group_mean) {
        preallocate_average_results(ngenes, ngroups, store.mean, output.mean);
    }
    if (opt.compute_group_detected) {
        preallocate_average_results(ngenes, ngroups, store.detected, output.detected);
    }
 
    const auto num_effect_sizes = sanisizer::product<typename std::vector<Stat_>::size_type>(ngenes, ngroups, ngroups);
 
    if (opt.compute_cohens_d) {
        store.cohens_d.resize(num_effect_sizes
#ifdef SCRAN_MARKERS_TEST_INIT
            , SCRAN_MARKERS_TEST_INIT
#endif
        );
        output.cohens_d = store.cohens_d.data();
    }
    if (opt.compute_auc) {
        store.auc.resize(num_effect_sizes
#ifdef SCRAN_MARKERS_TEST_INIT
            , SCRAN_MARKERS_TEST_INIT
#endif
        );
        output.auc = store.auc.data();
    }
    if (opt.compute_delta_mean) {
        store.delta_mean.resize(num_effect_sizes
#ifdef SCRAN_MARKERS_TEST_INIT
            , SCRAN_MARKERS_TEST_INIT
#endif
        );
        output.delta_mean = store.delta_mean.data();
    }
    if (opt.compute_delta_detected) {
        store.delta_detected.resize(num_effect_sizes
#ifdef SCRAN_MARKERS_TEST_INIT
            , SCRAN_MARKERS_TEST_INIT
#endif
        );
        output.delta_detected = store.delta_detected.data();
    }
 
    return output;
}
 
template<
    bool single_block_,
    typename Value_,
    typename Index_,
    typename Group_,
    typename Block_,
    typename Stat_
>
void score_markers_pairwise(
    const tatami::Matrix<Value_, Index_>& matrix, 
    const std::size_t ngroups,
    const Group_* const group, 
    const std::size_t nblocks,
    const Block_* const block,
    const std::size_t ncombos,
    const std::size_t* const combo,
    const std::vector<Index_>& combo_sizes,
    const ScoreMarkersPairwiseOptions& options,
    const ScoreMarkersPairwiseBuffers<Stat_>& output
) {
    const auto ngenes = matrix.nrow();
    const auto payload_size = sanisizer::product<typename std::vector<Stat_>::size_type>(ngenes, ncombos);
    std::vector<Stat_> combo_means, combo_vars, combo_detected;
    if (!output.mean.empty() || output.cohens_d != NULL || output.delta_mean != NULL) {
        combo_means.resize(payload_size);
    }
    if (output.cohens_d != NULL) {
        combo_vars.resize(payload_size);
    }
    if (!output.detected.empty() || output.delta_detected != NULL) {
        combo_detected.resize(payload_size);
    }
 
    // For a single block, this usually doesn't really matter, but we do it for consistency with the multi-block case,
    // and to account for variable weighting where non-zero block sizes get zero weight.
    BlockAverageInfo<Stat_> average_info;
    if (options.block_average_policy == BlockAveragePolicy::MEAN) {
        average_info = BlockAverageInfo<Stat_>(
            scran_blocks::compute_weights<Stat_>(
                combo_sizes,
                options.block_weight_policy,
                options.variable_block_weight_parameters
            )
        );
    } else {
        average_info = BlockAverageInfo<Stat_>(options.block_quantile);
    }
 
    if (output.auc != NULL || matrix.prefer_rows()) {
        scan_matrix_by_row_full_auc<single_block_>(
            matrix, 
            ngroups,
            group,
            nblocks,
            block,
            ncombos,
            combo,
            combo_sizes,
            average_info,
            combo_means,
            combo_vars,
            combo_detected,
            output.auc,
            options.threshold,
            options.num_threads
        );
 
    } else {
        scan_matrix_by_column(
            matrix,
            [&]{
                if constexpr(single_block_) {
                    return ngroups;
                } else {
                    return ncombos;
                }
            }(),
            [&]{
                if constexpr(single_block_) {
                    return group;
                } else {
                    return combo;
                }
            }(),
            combo_sizes,
            combo_means,
            combo_vars,
            combo_detected,
            options.num_threads
        );
    }
 
    process_simple_pairwise_effects(
        matrix.nrow(),
        ngroups,
        nblocks,
        ncombos,
        combo_means,
        combo_vars,
        combo_detected,
        options.threshold,
        average_info,
        output,
        options.num_threads
    );
}
 
}
template<typename Value_, typename Index_, typename Group_, typename Stat_>
void score_markers_pairwise(
    const tatami::Matrix<Value_, Index_>& matrix, 
    const Group_* const group, 
    const ScoreMarkersPairwiseOptions& options,
    const ScoreMarkersPairwiseBuffers<Stat_>& output
) {
    const Index_ NC = matrix.ncol();
    const auto group_sizes = tatami_stats::tabulate_groups(group, NC); 
    const auto ngroups = sanisizer::cast<std::size_t>(group_sizes.size());
 
    internal::score_markers_pairwise<true>(
        matrix,
        ngroups,
        group,
        1,
        static_cast<int*>(NULL),
        ngroups,
        static_cast<std::size_t*>(NULL),
        group_sizes,
        options,
        output
    );
}
 
template<typename Value_, typename Index_, typename Group_, typename Block_, typename Stat_>
void score_markers_pairwise_blocked(
    const tatami::Matrix<Value_, Index_>& matrix, 
    const Group_* const group, 
    const Block_* const block,
    const ScoreMarkersPairwiseOptions& options,
    const ScoreMarkersPairwiseBuffers<Stat_>& output
) {
    const Index_ NC = matrix.ncol();
    const auto ngroups = output.mean.size();
    const auto nblocks = tatami_stats::total_groups(block, NC); 
 
    const auto combinations = internal::create_combinations(ngroups, group, block, NC);
    const auto combo_sizes = internal::tabulate_combinations<Index_>(ngroups, nblocks, combinations);
    const auto ncombos = combo_sizes.size();
 
    internal::score_markers_pairwise<false>(
        matrix,
        sanisizer::cast<std::size_t>(ngroups),
        group,
        sanisizer::cast<std::size_t>(nblocks),
        block,
        sanisizer::cast<std::size_t>(ncombos),
        combinations.data(),
        combo_sizes,
        options,
        output
    );
}
 
template<typename Stat_ = double, typename Value_, typename Index_, typename Group_>
ScoreMarkersPairwiseResults<Stat_> score_markers_pairwise(const tatami::Matrix<Value_, Index_>& matrix, const Group_* const group, const ScoreMarkersPairwiseOptions& options) {
    const auto ngroups = tatami_stats::total_groups(group, matrix.ncol());
    ScoreMarkersPairwiseResults<Stat_> res;
    auto buffers = internal::preallocate_pairwise_results(matrix.nrow(), ngroups, res, options);
    score_markers_pairwise(matrix, group, options, buffers);
    return res; 
}
 
template<typename Stat_ = double, typename Value_, typename Index_, typename Group_, typename Block_>
ScoreMarkersPairwiseResults<Stat_> score_markers_pairwise_blocked(
    const tatami::Matrix<Value_, Index_>& matrix,
    const Group_* const group,
    const Block_* const block,
    const ScoreMarkersPairwiseOptions& options)
{
    const auto ngroups = tatami_stats::total_groups(group, matrix.ncol());
    ScoreMarkersPairwiseResults<Stat_> res;
    const auto buffers = internal::preallocate_pairwise_results(matrix.nrow(), ngroups, res, options);
    score_markers_pairwise_blocked(matrix, group, block, options, buffers);
    return res;
}
 
}
 
#endif