LogNormCounts.hpp

Go to the documentation of this file.

#ifndef SCRAN_LOG_NORM_COUNTS_H
#define SCRAN_LOG_NORM_COUNTS_H
 
#include "../utils/macros.hpp"
 
#include <algorithm>
#include <vector>
#include <numeric>
 
#include "tatami/tatami.hpp"
 
#include "CenterSizeFactors.hpp"
#include "SanitizeSizeFactors.hpp"
#include "ChoosePseudoCount.hpp"
 
namespace scran {
 
class LogNormCounts {
public:
    struct Defaults {
        static constexpr double pseudo_count = 1;
 
        static constexpr bool sparse_addition = true;
 
        static constexpr bool choose_pseudo_count = false;
 
        static constexpr bool center = true;
 
        static constexpr bool handle_zeros = false;
 
        static constexpr bool handle_non_finite = false;
 
        static constexpr int num_threads = 1;
    };
 
private:
    double pseudo_count = Defaults::pseudo_count;
    bool sparse_addition = Defaults::sparse_addition;
    bool handle_zeros = Defaults::handle_zeros;
    bool handle_non_finite = Defaults::handle_non_finite;
    int nthreads = Defaults::num_threads;
 
    bool center = Defaults::center;
    CenterSizeFactors centerer;
 
    bool choose_pseudo_count = Defaults::choose_pseudo_count;
    ChoosePseudoCount pseudo_chooser;
 
public:
    LogNormCounts& set_pseudo_count (double p = Defaults::pseudo_count) {
        pseudo_count = p;
        return *this;
    }
 
    LogNormCounts& set_sparse_addition(bool a = Defaults::sparse_addition) {
        sparse_addition = a;
        return *this;
    }
 
    LogNormCounts& set_center(bool c = Defaults::center) {
        center = c;
        return *this;
    }
 
    LogNormCounts& set_block_mode(CenterSizeFactors::BlockMode b = CenterSizeFactors::Defaults::block_mode) {
        centerer.set_block_mode(b);
        return *this;
    }
 
    LogNormCounts& set_handle_zeros(bool z = Defaults::handle_zeros) {
        handle_zeros = z;
        return *this;
    }
 
    LogNormCounts& set_handle_non_finite(bool n = Defaults::handle_non_finite) {
        handle_non_finite = n;
        return *this;
    }
 
    LogNormCounts& set_num_threads(int n = Defaults::num_threads) {
        nthreads = n;
        return *this;
    }
 
public:
    LogNormCounts& set_choose_pseudo_count(bool c = Defaults::choose_pseudo_count) {
        choose_pseudo_count = c;
        return *this;
    }
 
    LogNormCounts& set_max_bias(double m = ChoosePseudoCount::Defaults::max_bias) {
        pseudo_chooser.set_max_bias(m);
        return *this;
    }
 
    LogNormCounts& set_quantile(double q = ChoosePseudoCount::Defaults::quantile) {
        pseudo_chooser.set_quantile(q);
        return *this;
    }
 
    LogNormCounts& set_min_value(double m = ChoosePseudoCount::Defaults::min_value) {
        pseudo_chooser.set_min_value(m);
        return *this;
    }
 
public:
    template<class MAT, class V>
    std::shared_ptr<MAT> run(std::shared_ptr<MAT> mat, V size_factors) const {
        return run_blocked(std::move(mat), std::move(size_factors), static_cast<int*>(NULL));
    }
 
    template<class MAT, class V, typename B>
    std::shared_ptr<MAT> run_blocked(std::shared_ptr<MAT> mat, V size_factors, const B* block) const {
        // One might ask why we don't require a pointer for size_factors here.
        // It's because size_factors need to be moved into the Delayed operation
        // anyway, so we might as well ask the user to construct a vector for us.
        if (size_factors.size() != mat->ncol()) {
            throw std::runtime_error("number of size factors and columns are not equal");
        }
 
        SizeFactorValidity cresults;
        if (center) {
            cresults = centerer.run_blocked(size_factors.size(), size_factors.data(), block);
        } else {
            cresults = validate_size_factors(size_factors.size(), size_factors.data());
        }
 
        SanitizeSizeFactors sanitizer;
        sanitizer.set_handle_zero(handle_zeros ? SanitizeSizeFactors::HandlerAction::SANITIZE : SanitizeSizeFactors::HandlerAction::ERROR);
        sanitizer.set_handle_non_finite(handle_non_finite ? SanitizeSizeFactors::HandlerAction::SANITIZE : SanitizeSizeFactors::HandlerAction::ERROR);
        sanitizer.run(size_factors.size(), size_factors.data(), cresults);
 
        double current_pseudo = pseudo_count;
        if (choose_pseudo_count) {
            current_pseudo = pseudo_chooser.run(size_factors.size(), size_factors.data());
        }
 
        if (sparse_addition && current_pseudo != 1) {
            for (auto& d : size_factors) {
                d *= current_pseudo;
            }
            current_pseudo = 1; // effectively 1 now.
        }
 
        typedef typename MAT::value_type Value_;
        auto div = tatami::make_DelayedUnaryIsometricOp(std::move(mat), tatami::make_DelayedDivideVectorHelper<true, 1, Value_>(std::move(size_factors)));
 
        if (current_pseudo == 1) {
            return tatami::make_DelayedUnaryIsometricOp(std::move(div), tatami::DelayedLog1pHelper<Value_>(2.0));
        } else {
            auto add = tatami::make_DelayedUnaryIsometricOp(std::move(div), tatami::make_DelayedAddScalarHelper<Value_>(current_pseudo));
            return tatami::make_DelayedUnaryIsometricOp(std::move(add), tatami::DelayedLogHelper<Value_>(2.0));
        }
    }
 
public:
    template<class MAT>
    std::shared_ptr<MAT> run(std::shared_ptr<MAT> mat) const {
        auto size_factors = tatami::column_sums(mat.get(), nthreads);
        return run_blocked(std::move(mat), std::move(size_factors), static_cast<int*>(NULL));
    }
 
    template<class MAT, typename B>
    std::shared_ptr<MAT> run_blocked(std::shared_ptr<MAT> mat, const B* block) const {
        auto size_factors = tatami::column_sums(mat.get(), nthreads);
        return run_blocked(mat, std::move(size_factors), block);
    }
};
 
};
 
#endif