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 {
 
struct ComputeMeanSizeFactorOptions {
    bool ignore_invalid = true;
 
    SizeFactorDiagnostics* diagnostics = NULL;
};
 
template<typename SizeFactor_>
SizeFactor_ compute_mean_size_factor(const std::size_t num, const SizeFactor_* const size_factors, const ComputeMeanSizeFactorOptions& options) {
    static_assert(std::is_floating_point<SizeFactor_>::value);
    SizeFactor_ mean = 0;
    I<decltype(num)> denom = 0;
 
    if (options.ignore_invalid) {
        SizeFactorDiagnostics tmpdiag;
        auto& diag = (options.diagnostics == NULL ? tmpdiag : *(options.diagnostics));
        for (I<decltype(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_, typename Block_>
std::vector<SizeFactor_> compute_mean_size_factor_blocked(
    const std::size_t num,
    const SizeFactor_* const size_factors,
    const Block_* const block,
    const ComputeMeanSizeFactorOptions& 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<I<decltype(num)>> >(ngroups);
 
    if (options.ignore_invalid) {
        SizeFactorDiagnostics tmpdiag;
        auto& diag = (options.diagnostics == NULL ? tmpdiag : *(options.diagnostics));
        for (I<decltype(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 (I<decltype(num)> i = 0; i < num; ++i) {
            const auto b = block[i];
            group_mean[b] += size_factors[i];
            ++(group_num[b]);
        }
    }
 
    for (I<decltype(ngroups)> g = 0; g < ngroups; ++g) {
        if (group_num[g]) {
            group_mean[g] /= group_num[g];
        }
    }
 
    return group_mean;
}
 
struct CenterSizeFactorsOptions {
    bool ignore_invalid = true;
 
    double center = 1;
 
    SizeFactorDiagnostics* diagnostics = NULL;
 
    bool report_final = false;
};
 
template<typename SizeFactor_>
SizeFactor_ center_size_factors(const std::size_t num, SizeFactor_* const size_factors, const CenterSizeFactorsOptions& options) {
    ComputeMeanSizeFactorOptions copt;
    copt.ignore_invalid = options.ignore_invalid;
    copt.diagnostics = options.diagnostics;
    const auto mean = compute_mean_size_factor(num, size_factors, copt);
 
    if (mean == 0) {
        return 0;
    }
 
    const SizeFactor_ mult = options.center / mean;
    for (I<decltype(num)> i = 0; i < num; ++i){
        size_factors[i] *= mult;
    }
 
    if (options.report_final) {
        return options.center;
    } else {
        return mean;
    }
}
 
enum class CenterBlockMode : char { PER_BLOCK, LOWEST, CUSTOM };
 
struct CenterSizeFactorsBlockedOptions {
    bool ignore_invalid = true;
 
    SizeFactorDiagnostics* diagnostics = NULL;
 
    CenterBlockMode block_mode = CenterBlockMode::LOWEST;
 
    std::optional<std::vector<double> > custom_centers;
 
    bool report_final = false;
};
 
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,
    const CenterSizeFactorsBlockedOptions& options
) {
    ComputeMeanSizeFactorOptions copt;
    copt.ignore_invalid = options.ignore_invalid;
    copt.diagnostics = options.diagnostics;
    auto group_mean = compute_mean_size_factor_blocked(num, size_factors, block, copt);
    const auto ngroups = group_mean.size();
 
    if (options.block_mode == CenterBlockMode::PER_BLOCK) {
        std::vector<SizeFactor_> fac;
        fac.reserve(ngroups);
        for (I<decltype(ngroups)> g = 0; g < ngroups; ++g) {
            const auto gm = group_mean[g];
            if (gm) {
                fac.emplace_back(1 / gm);
            } else {
                fac.emplace_back(1); // i.e., no-op.
            }
        }
 
        for (I<decltype(num)> i = 0; i < num; ++i) {
            size_factors[i] *= fac[block[i]];
        }
 
        if (options.report_final) {
            for (auto& gm : group_mean) {
                if (gm) {
                    gm = 1;
                }
            }
        }
        return group_mean;
 
    } 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) {
            const SizeFactor_ mult = 1 / min;
            for (I<decltype(num)> i = 0; i < num; ++i) {
                size_factors[i] *= mult;
            }
 
            if (options.report_final) {
                for (auto& gm : group_mean) {
                    gm /= min;
                }
            }
        }
 
        return group_mean;
 
    } else { // i.e., options.block_mode == CenterBlockMode::CUSTOM
        if (!options.custom_centers.has_value()) {
            throw std::runtime_error("'custom_centers' should be set for custom block centers"); 
        }
        const auto& custom = *(options.custom_centers);
        const auto ngroups = group_mean.size();
        if (custom.size() != ngroups) {
            throw std::runtime_error("length of 'custom_centers' should be equal to the number of groups"); 
        }
 
        std::vector<SizeFactor_> fac;
        fac.reserve(ngroups);
        for (I<decltype(ngroups)> g = 0; g < ngroups; ++g) {
            const auto gm = group_mean[g];
            if (gm) {
                fac.emplace_back(custom[g] / gm);
            } else {
                fac.emplace_back(1);
            }
        }
 
        for (I<decltype(num)> i = 0; i < num; ++i) {
            size_factors[i] *= fac[block[i]];
        }
 
        if (options.report_final) {
            for (I<decltype(ngroups)> g = 0; g < ngroups; ++g) {
                auto& gm = group_mean[g];
                if (gm) {
                    gm = custom[g];
                }
            }
        }
        return group_mean;
    }
}
 
}
 
#endif