simple_pca.hpp

Go to the documentation of this file.

#ifndef SCRAN_PCA_SIMPLE_PCA_HPP
#define SCRAN_PCA_SIMPLE_PCA_HPP
 
#include <vector>
#include <type_traits>
#include <algorithm>
 
#include "tatami/tatami.hpp"
#include "tatami_stats/tatami_stats.hpp"
#include "irlba/irlba.hpp"
#include "irlba/parallel.hpp"
#include "Eigen/Dense"
#include "sanisizer/sanisizer.hpp"
 
#include "utils.hpp"
 
namespace scran_pca {
 
struct SimplePcaOptions {
    SimplePcaOptions() {
        // Avoid throwing an error if too many PCs are requested.
        irlba_options.cap_number = true;
    }
    int number = 25;
 
    bool scale = false;
 
    bool transpose = true;
 
    bool realize_matrix = true;
 
    int num_threads = 1;
 
    irlba::Options irlba_options;
};
 
namespace internal {
 
template<bool sparse_, typename Value_, typename Index_, class EigenVector_>
void compute_row_means_and_variances(const tatami::Matrix<Value_, Index_>& mat, const int num_threads, EigenVector_& center_v, EigenVector_& scale_v) {
    const auto ngenes = mat.nrow();
 
    if (mat.prefer_rows()) {
        tatami::parallelize([&](const int, const Index_ start, const Index_ length) -> void {
            auto ext = tatami::consecutive_extractor<sparse_>(mat, true, start, length, [&]{
                tatami::Options opt;
                opt.sparse_extract_index = false;
                return opt;
            }());
            const auto ncells = mat.ncol();
            auto vbuffer = tatami::create_container_of_Index_size<std::vector<Value_> >(ncells);
 
            for (Index_ g = start, end = start + length; g < end; ++g) {
                const auto results = [&]{
                    if constexpr(sparse_) {
                        auto range = ext->fetch(vbuffer.data(), NULL);
                        return tatami_stats::variances::direct(range.value, range.number, ncells, /* skip_nan = */ false);
                    } else {
                        auto ptr = ext->fetch(vbuffer.data());
                        return tatami_stats::variances::direct(ptr, ncells, /* skip_nan = */ false);
                    }
                }();
                center_v.coeffRef(g) = results.first;
                scale_v.coeffRef(g) = results.second;
            }
        }, ngenes, num_threads);
 
    } else {
        tatami::parallelize([&](int t, Index_ start, Index_ length) -> void {
            const auto ncells = mat.ncol();
            auto ext = tatami::consecutive_extractor<sparse_>(mat, false, static_cast<Index_>(0), ncells, start, length);
 
            typedef typename EigenVector_::Scalar Scalar;
            tatami_stats::LocalOutputBuffer<Scalar> cbuffer(t, start, length, center_v.data());
            tatami_stats::LocalOutputBuffer<Scalar> sbuffer(t, start, length, scale_v.data());
 
            auto running = [&]{
                if constexpr(sparse_) {
                    return tatami_stats::variances::RunningSparse<Scalar, Value_, Index_>(length, cbuffer.data(), sbuffer.data(), /* skip_nan = */ false, /* subtract = */ start);
                } else {
                    return tatami_stats::variances::RunningDense<Scalar, Value_, Index_>(length, cbuffer.data(), sbuffer.data(), /* skip_nan = */ false);
                }
            }();
 
            auto vbuffer = tatami::create_container_of_Index_size<std::vector<Value_> >(length);
            auto ibuffer = [&]{
                if constexpr(sparse_) {
                    return tatami::create_container_of_Index_size<std::vector<Index_> >(length);
                } else {
                    return false;
                }
            }();
 
            for (Index_ c = 0; c < ncells; ++c) {
                if constexpr(sparse_) {
                    const auto range = ext->fetch(vbuffer.data(), ibuffer.data());
                    running.add(range.value, range.index, range.number);
                } else {
                    const auto ptr = ext->fetch(vbuffer.data());
                    running.add(ptr);
                }
            }
 
            running.finish();
            cbuffer.transfer();
            sbuffer.transfer();
        }, ngenes, num_threads);
    }
}
 
template<class IrlbaMatrix_, class EigenMatrix_, class EigenVector_>
auto run_irlba_deferred(
    const IrlbaMatrix_& mat,
    const SimplePcaOptions& options,
    EigenMatrix_& components, 
    EigenMatrix_& rotation, 
    EigenVector_& variance_explained,
    EigenVector_& center_v,
    EigenVector_& scale_v)
{
    irlba::Centered<IrlbaMatrix_, EigenVector_> centered(mat, center_v);
    if (options.scale) {
        irlba::Scaled<true, decltype(I(centered)), EigenVector_> scaled(centered, scale_v, true);
        return irlba::compute(scaled, options.number, components, rotation, variance_explained, options.irlba_options);
    } else {
        return irlba::compute(centered, options.number, components, rotation, variance_explained, options.irlba_options);
    }
}
 
template<typename Value_, typename Index_, class EigenMatrix_, class EigenVector_>
void run_sparse(
    const tatami::Matrix<Value_, Index_>& mat, 
    const SimplePcaOptions& options,
    EigenMatrix_& components, 
    EigenMatrix_& rotation, 
    EigenVector_& variance_explained,
    EigenVector_& center_v,
    EigenVector_& scale_v,
    typename EigenVector_::Scalar& total_var,
    bool& converged)
{
    const auto ngenes = mat.nrow();
    sanisizer::resize(center_v, ngenes);
    sanisizer::resize(scale_v, ngenes);
 
    if (options.realize_matrix) {
        // 'extracted' contains row-major contents...
        auto extracted = tatami::retrieve_compressed_sparse_contents<Value_, Index_>(
            mat, 
            /* row = */ true, 
            [&]{
                tatami::RetrieveCompressedSparseContentsOptions opt;
                opt.two_pass = false;
                opt.num_threads = options.num_threads;
                return opt;
            }()
        );
 
        // But we effectively transpose it to CSC with genes in columns.
        const Index_ ncells = mat.ncol();
        irlba::ParallelSparseMatrix emat(
            ncells,
            ngenes,
            std::move(extracted.value),
            std::move(extracted.index),
            std::move(extracted.pointers), 
            true,
            options.num_threads
        ); 
 
        tatami::parallelize([&](const int, const Index_ start, const Index_ length) -> void {
            const auto& pointers = emat.get_pointers();
            const auto& values = emat.get_values();
            for (Index_ g = start, end = start + length; g < end; ++g) {
                const auto offset = pointers[g];
                const auto next_offset = pointers[g + 1]; // increment won't overflow as 'g + 1 <= end'.
                const Index_ num_nonzero = next_offset - offset;
                const auto results = tatami_stats::variances::direct(values.data() + offset, num_nonzero, ncells, /* skip_nan = */ false);
                center_v.coeffRef(g) = results.first;
                scale_v.coeffRef(g) = results.second;
            }
        }, ngenes, options.num_threads);
 
        total_var = internal::process_scale_vector(options.scale, scale_v);
        const auto out = run_irlba_deferred(emat, options, components, rotation, variance_explained, center_v, scale_v);
        converged = out.first;
 
    } else {
        compute_row_means_and_variances<true>(mat, options.num_threads, center_v, scale_v);
        total_var = internal::process_scale_vector(options.scale, scale_v);
        const auto out = run_irlba_deferred(
            internal::TransposedTatamiWrapper<EigenVector_, Value_, Index_>(mat, options.num_threads), 
            options, 
            components, 
            rotation, 
            variance_explained, 
            center_v, 
            scale_v
        );
        converged = out.first;
    }
}
 
template<typename Value_, typename Index_, class EigenMatrix_, class EigenVector_>
void run_dense(
    const tatami::Matrix<Value_, Index_>& mat, 
    const SimplePcaOptions& options,
    EigenMatrix_& components, 
    EigenMatrix_& rotation, 
    EigenVector_& variance_explained, 
    EigenVector_& center_v,
    EigenVector_& scale_v,
    typename EigenVector_::Scalar& total_var,
    bool& converged)
{
    const Index_ ngenes = mat.nrow();
    sanisizer::resize(center_v, ngenes);
    sanisizer::resize(scale_v, ngenes);
 
    if (options.realize_matrix) {
        // Create a matrix with genes in columns.
        const Index_ ncells = mat.ncol();
        EigenMatrix_ emat(
            sanisizer::cast<decltype(I(std::declval<EigenMatrix_>().rows()))>(ncells),
            sanisizer::cast<decltype(I(std::declval<EigenMatrix_>().cols()))>(ngenes)
        );
 
        // If emat is row-major, we want to fill it with columns of 'mat', so row_major = false.
        // If emat is column-major, we want to fill it with rows of 'mat', so row_major = true.
        tatami::convert_to_dense(
            mat,
            /* row_major = */ !emat.IsRowMajor,
            emat.data(),
            [&]{
                tatami::ConvertToDenseOptions opt;
                opt.num_threads = options.num_threads;
                return opt;
            }()
        );
 
        center_v.array() = emat.array().colwise().sum();
        if (ncells) {
            center_v /= ncells;
        } else {
            std::fill(center_v.begin(), center_v.end(), std::numeric_limits<typename EigenVector_::Scalar>::quiet_NaN());
        }
        emat.array().rowwise() -= center_v.adjoint().array(); // applying it to avoid wasting time with deferred operations inside IRLBA.
 
        scale_v.array() = emat.array().colwise().squaredNorm();
        if (ncells > 1) {
            scale_v /= ncells - 1;
        } else {
            std::fill(scale_v.begin(), scale_v.end(), std::numeric_limits<typename EigenVector_::Scalar>::quiet_NaN());
        }
 
        total_var = internal::process_scale_vector(options.scale, scale_v);
        if (options.scale) {
            emat.array().rowwise() /= scale_v.adjoint().array();
        }
 
        const auto out = irlba::compute(emat, options.number, components, rotation, variance_explained, options.irlba_options);
        converged = out.first;
 
    } else {
        compute_row_means_and_variances<false>(mat, options.num_threads, center_v, scale_v);
        total_var = internal::process_scale_vector(options.scale, scale_v);
        const auto out = run_irlba_deferred(
            internal::TransposedTatamiWrapper<EigenVector_, Value_, Index_>(mat, options.num_threads), 
            options, 
            components, 
            rotation, 
            variance_explained, 
            center_v, 
            scale_v
        );
        converged = out.first;
    }
}
 
}
template<typename EigenMatrix_, typename EigenVector_>
struct SimplePcaResults {
    EigenMatrix_ components;
 
    EigenVector_ variance_explained;
 
    typename EigenVector_::Scalar total_variance = 0;
 
    EigenMatrix_ rotation;
 
    EigenVector_ center;
 
    EigenVector_ scale;
 
    bool converged = false;
};
 
template<typename Value_, typename Index_, typename EigenMatrix_, class EigenVector_>
void simple_pca(const tatami::Matrix<Value_, Index_>& mat, const SimplePcaOptions& options, SimplePcaResults<EigenMatrix_, EigenVector_>& output) {
    irlba::EigenThreadScope t(options.num_threads);
 
    if (mat.sparse()) {
        internal::run_sparse(mat, options, output.components, output.rotation, output.variance_explained, output.center, output.scale, output.total_variance, output.converged);
    } else {
        internal::run_dense(mat, options, output.components, output.rotation, output.variance_explained, output.center, output.scale, output.total_variance, output.converged);
    }
 
    internal::clean_up(mat.ncol(), output.components, output.variance_explained);
    if (options.transpose) {
        output.components.adjointInPlace();
    }
 
    if (!options.scale) {
        output.scale = EigenVector_();
    }
}
 
template<typename EigenMatrix_ = Eigen::MatrixXd, class EigenVector_ = Eigen::VectorXd, typename Value_, typename Index_>
SimplePcaResults<EigenMatrix_, EigenVector_> simple_pca(const tatami::Matrix<Value_, Index_>& mat, const SimplePcaOptions& options) {
    SimplePcaResults<EigenMatrix_, EigenVector_> output;
    simple_pca(mat, options, output);
    return output;
}
 
}
 
#endif