quick_grouped_size_factors.hpp

Go to the documentation of this file.

#ifndef SCRAN_QUICK_GROUPED_SIZE_FACTORS_HPP
#define SCRAN_QUICK_GROUPED_SIZE_FACTORS_HPP
 
#include "../utils/macros.hpp"
 
#include <algorithm>
#include <vector>
#include <cmath>
#include <functional>
#include <memory>
 
#include "tatami/tatami.hpp"
#include "kmeans/Kmeans.hpp"
#include "kmeans/InitializePCAPartition.hpp"
 
#include "../utils/blocking.hpp"
#include "../dimensionality_reduction/SimplePca.hpp"
#include "LogNormCounts.hpp"
#include "GroupedSizeFactors.hpp"
 
namespace scran {
 
namespace quick_grouped_size_factors {
 
template<
    typename Block_ = int,
    typename SizeFactor_ = double
>
struct Options {
    int rank = 25;
 
    std::function<size_t(size_t)> clusters;
 
    const Block_* block = NULL;
 
    const SizeFactor_* initial_factors = NULL;
 
    int num_threads = 1;
};
 
namespace internal {
 
template<
    typename Value_, 
    typename Index_
>
auto cluster(const tatami::Matrix<Value_, Index_>* mat, int rank, size_t clusters, int num_threads) {
    SimplePca pca_runner;
    pca_runner.set_rank(rank);
    pca_runner.set_num_threads(num_threads);
    auto pc_out = pca_runner.run(mat);
    const auto& pcs = pc_out.pcs;
 
    kmeans::Kmeans kmeans_runner;
    kmeans_runner.set_num_threads(num_threads);
    kmeans::InitializePCAPartition<Value_, Index_, Index_> init;
    return kmeans_runner.run(
        pcs.rows(), 
        pcs.cols(), 
        pcs.data(), 
        clusters,
        &init
    );
}
 
}
template<
    typename Value_, 
    typename Index_, 
    typename OutputFactor_, 
    typename Block_,
    typename SizeFactor_
>
void run(const tatami::Matrix<Value_, Index_>* mat, OutputFactor_* output, const Options<Block_, SizeFactor_>& opt) {
    std::vector<Index_> clusters;
    Index_ NC = mat->ncol();
    auto ptr = tatami::wrap_shared_ptr(mat);
 
    LogNormCounts logger;
    logger.set_num_threads(opt.num_threads);
 
    auto fun = opt.clusters;
    if (!fun) {
        fun = [](size_t n) -> size_t {
            size_t candidate = std::sqrt(static_cast<double>(n));
            return std::min(candidate, static_cast<size_t>(50));
        };
    }
 
    if (opt.block) {
        auto nblocks = count_ids(NC, opt.block);
        std::vector<std::vector<Index_> > assignments(nblocks);
        for (Index_ c = 0; c < NC; ++c) {
            assignments[opt.block[c]].push_back(c);
        }
 
        clusters.resize(NC);
        Index_ last_cluster = 0;
 
        for (size_t b = 0; b < nblocks; ++b) {
            const auto& inblock = assignments[b];
            auto subptr = tatami::make_DelayedSubset<1>(ptr, tatami::ArrayView<Index_>(inblock.data(), inblock.size()));
 
            std::shared_ptr<tatami::Matrix<Value_, Index_> > normalized;
            if (opt.initial_factors) {
                std::vector<SizeFactor_> fac;
                fac.reserve(inblock.size());
                for (auto i : inblock) {
                    fac.push_back(opt.initial_factors[i]);
                }
                normalized = logger.run(std::move(subptr), std::move(fac));
            } else {
                normalized = logger.run(std::move(subptr));
            }
 
            auto res = internal::cluster(normalized.get(), opt.rank, fun(inblock.size()), opt.num_threads);
            auto cIt = res.clusters.begin();
            for (auto i : inblock) {
                clusters[i] = *cIt + last_cluster;
                ++cIt;
            }
            last_cluster += *std::max_element(res.clusters.begin(), res.clusters.end()) + 1;
        }
 
    } else {
        std::shared_ptr<const tatami::Matrix<Value_, Index_> > normalized; // TODO: avoid propagating const'ness from LogNormCounts.
        if (opt.initial_factors) {
            std::vector<SizeFactor_> fac(opt.initial_factors, opt.initial_factors + NC);
            normalized = logger.run(std::move(ptr), std::move(fac));
        } else {
            normalized = logger.run(std::move(ptr));
        }
 
        auto res = internal::cluster(normalized.get(), opt.rank, fun(NC), opt.num_threads);
        clusters = std::move(res.clusters);
    }
 
    GroupedSizeFactors group_runner;
    group_runner.set_num_threads(opt.num_threads);
    group_runner.run(mat, clusters.data(), output);
    return;
}
 
template<
    typename Value_, 
    typename Index_, 
    typename OutputFactor_
>
void run(const tatami::Matrix<Value_, Index_>* mat, OutputFactor_* output) {
    run(mat, output, Options<>());
}
 
template<
    typename OutputFactor_ = double, 
    typename Value_, 
    typename Index_, 
    typename Block_,
    typename SizeFactor_
>
std::vector<OutputFactor_> run(const tatami::Matrix<Value_, Index_>* mat, const Options<Block_, SizeFactor_>& opt) {
    std::vector<OutputFactor_> output(mat->ncol());
    run(mat, output.data(), opt);
    return output;
}
 
template<
    typename OutputFactor_ = double, 
    typename Value_, 
    typename Index_
>
std::vector<OutputFactor_> run(const tatami::Matrix<Value_, Index_>* mat) {
    std::vector<OutputFactor_> output(mat->ncol());
    run(mat, output.data(), Options<>());
    return output;
}
 
}
 
}
 
#endif