PerCellQcMetrics.hpp

Go to the documentation of this file.

#ifndef SCRAN_PER_CELL_QC_METRICS_HPP
#define SCRAN_PER_CELL_QC_METRICS_HPP
 
#include "../utils/macros.hpp"
 
#include <vector>
#include <algorithm>
#include <limits>
#include <cstdint>
 
#include "tatami/tatami.hpp"
#include "../utils/vector_to_pointers.hpp"
 
namespace scran {
 
class PerCellQcMetrics {
public:
    struct Results {
        Results() {}
 
        Results(size_t nsubsets) : subset_total(nsubsets), subset_detected(nsubsets) {}
        std::vector<double> total;
 
        std::vector<int> detected;
 
        std::vector<int> max_index;
 
        std::vector<double> max_count;
 
        std::vector<std::vector<double> > subset_total;
 
        std::vector<std::vector<int> > subset_detected;
    };
 
public:
    template<typename Float = double, typename Integer = int>
    struct Buffers {
        Buffers() {}
 
        Buffers(size_t nsubsets) : subset_total(nsubsets, NULL), subset_detected(nsubsets, NULL) {}
        Float* total = NULL;
 
        Integer* detected = NULL;
 
        Integer* max_index = NULL;
 
        Float* max_count = NULL;
 
        std::vector<Float*> subset_total;
 
        std::vector<Integer*> subset_detected;
 
        bool already_zeroed = false;
    };
 
public:
    struct Defaults {
        static constexpr bool compute_total = true;
 
        static constexpr bool compute_detected = true;
 
        static constexpr bool compute_max_count = true;
 
        static constexpr bool compute_max_index = true;
 
        static constexpr bool compute_subset_total = true;
 
        static constexpr bool compute_subset_detected = true;
 
        static constexpr int num_threads = 1;
    };
 
private:
    bool compute_total = Defaults::compute_total;
    bool compute_detected = Defaults::compute_detected;
    bool compute_max_count = Defaults::compute_max_count;
    bool compute_max_index = Defaults::compute_max_index;
    bool compute_subset_total = Defaults::compute_subset_total;
    bool compute_subset_detected = Defaults::compute_subset_detected;
    int num_threads = Defaults::num_threads;
 
public:
    PerCellQcMetrics& set_compute_total(bool s = Defaults::compute_total) {
        compute_total = s;
        return *this;        
    }
 
    PerCellQcMetrics& set_compute_detected(bool s = Defaults::compute_detected) {
        compute_detected = s;
        return *this;
    }
 
    PerCellQcMetrics& set_compute_max_count(bool s = Defaults::compute_max_count) {
        compute_max_count = s;
        return *this;        
    }
 
    PerCellQcMetrics& set_compute_max_index(bool s = Defaults::compute_max_index) {
        compute_max_index = s;
        return *this;        
    }
 
    PerCellQcMetrics& set_compute_subset_total(bool s = Defaults::compute_subset_total) {
        compute_subset_total = s;
        return *this;
    }
 
    PerCellQcMetrics& set_compute_subset_detected(bool s = Defaults::compute_subset_detected) {
        compute_subset_detected = s;
        return *this;
    }
 
    PerCellQcMetrics& set_num_threads(int n = Defaults::num_threads) {
        num_threads = n;
        return *this;
    }
 
public:
    template<class Matrix, typename Subset = const uint8_t*>
    Results run(const Matrix* mat, const std::vector<Subset>& subsets) const {
        Results output;
        Buffers<> buffers;
        auto ncells = mat->ncol();
 
        if (compute_total) {
            output.total.resize(ncells);
            buffers.total = output.total.data();
        }
        if (compute_detected) {
            output.detected.resize(ncells);
            buffers.detected = output.detected.data();
        }
        if (compute_max_index) {
            output.max_index.resize(ncells);
            buffers.max_index = output.max_index.data();
        }
        if (compute_max_count) {
            output.max_count.resize(ncells, pick_fill_value<double>());
            buffers.max_count = output.max_count.data();
        }
 
        size_t nsubsets = subsets.size();
 
        if (compute_subset_total) {
            output.subset_total.resize(nsubsets);
            buffers.subset_total.resize(nsubsets);
            for (size_t s = 0; s < nsubsets; ++s) {
                output.subset_total[s].resize(ncells);
                buffers.subset_total[s] = output.subset_total[s].data();
            }
        }
 
        if (compute_subset_detected) {
            output.subset_detected.resize(nsubsets);
            buffers.subset_detected.resize(nsubsets);
            for (size_t s = 0; s < nsubsets; ++s) {
                output.subset_detected[s].resize(ncells);
                buffers.subset_detected[s] = output.subset_detected[s].data();
            }
        }
 
        buffers.already_zeroed = true;
        run(mat, subsets, buffers);
        return output;
    }
 
private:
    template<typename T>
    static constexpr T pick_fill_value() {
        if constexpr(std::numeric_limits<T>::has_infinity) {
            return -std::numeric_limits<T>::infinity();
        } else {
            return std::numeric_limits<T>::lowest();
        }
    }
 
private:
    template<typename Data_, typename Index_, typename Subset_, typename Float_, typename Integer_>
    void compute_direct_dense(const tatami::Matrix<Data_, Index_>* mat, const std::vector<Subset_>& subsets, Buffers<Float_, Integer_>& output) const {
        std::vector<std::vector<int> > subset_indices;
        if (!output.subset_total.empty() || !output.subset_detected.empty()) {
            size_t nsubsets = subsets.size();
            subset_indices.resize(nsubsets);
            auto NR = mat->nrow();
 
            for (size_t s = 0; s < nsubsets; ++s) {
                auto& current = subset_indices[s];
                const auto& source = subsets[s];
 
                for (int i = 0; i < NR; ++i) {
                    if (source[i]) {
                        current.push_back(i);
                    }
                }
            }
        }
 
        tatami::parallelize([&](size_t, Index_ start, Index_ length) {
            auto NR = mat->nrow();
            auto ext = tatami::consecutive_extractor<false, false>(mat, start, length);
            std::vector<Data_> vbuffer(NR);
            bool do_max = output.max_index || output.max_count;
 
            for (Index_ c = start, end = start + length; c < end; ++c) {
                auto ptr = ext->fetch(c, vbuffer.data());
 
                if (output.total) {
                    output.total[c] = std::accumulate(ptr, ptr + NR, static_cast<Float_>(0));
                }
 
                if (output.detected) {
                    Integer_ count = 0;
                    for (Index_ r = 0; r < NR; ++r) {
                        count += (ptr[r] != 0);
                    }
                    output.detected[c] = count;
                }
 
                if (do_max) {
                    auto max_count = PerCellQcMetrics::pick_fill_value<Float_>();
                    Integer_ max_index = 0;
                    for (Index_ r = 0; r < NR; ++r) {
                        if (max_count < ptr[r]) {
                            max_count = ptr[r];
                            max_index = r;
                        }
                    }
 
                    if (output.max_index) {
                        output.max_index[c] = max_index;
                    }
                    if (output.max_count) {
                        output.max_count[c] = max_count;
                    }
                }
 
                size_t nsubsets = subset_indices.size();
                for (size_t s = 0; s < nsubsets; ++s) {
                    const auto& sub = subset_indices[s];
 
                    if (!output.subset_total.empty() && output.subset_total[s]) {
                        Float_ current = 0;
                        for (auto r : sub) {
                            current += ptr[r];
                        }
                        output.subset_total[s][c] = current;
                    }
 
                    if (!output.subset_detected.empty() && output.subset_detected[s]) {
                        Integer_ current = 0;
                        for (auto r : sub) {
                            current += ptr[r] != 0;
                        }
                        output.subset_detected[s][c] = current;
                    }
                }
            }
        }, mat->ncol(), num_threads);
    }
 
    template<typename Data_, typename Index_, typename Subset_, typename Float_, typename Integer_>
    void compute_direct_sparse(const tatami::Matrix<Data_, Index_>* mat, const std::vector<Subset_>& subsets, Buffers<Float_, Integer_>& output) const {
        tatami::Options opt;
        opt.sparse_ordered_index = false;
 
        tatami::parallelize([&](size_t, Index_ start, Index_ length) {
            auto NR = mat->nrow();
            auto ext = tatami::consecutive_extractor<false, true>(mat, start, length, opt);
            std::vector<Data_> vbuffer(NR);
            std::vector<Index_> ibuffer(NR);
 
            bool do_max = output.max_index || output.max_count;
            std::vector<unsigned char> internal_is_nonzero(output.max_index ? NR : 0);
 
            for (Index_ c = start, end = start + length; c < end; ++c) {
                auto range = ext->fetch(c, vbuffer.data(), ibuffer.data());
 
                if (output.total) {
                    output.total[c] = std::accumulate(range.value, range.value + range.number, static_cast<Float_>(0));
                }
 
                if (output.detected) {
                    Integer_ current = 0;
                    for (Index_ i = 0; i < range.number; ++i) {
                        current += (range.value[i] != 0);
                    }
                    output.detected[c] = current;
                }
 
                if (do_max) {
                    auto max_count = PerCellQcMetrics::pick_fill_value<Float_>();
                    Integer_ max_index = 0;
                    for (Index_ i = 0; i < range.number; ++i) {
                        if (max_count < range.value[i]) {
                            max_count = range.value[i];
                            max_index = range.index[i];
                        }
                    }
 
                    if (max_count <= 0 && range.number < NR) {
                        // Zero is the max.
                        max_count = 0;
 
                        // Finding the index of the first zero by tracking all
                        // indices with non-zero values. This isn't the fastest
                        // approach but it's simple and avoids assuming that
                        // indices are sorted. Hopefully we don't have to hit
                        // this section often.
                        if (output.max_index) {
                            for (Index_ i = 0; i < range.number; ++i) {
                                if (range.value[i]) {
                                    internal_is_nonzero[range.index[i]] = 1;
                                }
                            }
                            for (Index_ r = 0; r < NR; ++r) {
                                if (internal_is_nonzero[r] == 0) {
                                    max_index = r;
                                    break;
                                }
                            }
                            for (Index_ i = 0; i < range.number; ++i) { // setting back to zero.
                                internal_is_nonzero[range.index[i]] = 0;
                            }
                        }
                    }
 
                    if (output.max_index) {
                        output.max_index[c] = max_index;
                    }
                    if (output.max_count) {
                        output.max_count[c] = max_count;
                    }
                }
 
                size_t nsubsets = subsets.size();
                for (size_t s = 0; s < nsubsets; ++s) {
                    const auto& sub = subsets[s];
 
                    if (!output.subset_total.empty() && output.subset_total[s]) {
                        auto& current = output.subset_total[s][c];
                        for (Index_ i = 0; i < range.number; ++i) {
                            current += (sub[range.index[i]] != 0) * range.value[i];
                        }
                    }
 
                    if (!output.subset_detected.empty() && output.subset_detected[s]) {
                        auto& current = output.subset_detected[s][c];
                        for (Index_ i = 0; i < range.number; ++i) {
                            current += (sub[range.index[i]] != 0) * (range.value[i] != 0);
                        }
                    }
                }
            }
        }, mat->ncol(), num_threads);
    }
 
    template<typename Data_, typename Index_, typename Subset_, typename Float_, typename Integer_>
    void compute_running_dense(const tatami::Matrix<Data_, Index_>* mat, const std::vector<Subset_>& subsets, Buffers<Float_, Integer_>& output) const {
        tatami::parallelize([&](size_t, Index_ start, Index_ len) {
            auto NR = mat->nrow();
            auto ext = tatami::consecutive_extractor<true, false>(mat, 0, NR, start, len);
            std::vector<Data_> vbuffer(len);
            bool do_max = output.max_index || output.max_count;
            std::vector<Float_> internal_max_count(do_max && !output.max_count ? len : 0);
 
            for (Index_ r = 0; r < NR; ++r) {
                auto ptr = ext->fetch(r, vbuffer.data());
 
                if (output.total) {
                    auto outt = output.total + start;
                    for (Index_ i = 0; i < len; ++i) {
                        outt[i] += ptr[i];
                    }
                }
 
                if (output.detected) {
                    auto outd = output.detected + start;
                    for (Index_ i = 0; i < len; ++i) {
                        outd[i] += (ptr[i] != 0);
                    }
                }
 
                if (do_max) {
                    auto outmc = (output.max_count ? output.max_count + start : internal_max_count.data());
 
                    if (r == 0) {
                        std::copy(ptr, ptr + len, outmc);
                        if (output.max_index) {
                            auto outmi = output.max_index + start;
                            std::fill(outmi, outmi + len, 0);
                        }
 
                    } else {
                        for (Index_ i = 0; i < len; ++i) {
                            auto& curmax = outmc[i];
                            if (curmax < ptr[i]) {
                                curmax = ptr[i];
                                if (output.max_index) {
                                    output.max_index[i + start] = r;
                                }
                            }
                        }
                    }
                }
 
                size_t nsubsets = subsets.size();
                for (size_t s = 0; s < nsubsets; ++s) {
                    const auto& sub = subsets[s];
                    if (sub[r] == 0) {
                        continue;
                    }
 
                    if (!output.subset_total.empty() && output.subset_total[s]) {
                        auto current = output.subset_total[s] + start;
                        for (Index_ i = 0; i < len; ++i) {
                            current[i] += ptr[i];
                        }
                    }
 
                    if (!output.subset_detected.empty() && output.subset_detected[s]) {
                        auto current = output.subset_detected[s] + start;
                        for (Index_ i = 0; i < len; ++i) {
                            current[i] += (ptr[i] != 0);
                        }
                    }
                }
            }
        }, mat->ncol(), num_threads);
    }
 
    template<typename Data_, typename Index_, typename Subset_, typename Float_, typename Integer_>
    void compute_running_sparse(const tatami::Matrix<Data_, Index_>* mat, const std::vector<Subset_>& subsets, Buffers<Float_, Integer_>& output) const {
        tatami::Options opt;
        opt.sparse_ordered_index = false;
 
        Index_ NC = mat->ncol();
        bool do_max = output.max_index || output.max_count;
        std::vector<Float_> internal_max_count(do_max && !output.max_count ? NC : 0);
        std::vector<Integer_> internal_last_consecutive_nonzero(do_max ? NC : 0);
 
        tatami::parallelize([&](size_t t, Index_ start, Index_ len) {
            auto NR = mat->nrow();
            auto ext = tatami::consecutive_extractor<true, true>(mat, 0, NR, start, len, opt);
            std::vector<Data_> vbuffer(len);
            std::vector<Index_> ibuffer(len);
 
            for (Index_ r = 0; r < NR; ++r) {
                auto range = ext->fetch(r, vbuffer.data(), ibuffer.data());
 
                if (output.total) {
                    for (Index_ i = 0; i < range.number; ++i) {
                        output.total[range.index[i]] += range.value[i];
                    }
                }
 
                if (output.detected) {
                    for (Index_ i = 0; i < range.number; ++i) {
                        output.detected[range.index[i]] += (range.value[i] != 0);
                    }
                }
 
                if (do_max) {
                    auto outmc = (output.max_count ? output.max_count : internal_max_count.data());
 
                    for (Index_ i = 0; i < range.number; ++i) {
                        auto& curmax = outmc[range.index[i]];
                        if (curmax < range.value[i]) {
                            curmax = range.value[i];
                            if (output.max_index) {
                                output.max_index[range.index[i]] = r;
                            }
                        }
 
                        // Getting the index of the last consecutive non-zero entry, so that
                        // we can check if zero is the max and gets its first occurrence, if necessary.
                        auto& last = internal_last_consecutive_nonzero[range.index[i]];
                        if (static_cast<Index_>(last) == r) {
                            if (range.value[i] != 0) {
                                ++last;
                            }
                        }
                    }
                }
 
                size_t nsubsets = subsets.size();
                for (size_t s = 0; s < nsubsets; ++s) {
                    const auto& sub = subsets[s];
                    if (sub[r] == 0) {
                        continue;
                    }
 
                    if (!output.subset_total.empty() && output.subset_total[s]) {
                        auto& current = output.subset_total[s];
                        for (size_t i = 0; i < range.number; ++i) {
                            current[range.index[i]] += range.value[i];
                        }
                    }
 
                    if (!output.subset_detected.empty() && output.subset_detected[s]) {
                        auto& current = output.subset_detected[s];
                        for (size_t i = 0; i < range.number; ++i) {
                            current[range.index[i]] += (range.value[i] != 0);
                        }
                    }
                }
            }
        }, NC, num_threads);
 
        if (do_max) {
            auto outmc = (output.max_count ? output.max_count : internal_max_count.data());
            auto NR = mat->nrow();
 
            // Checking anything with non-positive maximum, and replacing it with zero
            // if there are any zeros (i.e., consecutive non-zeros is not equal to the number of rows).
            for (Index_ c = 0; c < NC; ++c) {
                auto& current = outmc[c];
                if (current > 0) {
                    continue;
                }
 
                auto last_nz = internal_last_consecutive_nonzero[c];
                if (last_nz == NR) {
                    continue;
                }
 
                current = 0;
                if (output.max_index) {
                    output.max_index[c] = last_nz;
                }
            }
        }
    }
 
public:
    template<class Matrix, typename Subset = const uint8_t*, typename Float, typename Integer>
    void run(const Matrix* mat, const std::vector<Subset>& subsets, Buffers<Float, Integer>& output) const {
        if (!output.already_zeroed) {
            size_t n = mat->ncol();
            auto check_and_fill = [&](auto* ptr, auto value) -> void {
                if (ptr) {
                    std::fill(ptr, ptr + n, value);
                }
            };
 
            check_and_fill(output.total, static_cast<Float>(0));
            check_and_fill(output.detected, static_cast<Integer>(0));
            check_and_fill(output.max_count, pick_fill_value<Float>());
            check_and_fill(output.max_index, static_cast<Integer>(0));
 
            for (size_t s = 0; s < subsets.size(); ++s) {
                if (!output.subset_total.empty()) {
                    check_and_fill(output.subset_total[s], static_cast<Float>(0));
                }
                if (!output.subset_detected.empty()) {
                    check_and_fill(output.subset_detected[s], static_cast<Integer>(0));
                }
            }
        }
 
        if (mat->sparse()) {
            if (mat->prefer_rows()) {
                compute_running_sparse(mat, subsets, output);
            } else {
                compute_direct_sparse(mat, subsets, output);
            }
        } else {
            if (mat->prefer_rows()) {
                compute_running_dense(mat, subsets, output);
            } else {
                compute_direct_dense(mat, subsets, output);
            }
        }
    }
};
 
}
 
#endif