per_cell_qc_metrics.hpp

Go to the documentation of this file.

#ifndef SCRAN_QC_PER_CELL_QC_METRICS_HPP
#define SCRAN_QC_PER_CELL_QC_METRICS_HPP
 
#include <vector>
#include <algorithm>
#include <limits>
#include <cstddef>
#include <optional>
 
#include "tatami/tatami.hpp"
#include "sanisizer/sanisizer.hpp"
 
#include "utils.hpp"
 
namespace scran_qc {
 
struct PerCellQcMetricsOptions {
    bool compute_sum = true;
 
    bool compute_detected = true;
 
    bool compute_max_value = true;
 
    bool compute_max_index = true;
 
    bool compute_subset_sum = true;
 
    bool compute_subset_detected = true;
 
    bool subset_containers_have_indices = true;
 
    int num_threads = 1;
};
 
template<typename Sum_, typename Detected_, typename Value_, typename Index_>
struct PerCellQcMetricsBuffers {
    PerCellQcMetricsBuffers() = default;
 
    PerCellQcMetricsBuffers(const std::size_t nsubsets) : 
        subset_sum(sanisizer::cast<I<decltype(subset_sum.size())> >(nsubsets), NULL),
        subset_detected(sanisizer::cast<I<decltype(subset_detected.size())> >(nsubsets), NULL)
    {}
    Sum_* sum = NULL;
 
    Detected_* detected = NULL;
 
    Value_* max_value = NULL;
 
    Index_* max_index = NULL;
 
    std::vector<Sum_*> subset_sum;
 
    std::vector<Detected_*> subset_detected;
};
 
template<typename Value_, typename Index_, typename Subset_, typename Sum_, typename Detected_>
void per_cell_qc_metrics_direct_dense(
    const tatami::Matrix<Value_, Index_>& mat,
    const std::vector<Subset_>& subsets,
    const PerCellQcMetricsBuffers<Sum_, Detected_, Value_, Index_>& output,
    const PerCellQcMetricsOptions& options
) {
    const auto NR = mat.nrow();
    const bool report_max = output.max_value || output.max_index;
    const auto nsubsets = subsets.size();
    const bool report_subsets = output.subset_sum.size() || output.subset_detected.size();
 
    std::optional<std::vector<std::vector<Index_> > > subset_indices;
    if (report_subsets) {
        [&](){ // use an IIFE for easy control via return.
            if constexpr(!std::is_pointer<Subset_>::value) {
                if (options.subset_containers_have_indices) {
                    return;
                }
            }
 
            subset_indices.emplace(sanisizer::cast<I<decltype(subset_indices->size())> >(nsubsets));
            for (I<decltype(nsubsets)> s = 0; s < nsubsets; ++s) {
                auto& current = (*subset_indices)[s];
                const auto& source = subsets[s];
                for (I<decltype(NR)> i = 0; i < NR; ++i) {
                    if (source[i]) {
                        current.push_back(i);
                    }
                }
            }
        }();
    }
 
    tatami::parallelize([&](const int, const Index_ start, const Index_ length) -> void {
        auto ext = tatami::consecutive_extractor<false>(mat, false, start, length);
        auto vbuffer = tatami::create_container_of_Index_size<std::vector<Value_> >(NR);
 
        for (Index_ c = start, end = start + length; c < end; ++c) {
            auto ptr = ext->fetch(c, vbuffer.data());
 
            if (output.sum) {
                output.sum[c] = std::accumulate(ptr, ptr + NR, static_cast<Sum_>(0));
            }
 
            if (output.detected) {
                Detected_ count = 0;
                for (I<decltype(NR)> r = 0; r < NR; ++r) {
                    count += (ptr[r] != 0);
                }
                output.detected[c] = count;
            }
 
            if (report_max) {
                if (NR) {
                    const auto it = std::max_element(ptr, ptr + NR);
                    if (output.max_value) {
                        output.max_value[c] = *it;
                    }
                    if (output.max_index) {
                        output.max_index[c] = it - ptr;
                    }
                } else {
                    if (output.max_value) {
                        output.max_value[c] = 0;
                    }
                    if (output.max_index) {
                        output.max_index[c] = 0;
                    }
                }
            }
 
            if (report_subsets) {
                for (I<decltype(nsubsets)> s = 0; s < nsubsets; ++s) {
                    if constexpr(!std::is_pointer<Subset_>::value) {
                        // Seeing if the subsets are already indices, in which case we use it directly.
                        // It is assumed that there are no duplicate values here.
                        if (options.subset_containers_have_indices) {
                            const auto& sub = subsets[s];
                            if (!output.subset_sum.empty() && output.subset_sum[s]) {
                                Sum_ current = 0;
                                for (const auto r : sub) {
                                    current += ptr[r];
                                }
                                output.subset_sum[s][c] = current;
                            }
                            if (!output.subset_detected.empty() && output.subset_detected[s]) {
                                Detected_ current = 0;
                                for (const auto r : sub) {
                                    current += ptr[r] != 0;
                                }
                                output.subset_detected[s][c] = current;
                            }
                            continue;
                        }
                    }
 
                    const auto& sub = (*subset_indices)[s];
                    if (!output.subset_sum.empty() && output.subset_sum[s]) {
                        Sum_ current = 0;
                        for (const auto r : sub) {
                            current += ptr[r];
                        }
                        output.subset_sum[s][c] = current;
                    }
                    if (!output.subset_detected.empty() && output.subset_detected[s]) {
                        Detected_ current = 0;
                        for (const auto r : sub) {
                            current += ptr[r] != 0;
                        }
                        output.subset_detected[s][c] = current;
                    }
                }
            }
        }
    }, mat.ncol(), options.num_threads);
}
 
template<typename Index_, typename Subset_>
std::vector<std::vector<unsigned char> > boolify_subsets(const Index_ NR, const std::vector<Subset_>& subsets) {
    const auto nsubsets = subsets.size();
    auto output = sanisizer::create<std::vector<std::vector<unsigned char> > >(nsubsets);
    for (I<decltype(nsubsets)> s = 0; s < nsubsets; ++s) {
        auto& current = output[s];
        tatami::resize_container_to_Index_size(current, NR);
        for (const auto i : subsets[s]) {
            current[i] = 1;
        }
    }
    return output;
}
 
template<typename Value_, typename Index_, typename Subset_, typename Sum_, typename Detected_>
void per_cell_qc_metrics_direct_sparse(
    const tatami::Matrix<Value_, Index_>& mat,
    const std::vector<Subset_>& subsets,
    const PerCellQcMetricsBuffers<Sum_, Detected_, Value_, Index_>& output,
    const PerCellQcMetricsOptions& options
) {
    const auto NR = mat.nrow();
    const bool report_max = output.max_value || output.max_index;
    const auto nsubsets = subsets.size();
    const bool report_subsets = output.subset_sum.size() || output.subset_detected.size();
 
    std::optional<std::vector<std::vector<unsigned char> > > is_in_subset;
    if (report_subsets) {
        if constexpr(!std::is_pointer<Subset_>::value) {
            if (options.subset_containers_have_indices) {
                is_in_subset = boolify_subsets(NR, subsets);
            }
        }
    }
 
    tatami::parallelize([&](const int, const Index_ start, const Index_ length) -> void {
        auto ext = tatami::consecutive_extractor<true>(mat, false, start, length);
        auto vbuffer = tatami::create_container_of_Index_size<std::vector<Value_> >(NR);
        auto ibuffer = tatami::create_container_of_Index_size<std::vector<Index_> >(NR);
 
        for (Index_ c = start, end = start + length; c < end; ++c) {
            auto range = ext->fetch(vbuffer.data(), ibuffer.data());
 
            if (output.sum) {
                output.sum[c] = std::accumulate(range.value, range.value + range.number, static_cast<Sum_>(0));
            }
 
            if (output.detected) {
                Detected_ current = 0;
                for (Index_ i = 0; i < range.number; ++i) {
                    current += (range.value[i] != 0);
                }
                output.detected[c] = current;
            }
 
            if (report_max) {
                if (range.number) {
                    const auto it = std::max_element(range.value, range.value + range.number);
                    if (*it > 0 || range.number == NR) {
                        if (output.max_value) {
                            output.max_value[c] = *it;
                        }
                        if (output.max_index) {
                            output.max_index[c] = range.index[it - range.value];
                        }
                    } else {
                        if (output.max_value) {
                            output.max_value[c] = 0;
                        }
                        if (output.max_index) {
                            if (*it < 0) {
                                // Find the first structural zero. 
                                output.max_index[c] = range.number;
                                for (Index_ i = 0; i < range.number; ++i) {
                                    if (range.index[i] != i) { 
                                        output.max_index[c] = i;
                                        break;
                                    }
                                }
                            } else {
                                // Find the first structural zero that occurs before the structural non-zero with a value of zero.
                                const Index_ candidate = it - range.value; 
                                output.max_index[c] = range.index[candidate];
                                for (Index_ i = 0; i <= candidate; ++i) {
                                    if (range.index[i] != i) {
                                        output.max_index[c] = i;
                                        break;
                                    }
                                }
                            }
                        }
                    }
                } else {
                    if (output.max_value) {
                        output.max_value[c] = 0;
                    }
                    if (output.max_index) {
                        output.max_index[c] = 0;
                    } 
                }
            }
 
            if (report_subsets) {
                for (I<decltype(nsubsets)> s = 0; s < nsubsets; ++s) {
                    if constexpr(!std::is_pointer<Subset_>::value) {
                        if (options.subset_containers_have_indices) {
                            const auto& sub = (*is_in_subset)[s];
                            if (!output.subset_sum.empty() && output.subset_sum[s]) {
                                Sum_ current = 0;
                                for (Index_ i = 0; i < range.number; ++i) {
                                    current += (sub[range.index[i]] != 0) * range.value[i];
                                }
                                output.subset_sum[s][c] = current;
                            }
                            if (!output.subset_detected.empty() && output.subset_detected[s]) {
                                Detected_ current = 0;
                                for (Index_ i = 0; i < range.number; ++i) {
                                    current += (sub[range.index[i]] != 0) * (range.value[i] != 0);
                                }
                                output.subset_detected[s][c] = current;
                            }
                            continue;
                        }
                    }
 
                    const auto& sub = subsets[s];
                    if (!output.subset_sum.empty() && output.subset_sum[s]) {
                        Sum_ current = 0;
                        for (Index_ i = 0; i < range.number; ++i) {
                            current += (sub[range.index[i]] != 0) * range.value[i];
                        }
                        output.subset_sum[s][c] = current;
                    }
                    if (!output.subset_detected.empty() && output.subset_detected[s]) {
                        Detected_ current = 0;
                        for (Index_ i = 0; i < range.number; ++i) {
                            current += (sub[range.index[i]] != 0) * (range.value[i] != 0);
                        }
                        output.subset_detected[s][c] = current;
                    }
                }
            }
        }
    }, mat.ncol(), options.num_threads);
}
 
template<typename Value_, typename Index_, typename Subset_, typename Sum_, typename Detected_>
void per_cell_qc_metrics_running(
    const tatami::Matrix<Value_, Index_>& mat,
    const std::vector<Subset_>& subsets,
    const PerCellQcMetricsBuffers<Sum_, Detected_, Value_, Index_>& output,
    const PerCellQcMetricsOptions& options
) {
    const auto NR = mat.nrow();
    const auto NC = mat.ncol();
    const bool is_sparse = mat.is_sparse();
 
    /************************************
     *** Setting up result containers ***
     ************************************/
 
    const bool report_max = output.max_value || output.max_index;
    std::optional<std::vector<Value_> > tmp_max_value;
    Value_* max_value_output_ptr;
    if (report_max) {
        if (!output.max_value) {
            tmp_max_value.emplace(tatami::cast_Index_to_container_size<std::vector<Index_> >(NC));
            max_value_output_ptr = tmp_max_value->data();
        } else {
            max_value_output_ptr = output.max_value;
        }
    }
 
    const bool do_parallel = options.num_threads > 1;
    std::optional<std::vector<std::optional<std::vector<Sum_> > > > partial_sum;
    std::optional<std::vector<std::optional<std::vector<Detected_> > > > partial_detected;
    std::optional<std::vector<std::optional<std::vector<Value_> > > > partial_max_value;
    std::optional<std::vector<std::optional<std::vector<Index_> > > > partial_max_index;
    std::optional<std::vector<std::optional<std::vector<std::vector<Sum_> > > > > partial_subset_sum;
    std::optional<std::vector<std::optional<std::vector<std::vector<Detected_> > > > > partial_subset_detected;
    if (do_parallel) {
        if (output.sum) {
            partial_sum.emplace(tatami::cast_Index_to_container_size<I<decltype(*partial_sum)> >(options.num_threads - 1));
        }
        if (output.detected) {
            partial_detected.emplace(tatami::cast_Index_to_container_size<I<decltype(*partial_detected)> >(options.num_threads - 1));
        }
        if (report_max) {
            partial_max_value.emplace(tatami::cast_Index_to_container_size<I<decltype(*partial_max_value)> >(options.num_threads - 1));
            if (output.max_index) {
                partial_max_index.emplace(tatami::cast_Index_to_container_size<I<decltype(*partial_max_index)> >(options.num_threads - 1));
            }
        }
        if (output.subset_sum.size()) {
            partial_subset_sum.emplace(tatami::cast_Index_to_container_size<I<decltype(*partial_subset_sum)> >(options.num_threads - 1));
        }
        if (output.subset_detected.size()) {
            partial_subset_detected.emplace(tatami::cast_Index_to_container_size<I<decltype(*partial_subset_detected)> >(options.num_threads - 1));
        }
    }
 
    // Zeroing the output arrays. 
    if (output.sum) {
        std::fill_n(output.sum, NC, 0);
    }
    if (output.detected) {
        std::fill_n(output.detected, NC, 0);
    }
    if (report_max && NR == 0) { // no need to zero if it's not empty, as it'll get filled by thread 0 upon encountering the first column.
        std::fill_n(output.max_value, NC, 0);
        if (output.max_index) {
            std::fill_n(output.max_index, NC, 0);
        }
    }
    for (const auto sptr : output.subset_sum) {
        if (sptr) {
            std::fill_n(sptr, NC, 0);
        }
    }
    for (const auto sptr : output.subset_detected) {
        if (sptr) {
            std::fill_n(sptr, NC, 0);
        }
    }
 
    const auto nsubsets = subsets.size();
    const bool report_subsets = output.subset_sum.size() || output.subset_detected.size();
    std::optional<std::vector<std::vector<unsigned char> > > is_in_subset;
    if (report_subsets) {
        if constexpr(!std::is_pointer<Subset_>::value) {
            if (options.subset_containers_have_indices) {
                is_in_subset = boolify_subsets(NR, subsets);
            }
        }
    }
 
    const auto num_used = tatami::parallelize([&](int thread, Index_ start, Index_ len) -> void {
        /*********************************************************
         *** Thread-local containers to mitigate false sharing ***
         *********************************************************/
 
        Sum_* sum_ptr = NULL;
        std::optional<std::vector<Sum_> > sum_buffer;
        Detected_* detected_ptr = NULL;
        std::optional<std::vector<Detected_> > detected_buffer;
        Value_* max_value_ptr = NULL;
        std::optional<std::vector<Value_> > max_value_buffer;
        Index_* max_index_ptr = NULL;
        std::optional<std::vector<Index_> > max_index_buffer;
 
        Sum_* const * subset_sum_ptr = NULL;
        std::optional<std::vector<Sum_*> > subset_sum_ptrs;
        std::optional<std::vector<std::vector<Sum_> > > subset_sum_buffers;
        Detected_* const * subset_detected_ptr = NULL;
        std::optional<std::vector<Detected_*> > subset_detected_ptrs;
        std::optional<std::vector<std::vector<Detected_> > > subset_detected_buffers;
 
        if (!do_parallel || thread == 0) {
            sum_ptr = output.sum;
            detected_ptr = output.detected;
            max_value_ptr = max_value_output_ptr;
            max_index_ptr = output.max_index;
            if (output.subset_sum.size()) {
                subset_sum_ptr = output.subset_sum.data();
            }
            if (output.subset_detected.size()) {
                subset_detected_ptr = output.subset_detected.data();
            }
 
        } else {
            if (output.sum) {
                sum_buffer.emplace(tatami::cast_Index_to_container_size<I<decltype(*sum_buffer)> >(NC));
                sum_ptr = sum_buffer->data();
            }
            if (output.detected) {
                detected_buffer.emplace(tatami::cast_Index_to_container_size<I<decltype(*detected_buffer)> >(NC));
                detected_ptr = detected_buffer->data();
            }
            if (report_max) {
                max_value_buffer.emplace(tatami::cast_Index_to_container_size<I<decltype(*max_value_buffer)> >(NC));
                max_value_ptr = max_value_buffer->data();
                if (output.max_index) {
                    max_index_buffer.emplace(tatami::cast_Index_to_container_size<I<decltype(*max_index_buffer)> >(NC));
                    max_index_ptr = max_index_buffer->data();
                }
            }
            if (output.subset_sum.size()) {
                subset_sum_ptrs.emplace(tatami::cast_Index_to_container_size<I<decltype(*subset_sum_ptrs)> >(nsubsets));
                subset_sum_buffers.emplace(tatami::cast_Index_to_container_size<I<decltype(*subset_sum_buffers)> >(nsubsets));
                for (I<decltype(nsubsets)> s = 0; s < nsubsets; ++s) {
                    if (output.subset_sum[s]) {
                        tatami::resize_container_to_Index_size((*subset_sum_buffers)[s], NC); 
                        (*subset_sum_ptrs)[s] = (*subset_sum_buffers)[s].data();
                    } else {
                        (*subset_sum_ptrs)[s] = NULL;
                    }
 
                }
            }
            if (output.subset_detected.size()) {
                subset_detected_ptrs.emplace(tatami::cast_Index_to_container_size<I<decltype(*subset_detected_ptrs)> >(nsubsets));
                subset_detected_buffers.emplace(tatami::cast_Index_to_container_size<I<decltype(*subset_detected_buffers)> >(nsubsets));
                for (I<decltype(nsubsets)> s = 0; s < nsubsets; ++s) {
                    if (output.subset_detected[s]) {
                        tatami::resize_container_to_Index_size((*subset_detected_buffers)[s], NC); 
                        (*subset_detected_ptrs)[s] = (*subset_detected_buffers)[s].data();
                    } else {
                        (*subset_detected_ptrs)[s] = NULL;
                    }
                }
            }
        }
 
        if (is_sparse) {
            /*****************************
             *** Sparse loop over rows ***
             *****************************/
 
            tatami::Options opt;
            opt.sparse_ordered_index = false;
            auto ext = tatami::consecutive_extractor<true>(mat, true, start, len, opt);
            auto vbuffer = tatami::create_container_of_Index_size<std::vector<Value_> >(NC);
            auto ibuffer = tatami::create_container_of_Index_size<std::vector<Index_> >(NC);
 
            // nonzeros_at_start contains the number of consecutive non-zero elements at the start, i.e., from r = 0.
            std::optional<std::vector<Index_> > nonzeros_at_start;
            if (max_value_ptr) {
                nonzeros_at_start.emplace(tatami::cast_Index_to_container_size<std::vector<Index_> >(NC));
            }
 
            for (Index_ r = 0; r < len; ++r) {
                auto range = ext->fetch(vbuffer.data(), ibuffer.data());
                if (sum_ptr) {
                    for (Index_ i = 0; i < range.number; ++i) {
                        sum_ptr[range.index[i]] += range.value[i];
                    }
                }
                if (detected_ptr) {
                    for (Index_ i = 0; i < range.number; ++i) {
                        detected_ptr[range.index[i]] += (range.value[i] != 0);
                    }
                }
 
                if (report_max) {
                    if (r == 0) {
                        std::fill_n(max_value_ptr, NC, 0);
                        for (Index_ i = 0; i < range.number; ++i) {
                            const auto j = range.index[i];
                            max_value_ptr[j] = range.value[i];
                            (*nonzeros_at_start)[j] = 1; 
                        }
                        if (max_index_ptr) {
                            std::fill_n(max_index_ptr, NC, start);
                        }
 
                    } else {
                        for (Index_ i = 0; i < range.number; ++i) {
                            const auto val = range.value[i];
                            const auto j = range.index[i]; 
                            auto& curmax = max_value_ptr[j];
                            if (curmax < val) {
                                curmax = val;
                                if (max_index_ptr) {
                                    max_index_ptr[j] = start + r;
                                }
                            }
                            auto& last = (*nonzeros_at_start)[j];
                            if (last == r) {
                                // If we don't have an unbroken run of structural non-zeros from the start, we stop incrementing.
                                ++last;
                            }
                        }
                    }
                }
 
                if (report_subsets) {
                    for (I<decltype(nsubsets)> s = 0; s < nsubsets; ++s) {
                        const auto in_subset = [&]() -> bool {
                            if constexpr(!std::is_pointer<Subset_>::value) {
                                if (options.subset_containers_have_indices) {
                                    return (*is_in_subset)[s][r];
                                }
                            }
                            return subsets[s][r];
                        }();
                        if (!in_subset) {
                            continue;
                        }
 
                        if (subset_sum_ptr && subset_sum_ptr[s]) {
                            const auto current = subset_sum_ptr[s];
                            for (Index_ i = 0; i < range.number; ++i) {
                                current[range.index[i]] += range.value[i];
                            }
                        }
 
                        if (subset_detected_ptr && subset_detected_ptr[s]) {
                            const auto current = subset_detected_ptr[s];
                            for (Index_ i = 0; i < range.number; ++i) {
                                current[range.index[i]] += (range.value[i] != 0);
                            }
                        }
                    }
                }
            }
 
            if (report_max) {
                // Checking anything with non-positive maximum, and replacing it with zero if there are any structural zeros.
                for (Index_ c = 0; c < NC; ++c) {
                    const auto last_nz = (*nonzeros_at_start)[c];
                    if (last_nz == len) { // i.e., no structural zeros.
                        continue;
                    }
                    auto& current = max_value_ptr[c];
                    if (current > 0) { // doesn't defeat the current maximum.
                        continue;
                    }
 
                    if (current < 0) {
                        current = 0;
                        if (max_index_ptr) {
                            max_index_ptr[c] = start + last_nz;
                        }
                    } else {
                        // Sometimes, structural non-zeros have a value of zero.
                        // This check ensures that we return the first occurrence of any zero.
                        if (max_index_ptr) {
                            const Index_ first_structural_zero = start + last_nz;
                            if (first_structural_zero < max_index_ptr[c]) {
                                max_index_ptr[c] = first_structural_zero;
                            }
                        }
                    }
                }
            }
 
        } else {
            /****************************
             *** Dense loop over rows ***
             ****************************/
 
            auto ext = tatami::consecutive_extractor<false>(mat, true, start, len);
            auto vbuffer = tatami::create_container_of_Index_size<std::vector<Value_> >(NC);
            for (Index_ r = 0; r < len; ++r) {
                auto ptr = ext->fetch(vbuffer.data());
 
                if (sum_ptr) {
                    for (Index_ i = 0; i < NC; ++i) {
                        sum_ptr[i] += ptr[i];
                    }
                }
 
                if (detected_ptr) {
                    for (Index_ i = 0; i < NC; ++i) {
                        detected_ptr[i] += (ptr[i] != 0);
                    }
                }
 
                if (report_max) {
                    if (r == 0) {
                        std::copy_n(ptr, NC, max_value_ptr);
                        if (max_index_ptr) {
                            std::fill_n(max_index_ptr, NC, start);
                        }
                    } else {
                        for (Index_ i = 0; i < NC; ++i) {
                            auto& curmax = max_value_ptr[i];
                            if (curmax < ptr[i]) {
                                curmax = ptr[i];
                                if (max_index_ptr) {
                                    max_index_ptr[i] = start + r;
                                }
                            }
                        }
                    }
                }
 
                if (report_subsets) {
                    for (I<decltype(nsubsets)> s = 0; s < nsubsets; ++s) {
                        const auto in_subset = [&]() -> bool {
                            if constexpr(!std::is_pointer<Subset_>::value) {
                                if (options.subset_containers_have_indices) {
                                    return (*is_in_subset)[s][r];
                                }
                            }
                            return subsets[s][r];
                        }();
                        if (!in_subset) {
                            continue;
                        }
 
                        if (subset_sum_ptr && subset_sum_ptr[s]) {
                            const auto current = subset_sum_ptr[s];
                            for (Index_ i = 0; i < NC; ++i) {
                                current[i] += ptr[i];
                            }
                        }
 
                        if (subset_detected_ptr && subset_detected_ptr[s]) {
                            const auto current = subset_detected_ptr[s];
                            for (Index_ i = 0; i < NC; ++i) {
                                current[i] += (ptr[i] != 0);
                            }
                        }
                    }
                }
            }
        }
 
        /********************************************
         *** Migrate results to serial containers ***
         ********************************************/
 
        if (do_parallel) {
            if (thread > 0) {
                if (output.sum) {
                    (*partial_sum)[thread - 1] = std::move(sum_buffer);
                }
                if (output.detected) {
                    (*partial_detected)[thread - 1] = std::move(detected_buffer);
                }
                if (report_max) {
                    (*partial_max_value)[thread - 1] = std::move(max_value_buffer);
                    if (output.max_index) {
                        (*partial_max_index)[thread - 1] = std::move(max_index_buffer);
                    }
                }
                if (output.subset_sum.size()) {
                    (*partial_subset_sum)[thread - 1] = std::move(subset_sum_buffers);
                }
                if (output.subset_detected.size()) {
                    (*partial_subset_detected)[thread - 1] = std::move(subset_detected_buffers);
                }
            }
        }
    }, NR, options.num_threads);
 
    /************************************
     *** Reduction into output arrays ***
     ************************************/
 
    if (do_parallel) {
        if (output.sum) {
            for (int u = 1; u < num_used; ++u) {
                const auto& cursum = *((*partial_sum)[u - 1]);
                for (Index_ c = 0; c < NC; ++c) {
                    output.sum[c] += cursum[c];
                }
            }
        }
 
        if (output.detected) {
            for (int u = 1; u < num_used; ++u) {
                const auto& curdetected = *((*partial_detected)[u - 1]);
                for (Index_ c = 0; c < NC; ++c) {
                    output.detected[c] += curdetected[c];
                }
            }
        }
 
        // All used threads will have processed non-empty ranges, so we don't need to worry about the validity of the maxima from each thread.
        if (output.max_index) {
            for (int u = 1; u < num_used; ++u) {
                const auto& curmaxval = *((*partial_max_value)[u - 1]);
                const auto& curmaxidx = *((*partial_max_index)[u - 1]);
                for (Index_ c = 0; c < NC; ++c) {
                    if (curmaxval[c] > max_value_output_ptr[c]) {
                        max_value_output_ptr[c] = curmaxval[c];
                        output.max_index[c] = curmaxidx[c];
                    }
                }
            }
        } else if (output.max_value) {
            for (int u = 1; u < num_used; ++u) {
                const auto& curmaxval = *((*partial_max_value)[u - 1]);
                for (Index_ c = 0; c < NC; ++c) {
                    if (curmaxval[c] > output.max_value[c]) {
                        output.max_value[c] = curmaxval[c];
                    }
                }
            }
        }
 
        if (output.subset_sum.size()) {
            for (I<decltype(nsubsets)> s = 0; s < nsubsets; ++s) {
                const auto outptr = output.subset_sum[s];
                if (outptr == NULL) {
                    continue;
                }
                for (int u = 1; u < num_used; ++u) {
                    const auto& cursubset = (*((*partial_subset_sum)[u - 1]))[s];
                    for (Index_ c = 0; c < NC; ++c) {
                        outptr[c] += cursubset[c];
                    }
                }
            }
        }
 
        if (output.subset_detected.size()) {
            for (I<decltype(nsubsets)> s = 0; s < nsubsets; ++s) {
                const auto outptr = output.subset_detected[s];
                if (outptr == NULL) {
                    continue;
                }
                for (int u = 1; u < num_used; ++u) {
                    const auto& cursubset = (*((*partial_subset_detected)[u - 1]))[s];
                    for (Index_ c = 0; c < NC; ++c) {
                        outptr[c] += cursubset[c];
                    }
                }
            }
        }
    }
}
template<typename Sum_, typename Detected_, typename Value_, typename Index_>
struct PerCellQcMetricsResults {
    PerCellQcMetricsResults() = default;
 
    PerCellQcMetricsResults(const std::size_t nsubsets) : 
        subset_sum(sanisizer::cast<I<decltype(subset_sum.size())> >(nsubsets)),
        subset_detected(sanisizer::cast<I<decltype(subset_detected.size())> >(nsubsets))
    {}
    std::vector<Sum_> sum;
 
    std::vector<Detected_> detected;
 
    std::vector<Value_> max_value;
 
    std::vector<Index_> max_index;
 
    std::vector<std::vector<Sum_> > subset_sum;
 
    std::vector<std::vector<Detected_> > subset_detected;
};
 
template<typename Value_, typename Index_, typename Subset_, typename Sum_, typename Detected_>
void per_cell_qc_metrics(
    const tatami::Matrix<Value_, Index_>& mat,
    const std::vector<Subset_>& subsets, 
    const PerCellQcMetricsBuffers<Sum_, Detected_, Value_, Index_>& output,
    const PerCellQcMetricsOptions& options)
{
    if (mat.prefer_rows()) {
        per_cell_qc_metrics_running(mat, subsets, output, options);
    } else {
        if (mat.sparse()) {
            per_cell_qc_metrics_direct_sparse(mat, subsets, output, options);
        } else {
            per_cell_qc_metrics_direct_dense(mat, subsets, output, options);
        }
    }
}
 
template<typename Sum_ = double, typename Detected_ = int, typename Value_, typename Index_, typename Subset_>
PerCellQcMetricsResults<Sum_, Detected_, Value_, Index_> per_cell_qc_metrics(
    const tatami::Matrix<Value_, Index_>& mat,
    const std::vector<Subset_>& subsets,
    const PerCellQcMetricsOptions& options)
{
    PerCellQcMetricsResults<Sum_, Detected_, Value_, Index_> output;
    PerCellQcMetricsBuffers<Sum_, Detected_, Value_, Index_> buffers;
    const auto ncells = mat.ncol();
 
    if (options.compute_sum) {
        tatami::resize_container_to_Index_size(output.sum, ncells
#ifdef SCRAN_QC_TEST_INIT
            , SCRAN_QC_TEST_INIT
#endif
        );
        buffers.sum = output.sum.data();
    }
 
    if (options.compute_detected) {
        tatami::resize_container_to_Index_size(output.detected, ncells
#ifdef SCRAN_QC_TEST_INIT
            , SCRAN_QC_TEST_INIT
#endif
        );
        buffers.detected = output.detected.data();
    }
 
    if (options.compute_max_value) {
        tatami::resize_container_to_Index_size(output.max_value, ncells
#ifdef SCRAN_QC_TEST_INIT
            , SCRAN_QC_TEST_INIT
#endif
        );
        buffers.max_value = output.max_value.data();
    }
    if (options.compute_max_index) {
        tatami::resize_container_to_Index_size(output.max_index, ncells
#ifdef SCRAN_QC_TEST_INIT
            , SCRAN_QC_TEST_INIT
#endif
        );
        buffers.max_index = output.max_index.data();
    }
 
    const auto nsubsets = subsets.size();
 
    if (options.compute_subset_sum) {
        sanisizer::resize(output.subset_sum, nsubsets);
        sanisizer::resize(buffers.subset_sum, nsubsets);
        for (I<decltype(nsubsets)> s = 0; s < nsubsets; ++s) {
            tatami::resize_container_to_Index_size(output.subset_sum[s], ncells
#ifdef SCRAN_QC_TEST_INIT
                , SCRAN_QC_TEST_INIT
#endif
            );
            buffers.subset_sum[s] = output.subset_sum[s].data();
        }
    }
 
    if (options.compute_subset_detected) {
        sanisizer::resize(output.subset_detected, nsubsets);
        sanisizer::resize(buffers.subset_detected, nsubsets);
        for (I<decltype(nsubsets)> s = 0; s < nsubsets; ++s) {
            tatami::resize_container_to_Index_size(output.subset_detected[s], ncells
#ifdef SCRAN_QC_TEST_INIT
                , SCRAN_QC_TEST_INIT
#endif
            );
            buffers.subset_detected[s] = output.subset_detected[s].data();
        }
    }
 
    per_cell_qc_metrics(mat, subsets, buffers, options);
    return output;
}
 
}
 
#endif