CenterSizeFactors.hpp

Go to the documentation of this file.

#ifndef SCRAN_CENTER_SIZE_FACTORS_HPP
#define SCRAN_CENTER_SIZE_FACTORS_HPP
 
#include "../utils/macros.hpp"
 
#include "../utils/blocking.hpp"
 
#include "SanitizeSizeFactors.hpp"
#include <stdexcept>
#include <vector>
#include <algorithm>
 
namespace scran {
 
class CenterSizeFactors {
public:
    enum BlockMode { PER_BLOCK, LOWEST };
 
    struct Defaults {
        static constexpr BlockMode block_mode = LOWEST;
 
        static constexpr bool ignore_zeros = true;
    };
 
private:
    BlockMode block_mode = Defaults::block_mode;
 
    bool ignore_zeros = Defaults::ignore_zeros;
 
public:
    CenterSizeFactors& set_block_mode(BlockMode b = Defaults::block_mode) {
        block_mode = b;
        return *this;
    }
 
    CenterSizeFactors& set_ignore_zeros(bool i = Defaults::ignore_zeros) {
        ignore_zeros = i;
        return *this;
    }
 
public:
    template<typename T>
    SizeFactorValidity run(size_t n, T* size_factors) const {
        size_t num_used = 0;
        SizeFactorValidity output;
 
        double mean = 0;
        for (size_t i = 0; i < n; ++i) {
            const auto& current = size_factors[i];
            if (current < 0) {
                output.has_negative = true;
                continue;
            } else if (current == 0) {
                output.has_zero = true;
                if (ignore_zeros) {
                    continue;
                }
            } else if (std::isnan(current)) {
                output.has_nan = true;
                continue;
            } else if (std::isinf(current)) {
                output.has_infinite = true;
                continue;
            }
 
            ++num_used;
            mean += current;
        }
 
        if (mean) {
            mean /= num_used;
            for (size_t i = 0; i < n; ++i){
                size_factors[i] /= mean;
            }
        }
 
        return output;
    }
 
    template<typename T, typename B>
    SizeFactorValidity run_blocked(size_t n, T* size_factors, const B* block) const {
        if (block == NULL) {
            return run(n, size_factors);
        } 
 
        size_t ngroups = count_ids(n, block);
        std::vector<double> group_mean(ngroups);
        std::vector<double> group_num(ngroups);
 
        SizeFactorValidity output;
        for (size_t i = 0; i < n; ++i) {
            const auto& current = size_factors[i];
            if (current < 0) {
                output.has_negative = true;
                continue;
            } else if (current == 0) {
                output.has_zero = true;
                if (ignore_zeros) {
                    continue;
                }
            } else if (std::isinf(current)) {
                output.has_infinite = true;
                continue;
            } else if (std::isnan(current)) {
                output.has_nan = true;
                continue;
            }
 
            const auto& b = block[i];
            group_mean[b] += size_factors[i];
            ++group_num[b];
        }
 
        for (size_t g = 0; g < ngroups; ++g) {
            if (group_num[g]) {
                group_mean[g] /= group_num[g];
            }
        }
 
        if (block_mode == PER_BLOCK) {
            for (size_t i = 0; i < n; ++i) {
                const auto& div = group_mean[block[i]];
                if (div) {
                    size_factors[i] /= div;
                }
            }
 
        } else if (block_mode == LOWEST) {
            // Ignore groups with means of zeros, either because they're full of zeros themselves
            // or they have no cells associated with them.
            double min = std::numeric_limits<double>::infinity();
            for (size_t b = 0; b < ngroups; ++b) {
                const auto& div = group_mean[b];
                if (div && div < min) {
                    min = div;
                }
            }
 
            if (std::isfinite(min)) {
                for (size_t i = 0; i < n; ++i) {
                    size_factors[i] /= min;
                }
            }
        }
 
        return output;
    }
};
 
}
 
#endif