ScoreFeatureSet.hpp

Go to the documentation of this file.

#ifndef SCRAN_SCORE_FEATURE_SET_HPP
#define SCRAN_SCORE_FEATURE_SET_HPP
 
#include "../utils/macros.hpp"
 
#include <algorithm>
#include <vector>
#include "tatami/tatami.hpp"
 
#include "irlba/irlba.hpp"
#include "Eigen/Dense"
 
#include "../dimensionality_reduction/SimplePca.hpp"
#include "../dimensionality_reduction/ResidualPca.hpp"
 
namespace scran {
 
class ScoreFeatureSet {
public:
    struct Defaults {
        static constexpr WeightPolicy block_weight_policy = WeightPolicy::VARIABLE;
 
        static constexpr int num_threads = 1;
 
        static constexpr bool scale = false;
    };
 
private:
    WeightPolicy block_weight_policy = Defaults::block_weight_policy;
    int nthreads = Defaults::num_threads;
    bool scale = Defaults::scale;
 
public:
    ScoreFeatureSet& set_block_weight_policy(WeightPolicy b = Defaults::block_weight_policy) {
        block_weight_policy = b;
        return *this;
    }
 
    ScoreFeatureSet& set_num_threads(int n = Defaults::num_threads) {
        nthreads = n;
        return *this;
    }
 
    ScoreFeatureSet& set_scale(bool s = Defaults::scale) {
        scale = s;
        return *this;
    }
 
public:
    struct Results {
        std::vector<double> scores;
 
        std::vector<double> weights;
    };
 
private:
    /*
     * We have the first PC 'P' and the first rotation vector 'R', plus a centering vector 'C' 
     * and scaling vector 'S'. The low-rank approximation is defined as (using R syntax):
     *
     *     L = outer(R, P) * S + C 
     *       = outer(R * S, P) + C
     *
     * Remember that we want the column means of the rank-1 approximation, so:
     *
     *     colMeans(L) = mean(R * S) * P + colMeans(C)
     *
     * If scale = false, then S can be dropped from the above expression.
     */
    double compute_multiplier(const Eigen::MatrixXd& rotation, const Eigen::VectorXd& scale_v) const {
        auto first_rot = rotation.col(0);
 
        double multiplier = 0;
        if (scale) {
            for (Eigen::Index i = 0, end = first_rot.size(); i < end; ++i) {
                multiplier += scale_v.coeff(i) * first_rot.coeff(i);
            }
        } else {
            multiplier = std::accumulate(first_rot.begin(), first_rot.end(), 0.0);
        }
 
        // no need to protect against zero rows, as that should already be caught.
        return multiplier / first_rot.size();
    }
 
    void transfer_rotation(const Eigen::MatrixXd& rotation, std::vector<double>& weights) const {
        auto first_rot = rotation.col(0);
        weights.insert(weights.end(), first_rot.begin(), first_rot.end());
        return;
    }
 
public:
    template<typename T, typename IDX, typename X, typename Block>
    Results run_blocked(const tatami::Matrix<T, IDX>* mat, const X* features, const Block* block) const {
        std::shared_ptr<const tatami::Matrix<T, IDX> > subsetted = pca_utils::subset_matrix_by_features(mat, features);
        auto NR = subsetted->nrow();
        auto NC = subsetted->ncol();
 
        // Catching edge cases.
        if (NR == 0) {
            Results output;
            output.scores.resize(NC);
            return output;
        } else if (NR == 1) {
            Results output;
            output.weights.push_back(1);
            output.scores = subsetted->dense_row()->fetch(0);
            return output;
        } else if (NC == 0) {
            Results output;
            output.weights.resize(NR);
            return output;
        }
 
        Results output;
 
        if (block == NULL) {
            SimplePca runner;
            runner.set_rank(1);
            runner.set_scale(scale);
            runner.set_num_threads(nthreads);
            runner.set_return_rotation(true).set_return_scale(scale).set_return_center(true);
 
            auto temp = runner.run(subsetted.get());
            transfer_rotation(temp.rotation, output.weights);
            double multiplier = compute_multiplier(temp.rotation, temp.scale);
            double shift = std::accumulate(temp.center.begin(), temp.center.end(), 0.0) / temp.center.size();
 
            output.scores.resize(temp.pcs.cols());
            for (Eigen::Index c = 0, end = temp.pcs.cols(); c < end; ++c) {
                output.scores[c] = temp.pcs.coeff(0, c) * multiplier + shift;
            }
 
        } else {
            ResidualPca runner;
            runner.set_rank(1);
            runner.set_scale(scale);
            runner.set_num_threads(nthreads);
            runner.set_return_rotation(true).set_return_scale(scale).set_return_center(true);
            runner.set_block_weight_policy(block_weight_policy);
 
            auto temp = runner.run(subsetted.get(), block);
            transfer_rotation(temp.rotation, output.weights);
            double multiplier = compute_multiplier(temp.rotation, temp.scale);
 
            // 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.
            Eigen::VectorXd shift = temp.center.colwise().sum() / temp.center.rows();
            output.scores.resize(temp.pcs.cols());
            for (Eigen::Index c = 0, end = temp.pcs.cols(); c < end; ++c) {
                output.scores[c] = temp.pcs.coeff(0, c) * multiplier + shift.coeff(block[c]);
            }
        }
 
        return output;
    }
 
    template<typename T, typename IDX, typename X>
    Results run(const tatami::Matrix<T, IDX>* mat, const X* features) const {
        return run_blocked(mat, features, static_cast<unsigned char*>(NULL));
    }
};
 
}
 
#endif