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 "sanisizer/sanisizer.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_cells, const SizeFactor_* const size_factors, const ComputeMeanSizeFactorOptions& options) {
    static_assert(std::is_floating_point<SizeFactor_>::value);
    SizeFactor_ mean = 0;
    I<decltype(num_cells)> denom = 0;
 
    if (options.ignore_invalid) {
        SizeFactorDiagnostics tmpdiag;
        for (I<decltype(num_cells)> i = 0; i < num_cells; ++i) {
            const auto val = size_factors[i];
            if (!internal::is_invalid(val, tmpdiag)) {
                mean += val;
                ++denom;
            }
        }
        if (options.diagnostics != NULL) {
            *(options.diagnostics) = tmpdiag;
        }
 
    } else {
        mean = std::accumulate(size_factors, size_factors + num_cells, static_cast<SizeFactor_>(0));
        denom = num_cells;
    }
 
    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_cells,
    const SizeFactor_* const size_factors,
    const Block_* const block,
    const std::size_t num_blocks,
    const ComputeMeanSizeFactorOptions& options
) {
    static_assert(std::is_floating_point<SizeFactor_>::value);
    auto block_mean = sanisizer::create<std::vector<SizeFactor_> >(num_blocks);
    auto block_num = sanisizer::create<std::vector<I<decltype(num_cells)> > >(num_blocks);
 
    if (options.ignore_invalid) {
        SizeFactorDiagnostics tmpdiag;
        for (I<decltype(num_cells)> i = 0; i < num_cells; ++i) {
            const auto val = size_factors[i];
            if (!internal::is_invalid(val, tmpdiag)) {
                const auto b = block[i];
                block_mean[b] += val;
                ++(block_num[b]);
            }
        }
        if (options.diagnostics != NULL) {
            *(options.diagnostics) = tmpdiag;
        }
 
    } else {
        for (I<decltype(num_cells)> i = 0; i < num_cells; ++i) {
            const auto b = block[i];
            block_mean[b] += size_factors[i];
            ++(block_num[b]);
        }
    }
 
    for (I<decltype(num_blocks)> g = 0; g < num_blocks; ++g) {
        if (block_num[g]) {
            block_mean[g] /= block_num[g];
        }
    }
 
    return block_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_cells, 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_cells, size_factors, copt);
 
    if (mean == 0) {
        return 0;
    }
 
    const SizeFactor_ mult = options.center / mean;
    for (I<decltype(num_cells)> i = 0; i < num_cells; ++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_cells,
    SizeFactor_* const size_factors,
    const Block_* const block,
    const std::size_t num_blocks,
    const CenterSizeFactorsBlockedOptions& options
) {
    ComputeMeanSizeFactorOptions copt;
    copt.ignore_invalid = options.ignore_invalid;
    copt.diagnostics = options.diagnostics;
    auto block_mean = compute_mean_size_factor_blocked(num_cells, size_factors, block, num_blocks, copt);
 
    if (options.block_mode == CenterBlockMode::PER_BLOCK) {
        std::vector<SizeFactor_> fac;
        fac.reserve(num_blocks);
        for (I<decltype(num_blocks)> g = 0; g < num_blocks; ++g) {
            const auto gm = block_mean[g];
            if (gm) {
                fac.emplace_back(1 / gm);
            } else {
                fac.emplace_back(1); // i.e., no-op.
            }
        }
 
        for (I<decltype(num_cells)> i = 0; i < num_cells; ++i) {
            size_factors[i] *= fac[block[i]];
        }
 
        if (options.report_final) {
            for (auto& gm : block_mean) {
                if (gm) {
                    gm = 1;
                }
            }
        }
 
        return block_mean;
 
    } else if (options.block_mode == CenterBlockMode::LOWEST) {
        SizeFactor_ min = 0;
        bool found = false;
        for (const auto m : block_mean) {
            // Ignore blocks 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_cells)> i = 0; i < num_cells; ++i) {
                size_factors[i] *= mult;
            }
 
            if (options.report_final) {
                for (auto& gm : block_mean) {
                    gm /= min;
                }
            }
        }
 
        return block_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);
        if (custom.size() != num_blocks) {
            throw std::runtime_error("length of 'custom_centers' should be equal to the number of blocks"); 
        }
 
        std::vector<SizeFactor_> fac;
        fac.reserve(num_blocks);
        for (I<decltype(num_blocks)> g = 0; g < num_blocks; ++g) {
            const auto gm = block_mean[g];
            if (gm) {
                fac.emplace_back(custom[g] / gm);
            } else {
                fac.emplace_back(1);
            }
        }
 
        for (I<decltype(num_cells)> i = 0; i < num_cells; ++i) {
            size_factors[i] *= fac[block[i]];
        }
 
        if (options.report_final) {
            for (I<decltype(num_blocks)> g = 0; g < num_blocks; ++g) {
                auto& gm = block_mean[g];
                if (gm) {
                    gm = custom[g];
                }
            }
        }
 
        return block_mean;
    }
}
 
}
 
#endif