compute.hpp

Go to the documentation of this file.

#ifndef GSDECON_COMPUTE_HPP
#define GSDECON_COMPUTE_HPP
 
#include <algorithm>
#include <numeric>
 
#include "tatami/tatami.hpp"
#include "irlba/irlba.hpp"
#include "scran_pca/scran_pca.hpp"
#include "sanisizer/sanisizer.hpp"
#include "quickstats/quickstats.hpp"
 
#include "Options.hpp"
#include "Results.hpp"
#include "utils.hpp"
 
namespace gsdecon {
 
template<typename Value_, typename Index_, typename Float_>
irlba::Metrics compute(const tatami::Matrix<Value_, Index_>& matrix, 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::SimplePcaOptions sopt;
    sopt.number = options.rank;
    sopt.scale = options.scale;
    sopt.realize_matrix = options.realize_matrix;
    sopt.num_threads = options.num_threads;
    sopt.irlba_options = options.irlba_options;
    const auto res = scran_pca::simple_pca(matrix, sopt);
 
    quickstats::PairwiseSumWorkspace<Float_> pswrk;
    quickstats::PairwiseSumOptions psopt;
    const Float_ shift = quickstats::pairwise_sum(res.center.size(), res.center.data(), pswrk, psopt) / matrix.nrow();
    std::fill_n(output.scores, matrix.ncol(), shift);
 
    process_output(res.rotation, res.components, options.scale, res.scale, output);
 
    return res.metrics;
}
 
template<typename Float_ = double, typename Value_, typename Index_>
Results<Float_> compute(const tatami::Matrix<Value_, Index_>& matrix, 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(matrix, options, buffers);
    return output;
}
 
}
 
#endif