read_layered_sparse_from_matrix_market.hpp

Go to the documentation of this file.

#ifndef TATAMI_LAYERED_READ_LAYERED_SPARSE_FROM_MATRIX_MARKET_HPP
#define TATAMI_LAYERED_READ_LAYERED_SPARSE_FROM_MATRIX_MARKET_HPP
 
#include <vector>
#include <algorithm>
#include <cstddef>
 
#include "byteme/byteme.hpp"
#include "eminem/eminem.hpp"
#include "tatami/tatami.hpp"
#include "sanisizer/sanisizer.hpp"
 
#include "utils.hpp"
 
namespace tatami_layered {
 
template<typename Value_, typename Index_, typename ColumnIndex_, class Creator_>
std::shared_ptr<tatami::Matrix<Value_, Index_> > read_layered_sparse_from_matrix_market(Creator_ create, const Index_ chunk_size, const int num_threads) {
    Index_ NR, NC, nchunks, leftovers;
 
    std::vector<Holder< std::uint8_t, Index_, ColumnIndex_> > store8;
    std::vector<Holder<std::uint16_t, Index_, ColumnIndex_> > store16;
    std::vector<Holder<std::uint32_t, Index_, ColumnIndex_> > store32;
 
    std::vector<std::vector<Index_> > identities8, identities16, identities32;
    std::vector<std::vector<Index_> > assigned_position;
    std::vector<std::vector<Category> > assigned_category;
 
    eminem::ParserOptions eopt;
    eopt.num_threads = num_threads;
 
    // First pass, scanning for the max and number.
    {
        auto reader = create();
        byteme::PerByteSerial<char, byteme::Reader*> pb(&reader);
        eminem::Parser<I<decltype(&pb)>, Index_> parser(&pb, eopt);
 
        parser.scan_preamble();
        NR = parser.get_nrows();
        NC = parser.get_ncols();
        leftovers = NC % chunk_size;
        nchunks = sanisizer::max(1, NC / chunk_size + (leftovers != 0));
 
        tatami::resize_container_to_Index_size(store8, nchunks);
        tatami::resize_container_to_Index_size(store16, nchunks);
        tatami::resize_container_to_Index_size(store32, nchunks);
        tatami::resize_container_to_Index_size(identities8, nchunks);
        tatami::resize_container_to_Index_size(identities16, nchunks);
        tatami::resize_container_to_Index_size(identities32, nchunks);
        tatami::resize_container_to_Index_size(assigned_position, nchunks);
        tatami::resize_container_to_Index_size(assigned_category, nchunks);
 
        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<Index_> > >(nchunks);
        for (auto& x : num_per_chunk) {
            tatami::resize_container_to_Index_size(x, NR);
        }
 
        auto handler = [&](const Index_ r, const Index_ c, const Category cat) -> void {
            const auto chunk = (c - 1) / chunk_size;
            auto& maxcat = max_per_chunk[chunk][r - 1];
            maxcat = std::max(maxcat, cat);
            ++num_per_chunk[chunk][r - 1];
        };
 
        const auto& banner = parser.get_banner();
        if (banner.field == eminem::Field::INTEGER) {
            parser.template scan_integer<std::uint32_t>([&](const Index_ r, const Index_ c, const std::uint32_t val) -> void {
                handler(r, c, categorize(val));
            });
        } else if (banner.field == eminem::Field::DOUBLE || banner.field == eminem::Field::REAL) {
            parser.scan_real([&](const Index_ r, const Index_ c, const double val) -> void {
                handler(r, c, categorize(val));
            });
        } else {
            throw std::runtime_error("expected a numeric field in the Matrix Market file");
        }
 
        allocate_rows(
            max_per_chunk, 
            num_per_chunk, 
            identities8, 
            identities16, 
            identities32, 
            store8, 
            store16, 
            store32, 
            assigned_category, 
            assigned_position
        );
    }
 
    // Now allocating.
    {
        std::vector<std::vector<std::size_t> > output_positions(nchunks);
        for (I<decltype(nchunks)> chunk = 0; chunk < nchunks; ++chunk) {
            tatami::resize_container_to_Index_size(output_positions[chunk], NR);
            for (I<decltype(NR)> r = 0; r < NR; ++r) {
                output_positions[chunk][r] = get_sparse_ptr(store8, store16, store32, assigned_category, assigned_position, chunk, r);
            }
        }
 
        auto reader = create();
        byteme::PerByteSerial<char, byteme::Reader*> pb(&reader);
        eminem::Parser<I<decltype(&pb)>, Index_> parser(&pb, eopt);
 
        auto handler = [&](Index_ r, Index_ c, const auto val) -> void {
            --c;
            const Index_ chunk = c / chunk_size;
            const Index_ offset = c % chunk_size;
            --r;
            fill_sparse_value(store8, store16, store32, assigned_category[chunk][r], chunk, offset, val, output_positions[chunk][r]++);
        };
 
        parser.scan_preamble();
        const auto& banner = parser.get_banner();
        if (banner.field == eminem::Field::INTEGER) {
            parser.template scan_integer<std::uint32_t>([&](const Index_ r, const Index_ c, const std::uint32_t val) -> void {
                handler(r, c, val);
            });
        } else if (banner.field == eminem::Field::DOUBLE || banner.field == eminem::Field::REAL) {
            parser.scan_real([&](const Index_ r, const Index_ c, const double val) -> void {
                handler(r, c, val);
            });
        }
 
        // Checking that the column indices are sorted properly.
        auto sorter = [&](auto& store) -> void {
            std::vector<std::pair<I<decltype(store[0].index[0])>, I<decltype(store[0].value[0])>> > buffer;
            buffer.reserve(chunk_size);
 
            for (auto& st : store) {
                const auto num_ptr = st.ptr.size();
                for (I<decltype(num_ptr)> r = 1; r < num_ptr; ++r) {
                    const auto start = st.ptr[r - 1], end = st.ptr[r];
 
                    if (!std::is_sorted(st.index.begin() + start, st.index.begin() + end)) {
                        buffer.clear();
                        for (auto i = start; i < end; ++i) {
                            buffer.emplace_back(st.index[i], st.value[i]);
                        }
 
                        std::sort(buffer.begin(), buffer.end());
                        auto bIt = buffer.begin();
                        for (auto i = start; i < end; ++i, ++bIt) {
                            st.index[i] = bIt->first;
                            st.value[i] = bIt->second;
                        }
                    }
                }
            }
        };
 
        sorter(store8);
        sorter(store16);
        sorter(store32);
    }
 
    return consolidate_matrices<Value_, Index_>(
        identities8, 
        identities16, 
        identities32, 
        std::move(store8), 
        std::move(store16), 
        std::move(store32),
        NR,
        chunk_size,
        leftovers
    );
}
struct ReadLayeredSparseFromMatrixMarketOptions {
    std::size_t chunk_size = sanisizer::cap<std::size_t>(65536);
 
    std::size_t buffer_size = sanisizer::cap<std::size_t>(65536);
 
    int num_threads = 1;
};
 
template<typename Value_ = double, typename Index_ = int, typename ColumnIndex_ = std::uint16_t>
std::shared_ptr<tatami::Matrix<Value_, Index_> > read_layered_sparse_from_matrix_market_text_file(const char* filepath, const ReadLayeredSparseFromMatrixMarketOptions& options) {
    return read_layered_sparse_from_matrix_market<Value_, Index_, ColumnIndex_>(
        [&]() -> auto {
            return byteme::RawFileReader(filepath, [&]{
                byteme::RawFileReaderOptions opt;
                opt.buffer_size = options.buffer_size;
                return opt;
            }());
        },
        check_chunk_size<Index_, ColumnIndex_>(options.chunk_size),
        options.num_threads
    );
}
 
// Back-compatibility.
template<typename Value_ = double, typename Index_ = int, typename ColumnIndex_ = std::uint16_t>
std::shared_ptr<tatami::Matrix<Value_, Index_> > read_layered_sparse_from_matrix_market_text_file(const char* filepath, Index_ chunk_size = 65536, std::size_t buffer_size = 65536) {
    return read_layered_sparse_from_matrix_market_text_file<Value_, Index_, ColumnIndex_>(filepath, [&]{
        ReadLayeredSparseFromMatrixMarketOptions opt;
        opt.chunk_size = chunk_size;
        opt.buffer_size = buffer_size;
        return opt;
    }());
}
#if __has_include("zlib.h")
 
template<typename Value_ = double, typename Index_ = int, typename ColumnIndex_ = std::uint16_t>
std::shared_ptr<tatami::Matrix<Value_, Index_> > read_layered_sparse_from_matrix_market_some_file(const char* filepath, const ReadLayeredSparseFromMatrixMarketOptions& options) {
    return read_layered_sparse_from_matrix_market<Value_, Index_, ColumnIndex_>(
        [&]() -> auto {
            return byteme::SomeFileReader(filepath, [&]{
                byteme::SomeFileReaderOptions opt;
                opt.buffer_size = options.buffer_size;
                return opt;
            }());
        },
        check_chunk_size<Index_, ColumnIndex_>(options.chunk_size),
        options.num_threads
    );
}
 
template<typename Value_ = double, typename Index_ = int, typename ColumnIndex_ = std::uint16_t>
std::shared_ptr<tatami::Matrix<Value_, Index_> > read_layered_sparse_from_matrix_market_gzip_file(const char* filepath, const ReadLayeredSparseFromMatrixMarketOptions& options) {
    return read_layered_sparse_from_matrix_market<Value_, Index_, ColumnIndex_>(
        [&]() -> auto {
            return byteme::GzipFileReader(filepath, [&]{
                byteme::GzipFileReaderOptions opt;
                opt.buffer_size = options.buffer_size;
                return opt;
            }());
        },
        check_chunk_size<Index_, ColumnIndex_>(options.chunk_size),
        options.num_threads
    );
}
 
// Back-compatibility.
template<typename Value_ = double, typename Index_ = int, typename ColumnIndex_ = std::uint16_t>
std::shared_ptr<tatami::Matrix<Value_, Index_> > read_layered_sparse_from_matrix_market_some_file(const char* filepath, Index_ chunk_size = 65536, std::size_t buffer_size = 65536) {
    return read_layered_sparse_from_matrix_market_some_file<Value_, Index_, ColumnIndex_>(filepath, [&]{
        ReadLayeredSparseFromMatrixMarketOptions opt;
        opt.chunk_size = chunk_size;
        opt.buffer_size = buffer_size;
        return opt;
    }());
}
 
template<typename Value_ = double, typename Index_ = int, typename ColumnIndex_ = std::uint16_t>
std::shared_ptr<tatami::Matrix<Value_, Index_> > read_layered_sparse_from_matrix_market_gzip_file(const char* filepath, Index_ chunk_size = 65536, std::size_t buffer_size = 65536) {
    return read_layered_sparse_from_matrix_market_gzip_file<Value_, Index_, ColumnIndex_>(filepath, [&]{
        ReadLayeredSparseFromMatrixMarketOptions opt;
        opt.chunk_size = chunk_size;
        opt.buffer_size = buffer_size;
        return opt;
    }());
}
#endif
 
template<typename Value_ = double, typename Index_ = int, typename ColumnIndex_ = std::uint16_t>
std::shared_ptr<tatami::Matrix<Value_, Index_> > read_layered_sparse_from_matrix_market_text_buffer(
    const unsigned char* contents,
    std::size_t length,
    const ReadLayeredSparseFromMatrixMarketOptions& options)
{
    return read_layered_sparse_from_matrix_market<Value_, Index_, ColumnIndex_>(
        [&]() -> auto {
            return byteme::RawBufferReader(contents, length);
        },
        check_chunk_size<Index_, ColumnIndex_>(options.chunk_size),
        options.num_threads
    );
}
 
// Back-compatibility.
template<typename Value_ = double, typename Index_ = int, typename ColumnIndex_ = std::uint16_t>
std::shared_ptr<tatami::Matrix<Value_, Index_> > read_layered_sparse_from_matrix_market_text_buffer(const unsigned char* contents, std::size_t length, Index_ chunk_size = 65536) {
    return read_layered_sparse_from_matrix_market_text_buffer<Value_, Index_, ColumnIndex_>(contents, length, [&]{
        ReadLayeredSparseFromMatrixMarketOptions opt;
        opt.chunk_size = chunk_size;
        return opt;
    }());
}
#if __has_include("zlib.h")
 
template<typename Value_ = double, typename Index_ = int, typename ColumnIndex_ = std::uint16_t>
std::shared_ptr<tatami::Matrix<Value_, Index_> > read_layered_sparse_from_matrix_market_some_buffer(
    const unsigned char* contents,
    std::size_t length,
    const ReadLayeredSparseFromMatrixMarketOptions& options)
{
    return read_layered_sparse_from_matrix_market<Value_, Index_, ColumnIndex_>(
        [&]() -> auto {
            return byteme::SomeBufferReader(contents, length, [&]{
                byteme::SomeBufferReaderOptions opt;
                opt.buffer_size = options.buffer_size;
                return opt;
            }());
        },
        check_chunk_size<Index_, ColumnIndex_>(options.chunk_size),
        options.num_threads
    );
}
 
template<typename Value_ = double, typename Index_ = int, typename ColumnIndex_ = std::uint16_t>
std::shared_ptr<tatami::Matrix<Value_, Index_> > read_layered_sparse_from_matrix_market_zlib_buffer(
    const unsigned char* contents,
    std::size_t length,
    const ReadLayeredSparseFromMatrixMarketOptions& options)
{
    return read_layered_sparse_from_matrix_market<Value_, Index_, ColumnIndex_>(
        [&]() -> auto {
            return byteme::ZlibBufferReader(contents, length, [&]{
                byteme::ZlibBufferReaderOptions opt;
                opt.buffer_size = options.buffer_size;
                return opt;
            }());
        },
        check_chunk_size<Index_, ColumnIndex_>(options.chunk_size),
        options.num_threads
    );
}
 
// Back-compatibility.
template<typename Value_ = double, typename Index_ = int, typename ColumnIndex_ = std::uint16_t>
std::shared_ptr<tatami::Matrix<Value_, Index_> > read_layered_sparse_from_matrix_market_some_buffer(
    const unsigned char* contents,
    std::size_t length,
    Index_ chunk_size = 65536,
    std::size_t buffer_size = 65536)
{
    return read_layered_sparse_from_matrix_market_some_buffer<Value_, Index_, ColumnIndex_>(contents, length, [&]{
        ReadLayeredSparseFromMatrixMarketOptions opt;
        opt.chunk_size = chunk_size;
        opt.buffer_size = buffer_size;
        return opt;
    }());
}
 
template<typename Value_ = double, typename Index_ = int, typename ColumnIndex_ = std::uint16_t>
std::shared_ptr<tatami::Matrix<Value_, Index_> > read_layered_sparse_from_matrix_market_zlib_buffer(
    const unsigned char* contents,
    std::size_t length,
    Index_ chunk_size = 65536,
    std::size_t buffer_size = 65536)
{
    return read_layered_sparse_from_matrix_market_zlib_buffer<Value_, Index_, ColumnIndex_>(contents, length, [&]{
        ReadLayeredSparseFromMatrixMarketOptions opt;
        opt.chunk_size = chunk_size;
        opt.buffer_size = buffer_size;
        return opt;
    }());
}
#endif
 
}
 
#endif