center_size_factors.hpp

Go to the documentation of this file.

#ifndef SCRAN_NORM_CENTER_SIZE_FACTORS_HPP
#define SCRAN_NORM_CENTER_SIZE_FACTORS_HPP
 
#include <vector>
#include <numeric>
#include <algorithm>
#include <type_traits>
#include <cstddef>
 
#include "tatami_stats/tatami_stats.hpp"
 
#include "sanitize_size_factors.hpp"
#include "utils.hpp"
 
namespace scran_norm {
 
enum class CenterBlockMode : char { PER_BLOCK, LOWEST };
 
struct CenterSizeFactorsOptions {
    CenterBlockMode block_mode = CenterBlockMode::LOWEST;
 
    bool ignore_invalid = true;
};
 
template<typename SizeFactor_>
SizeFactor_ center_size_factors_mean(const std::size_t num, const SizeFactor_* const size_factors, SizeFactorDiagnostics* const diagnostics, const CenterSizeFactorsOptions& options) {
    static_assert(std::is_floating_point<SizeFactor_>::value);
    SizeFactor_ mean = 0;
    decltype(I(num)) denom = 0;
 
    if (options.ignore_invalid) {
        SizeFactorDiagnostics tmpdiag;
        auto& diag = (diagnostics == NULL ? tmpdiag : *diagnostics);
        for (decltype(I(num)) i = 0; i < num; ++i) {
            const auto val = size_factors[i];
            if (!internal::is_invalid(val, diag)) {
                mean += val;
                ++denom;
            }
        }
    } else {
        mean = std::accumulate(size_factors, size_factors + num, static_cast<SizeFactor_>(0));
        denom = num;
    }
 
    if (denom) {
        return mean/denom;
    } else {
        return 0;
    }
}
 
template<typename SizeFactor_>
SizeFactor_ center_size_factors(const std::size_t num, SizeFactor_* const size_factors, SizeFactorDiagnostics* const diagnostics, const CenterSizeFactorsOptions& options) {
    const auto mean = center_size_factors_mean(num, size_factors, diagnostics, options);
    if (mean) {
        for (decltype(I(num)) i = 0; i < num; ++i){
            size_factors[i] /= mean;
        }
    }
    return mean;
}
 
template<typename SizeFactor_, typename Block_>
std::vector<SizeFactor_> center_size_factors_blocked_mean(
    const std::size_t num,
    const SizeFactor_* const size_factors,
    const Block_* const block,
    SizeFactorDiagnostics* const diagnostics,
    const CenterSizeFactorsOptions& options)
{
    static_assert(std::is_floating_point<SizeFactor_>::value);
    const auto ngroups = tatami_stats::total_groups(block, num);
    auto group_mean = sanisizer::create<std::vector<SizeFactor_> >(ngroups);
    auto group_num = sanisizer::create<std::vector<decltype(I(num))> >(ngroups);
 
    if (options.ignore_invalid) {
        SizeFactorDiagnostics tmpdiag;
        auto& diag = (diagnostics == NULL ? tmpdiag : *diagnostics);
        for (decltype(I(num)) i = 0; i < num; ++i) {
            const auto val = size_factors[i];
            if (!internal::is_invalid(val, diag)) {
                const auto b = block[i];
                group_mean[b] += val;
                ++(group_num[b]);
            }
        }
    } else {
        for (decltype(I(num)) i = 0; i < num; ++i) {
            const auto b = block[i];
            group_mean[b] += size_factors[i];
            ++(group_num[b]);
        }
    }
 
    for (decltype(I(ngroups)) g = 0; g < ngroups; ++g) {
        if (group_num[g]) {
            group_mean[g] /= group_num[g];
        }
    }
 
    return group_mean;
}
 
template<typename SizeFactor_, typename Block_>
std::vector<SizeFactor_> center_size_factors_blocked(
    const std::size_t num,
    SizeFactor_* const size_factors,
    const Block_* const block,
    SizeFactorDiagnostics* const diagnostics,
    const CenterSizeFactorsOptions& options)
{
    const auto group_mean = center_size_factors_blocked_mean(num, size_factors, block, diagnostics, options);
 
    if (options.block_mode == CenterBlockMode::PER_BLOCK) {
        for (decltype(I(num)) i = 0; i < num; ++i) {
            const auto& div = group_mean[block[i]];
            if (div) {
                size_factors[i] /= div;
            }
        }
 
    } else if (options.block_mode == CenterBlockMode::LOWEST) {
        SizeFactor_ min = 0;
        bool found = false;
        for (const auto m : group_mean) {
            // Ignore groups with means of zeros, either because they're full
            // of zeros themselves or they have no cells associated with them.
            if (m) {
                if (!found || m < min) {
                    min = m;
                    found = true;
                }
            }
        }
 
        if (min > 0) {
            for (decltype(I(num)) i = 0; i < num; ++i) {
                size_factors[i] /= min;
            }
        }
    }
 
    return group_mean;
}
 
}
 
#endif