convert_to_layered_sparse.hpp

Go to the documentation of this file.

#ifndef TATAMI_CONVERT_TO_LAYERED_SPARSE_HPP
#define TATAMI_CONVERT_TO_LAYERED_SPARSE_HPP
 
#include <cstdint>
#include <vector>
#include <memory>
#include <limits>
#include <algorithm>
 
#include "tatami/tatami.hpp"
#include "sanisizer/sanisizer.hpp"
 
#include "utils.hpp"
 
namespace tatami_layered {
 
template<typename ColIndex_, typename ValueOut_ = double, typename IndexOut_ = int, typename ValueIn_, typename IndexIn_>
std::shared_ptr<tatami::Matrix<ValueOut_, IndexOut_> > convert_by_row(const tatami::Matrix<ValueIn_, IndexIn_>& mat, const IndexIn_ chunk_size, const int nthreads) {
    const auto NR = mat.nrow(), NC = mat.ncol();
    const IndexIn_ leftovers = NC % chunk_size;
    const IndexIn_ nchunks = sanisizer::max(1, NC / chunk_size + (leftovers != 0));
 
    auto store8  = tatami::create_container_of_Index_size<std::vector<Holder< std::uint8_t, IndexOut_, ColIndex_> > >(nchunks);
    auto store16 = tatami::create_container_of_Index_size<std::vector<Holder<std::uint16_t, IndexOut_, ColIndex_> > >(nchunks);
    auto store32 = tatami::create_container_of_Index_size<std::vector<Holder<std::uint32_t, IndexOut_, ColIndex_> > >(nchunks);
 
    auto identities8  = tatami::create_container_of_Index_size<std::vector<std::vector<IndexOut_> > >(nchunks);
    auto identities16 = tatami::create_container_of_Index_size<std::vector<std::vector<IndexOut_> > >(nchunks);
    auto identities32 = tatami::create_container_of_Index_size<std::vector<std::vector<IndexOut_> > >(nchunks);
 
    auto assigned_position = tatami::create_container_of_Index_size<std::vector<std::vector<IndexOut_> > >(nchunks);
    auto assigned_category = tatami::create_container_of_Index_size<std::vector<std::vector<Category> > >(nchunks);
 
    // First pass to define the allocations.
    {
        auto max_per_chunk = tatami::create_container_of_Index_size<std::vector<std::vector<Category> > >(nchunks);
        for (auto& x : max_per_chunk) {
            tatami::resize_container_to_Index_size(x, NR);
        }
 
        auto num_per_chunk = tatami::create_container_of_Index_size<std::vector<std::vector<IndexIn_> > >(nchunks);
        for (auto& x : num_per_chunk) {
            tatami::resize_container_to_Index_size(x, NR);
        }
 
        if (mat.sparse()) {
            tatami::parallelize([&](const int, const IndexIn_ start, const IndexIn_ length) -> void {
                auto ext = tatami::consecutive_extractor<true>(mat, true, start, length, [&]{
                    tatami::Options opt;
                    opt.sparse_ordered_index = false;
                    return opt;
                }());
                auto dbuffer = tatami::create_container_of_Index_size<std::vector<ValueIn_> >(NC);
                auto ibuffer = tatami::create_container_of_Index_size<std::vector<IndexIn_> >(NC);
 
                for (IndexIn_ r = start, end = start + length; r < end; ++r) {
                    const auto range = ext->fetch(r, dbuffer.data(), ibuffer.data());
                    for (IndexIn_ i = 0; i < range.number; ++i) {
                        if (range.value[i]) {
                            const auto chunk = range.index[i] / chunk_size;
                            const auto cat = categorize(range.value[i]);
                            max_per_chunk[chunk][r] = std::max(max_per_chunk[chunk][r], cat);
                            ++num_per_chunk[chunk][r];
                        }
                    }
                }
            }, NR, nthreads);
 
        } else {
            tatami::parallelize([&](const int, const IndexIn_ start, const IndexIn_ length) -> void {
                auto ext = tatami::consecutive_extractor<false>(mat, true, start, length);
                auto dbuffer = tatami::create_container_of_Index_size<std::vector<ValueIn_> >(NC);
 
                for (IndexIn_ r = start, end = start + length; r < end; ++r) {
                    auto ptr = ext->fetch(r, dbuffer.data());
                    for (IndexIn_ c = 0; c < NC; ++c) {
                        if (ptr[c]) {
                            const auto chunk = c / chunk_size;
                            const auto cat = categorize(ptr[c]);
                            max_per_chunk[chunk][r] = std::max(max_per_chunk[chunk][r], cat);
                            ++num_per_chunk[chunk][r];
                        }
                    }
                }
            }, NR, nthreads);
        }
 
        allocate_rows(
            max_per_chunk, 
            num_per_chunk, 
            identities8, 
            identities16, 
            identities32, 
            store8, 
            store16, 
            store32, 
            assigned_category, 
            assigned_position
        );
    }
 
    // Second pass to actually fill the vectors.
    {
        tatami::parallelize([&](const int, const IndexIn_ start, const IndexIn_ length) -> void {
            auto output_positions = tatami::create_container_of_Index_size<std::vector<std::size_t> >(nchunks);
            auto dbuffer = tatami::create_container_of_Index_size<std::vector<ValueIn_> >(NC);
 
            if (mat.sparse()) {
                auto ibuffer = tatami::create_container_of_Index_size<std::vector<IndexIn_> >(NC);
                auto ext = tatami::consecutive_extractor<true>(mat, true, start, length);
 
                for (IndexIn_ r = start, end = start + length; r < end; ++r) {
                    for (I<decltype(nchunks)> chunk = 0; chunk < nchunks; ++chunk) {
                        output_positions[chunk] = get_sparse_ptr(store8, store16, store32, assigned_category, assigned_position, chunk, r);
                    }
 
                    auto range = ext->fetch(r, dbuffer.data(), ibuffer.data());
                    for (IndexIn_ i = 0; i < range.number; ++i) {
                        if (range.value[i]) {
                            const IndexIn_ chunk = range.index[i] / chunk_size;
                            const IndexIn_ col = range.index[i] % chunk_size;
                            fill_sparse_value(store8, store16, store32, assigned_category[chunk][r], chunk, col, range.value[i], output_positions[chunk]++);
                        }
                    }
                }
 
            } else {
                auto ext = tatami::consecutive_extractor<false>(mat, true, start, length);
 
                for (IndexIn_ r = start, end = start + length; r < end; ++r) {
                    for (I<decltype(nchunks)> chunk = 0; chunk < nchunks; ++chunk) {
                        output_positions[chunk] = get_sparse_ptr(store8, store16, store32, assigned_category, assigned_position, chunk, r);
                    }
 
                    auto ptr = ext->fetch(r, dbuffer.data());
                    for (IndexIn_ c = 0; c < NC; ++c) {
                        if (ptr[c]) {
                            const IndexIn_ chunk = c / chunk_size;
                            const IndexIn_ col = c % chunk_size;
                            fill_sparse_value(store8, store16, store32, assigned_category[chunk][r], chunk, col, ptr[c], output_positions[chunk]++);
                        }
                    }
                }
            }
 
        }, NR, nthreads);
    }
 
    return consolidate_matrices<ValueOut_, IndexOut_>(
        identities8, 
        identities16, 
        identities32, 
        std::move(store8), 
        std::move(store16), 
        std::move(store32),
        NR,
        chunk_size,
        leftovers
    );
}
 
template<typename Value_, typename IndexIn_>
std::vector<std::vector<Value_> > convert_by_column_allocate_store_per_chunk(const IndexIn_ nchunks, const IndexIn_ NR) {
    auto output = tatami::create_container_of_Index_size<std::vector<std::vector<Value_> > >(nchunks);
    for (auto& x : output) {
        tatami::resize_container_to_Index_size(x, NR);
    }
    return output;
};
 
template<typename Value_, typename IndexIn_>
std::vector<std::vector<Value_> >& convert_by_column_get_store_per_chunk(
    int thread,
    std::vector<std::vector<Value_> >& base,
    std::optional<std::vector<std::vector<Value_> > >& tmp,
    const IndexIn_ nchunks,
    const IndexIn_ NR
) {
    if (thread) {
        tmp = convert_by_column_allocate_store_per_chunk<Value_>(nchunks, NR);
        return *tmp;
    } else {
        return base;
    }
}
 
template<typename Value_>
void convert_by_column_save_store_per_chunk(
    int thread,
    std::optional<std::vector<std::vector<Value_> > >& tmp,
    std::vector<std::vector<std::vector<Value_> > >& collected
) { 
    if (thread) {
        collected[thread - 1] = std::move(*tmp);
    }
};
 
template<typename ColIndex_, typename ValueOut_ = double, typename IndexOut_ = int, typename ValueIn_, typename IndexIn_>
std::shared_ptr<tatami::Matrix<ValueOut_, IndexOut_> > convert_by_column(const tatami::Matrix<ValueIn_, IndexIn_>& mat, const IndexIn_ chunk_size, const int nthreads) {
    const auto NR = mat.nrow(), NC = mat.ncol();
    const IndexIn_ leftovers = NC % chunk_size;
    const IndexIn_ nchunks = sanisizer::max(1, NC / chunk_size + (leftovers != 0));
 
    auto store8  = tatami::create_container_of_Index_size<std::vector<Holder< std::uint8_t, IndexOut_, ColIndex_> > >(nchunks);
    auto store16 = tatami::create_container_of_Index_size<std::vector<Holder<std::uint16_t, IndexOut_, ColIndex_> > >(nchunks);
    auto store32 = tatami::create_container_of_Index_size<std::vector<Holder<std::uint32_t, IndexOut_, ColIndex_> > >(nchunks);
 
    auto identities8  = tatami::create_container_of_Index_size<std::vector<std::vector<IndexOut_> > >(nchunks);
    auto identities16 = tatami::create_container_of_Index_size<std::vector<std::vector<IndexOut_> > >(nchunks);
    auto identities32 = tatami::create_container_of_Index_size<std::vector<std::vector<IndexOut_> > >(nchunks);
 
    auto assigned_position = tatami::create_container_of_Index_size<std::vector<std::vector<IndexOut_> > >(nchunks);
    auto assigned_category = tatami::create_container_of_Index_size<std::vector<std::vector<Category> > >(nchunks);
 
    // First pass to define the allocations.
    {
        auto max_per_chunk = convert_by_column_allocate_store_per_chunk<Category>(nchunks, NR);
        auto max_per_chunk_threaded = sanisizer::create<std::vector<std::vector<std::vector<Category> > > >(nthreads - 1);
        auto num_per_chunk = convert_by_column_allocate_store_per_chunk<IndexIn_>(nchunks, NR);
        auto num_per_chunk_threaded = sanisizer::create<std::vector<std::vector<std::vector<IndexIn_> > > >(nthreads - 1);
        int num_used;
 
        if (mat.sparse()) {
            num_used = tatami::parallelize([&](const int t, const IndexIn_ start, const IndexIn_ length) -> void {
                auto ext = tatami::consecutive_extractor<true>(mat, false, start, length, [&]{
                    tatami::Options opt;
                    opt.sparse_ordered_index = false;
                    return opt;
                }());
                auto dbuffer = tatami::create_container_of_Index_size<std::vector<ValueIn_> >(NR);
                auto ibuffer = tatami::create_container_of_Index_size<std::vector<IndexIn_> >(NR);
 
                std::optional<std::vector<std::vector<Category> > > max_tmp;
                auto& cur_max_per_chunk = convert_by_column_get_store_per_chunk(t, max_per_chunk, max_tmp, nchunks, NR);
                std::optional<std::vector<std::vector<IndexIn_> > > num_tmp;
                auto& cur_num_per_chunk = convert_by_column_get_store_per_chunk(t, num_per_chunk, num_tmp, nchunks, NR);
 
                for (IndexIn_ c = start, end = start + length; c < end; ++c) {
                    const auto range = ext->fetch(c, dbuffer.data(), ibuffer.data());
                    const auto chunk = c / chunk_size;
                    auto& max_vec = cur_max_per_chunk[chunk];
                    auto& num_vec = cur_num_per_chunk[chunk];
 
                    for (IndexIn_ i = 0; i < range.number; ++i) {
                        if (range.value[i]) {
                            const auto cat = categorize(range.value[i]);
                            const auto r = range.index[i];
                            max_vec[r] = std::max(max_vec[r], cat);
                            ++num_vec[r];
                        }
                    }
                }
 
                convert_by_column_save_store_per_chunk(t, max_tmp, max_per_chunk_threaded);
                convert_by_column_save_store_per_chunk(t, num_tmp, num_per_chunk_threaded);
            }, NC, nthreads);
 
        } else {
            num_used = tatami::parallelize([&](const int t, const IndexIn_ start, const IndexIn_ length) -> void {
                auto ext = tatami::consecutive_extractor<false>(mat, false, start, length);
                auto dbuffer = tatami::create_container_of_Index_size<std::vector<ValueIn_> >(NR);
 
                std::optional<std::vector<std::vector<Category> > > max_tmp;
                auto& cur_max_per_chunk = convert_by_column_get_store_per_chunk(t, max_per_chunk, max_tmp, nchunks, NR);
                std::optional<std::vector<std::vector<IndexIn_> > > num_tmp;
                auto& cur_num_per_chunk = convert_by_column_get_store_per_chunk(t, num_per_chunk, num_tmp, nchunks, NR);
 
                for (IndexIn_ c = start, end = start + length; c < end; ++c) {
                    const auto ptr = ext->fetch(c, dbuffer.data());
                    const auto chunk = c / chunk_size;
                    auto& max_vec = cur_max_per_chunk[chunk];
                    auto& num_vec = cur_num_per_chunk[chunk];
 
                    for (IndexIn_ r = 0; r < NR; ++r) {
                        if (ptr[r]) {
                            auto cat = categorize(ptr[r]);
                            max_vec[r] = std::max(max_vec[r], cat);
                            ++num_vec[r];
                        }
                    }
                }
 
                convert_by_column_save_store_per_chunk(t, max_tmp, max_per_chunk_threaded);
                convert_by_column_save_store_per_chunk(t, num_tmp, num_per_chunk_threaded);
            }, NC, nthreads);
        }
 
        for (int t = 1; t < num_used; ++t) {
            const auto& cur_max_per_chunk = max_per_chunk_threaded[t - 1];
            const auto& cur_num_per_chunk = num_per_chunk_threaded[t - 1];
            for (I<decltype(nchunks)> chunk = 0; chunk < nchunks; ++chunk) {
                for (IndexIn_ r = 0; r < NR; ++r) {
                    max_per_chunk[chunk][r] = std::max(max_per_chunk[chunk][r], cur_max_per_chunk[chunk][r]);
                    num_per_chunk[chunk][r] += cur_num_per_chunk[chunk][r];
                }
            }
        }
 
        allocate_rows(
            max_per_chunk, 
            num_per_chunk, 
            identities8, 
            identities16, 
            identities32, 
            store8, 
            store16, 
            store32, 
            assigned_category, 
            assigned_position
        );
    }
 
    // Second pass to actually fill the vectors.
    {
        tatami::parallelize([&](const int, const IndexIn_ start, const IndexIn_ length) -> void {
            auto output_positions = tatami::create_container_of_Index_size<std::vector<std::vector<std::size_t> > >(nchunks);
            for (I<decltype(nchunks)> chunk = 0; chunk < nchunks; ++chunk) {
                tatami::resize_container_to_Index_size(output_positions[chunk], length);
                for (IndexIn_ r = 0; r < length; ++r) {
                    output_positions[chunk][r] = get_sparse_ptr(store8, store16, store32, assigned_category, assigned_position, chunk, r + start);
                }
            }
 
            auto dbuffer = tatami::create_container_of_Index_size<std::vector<ValueIn_> >(length);
 
            if (mat.sparse()) {
                auto ibuffer = tatami::create_container_of_Index_size<std::vector<IndexIn_> >(length);
                auto ext = tatami::consecutive_extractor<true>(mat, false, static_cast<IndexIn_>(0), NC, start, length);
 
                for (IndexIn_ c = 0; c < NC; ++c) {
                    const auto range = ext->fetch(c, dbuffer.data(), ibuffer.data());
                    const auto chunk = c / chunk_size;
                    const IndexIn_ col = c % chunk_size;
                    auto& outpos = output_positions[chunk];
 
                    for (IndexIn_ i = 0; i < range.number; ++i) {
                        if (range.value[i]) {
                            const auto r = range.index[i];
                            fill_sparse_value(store8, store16, store32, assigned_category[chunk][r], chunk, col, range.value[i], outpos[r - start]++);
                        }
                    }
                }
 
            } else {
                auto ext = tatami::consecutive_extractor<false>(mat, false, static_cast<IndexIn_>(0), NC, start, length);
 
                for (IndexIn_ c = 0; c < NC; ++c) {
                    const auto ptr = ext->fetch(c, dbuffer.data());
                    const auto chunk = c / chunk_size;
                    const IndexIn_ col = c % chunk_size;
                    auto& outpos = output_positions[chunk];
 
                    for (IndexIn_ r = 0; r < NR; ++r) {
                        if (ptr[r]) {
                            fill_sparse_value(store8, store16, store32, assigned_category[chunk][r], chunk, col, ptr[r], outpos[r - start]++);
                        }
                    }
                }
            }
 
        }, NR, nthreads);
    }
 
    return consolidate_matrices<ValueOut_, IndexOut_>(
        identities8, 
        identities16, 
        identities32, 
        std::move(store8), 
        std::move(store16), 
        std::move(store32),
        NR,
        chunk_size,
        leftovers
    );
}
struct ConvertToLayeredSparseOptions {
    std::size_t chunk_size = sanisizer::cap<std::size_t>(65536);
 
    int num_threads = 1;
};
 
template<typename ValueOut_ = double, typename IndexOut_ = int, typename ColumnIndex_ = std::uint16_t, typename ValueIn_, typename IndexIn_>
std::shared_ptr<tatami::Matrix<ValueOut_, IndexOut_> > convert_to_layered_sparse(const tatami::Matrix<ValueIn_, IndexIn_>& mat, const ConvertToLayeredSparseOptions& options) {
    const IndexIn_ chunk_size = check_chunk_size<IndexIn_, ColumnIndex_>(options.chunk_size);
    if (mat.prefer_rows()) {
        return convert_by_row<ColumnIndex_, ValueOut_, IndexOut_>(mat, chunk_size, options.num_threads);
    } else {
        return convert_by_column<ColumnIndex_, ValueOut_, IndexOut_>(mat, chunk_size, options.num_threads);
    }
}
 
// Provided for back-compatibility.
template<typename ValueOut_ = double, typename IndexOut_ = int, typename ColumnIndex_ = std::uint16_t, typename ValueIn_, typename IndexIn_>
std::shared_ptr<tatami::Matrix<ValueOut_, IndexOut_> > convert_to_layered_sparse(const tatami::Matrix<ValueIn_, IndexIn_>& mat, IndexIn_ chunk_size = 65536, int num_threads = 1) {
    return convert_to_layered_sparse(mat, [&]{
        ConvertToLayeredSparseOptions opt;
        opt.chunk_size = chunk_size;
        opt.num_threads = num_threads;
        return opt;
    }());
}
 
template<typename ValueOut_ = double, typename IndexOut_ = int, typename ColumnIndex_ = std::uint16_t, typename ValueIn_, typename IndexIn_>
std::shared_ptr<tatami::Matrix<ValueOut_, IndexOut_> > convert_to_layered_sparse(const tatami::Matrix<ValueIn_, IndexIn_>* mat, IndexIn_ chunk_size = 65536, int num_threads = 1) {
    return convert_to_layered_sparse<ValueOut_, IndexOut_, ColumnIndex_, ValueIn_, IndexIn_>(*mat, chunk_size, num_threads);
}
}
 
#endif