blocked.hpp

Go to the documentation of this file.

#ifndef GSDECON_BLOCKED_HPP
#define GSDECON_BLOCKED_HPP
 
#include <vector>
 
#include "Eigen/Dense"
#include "tatami/tatami.hpp"
#include "irlba/irlba.hpp"
#include "scran_pca/scran_pca.hpp"
#include "sanisizer/sanisizer.hpp"
 
#include "Options.hpp"
#include "Results.hpp"
#include "utils.hpp"
 
namespace gsdecon {
 
template<typename Value_, typename Index_, typename Block_, typename Float_>
irlba::Metrics compute_blocked(
    const tatami::Matrix<Value_, Index_>& matrix,
    const Block_* const block,
    const std::size_t num_blocks,
    const Options& options,
    const Buffers<Float_>& output
) {
    if (check_edge_cases(matrix, options.rank, output)) {
        irlba::Metrics metrics;
        metrics.converged = true;
        return metrics;
    }
 
    scran_pca::BlockedPcaOptions bopt;
    bopt.number = options.rank;
    bopt.scale = options.scale;
    bopt.block_weight_policy = options.block_weight_policy;
    bopt.variable_block_weight_parameters = options.variable_block_weight_parameters;
    bopt.realize_matrix = options.realize_matrix;
    bopt.num_threads = options.num_threads;
    bopt.irlba_options = options.irlba_options;
    bopt.center_scores_by_block = true; // we'll add them back in afterwards.
    const auto res = scran_pca::blocked_pca(matrix, block, num_blocks, bopt);
 
    // Here, we restore the block-specific centers.
    static_assert(!Eigen::MatrixXd::IsRowMajor); // just double-checking...
    const auto nfeat = res.center.cols();
    const auto nblocks = res.center.rows();
    auto block_means = sanisizer::create<std::vector<Float_> >(nblocks);
 
    for (I<decltype(nfeat)> f = 0; f < nfeat; ++f) {
        for (I<decltype(nblocks)> b = 0; b < nblocks; ++b) {
            block_means[b] += res.center.coeff(b, f);
        }
    }
    for (auto& b : block_means) {
        b /= nfeat;
    }
 
    const auto ncells = res.components.cols();
    for (I<decltype(ncells)> c = 0; c < ncells; ++c) {
        output.scores[c] = block_means[block[c]];
    }
    process_output(res.rotation, res.components, options.scale, res.scale, output);
 
    return res.metrics;
}
 
template<typename Float_ = double, typename Value_, typename Index_, typename Block_>
Results<Float_> compute_blocked(
    const tatami::Matrix<Value_, Index_>& matrix,
    const Block_* const block,
    const std::size_t num_blocks,
    const Options& options
) {
    Results<Float_> output;
    sanisizer::resize(output.weights, matrix.nrow()
#ifdef SCRAN_QC_TEST_INIT
        , SCRAN_QC_TEST_INIT
#endif
    );
    sanisizer::resize(output.scores, matrix.ncol()
#ifdef SCRAN_QC_TEST_INIT
        , SCRAN_QC_TEST_INIT
#endif
    );
 
    Buffers<Float_> buffers;
    buffers.weights = output.weights.data();
    buffers.scores = output.scores.data();
 
    output.metrics = compute_blocked(matrix, block, num_blocks, options, buffers);
    return output;
}
 
}
 
#endif