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 "Options.hpp"
#include "Results.hpp"
#include "utils.hpp"
 
namespace gsdecon {
 
template<typename Value_, typename Index_, typename Block_, typename Float_>
void compute_blocked(const tatami::Matrix<Value_, Index_>& matrix, const Block_* block, const Options& options, const Buffers<Float_>& output) {
    if (internal::check_edge_cases(matrix, options.rank, output)) {
        return;
    }
 
    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;
    auto res = scran_pca::blocked_pca(matrix, block, bopt);
 
    // Here, we restore the block-specific centers. Don't be tempted into
    // using MultiBatchPca, as that doesn't yield a rank-1 approximation
    // that preserves global shifts between blocks.
    static_assert(!Eigen::MatrixXd::IsRowMajor); // just double-checking...
    const double* cptr = res.center.data();
    size_t nfeat = res.center.cols();
    size_t nblocks = res.center.rows();
    std::vector<double> block_means(nblocks);
 
    for (size_t f = 0; f < nfeat; ++f, cptr += nblocks) {
#ifdef _OPENMP
        #pragma omp simd
#endif
        for (size_t b = 0; b < nblocks; ++b) {
            block_means[b] += cptr[b];
        }
    }
    double denom = nfeat;
    for (auto& b : block_means) {
        b /= denom;
    }
 
    size_t ncells = res.components.cols();
#ifdef _OPENMP
    #pragma omp simd
#endif
    for (size_t c = 0; c < ncells; ++c) {
        output.scores[c] = block_means[block[c]];
    }
    internal::process_output(res.rotation, res.components, options.scale, res.scale, output);
}
 
template<typename Float_ = double, typename Value_, typename Index_, typename Block_>
Results<Float_> compute_blocked(const tatami::Matrix<Value_, Index_>& matrix, const Block_* block, const Options& options) {
    Results<Float_> output;
    output.weights.resize(matrix.nrow());
    output.scores.resize(matrix.ncol());
 
    Buffers<Float_> buffers;
    buffers.weights = output.weights.data();
    buffers.scores = output.scores.data();
 
    compute_blocked(matrix, block, options, buffers);
    return output;
}
 
}
 
#endif