PairwiseEffects.hpp

Go to the documentation of this file.

#ifndef SCRAN_PAIRWISE_EFFECTS_HPP
#define SCRAN_PAIRWISE_EFFECTS_HPP
 
#include "../utils/macros.hpp"
 
#include "MatrixCalculator.hpp"
#include "cohens_d.hpp"
#include "simple_diff.hpp"
 
#include "../utils/vector_to_pointers.hpp"
#include "../utils/blocking.hpp"
#include "../utils/average_vectors.hpp"
 
#include "tatami/tatami.hpp"
 
namespace scran {
 
class PairwiseEffects {
public:
    struct Defaults {
        static constexpr double threshold = 0;
 
        static constexpr int num_threads = 1;
 
        static constexpr bool compute_cohen = true;
 
        static constexpr bool compute_auc = true;
 
        static constexpr bool compute_lfc = true;
 
        static constexpr bool compute_delta_detected = true;
 
        static constexpr WeightPolicy block_weight_policy = WeightPolicy::VARIABLE;
 
        static constexpr VariableBlockWeightParameters variable_block_weight_parameters = VariableBlockWeightParameters();
    };
 
private:
    double threshold = Defaults::threshold;
    int nthreads = Defaults::num_threads;
    bool do_cohen = Defaults::compute_cohen;
    bool do_auc = Defaults::compute_auc;
    bool do_lfc = Defaults::compute_lfc;
    bool do_delta_detected = Defaults::compute_delta_detected;
 
    WeightPolicy block_weight_policy = Defaults::block_weight_policy;
    VariableBlockWeightParameters variable_block_weight_parameters = Defaults::variable_block_weight_parameters;
 
public:
    PairwiseEffects& set_threshold(double t = Defaults::threshold) {
        threshold = t;
        return *this;
    }
 
    PairwiseEffects& set_num_threads(int n = Defaults::num_threads) {
        nthreads = n;
        return *this;
    }
 
    PairwiseEffects& set_compute_cohen(bool c = Defaults::compute_cohen) {
        do_cohen = c;
        return *this;
    }
 
    PairwiseEffects& set_compute_auc(bool c = Defaults::compute_auc) {
        do_auc = c;
        return *this;
    }
 
    PairwiseEffects& set_compute_lfc(bool c = Defaults::compute_lfc) {
        do_lfc = c;
        return *this;
    }
 
    PairwiseEffects& set_compute_delta_detected(bool c = Defaults::compute_delta_detected) {
        do_delta_detected = c;
        return *this;
    }
 
    PairwiseEffects& set_block_weight_policy(WeightPolicy w = Defaults::block_weight_policy) {
        block_weight_policy = w;
        return *this;
    }
 
    PairwiseEffects& set_variable_block_weight_parameters(VariableBlockWeightParameters v = Defaults::variable_block_weight_parameters) {
        variable_block_weight_parameters = v;
        return *this;
    }
 
public:
    template<typename Data_, typename Index_, typename Group_, typename Stat_>
    void run(
        const tatami::Matrix<Data_, Index_>* p, 
        const Group_* group, 
        std::vector<Stat_*> means, 
        std::vector<Stat_*> detected, 
        Stat_* cohen,
        Stat_* auc,
        Stat_* lfc,
        Stat_* delta_detected) 
    const {
        int ngroups = means.size();
        differential_analysis::MatrixCalculator runner(nthreads, threshold, block_weight_policy, variable_block_weight_parameters);
        Overlord overlord(auc);
        auto state = runner.run(p, group, ngroups, overlord);
        process_simple_effects(p->nrow(), ngroups, 1, state, means, detected, cohen, lfc, delta_detected);
    }
 
    template<typename Data_, typename Index_, typename Group_, typename Block_, typename Stat_>
    void run_blocked(
        const tatami::Matrix<Data_, Index_>* p, 
        const Group_* group, 
        const Block_* block, 
        std::vector<Stat_*> means, 
        std::vector<Stat_*> detected, 
        Stat_* cohen,
        Stat_* auc,
        Stat_* lfc,
        Stat_* delta_detected)
    const {
        if (block == NULL) {
            run(p, group, std::move(means), std::move(detected), cohen, auc, lfc, delta_detected);
 
        } else {
            auto ngenes = p->nrow();
            int ngroups = means.size();
            int nblocks = count_ids(p->ncol(), block);
 
            int ncombos = ngroups * nblocks;
            std::vector<std::vector<Stat_> > means_store(ncombos), detected_store(ncombos);
            std::vector<Stat_*> means2(ncombos), detected2(ncombos);
            for (int c = 0; c < ncombos; ++c) {
                means_store[c].resize(ngenes);
                detected_store[c].resize(ngenes);
                means2[c] = means_store[c].data();
                detected2[c] = detected_store[c].data();
            }
 
            differential_analysis::MatrixCalculator runner(nthreads, threshold, block_weight_policy, variable_block_weight_parameters);
            Overlord overlord(auc);
            auto state = runner.run_blocked(p, group, ngroups, block, nblocks, overlord);
            process_simple_effects(p->nrow(), ngroups, nblocks, state, means2, detected2, cohen, lfc, delta_detected);
 
            // Averaging the remaining statistics.
            std::vector<double> weights(nblocks);
            std::vector<Stat_*> mstats(nblocks), dstats(nblocks);
 
            for (int gr = 0; gr < ngroups; ++gr) {
                for (int b = 0; b < nblocks; ++b) {
                    size_t offset = gr * static_cast<size_t>(nblocks) + b;
                    weights[b] = state.level_weight[offset];
                    mstats[b] = means2[offset];
                    dstats[b] = detected2[offset];
                }
 
                average_vectors_weighted(ngenes, mstats, weights.data(), means[gr]);
                average_vectors_weighted(ngenes, dstats, weights.data(), detected[gr]);
            }
        }
    }
 
private:
    template<typename Stat_>
    struct Overlord {
        Overlord(Stat_* auc_) : auc(auc_) {}
 
        bool needs_auc() const { 
            return auc != NULL;
        }
 
        Stat_* auc;
 
        Stat_* prepare_auc_buffer(size_t gene, size_t ngroups) { 
            return auc + gene * ngroups * ngroups;
        }
    };
 
    template<typename Stat_>
    void process_simple_effects(
        size_t ngenes, // using size_t consistently here, to avoid integer overflow bugs when computing products.
        size_t ngroups,
        size_t nblocks,
        const differential_analysis::MatrixCalculator::State& state, 
        std::vector<Stat_*>& means, 
        std::vector<Stat_*>& detected, 
        Stat_* cohen, 
        Stat_* lfc, 
        Stat_* delta_detected) 
    const {
        const auto& level_size = state.level_size;
        size_t nlevels = level_size.size();
 
        tatami::parallelize([&](size_t, size_t start, size_t length) -> void {
            auto in_offset = nlevels * start;
            const auto* tmp_means = state.means.data() + in_offset;
            const auto* tmp_variances = state.variances.data() + in_offset;
            const auto* tmp_detected = state.detected.data() + in_offset;
 
            size_t squared = ngroups * ngroups;
            size_t out_offset = start * squared;
            for (size_t gene = start, end = start + length; gene < end; ++gene, out_offset += squared) {
                for (size_t l = 0; l < nlevels; ++l) {
                    means[l][gene] = tmp_means[l];
                    detected[l][gene] = tmp_detected[l];
                }
 
                if (cohen != NULL) {
                    differential_analysis::compute_pairwise_cohens_d(tmp_means, tmp_variances, state.level_weight, ngroups, nblocks, threshold, cohen + out_offset);
                }
 
                if (delta_detected != NULL) {
                    differential_analysis::compute_pairwise_simple_diff(tmp_detected, state.level_weight, ngroups, nblocks, delta_detected + out_offset);
                }
 
                if (lfc != NULL) {
                    differential_analysis::compute_pairwise_simple_diff(tmp_means, state.level_weight, ngroups, nblocks, lfc + out_offset);
                }
 
                tmp_means += nlevels;
                tmp_variances += nlevels;
                tmp_detected += nlevels;
            }
        }, ngenes, nthreads);
    }
 
public:
    template<typename Stat_>
    struct Results {
        Results(size_t ngenes, size_t ngroups, bool do_cohen, bool do_auc, bool do_lfc, bool do_delta) :
            means(ngroups, std::vector<Stat_>(ngenes)),
            detected(ngroups, std::vector<Stat_>(ngenes))
        {
            size_t nelements = ngenes * ngroups * ngroups;
            if (do_cohen) {
                cohen.resize(nelements);
            }
            if (do_auc) {
                auc.resize(nelements);
            }
            if (do_lfc) {
                lfc.resize(nelements);
            }
            if (do_delta) {
                delta_detected.resize(nelements);
            }
        }
        std::vector<Stat_> cohen;
 
        std::vector<Stat_> auc;
 
        std::vector<Stat_> lfc;
 
        std::vector<Stat_> delta_detected;
 
        std::vector<std::vector<Stat_> > means;
 
        std::vector<std::vector<Stat_> > detected;
    };
 
    template<typename Stat_ = double, typename Data_ = double, typename Index_ = int, typename Group_ = int>
    Results<Stat_> run(const tatami::Matrix<Data_, Index_>* p, const Group_* group) {
        size_t ngroups = count_ids(p->ncol(), group);
        Results<Stat_> res(p->nrow(), ngroups, do_cohen, do_auc, do_lfc, do_delta_detected); 
        run(
            p, 
            group, 
            vector_to_pointers(res.means), 
            vector_to_pointers(res.detected), 
            harvest_pointer(res.cohen, do_cohen),
            harvest_pointer(res.auc, do_auc),
            harvest_pointer(res.lfc, do_lfc),
            harvest_pointer(res.delta_detected, do_delta_detected)
        );
        return res; 
    }
 
    template<typename Stat_ = double, typename Data_ = double, typename Index_ = int, typename Group_ = int, typename Block_ = int>
    Results<Stat_> run_blocked(const tatami::Matrix<Data_, Index_>* p, const Group_* group, const Block_* block) {
        size_t ngroups = count_ids(p->ncol(), group);
        Results<Stat_> res(p->nrow(), ngroups, do_cohen, do_auc, do_lfc, do_delta_detected); 
        run_blocked(
            p,
            group,
            block,
            vector_to_pointers(res.means),
            vector_to_pointers(res.detected),
            harvest_pointer(res.cohen, do_cohen),
            harvest_pointer(res.auc, do_auc),
            harvest_pointer(res.lfc, do_lfc),
            harvest_pointer(res.delta_detected, do_delta_detected)
        );
        return res;
    }
 
private:
    template<typename Ptr>
    static std::vector<Ptr> fetch_first(const std::vector<std::vector<Ptr> >& input) {
        std::vector<Ptr> output;
        output.reserve(input.size());
        for (const auto& i : input) {
            output.push_back(i.front());
        }
        return output;
    }
 
    template<typename Stat> 
    static Stat* harvest_pointer(std::vector<Stat>& source, bool use) {
        if (use) {
            return source.data();
        } else {
            return static_cast<Stat*>(NULL);
        }
    }
};
 
}
 
#endif