combine_factors.hpp

Go to the documentation of this file.

#ifndef SCRAN_AGGREGATE_COMBINE_FACTORS_HPP
#define SCRAN_AGGREGATE_COMBINE_FACTORS_HPP
 
#include <algorithm>
#include <vector>
#include <map>
#include <unordered_map>
#include <typeindex>
 
#include "clean_factor.hpp"
 
namespace scran_aggregate {
 
template<typename Factor_, typename Combined_>
std::vector<std::vector<Factor_> > combine_factors(size_t n, const std::vector<const Factor_*>& factors, Combined_* combined) {
    size_t nfac = factors.size();
    std::vector<std::vector<Factor_> > output(nfac);
 
    // Handling the special cases.
    if (nfac == 0) {
        std::fill_n(combined, n, 0);
        return output;
    }
    if (nfac == 1) {
        output[0] = clean_factor(n, factors.front(), combined);
        return output;
    }
 
    // Creating a hashmap on the combinations of each factor.
    struct Combination {
        Combination(size_t i) : index(i) {}
        size_t index;
    };
 
    auto unique = [&]{ // scoping this in an IIFE to release map memory sooner.
        // Using a map with a custom comparator that uses the index
        // of first occurrence of each factor as the key. Currently using a map
        // to (i) avoid issues with collisions of combined hashes and (ii)
        // avoid having to write more code for sorting a vector of arrays.
        auto cmp = [&](const Combination& left, const Combination& right) -> bool {
            for (auto curf : factors) {
                if (curf[left.index] < curf[right.index]) {
                    return true;
                } else if (curf[left.index] > curf[right.index]) {
                    return false;
                }
            }
            return false;
        };
 
        auto eq = [&](const Combination& left, const Combination& right) -> bool {
            for (auto curf : factors) {
                if (curf[left.index] != curf[right.index]) {
                    return false;
                }
            }
            return true;
        };
 
        std::map<Combination, Combined_, decltype(cmp)> mapping(std::move(cmp));
        for (size_t i = 0; i < n; ++i) {
            Combination current(i);
            auto mIt = mapping.find(current);
            if (mIt == mapping.end() || !eq(mIt->first, current)) {
                Combined_ alt = mapping.size();
                mapping.insert(mIt, std::make_pair(current, alt));
                combined[i] = alt;
            } else {
                combined[i] = mIt->second;
            }
        }
 
        return std::vector<std::pair<Combination, Combined_> >(mapping.begin(), mapping.end());
    }();
 
    // Remapping to a sorted set.
    size_t nuniq = unique.size();
    for (auto& ofac : output) {
        ofac.reserve(nuniq);
    }
    std::vector<Combined_> remapping(nuniq);
    for (size_t u = 0; u < nuniq; ++u) {
        auto ix = unique[u].first.index;
        for (size_t f = 0; f < nfac; ++f) {
            output[f].push_back(factors[f][ix]);
        }
        remapping[unique[u].second] = u;
    }
 
    // Mapping each cell to its sorted combination.
    for (size_t i = 0; i < n; ++i) {
        combined[i] = remapping[combined[i]];
    }
 
    return output;
}
 
template<typename Factor_, typename Number_, typename Combined_>
std::vector<std::vector<Factor_> > combine_factors_unused(size_t n, const std::vector<std::pair<const Factor_*, Number_> >& factors, Combined_* combined) {
    size_t nfac = factors.size();
    std::vector<std::vector<Factor_> > output(nfac);
 
    // Handling the special cases.
    if (nfac == 0) {
        std::fill_n(combined, n, 0);
        return output;
    }
    if (nfac == 1) {
        output[0].resize(factors[0].second);
        std::iota(output[0].begin(), output[0].end(), static_cast<Combined_>(0));
        std::copy_n(factors[0].first, n, combined);
        return output;
    }
 
    // We iterate from back to front, where the first factor is the slowest changing.
    std::copy_n(factors[nfac - 1].first, n, combined); 
    Combined_ mult = factors[nfac - 1].second;
    for (size_t f = nfac - 1; f > 0; --f) {
        const auto& finfo = factors[f - 1];
        auto ff = finfo.first;
        for (size_t i = 0; i < n; ++i) {
            combined[i] += mult * ff[i];
        }
        mult *= finfo.second;
    }
 
    auto ncombos = mult;
    Combined_ outer_repeats = mult;
    Combined_ inner_repeats = 1;
    for (size_t f = nfac; f > 0; --f) {
        auto& out = output[f - 1];
        out.reserve(ncombos);
 
        const auto& finfo = factors[f - 1];
        size_t initial_size = inner_repeats * finfo.second;
        out.resize(initial_size);
 
        if (inner_repeats == 1) {
            std::iota(out.begin(), out.end(), static_cast<Combined_>(0));
        } else {
            auto oIt = out.begin();
            for (Number_ l = 0; l < finfo.second; ++l) {
                std::fill_n(oIt, inner_repeats, l);
                oIt += inner_repeats;
            }
        }
        inner_repeats = initial_size;
 
        outer_repeats /= finfo.second;
        for (Combined_ r = 1; r < outer_repeats; ++r) {
            out.insert(out.end(), out.begin(), out.begin() + initial_size);
        }
    }
 
    return output;
}
 
}
 
#endif