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 "tatami_stats/tatami_stats.hpp"
 
#include <vector>
#include <numeric>
#include <algorithm>
#include <type_traits>
#include <cstddef>
 
#include "sanitize_size_factors.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(std::size_t num, const SizeFactor_* size_factors, SizeFactorDiagnostics* diagnostics, const CenterSizeFactorsOptions& options) {
    static_assert(std::is_floating_point<SizeFactor_>::value);
    SizeFactor_ mean = 0;
    decltype(num) denom = 0;
 
    if (options.ignore_invalid) {
        SizeFactorDiagnostics tmpdiag;
        auto& diag = (diagnostics == NULL ? tmpdiag : *diagnostics);
        for (decltype(num) i = 0; i < num; ++i) {
            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(std::size_t num, SizeFactor_* size_factors, SizeFactorDiagnostics* diagnostics, const CenterSizeFactorsOptions& options) {
    auto mean = center_size_factors_mean(num, size_factors, diagnostics, options);
    if (mean) {
        for (decltype(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(std::size_t num, const SizeFactor_* size_factors, const Block_* block, SizeFactorDiagnostics* diagnostics, const CenterSizeFactorsOptions& options) {
    static_assert(std::is_floating_point<SizeFactor_>::value);
    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(num)> >(ngroups);
 
    if (options.ignore_invalid) {
        SizeFactorDiagnostics tmpdiag;
        auto& diag = (diagnostics == NULL ? tmpdiag : *diagnostics);
        for (decltype(num) i = 0; i < num; ++i) {
            auto val = size_factors[i];
            if (!internal::is_invalid(val, diag)) {
                auto b = block[i];
                group_mean[b] += val;
                ++(group_num[b]);
            }
        }
    } else {
        for (decltype(num) i = 0; i < num; ++i) {
            auto b = block[i];
            group_mean[b] += size_factors[i];
            ++(group_num[b]);
        }
    }
 
    for (decltype(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(std::size_t num, SizeFactor_* size_factors, const Block_* block, SizeFactorDiagnostics* diagnostics, const CenterSizeFactorsOptions& options) {
    auto group_mean = center_size_factors_blocked_mean(num, size_factors, block, diagnostics, options);
 
    if (options.block_mode == CenterBlockMode::PER_BLOCK) {
        for (decltype(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 (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(num) i = 0; i < num; ++i) {
                size_factors[i] /= min;
            }
        }
    }
 
    return group_mean;
}
 
}
 
#endif