fastx_sequence_file.hpp

Go to the documentation of this file.

#pragma once
 
#include <libdeflate.h>
 
#include <algorithm>
#include <array>
#include <bit>
#include <condition_variable>
#include <cstdint>
#include <cstdlib>
#include <cstring>
#include <deque>
#include <filesystem>
#include <format>
#include <fstream>
#include <memory>
#include <mutex>
#include <optional>
#include <span>
#include <string>
#include <string_view>
#include <thread>
#include <utility>
#include <vector>
 
#include <cuddl/detail/ablation.hpp>
#include <cuddl/detail/dna.hpp>
#include <cuddl/error.hpp>
 
#if defined(__aarch64__) || defined(_M_ARM64)
    #include <arm_neon.h>
#endif
#if defined(__ARM_FEATURE_SVE)
    #include <arm_sve.h>
#elif defined(__x86_64__) || defined(_M_X64)
    #include <immintrin.h>
#endif
 
#if defined(__linux__)
    #include <fcntl.h>
    #include <sys/mman.h>
    #include <sys/stat.h>
    #include <unistd.h>
#endif
 
namespace cuddl::detail {
 
struct fastx_sequence_extent {
    char const* begin;
    char const* end;
};
 
struct fastx_packed_record {
    uint64_t word;
    uint64_t bases;
};
 
[[nodiscard]] inline constexpr bool fastx_is_sequence_whitespace(char ch) noexcept {
    return ch == '\n' || ch == '\r' || ch == ' ' || ch == '\t';
}
 
// SIMD line scanners: https://github.com/tdortman/cuSBF
// Source revision: ec678b0750b0067d4fc66e066c8cd7a67cc5e792
/*
Boost Software License - Version 1.0 - August 17th, 2003
 
Permission is hereby granted, free of charge, to any person or organization
obtaining a copy of the software and accompanying documentation covered by
this license (the "Software") to use, reproduce, display, distribute,
execute, and transmit the Software, and to prepare derivative works of the
Software, and to permit third-parties to whom the Software is furnished to
do so, all subject to the following:
 
The copyright notices in the Software and this entire statement, including
the above license grant, this restriction and the following disclaimer,
must be included in all copies of the Software, in whole or in part, and
all derivative works of the Software, unless such copies or derivative
works are solely in the form of machine-executable object code generated by
a source language processor.
 
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE, TITLE AND NON-INFRINGEMENT. IN NO EVENT
SHALL THE COPYRIGHT HOLDERS OR ANYONE DISTRIBUTING THE SOFTWARE BE LIABLE
FOR ANY DAMAGES OR OTHER LIABILITY, WHETHER IN CONTRACT, TORT OR OTHERWISE,
ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
DEALINGS IN THE SOFTWARE.
*/
[[nodiscard]] inline size_t fastx_line_end_scalar(std::string_view data, size_t position) {
    while (position < data.size() && data[position] != '\n' && data[position] != '\r') {
        ++position;
    }
    return position;
}
 
#if (defined(__GNUC__) || defined(__clang__)) && defined(__x86_64__)
[[nodiscard]] __attribute__((target("avx2"))) inline size_t
fastx_line_end_avx2(std::string_view data, size_t position) {
    const __m256i newline = _mm256_set1_epi8('\n');
    const __m256i carriage_return = _mm256_set1_epi8('\r');
    while (data.size() - position >= 32) {
        const __m256i bytes =
            _mm256_loadu_si256(reinterpret_cast<const __m256i*>(data.data() + position));
        const auto mask = static_cast<uint32_t>(_mm256_movemask_epi8(_mm256_or_si256(
            _mm256_cmpeq_epi8(bytes, newline), _mm256_cmpeq_epi8(bytes, carriage_return)
        )));
        if (mask != 0) {
            return position + static_cast<size_t>(__builtin_ctz(mask));
        }
        position += 32;
    }
    return fastx_line_end_scalar(data, position);
}
#endif
 
#if (defined(__GNUC__) || defined(__clang__)) && defined(__x86_64__)
[[nodiscard]] __attribute__((target("avx512f,avx512bw"))) inline size_t
fastx_line_end_avx512(std::string_view data, size_t position) {
    while (data.size() - position >= 64) {
        auto const bytes = _mm512_loadu_si512(data.data() + position);
        auto const mask = _mm512_cmpeq_epi8_mask(bytes, _mm512_set1_epi8('\n')) |
                          _mm512_cmpeq_epi8_mask(bytes, _mm512_set1_epi8('\r'));
        if (mask != 0) return position + std::countr_zero(mask);
        position += 64;
    }
    return fastx_line_end_scalar(data, position);
}
#endif
 
#if defined(__ARM_FEATURE_SVE)
[[nodiscard]] inline size_t fastx_line_end_sve(std::string_view data, size_t position) {
    while (position < data.size()) {
        const svbool_t active = svwhilelt_b8(position, data.size());
        const svuint8_t bytes =
            svld1_u8(active, reinterpret_cast<const uint8_t*>(data.data() + position));
        const svbool_t matches =
            svorr_b_z(active, svcmpeq_n_u8(active, bytes, '\n'), svcmpeq_n_u8(active, bytes, '\r'));
        if (svptest_any(active, matches)) {
            return position + svcntp_b8(active, svbrkb_z(active, matches));
        }
        position += svcntb();
    }
    return position;
}
#endif
 
#if defined(__aarch64__) || defined(_M_ARM64)
[[nodiscard]] inline size_t fastx_line_end_neon(std::string_view data, size_t position) {
    const uint8x16_t newline = vdupq_n_u8('\n');
    const uint8x16_t carriage_return = vdupq_n_u8('\r');
    while (data.size() - position >= 16) {
        const uint8x16_t bytes = vld1q_u8(reinterpret_cast<const uint8_t*>(data.data() + position));
        const uint8x16_t matches =
            vorrq_u8(vceqq_u8(bytes, newline), vceqq_u8(bytes, carriage_return));
        const uint64x2_t words = vreinterpretq_u64_u8(matches);
        const uint64_t low = vgetq_lane_u64(words, 0);
        if (low != 0) {
            return position + (std::countr_zero(low) >> 3);
        }
        const uint64_t high = vgetq_lane_u64(words, 1);
        if (high != 0) {
            return position + 8 + (std::countr_zero(high) >> 3);
        }
        position += 16;
    }
    return fastx_line_end_scalar(data, position);
}
#endif
 
using fastx_line_end_fn = decltype(&fastx_line_end_scalar);
 
[[nodiscard]] inline fastx_line_end_fn resolve_fastx_line_end() {
    if (ablated("simd-parsing")) return fastx_line_end_scalar;
#if defined(__ARM_FEATURE_SVE)
    return fastx_line_end_sve;
#elif defined(__aarch64__) || defined(_M_ARM64)
    return fastx_line_end_neon;
#elif (defined(__GNUC__) || defined(__clang__)) && defined(__x86_64__)
    if (__builtin_cpu_supports("avx512bw")) return fastx_line_end_avx512;
    return __builtin_cpu_supports("avx2") ? fastx_line_end_avx2 : fastx_line_end_scalar;
#else
    return fastx_line_end_scalar;
#endif
}
 
static const fastx_line_end_fn line_end = resolve_fastx_line_end();
 
[[nodiscard]] inline size_t fastx_line_end(std::string_view data, size_t position) {
    return line_end(data, position);
}
 
#if defined(__linux__)
[[nodiscard]] inline bool
read_whole_file(int fd, char* out, size_t size, bool direct = false, size_t capacity = 0) noexcept {
    constexpr size_t page = 4096;
    size_t done = 0;
    int const flags = direct && size > 64 * page ? ::fcntl(fd, F_GETFL) : -1;
    if (flags != -1 && ::fcntl(fd, F_SETFL, flags | O_DIRECT) == 0) {
        auto const whole = (size + page - 1) / page * page;
        if (reinterpret_cast<uintptr_t>(out) % page == 0) {
            // Whole pages straight into out; the last partial page too if out has room for it.
            auto const direct_end = capacity >= whole ? whole : size / page * page;
            while (done < direct_end && done < size) {
                auto const got =
                    ::pread(fd, out + done, direct_end - done, static_cast<off_t>(done));
                if (got <= 0) break;
                done = std::min(size, done + static_cast<size_t>(got));
            }
        }
        constexpr size_t bounce_bytes = size_t{8} << 20;
        static thread_local std::unique_ptr<char, decltype(&std::free)> bounce{
            static_cast<char*>(std::aligned_alloc(page, bounce_bytes)), &std::free
        };
        while (bounce && done < size && done % page == 0) {
            auto const wanted = std::min(bounce_bytes, (size - done + page - 1) / page * page);
            auto const got = ::pread(fd, bounce.get(), wanted, static_cast<off_t>(done));
            if (got <= 0) break;
            auto const used = std::min(size - done, static_cast<size_t>(got));
            std::memcpy(out + done, bounce.get(), used);
            done += used;
        }
        if (done < size) ::fcntl(fd, F_SETFL, flags);
    }
    while (done < size) {
        auto const got = ::pread(fd, out + done, size - done, static_cast<off_t>(done));
        if (got <= 0) break;
        done += static_cast<size_t>(got);
    }
    ::close(fd);
    return done == size;
}
#endif
 
class fastx_mapped_file {
   public:
    fastx_mapped_file() = default;
    fastx_mapped_file(fastx_mapped_file const&) = delete;
    fastx_mapped_file& operator=(fastx_mapped_file const&) = delete;
    fastx_mapped_file(fastx_mapped_file&& other) noexcept {
        *this = std::move(other);
    }
 
    fastx_mapped_file& operator=(fastx_mapped_file&& other) noexcept {
        if (this != &other) {
            release();
            owned_storage_ = std::move(other.owned_storage_);
            data_ = other.data_;
            size_ = other.size_;
            other.data_ = nullptr;
            other.size_ = 0;
#if defined(__linux__)
            mapped_ = other.mapped_;
            mapped_size_ = other.mapped_size_;
            other.mapped_ = nullptr;
            other.mapped_size_ = 0;
#endif
        }
        return *this;
    }
 
    ~fastx_mapped_file() {
        release();
    }
 
    [[nodiscard]] static Result<std::unique_ptr<fastx_mapped_file>>
    load(std::string const& path, bool read = false, bool direct = false) {
        auto buffer = std::make_unique<fastx_mapped_file>();
        auto loaded = buffer->load_from_path(path, read, direct);
        if (!loaded) return Err(loaded.error());
        return buffer;
    }
 
#if defined(__linux__)
    [[nodiscard]] static Result<std::unique_ptr<fastx_mapped_file>>
    adopt(int fd, size_t size, std::string const& path, bool read, bool direct = false) {
        auto buffer = std::make_unique<fastx_mapped_file>();
        auto loaded = buffer->load_from_descriptor(fd, size, path, read, direct);
        if (!loaded) return Err(loaded.error());
        return buffer;
    }
#endif
 
    [[nodiscard]] std::string_view data() const noexcept {
        return std::string_view{data_, size_};
    }
 
    [[nodiscard]] bool empty() const noexcept {
        return size_ == 0;
    }
 
   private:
    struct free_deleter {
        void operator()(char* data) const noexcept {
            std::free(data);
        }
    };
    std::unique_ptr<char, free_deleter> owned_storage_;
 
    [[nodiscard]] static size_t page_rounded(size_t size) noexcept {
        return (std::max<size_t>(size, 1) + 4095) / 4096 * 4096;
    }
 
    void allocate(size_t size) {
        owned_storage_.reset(static_cast<char*>(std::aligned_alloc(4096, page_rounded(size))));
        if (!owned_storage_) throw std::bad_alloc();
    }
    char const* data_{nullptr};
    size_t size_{0};
 
#if defined(__linux__)
    void* mapped_{nullptr};
    size_t mapped_size_{0};
#endif
 
    void release() {
#if defined(__linux__)
        if (mapped_ != nullptr) {
            munmap(mapped_, mapped_size_);
            mapped_ = nullptr;
            mapped_size_ = 0;
        }
#endif
        owned_storage_.reset();
        data_ = nullptr;
        size_ = 0;
    }
 
    [[nodiscard]] Result<void> load_from_path(
        std::string const& path,
        [[maybe_unused]] bool read,
        [[maybe_unused]] bool direct
    ) {
#if defined(__linux__)
        int const fd = ::open(path.c_str(), O_RDONLY);
        if (fd == -1) {
            return Err(Error::invalid_argument("Failed to open FASTA/FASTQ file: " + path));
        }
 
        struct stat file_status{};
        if (::fstat(fd, &file_status) != 0 || file_status.st_size < 0) {
            ::close(fd);
            return Err(Error::invalid_argument("Failed to stat FASTA/FASTQ file: " + path));
        }
        return load_from_descriptor(
            fd, static_cast<size_t>(file_status.st_size), path, read, direct
        );
#else
        std::ifstream input(path, std::ios::binary);
        if (!input.is_open()) {
            return Err(Error::invalid_argument("Failed to open FASTA/FASTQ file: " + path));
        }
        input.seekg(0, std::ios::end);
        auto const file_size = input.tellg();
        if (file_size < 0) {
            return Err(Error::invalid_argument("Failed to size FASTA/FASTQ file: " + path));
        }
        allocate(static_cast<size_t>(file_size));
        input.seekg(0, std::ios::beg);
        if (file_size > 0) {
            input.read(owned_storage_.get(), static_cast<std::streamsize>(file_size));
            if (!input) {
                return Err(Error::invalid_argument("Failed to read FASTA/FASTQ file: " + path));
            }
        }
        data_ = owned_storage_.get();
        size_ = static_cast<size_t>(file_size);
        return Ok();
#endif
    }
 
#if defined(__linux__)
    [[nodiscard]] Result<void>
    load_from_descriptor(int fd, size_t size, std::string const& path, bool read, bool direct) {
        if (size == 0) {
            ::close(fd);
            return Ok();
        }
        if (!read) {
            mapped_ = ::mmap(nullptr, size, PROT_READ, MAP_PRIVATE, fd, 0);
            ::close(fd);
            if (mapped_ == MAP_FAILED) {
                mapped_ = nullptr;
                return Err(Error::invalid_argument("Failed to mmap FASTA/FASTQ file: " + path));
            }
            mapped_size_ = size;
            data_ = static_cast<char const*>(mapped_);
            size_ = size;
            return Ok();
        }
        allocate(size);
        if (!read_whole_file(fd, owned_storage_.get(), size, direct, page_rounded(size))) {
            owned_storage_.reset();
            return Err(Error::invalid_argument("Failed to read FASTA/FASTQ file: " + path));
        }
        data_ = owned_storage_.get();
        size_ = size;
        return Ok();
    }
#endif
};
 
enum class fastx_format : uint8_t {
    unknown,
    fasta,
    fastq,
};
 
class fastx_fastq_reader {
   public:
    explicit fastx_fastq_reader(std::string_view data, std::string source_name = "<buffer>")
        : data_(data), source_name_(std::move(source_name)) {}
 
    [[nodiscard]] Result<bool> append_next_record(std::string& sequence) {
        auto header = CUDDL_TRY(read_header_line());
        if (header.empty()) {
            return false;
        }
 
        char const tag = header.front();
        if (format_ == fastx_format::unknown) {
            if (tag == '>') {
                format_ = fastx_format::fasta;
            } else if (tag == '@') {
                format_ = fastx_format::fastq;
            } else {
                return Err(parse_error("expected FASTA or FASTQ header"));
            }
        }
 
        char const expected = format_ == fastx_format::fasta ? '>' : '@';
        if (tag != expected) {
            return Err(parse_error("mixed FASTA and FASTQ records are not supported"));
        }
        if (format_ != fastx_format::fastq) {
            return Err(parse_error("mixed FASTA and FASTQ records are not supported"));
        }
        auto const base = sequence.size();
        CUDDL_TRY(read_fastq_sequence(sequence, base));
        return true;
    }
 
   private:
    std::string_view data_;
    std::string source_name_;
    size_t position_{0};
    fastx_format format_{fastx_format::unknown};
    uint64_t line_number_{};
 
    [[nodiscard]] Error parse_error(std::string_view message) const {
        return Error::invalid_argument(
            std::format("{}:{}: {}", source_name_, line_number_, message)
        );
    }
 
    [[nodiscard]] std::string_view read_line() {
        if (position_ >= data_.size()) {
            return {};
        }
 
        size_t const end = fastx_line_end(data_, position_);
 
        std::string_view const line = data_.substr(position_, end - position_);
        position_ = end;
        if (position_ < data_.size() && data_[position_] == '\r') {
            ++position_;
        }
        if (position_ < data_.size() && data_[position_] == '\n') {
            ++position_;
        }
        ++line_number_;
        return line;
    }
 
    [[nodiscard]] Result<std::string_view> read_header_line() {
        while (position_ < data_.size()) {
            std::string_view const line = read_line();
            if (!line.empty()) {
                return line;
            }
        }
        return std::string_view{};
    }
 
    [[nodiscard]] Result<void> read_fastq_sequence(std::string& sequence, size_t base) {
        while (position_ < data_.size()) {
            std::string_view const line = read_line();
            if (!line.empty() && line.front() == '+') {
                return read_fastq_qualities(sequence.size() - base);
            }
            sequence.append(line.data(), line.size());
        }
        return Err(parse_error("unterminated FASTQ record: missing '+' separator"));
    }
 
    [[nodiscard]] Result<void> read_fastq_qualities(uint64_t expected_length) {
        uint64_t quality_length = 0;
        while (quality_length < expected_length && position_ < data_.size()) {
            std::string_view const line = read_line();
            quality_length += line.size();
            if (quality_length > expected_length) {
                return Err(parse_error("FASTQ quality length exceeds sequence length"));
            }
        }
        if (quality_length != expected_length) {
            return Err(parse_error("FASTQ quality length does not match sequence length"));
        }
        return Ok();
    }
};
 
// Heap ownership keeps extent addresses stable, including short decompressed strings.
struct fastx_sequence_file {
    std::unique_ptr<fastx_mapped_file> file;
    std::string_view input;
    std::string decompressed;
    std::string sequence;
    std::vector<fastx_sequence_extent> extents;
    // Set when the caller supplied memory for the decompressed bytes, which then live outside
    // this object and stay valid as long as the caller keeps them. `decompressed` stays empty.
    char* decompressed_target = nullptr;
    size_t decompressed_size = 0;
    std::shared_ptr<void> storage_owner;
    // Set by pack_fastx_sequence_file, which consumes `extents`: runs of at least k bases as
    // two-bit words, sixteen bases each, every run starting on a word of its own.
    uint32_t const* words = nullptr;
    std::vector<fastx_packed_record> packed;
};
 
[[nodiscard]] inline uint32_t pack_bases16_portable(char const* bases, size_t count, uint32_t& word) {
    if (count < 16) {
        uint32_t bad = 0;
        word = 0;
        for (uint32_t j = 0; j < 16; ++j) {
            auto const symbol = j < count ? encode_base(bases[j]) : uint8_t{0xFF};
            word |= static_cast<uint32_t>(symbol & 3U) << (2 * j);
            bad |= static_cast<uint32_t>(symbol == 0xFFU) << j;
        }
        return bad;
    }
    uint32_t bad = 0;
    word = 0;
    for (uint32_t half = 0; half < 2; ++half) {
        uint64_t x;
        std::memcpy(&x, bases + 8 * half, sizeof(x));
        constexpr uint64_t low7 = 0x7F7F7F7F7F7F7F7FULL;
        // High bit of each byte of v that is zero, without carries between bytes.
        auto const zero = [](uint64_t v) { return ~(((v & low7) + low7) | v | low7); };
        auto const upper = x & 0xDFDFDFDFDFDFDFDFULL;
        auto const valid = zero(upper ^ 0x4141414141414141ULL) |
                           zero(upper ^ 0x4343434343434343ULL) |
                           zero(upper ^ 0x4747474747474747ULL) | zero(upper ^ 0x5454545454545454ULL);
        auto codes = (x >> 1) & 0x0303030303030303ULL;
        codes = (codes | (codes >> 6)) & 0x000F000F000F000FULL;
        codes = (codes | (codes >> 12)) & 0x000000FF000000FFULL;
        codes = (codes | (codes >> 24)) & 0xFFFFULL;
        word |= static_cast<uint32_t>(codes) << (16 * half);
        // Byte j's flag lands on bit j.
        bad |= static_cast<uint32_t>(
                   ((~valid & 0x8080808080808080ULL) >> 7) * 0x0102040810204080ULL >> 56
               )
               << (8 * half);
    }
    return bad;
}
 
inline uint64_t pack_bases_pieces(char const* bases, size_t count, uint32_t* words, size_t width) {
    uint64_t bad = 0;
    for (size_t piece = 0; piece < width / 16; ++piece) {
        auto const left = count > 16 * piece ? count - 16 * piece : 0;
        bad |= uint64_t{pack_bases16_portable(bases + 16 * piece, left, words[piece])} << (16 * piece);
    }
    return bad;
}
 
inline uint64_t pack_bases_portable(char const* bases, size_t count, uint32_t* words) {
    return pack_bases16_portable(bases, count, words[0]);
}
 
#if (defined(__x86_64__) || defined(__i386__)) && (defined(__GNUC__) || defined(__clang__))
__attribute__((target("avx2"))) inline uint64_t
pack_bases_avx2(char const* bases, size_t count, uint32_t* words) {
    if (count < 32) return pack_bases_pieces(bases, count, words, 32);
    auto const x = _mm256_loadu_si256(reinterpret_cast<__m256i const*>(bases));
    auto const upper = _mm256_and_si256(x, _mm256_set1_epi8(static_cast<char>(0xDF)));
    auto const valid = _mm256_or_si256(
        _mm256_or_si256(
            _mm256_cmpeq_epi8(upper, _mm256_set1_epi8('A')),
            _mm256_cmpeq_epi8(upper, _mm256_set1_epi8('C'))
        ),
        _mm256_or_si256(
            _mm256_cmpeq_epi8(upper, _mm256_set1_epi8('G')),
            _mm256_cmpeq_epi8(upper, _mm256_set1_epi8('T'))
        )
    );
    // Codes weighted 1, 4, 16, 64 and summed: every four bases make one byte.
    auto const codes = _mm256_and_si256(_mm256_srli_epi16(x, 1), _mm256_set1_epi8(3));
    auto const pairs = _mm256_maddubs_epi16(codes, _mm256_set1_epi16(0x0401));
    auto const quads = _mm256_madd_epi16(pairs, _mm256_set1_epi32(0x00100001));
    // Each lane's four bytes, then both lanes into one quadword.
    auto const gathered = _mm256_shuffle_epi8(
        quads, _mm256_setr_epi8(0, 4, 8, 12, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
                                0, 4, 8, 12, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1)
    );
    auto const both = _mm256_permutevar8x32_epi32(gathered, _mm256_setr_epi32(0, 4, 0, 0, 0, 0, 0, 0));
    _mm_storel_epi64(reinterpret_cast<__m128i*>(words), _mm256_castsi256_si128(both));
    return ~static_cast<uint64_t>(static_cast<uint32_t>(_mm256_movemask_epi8(valid))) &
           0xFFFFFFFFULL;
}
 
__attribute__((target("avx2,avx512f,avx512bw"))) inline uint64_t
pack_bases_avx512(char const* bases, size_t count, uint32_t* words) {
    if (count < 64) {
        // A short tail packs as two AVX2 halves.
        auto const low = pack_bases_avx2(bases, count, words);
        auto const high = pack_bases_avx2(bases + 32, count > 32 ? count - 32 : 0, words + 2);
        return low | (high << 32);
    }
    auto const x = _mm512_loadu_si512(bases);
    auto const upper = _mm512_and_si512(x, _mm512_set1_epi8(static_cast<char>(0xDF)));
    auto const valid = _mm512_cmpeq_epi8_mask(upper, _mm512_set1_epi8('A')) |
                       _mm512_cmpeq_epi8_mask(upper, _mm512_set1_epi8('C')) |
                       _mm512_cmpeq_epi8_mask(upper, _mm512_set1_epi8('G')) |
                       _mm512_cmpeq_epi8_mask(upper, _mm512_set1_epi8('T'));
    // Codes weighted 1, 4, 16, 64 and summed: every four bases make one byte.
    auto const codes = _mm512_and_si512(_mm512_srli_epi16(x, 1), _mm512_set1_epi8(3));
    auto const pairs = _mm512_maddubs_epi16(codes, _mm512_set1_epi16(0x0401));
    auto const quads = _mm512_madd_epi16(pairs, _mm512_set1_epi32(0x00100001));
    _mm_storeu_si128(reinterpret_cast<__m128i*>(words), _mm512_cvtepi32_epi8(quads));
    return ~static_cast<uint64_t>(valid);
}
#endif
 
#if defined(__aarch64__) || defined(_M_ARM64)
inline uint64_t pack_bases_neon(char const* bases, size_t count, uint32_t* words) {
    if (count < 64) return pack_bases_pieces(bases, count, words, 64);
    static constexpr int8_t shift_table[16] = {0, 2, 4, 6, 0, 2, 4, 6, 0, 2, 4, 6, 0, 2, 4, 6};
    auto const shifts = vld1q_s8(shift_table);
    for (size_t piece = 0; piece < 4; ++piece) {
        auto const x = vld1q_u8(reinterpret_cast<uint8_t const*>(bases) + 16 * piece);
        auto const upper = vandq_u8(x, vdupq_n_u8(0xDF));
        auto const valid = vorrq_u8(
            vorrq_u8(vceqq_u8(upper, vdupq_n_u8('A')), vceqq_u8(upper, vdupq_n_u8('C'))),
            vorrq_u8(vceqq_u8(upper, vdupq_n_u8('G')), vceqq_u8(upper, vdupq_n_u8('T')))
        );
        // One nibble per byte: set when the byte is not a base.
        auto const invalid = vget_lane_u64(
            vreinterpret_u64_u8(vshrn_n_u16(vreinterpretq_u16_u8(vmvnq_u8(valid)), 4)), 0
        );
        // Shift each code into its slot, then add four neighbours: their bits are disjoint.
        auto const codes = vshlq_u8(vandq_u8(vshrq_n_u8(x, 1), vdupq_n_u8(3)), shifts);
        auto const pairs = vpaddq_u8(codes, codes);
        words[piece] = vgetq_lane_u32(vreinterpretq_u32_u8(vpaddq_u8(pairs, pairs)), 0);
        if (invalid != 0) {
            return ~uint64_t{0} << (16 * piece + static_cast<size_t>(std::countr_zero(invalid)) / 4);
        }
    }
    return 0;
}
#endif
 
struct bases_packer {
    uint64_t (*pack)(char const* bases, size_t count, uint32_t* words);
    size_t width;
};
 
[[nodiscard]] inline bases_packer resolve_bases_packer() {
    if (ablated("simd-parsing")) return {pack_bases_portable, 16};
#if (defined(__x86_64__) || defined(__i386__)) && (defined(__GNUC__) || defined(__clang__))
    if (__builtin_cpu_supports("avx512bw")) return {pack_bases_avx512, 64};
    if (__builtin_cpu_supports("avx2")) return {pack_bases_avx2, 32};
#endif
#if defined(__aarch64__) || defined(_M_ARM64)
    return {pack_bases_neon, 64};
#else
    return {pack_bases_portable, 16};
#endif
}
 
static const bases_packer packer = resolve_bases_packer();
 
inline void pack_fastx_sequence_file(fastx_sequence_file& result, uint32_t k) {
    if (result.extents.empty() || k < 16) return;
    auto* const out = const_cast<char*>(result.extents.front().begin);
    uint64_t written = 0;
    for (auto const& extent : result.extents) {
        auto const* next = extent.begin;
        auto const* const end = extent.end;
        while (next < end) {
            while (next < end && encode_base(*next) == 0xFF) ++next;
            auto const* const start = next;
            auto const first = written;
            while (next < end) {
                uint32_t words[4];
                auto const bad = packer.pack(next, static_cast<size_t>(end - next), words);
                auto const bases =
                    bad == 0 ? packer.width : static_cast<size_t>(std::countr_zero(bad));
                // A partial step ends the run; a short run's partial words are never written.
                if (bases == packer.width || static_cast<size_t>(next + bases - start) >= k) {
                    auto const count = (bases + 15) / 16;
                    std::memcpy(out + written * sizeof(uint32_t), words, count * sizeof(uint32_t));
                    written += count;
                }
                next += bases;
                if (bases != packer.width) break;
            }
            auto const bases = static_cast<uint64_t>(next - start);
            if (bases >= k) {
                result.packed.push_back({first, bases});
            } else {
                written = first;
            }
        }
    }
    result.words = reinterpret_cast<uint32_t const*>(out);
    result.extents.clear();
}
 
struct decompression_target {
    char* data = nullptr;
    size_t capacity = 0;
    std::shared_ptr<void> owner;
};
 
struct decompression_source {
    decompression_target (*acquire)(void* context, size_t bytes) = nullptr;
    void* context = nullptr;
    bool direct_reads = false;
 
    [[nodiscard]] decompression_target request(size_t bytes) const {
        return acquire != nullptr ? acquire(context, bytes) : decompression_target{};
    }
};
 
struct shared_gzip_decompressor {
    libdeflate_decompressor* handle{libdeflate_alloc_decompressor()};
    shared_gzip_decompressor() = default;
    shared_gzip_decompressor(shared_gzip_decompressor const&) = delete;
    shared_gzip_decompressor& operator=(shared_gzip_decompressor const&) = delete;
    ~shared_gzip_decompressor() {
        libdeflate_free_decompressor(handle);
    }
};
 
[[nodiscard]] inline libdeflate_result decompress_gzip_member(
    libdeflate_decompressor* decoder,
    std::string_view input,
    char* output,
    size_t capacity,
    size_t& consumed,
    size_t& produced
) {
    auto const byte = [&](size_t i) {
        return static_cast<uint32_t>(static_cast<unsigned char>(input[i]));
    };
    if (input.size() < 18 || byte(0) != 0x1f || byte(1) != 0x8b || byte(2) != 8 ||
        (byte(3) & 0xe0) != 0) {
        return LIBDEFLATE_BAD_DATA;
    }
    auto const flags = byte(3);
    auto const end = input.size() - 8;
    size_t at = 10;
    if ((flags & 4) != 0) {
        if (end - at < 2) return LIBDEFLATE_BAD_DATA;
        auto const length = byte(at) | (byte(at + 1) << 8);
        at += 2;
        if (length > end - at) return LIBDEFLATE_BAD_DATA;
        at += length;
    }
    for (auto const flag : {8U, 16U}) {
        if ((flags & flag) != 0) {
            auto const terminator = input.find('\0', at);
            if (terminator == std::string_view::npos || terminator >= end) {
                return LIBDEFLATE_BAD_DATA;
            }
            at = terminator + 1;
        }
    }
    if ((flags & 2) != 0) {
        if (end - at < 2) return LIBDEFLATE_BAD_DATA;
        at += 2;
    }
    size_t deflate_bytes = 0;
    auto const status = libdeflate_deflate_decompress_ex(
        decoder, input.data() + at, end - at, output, capacity, &deflate_bytes, &produced
    );
    if (status != LIBDEFLATE_SUCCESS) return status;
    auto const trailer = at + deflate_bytes;
    auto const size = byte(trailer + 4) | (byte(trailer + 5) << 8) | (byte(trailer + 6) << 16) |
                      (byte(trailer + 7) << 24);
    if (static_cast<uint32_t>(produced) != size) return LIBDEFLATE_BAD_DATA;
    consumed = trailer + 8;
    return LIBDEFLATE_SUCCESS;
}
 
[[nodiscard]] inline Result<void>
gunzip_members_into(std::string_view input, std::string& output, std::string const& path) {
    static thread_local shared_gzip_decompressor shared;
    if (shared.handle == nullptr) {
        return Err(Error::resource("cannot allocate gzip decompressor"));
    }
    // Genomes compress around 3x; guessing low costs a realloc per member, not correctness.
    if (output.empty()) output.resize(std::max<size_t>(input.size() * 4, 1U << 16));
    size_t consumed = 0;
    size_t produced = 0;
    while (consumed < input.size()) {
        size_t member_in = 0;
        size_t member_out = 0;
        auto const status = decompress_gzip_member(
            shared.handle,
            input.substr(consumed),
            output.data() + produced,
            output.size() - produced,
            member_in,
            member_out
        );
        if (status == LIBDEFLATE_INSUFFICIENT_SPACE) {
            output.resize(output.size() * 2 + 1);
            continue;
        }
        if (status != LIBDEFLATE_SUCCESS) {
            return Err(Error::invalid_argument("gzip FASTX error: invalid stream: " + path));
        }
        consumed += member_in;
        produced += member_out;
    }
    output.resize(produced);
    return Ok();
}
 
[[nodiscard]] inline Result<size_t> gunzip_members_direct(
    std::string_view input,
    decompression_target target,
    bool& insufficient,
    std::string const& path
) {
    insufficient = false;
    static thread_local shared_gzip_decompressor shared;
    if (shared.handle == nullptr) {
        return Err(Error::resource("cannot allocate gzip decompressor"));
    }
    size_t consumed = 0;
    size_t produced = 0;
    while (consumed < input.size()) {
        size_t member_in = 0;
        size_t member_out = 0;
        auto const status = decompress_gzip_member(
            shared.handle,
            input.substr(consumed),
            target.data + produced,
            target.capacity - produced,
            member_in,
            member_out
        );
        if (status == LIBDEFLATE_INSUFFICIENT_SPACE) {
            // A caller-supplied buffer cannot grow, so the loader falls back instead.
            insufficient = true;
            return Result<size_t>::ok(0);
        }
        if (status != LIBDEFLATE_SUCCESS) {
            return Err(Error::invalid_argument("gzip FASTX error: invalid stream: " + path));
        }
        consumed += member_in;
        produced += member_out;
    }
    return produced;
}
 
[[nodiscard]] inline uint64_t gzip_decompressed_size(std::string const& path) {
    auto file = std::fopen(path.c_str(), "rb");
    if (file == nullptr) return 0;
    uint64_t size = 0;
    if (std::fseek(file, -4, SEEK_END) == 0) {
        unsigned char trailer[4] = {};
        if (std::fread(trailer, 1, sizeof(trailer), file) == sizeof(trailer)) {
            size = static_cast<uint64_t>(trailer[0]) | (static_cast<uint64_t>(trailer[1]) << 8) |
                   (static_cast<uint64_t>(trailer[2]) << 16) |
                   (static_cast<uint64_t>(trailer[3]) << 24);
        }
    }
    std::fclose(file);
    return size;
}
 
[[nodiscard]] inline size_t gzip_size_hint(std::string_view input) {
    if (input.size() < 8) return 0;
    auto const byte = [&](size_t index) {
        return static_cast<uint32_t>(static_cast<unsigned char>(input[index]));
    };
    auto const tail = input.size() - 4;
    return byte(tail) | (byte(tail + 1) << 8) | (byte(tail + 2) << 16) | (byte(tail + 3) << 24);
}
 
#if defined(__x86_64__) || defined(__i386__)
    #include <immintrin.h>
    #include <tmmintrin.h>
#elif defined(__ARM_NEON)
    #include <arm_neon.h>
#endif
 
namespace compact_impl {
 
struct compact_shuffle_table {
    alignas(64) uint8_t control[256][16];
    uint8_t kept[256];
    compact_shuffle_table() noexcept {
        for (int mask = 0; mask < 256; ++mask) {
            int at = 0;
            for (int bit = 0; bit < 8; ++bit) {
                if ((mask & (1 << bit)) == 0) control[mask][at++] = static_cast<uint8_t>(bit);
            }
            for (; at < 16; ++at) control[mask][at] = 0x80U;
            kept[mask] = static_cast<uint8_t>(8 - __builtin_popcount(static_cast<unsigned>(mask)));
        }
    }
};
 
[[nodiscard]] inline compact_shuffle_table const& compact_shuffle() noexcept {
    static compact_shuffle_table const table;
    return table;
}
 
#if defined(__ARM_NEON)
[[nodiscard]] inline unsigned compact_mask(uint8x16_t flags) noexcept {
    constexpr std::uint8_t weights[16] = {1, 2, 4, 8, 16, 32, 64, 128, 1, 2, 4, 8, 16, 32, 64, 128};
    auto const weighted = vandq_u8(flags, vld1q_u8(weights));
    return static_cast<unsigned>(vaddv_u8(vget_low_u8(weighted))) |
           (static_cast<unsigned>(vaddv_u8(vget_high_u8(weighted))) << 8U);
}
#endif
 
}  // namespace compact_impl
 
#if defined(__ARM_NEON)
inline char* compact_block_neon(uint8x16_t block, uint8x16_t whitespace, char* out) {
    if (vmaxvq_u8(whitespace) == 0) {
        vst1q_u8(reinterpret_cast<uint8_t*>(out), block);
        return out + 16;
    }
    auto const& table = compact_impl::compact_shuffle();
    auto const flags = compact_impl::compact_mask(whitespace);
    auto const low = flags & 0xFFU;
    auto const high = (flags >> 8U) & 0xFFU;
    auto const packed_low = vqtbl1q_u8(block, vld1q_u8(table.control[low]));
    auto const packed_high = vqtbl1q_u8(vextq_u8(block, block, 8), vld1q_u8(table.control[high]));
    // Half-width stores cannot overwrite input beyond the loaded block when compacting in place.
    vst1_u8(reinterpret_cast<uint8_t*>(out), vget_low_u8(packed_low));
    out += table.kept[low];
    vst1_u8(reinterpret_cast<uint8_t*>(out), vget_low_u8(packed_high));
    return out + table.kept[high];
}
 
inline char* compact_sequence_whitespace_neon(char const* first, char const* last, char* out) {
    auto const newline = vdupq_n_u8('\n');
    auto const carriage = vdupq_n_u8('\r');
    auto const space = vdupq_n_u8(' ');
    auto const tab = vdupq_n_u8('\t');
    auto const* at = first;
    while (static_cast<size_t>(last - at) >= 16) {
        auto const block = vld1q_u8(reinterpret_cast<std::uint8_t const*>(at));
        auto const whitespace = vorrq_u8(
            vorrq_u8(vceqq_u8(block, newline), vceqq_u8(block, carriage)),
            vorrq_u8(vceqq_u8(block, space), vceqq_u8(block, tab))
        );
        out = compact_block_neon(block, whitespace, out);
        at += 16;
    }
    for (; at != last; ++at) {
        if (!fastx_is_sequence_whitespace(*at)) *out++ = *at;
    }
    return out;
}
#endif
 
#if defined(__x86_64__) || defined(__i386__)
__attribute__((target("avx2"))) inline char*
compact_block_avx2(__m256i block, unsigned flags, char* out) {
    auto const& table = compact_impl::compact_shuffle();
    auto const low_first = flags & 0xFFU;
    auto const low_second = (flags >> 8U) & 0xFFU;
    auto const high_first = (flags >> 16U) & 0xFFU;
    auto const high_second = (flags >> 24U) & 0xFFU;
    auto const control_first = _mm256_set_m128i(
        _mm_load_si128(reinterpret_cast<__m128i const*>(table.control[high_first])),
        _mm_load_si128(reinterpret_cast<__m128i const*>(table.control[low_first]))
    );
    auto const control_second = _mm256_set_m128i(
        _mm_load_si128(reinterpret_cast<__m128i const*>(table.control[high_second])),
        _mm_load_si128(reinterpret_cast<__m128i const*>(table.control[low_second]))
    );
    auto const packed_first = _mm256_shuffle_epi8(block, control_first);
    auto const packed_second = _mm256_shuffle_epi8(_mm256_bsrli_epi128(block, 8), control_second);
    // A lane holds the kept bytes of two halves, which are eight bytes apart in the vector and
    // adjacent in the output. The halves leave the vector in the order 0-7, 16-23, 8-15, 24-31,
    // and the output wants them in byte order, so each half is stored where its own count puts it.
    _mm_storel_epi64(reinterpret_cast<__m128i*>(out), _mm256_castsi256_si128(packed_first));
    out += table.kept[low_first];
    _mm_storel_epi64(reinterpret_cast<__m128i*>(out), _mm256_castsi256_si128(packed_second));
    out += table.kept[low_second];
    _mm_storel_epi64(reinterpret_cast<__m128i*>(out), _mm256_extracti128_si256(packed_first, 1));
    out += table.kept[high_first];
    _mm_storel_epi64(reinterpret_cast<__m128i*>(out), _mm256_extracti128_si256(packed_second, 1));
    out += table.kept[high_second];
    return out;
}
 
__attribute__((target("avx2"))) inline char*
compact_sequence_whitespace_avx2(char const* first, char const* last, char* out) {
    auto const newline = _mm256_set1_epi8('\n');
    auto const carriage = _mm256_set1_epi8('\r');
    auto const space = _mm256_set1_epi8(' ');
    auto const tab = _mm256_set1_epi8('\t');
    auto const* at = first;
    while (static_cast<size_t>(last - at) >= 32) {
        auto const block = _mm256_loadu_si256(reinterpret_cast<__m256i const*>(at));
        auto const whitespace = _mm256_or_si256(
            _mm256_or_si256(_mm256_cmpeq_epi8(block, newline), _mm256_cmpeq_epi8(block, carriage)),
            _mm256_or_si256(_mm256_cmpeq_epi8(block, space), _mm256_cmpeq_epi8(block, tab))
        );
        out =
            compact_block_avx2(block, static_cast<unsigned>(_mm256_movemask_epi8(whitespace)), out);
        at += 32;
    }
    for (; at != last; ++at) {
        if (!fastx_is_sequence_whitespace(*at)) *out++ = *at;
    }
    return out;
}
#endif
 
static bool const simd_compaction = !ablated("simd-parsing");
 
#if (defined(__x86_64__) || defined(__i386__)) && (defined(__GNUC__) || defined(__clang__))
static bool const compact_avx512 = simd_compaction && __builtin_cpu_supports("avx512bw") &&
                                   __builtin_cpu_supports("avx512vbmi2");
static bool const compact_avx2 = simd_compaction && __builtin_cpu_supports("avx2");
 
__attribute__((target("avx512f,avx512bw,avx512vbmi2,popcnt"))) inline char*
compact_sequence_whitespace_avx512(char const* first, char const* last, char* out) {
    auto const* at = first;
    while (static_cast<size_t>(last - at) >= 64) {
        auto const block = _mm512_loadu_si512(at);
        auto const keep =
            ~(_mm512_cmpeq_epi8_mask(block, _mm512_set1_epi8('\n')) |
              _mm512_cmpeq_epi8_mask(block, _mm512_set1_epi8('\r')) |
              _mm512_cmpeq_epi8_mask(block, _mm512_set1_epi8(' ')) |
              _mm512_cmpeq_epi8_mask(block, _mm512_set1_epi8('\t')));
        _mm512_storeu_si512(out, _mm512_maskz_compress_epi8(keep, block));
        out += std::popcount(keep);
        at += 64;
    }
    for (; at != last; ++at) {
        if (!fastx_is_sequence_whitespace(*at)) *out++ = *at;
    }
    return out;
}
#endif
 
inline char* compact_sequence_whitespace(char const* first, char const* last, char* out) {
    if (static_cast<size_t>(last - first) >= 64 && simd_compaction) {
#if defined(__ARM_NEON)
        return compact_sequence_whitespace_neon(first, last, out);
#elif defined(__x86_64__) || defined(__i386__)
        if (compact_avx512) return compact_sequence_whitespace_avx512(first, last, out);
        if (compact_avx2) {
            return compact_sequence_whitespace_avx2(first, last, out);
        }
#endif
    }
    constexpr uint64_t ones = 0x0101010101010101ULL;
    auto const has_zero_byte = [](uint64_t word) {
        return (word - ones) & ~word & 0x8080808080808080ULL;
    };
    auto const has_byte = [&](uint64_t word, uint64_t byte) {
        return has_zero_byte(word ^ (ones * byte));
    };
    auto const* word = first;
    while (static_cast<size_t>(last - word) >= sizeof(uint64_t)) {
        uint64_t value = 0;
        std::memcpy(&value, word, sizeof(value));
        auto const whitespace = has_byte(value, ' ') | has_byte(value, '\t') |
                                has_byte(value, '\n') | has_byte(value, '\r');
        if (whitespace == 0) {
            std::memcpy(out, word, sizeof(value));
            out += sizeof(value);
        } else {
            for (size_t index = 0; index < sizeof(value); ++index) {
                if (!fastx_is_sequence_whitespace(word[index])) *out++ = word[index];
            }
        }
        word += sizeof(value);
    }
    for (auto const* at = word; at != last; ++at) {
        if (!fastx_is_sequence_whitespace(*at)) *out++ = *at;
    }
    return out;
}
 
[[nodiscard]] inline char const* fasta_next_header(char const* first, char const* last) {
    for (auto const* at = first; at != last; ++at) {
        at = static_cast<char const*>(std::memchr(at, '>', static_cast<size_t>(last - at)));
        if (at == nullptr) return last;
        if (at == first || at[-1] == '\n' || at[-1] == '\r') return at;
    }
    return last;
}
 
#if defined(__x86_64__) || defined(__i386__)
__attribute__((target("avx512f,avx512bw,avx512vbmi2,popcnt"))) inline char*
fasta_sequence_compact_avx512(char const* first, char const* last, char* out, char const*& stop) {
    auto const newline = _mm512_set1_epi8('\n');
    auto const carriage = _mm512_set1_epi8('\r');
    auto const space = _mm512_set1_epi8(' ');
    auto const tab = _mm512_set1_epi8('\t');
    auto const marker = _mm512_set1_epi8('>');
    // Bit set when the byte before the block ends a line; @p first follows a line ending.
    uint64_t carry = 1;
    auto const* at = first;
    while (static_cast<size_t>(last - at) >= 64) {
        auto const block = _mm512_loadu_si512(at);
        auto const ends =
            _mm512_cmpeq_epi8_mask(block, newline) | _mm512_cmpeq_epi8_mask(block, carriage);
        auto const headers = _mm512_cmpeq_epi8_mask(block, marker) & ((ends << 1U) | carry);
        if (headers != 0) {
            stop = at + std::countr_zero(headers);
            return compact_sequence_whitespace(at, stop, out);
        }
        carry = ends >> 63U;
        auto const keep =
            ~(ends | _mm512_cmpeq_epi8_mask(block, space) | _mm512_cmpeq_epi8_mask(block, tab));
        _mm512_storeu_si512(out, _mm512_maskz_compress_epi8(keep, block));
        out += std::popcount(keep);
        at += 64;
    }
    bool line_start = carry != 0;
    for (; at != last; ++at) {
        char const ch = *at;
        if (ch == '>' && line_start) break;
        line_start = ch == '\n' || ch == '\r';
        if (!fastx_is_sequence_whitespace(ch)) *out++ = ch;
    }
    stop = at;
    return out;
}
 
__attribute__((target("avx2"))) inline char*
fasta_sequence_compact_avx2(char const* first, char const* last, char* out, char const*& stop) {
    auto const newline = _mm256_set1_epi8('\n');
    auto const carriage = _mm256_set1_epi8('\r');
    auto const space = _mm256_set1_epi8(' ');
    auto const tab = _mm256_set1_epi8('\t');
    auto const marker = _mm256_set1_epi8('>');
    // Bit set when the byte before the block ends a line; @p first follows a line ending.
    unsigned carry = 1;
    auto const* at = first;
    while (static_cast<size_t>(last - at) >= 32) {
        auto const block = _mm256_loadu_si256(reinterpret_cast<__m256i const*>(at));
        auto const line_end =
            _mm256_or_si256(_mm256_cmpeq_epi8(block, newline), _mm256_cmpeq_epi8(block, carriage));
        auto const ends = static_cast<unsigned>(_mm256_movemask_epi8(line_end));
        auto const markers =
            static_cast<unsigned>(_mm256_movemask_epi8(_mm256_cmpeq_epi8(block, marker)));
        auto const headers = markers & ((ends << 1U) | carry);
        if (headers != 0) {
            stop = at + std::countr_zero(headers);
            return compact_sequence_whitespace(at, stop, out);
        }
        carry = ends >> 31U;
        auto const whitespace = _mm256_or_si256(
            line_end,
            _mm256_or_si256(_mm256_cmpeq_epi8(block, space), _mm256_cmpeq_epi8(block, tab))
        );
        out =
            compact_block_avx2(block, static_cast<unsigned>(_mm256_movemask_epi8(whitespace)), out);
        at += 32;
    }
    bool line_start = carry != 0;
    for (; at != last; ++at) {
        char const ch = *at;
        if (ch == '>' && line_start) break;
        line_start = ch == '\n' || ch == '\r';
        if (!fastx_is_sequence_whitespace(ch)) *out++ = ch;
    }
    stop = at;
    return out;
}
#endif
 
#if defined(__ARM_NEON)
inline char*
fasta_sequence_compact_neon(char const* first, char const* last, char* out, char const*& stop) {
    unsigned carry = 1;
    auto const* at = first;
    while (static_cast<size_t>(last - at) >= 16) {
        auto const block = vld1q_u8(reinterpret_cast<uint8_t const*>(at));
        auto const line_end =
            vorrq_u8(vceqq_u8(block, vdupq_n_u8('\n')), vceqq_u8(block, vdupq_n_u8('\r')));
        auto const ends = compact_impl::compact_mask(line_end);
        auto const markers = compact_impl::compact_mask(vceqq_u8(block, vdupq_n_u8('>')));
        auto const headers = markers & ((ends << 1U) | carry);
        if (headers != 0) {
            stop = at + std::countr_zero(headers);
            return compact_sequence_whitespace(at, stop, out);
        }
        carry = ends >> 15U;
        auto const whitespace = vorrq_u8(
            line_end, vorrq_u8(vceqq_u8(block, vdupq_n_u8(' ')), vceqq_u8(block, vdupq_n_u8('\t')))
        );
        out = compact_block_neon(block, whitespace, out);
        at += 16;
    }
    bool line_start = carry != 0;
    for (; at != last; ++at) {
        char const ch = *at;
        if (ch == '>' && line_start) break;
        line_start = ch == '\n' || ch == '\r';
        if (!fastx_is_sequence_whitespace(ch)) *out++ = ch;
    }
    stop = at;
    return out;
}
#endif
 
inline char*
fasta_sequence_compact(char const* first, char const* last, char* out, char const*& stop) {
#if (defined(__x86_64__) || defined(__i386__)) && (defined(__GNUC__) || defined(__clang__))
    if (compact_avx512) return fasta_sequence_compact_avx512(first, last, out, stop);
    if (compact_avx2) return fasta_sequence_compact_avx2(first, last, out, stop);
#endif
#if defined(__ARM_NEON)
    if (simd_compaction) return fasta_sequence_compact_neon(first, last, out, stop);
#endif
    stop = fasta_next_header(first, last);
    return compact_sequence_whitespace(first, stop, out);
}
 
[[nodiscard]] inline std::vector<fastx_sequence_extent>
fastx_fasta_compact_extents(std::string_view data, char* out) {
    std::vector<fastx_sequence_extent> extents;
    size_t header = data.find('>');
    while (header != std::string_view::npos && header != 0 && data[header - 1] != '\n' &&
           data[header - 1] != '\r') {
        header = data.find('>', header + 1);
    }
    if (header == std::string_view::npos) return extents;
    auto const* const last = data.data() + data.size();
    while (true) {
        auto const end = fastx_line_end(data, header);
        if (end + 1 >= data.size()) break;
        auto const* const sequence = data.data() + end + 1;
        char const* stop = nullptr;
        auto* const begin = out;
        out = fasta_sequence_compact(sequence, last, out, stop);
        if (stop != sequence) extents.push_back({begin, out});
        if (stop == last) break;
        header = static_cast<size_t>(stop - data.data());
    }
    return extents;
}
 
inline void compact_fastx_sequence_extents(fastx_sequence_file& result) {
    if (result.extents.empty()) return;
    // Extents point into whichever buffer holds the bytes: a page-locked target, the heap copy of
    // a compressed file, or the FASTQ sequence buffer. Those are writable; a memory-mapped plain
    // file is not, so its bases are compacted into `decompressed` instead.
    struct writable_buffer {
        char* data;
        size_t size;
    };
    auto const buffers = std::array{
        writable_buffer{result.decompressed_target, result.decompressed_size},
        writable_buffer{result.decompressed.data(), result.decompressed.size()},
        writable_buffer{result.sequence.data(), result.sequence.size()},
    };
    auto const writable = [&](char const* first, size_t size) {
        return std::any_of(buffers.begin(), buffers.end(), [&](writable_buffer const& buffer) {
            return buffer.data != nullptr && buffer.size >= size && first >= buffer.data &&
                   first + size <= buffer.data + buffer.size;
        });
    };
    size_t raw_total = 0;
    for (auto const& extent : result.extents)
        raw_total += static_cast<size_t>(extent.end - extent.begin);
    // Reserved up front so appending cannot move the buffer the extents are about to point at.
    std::string copied;
    copied.reserve(raw_total);
    std::vector<std::pair<size_t, std::pair<size_t, size_t>>> offsets;
    for (size_t index = 0; index < result.extents.size(); ++index) {
        auto const& extent = result.extents[index];
        auto const size = static_cast<size_t>(extent.end - extent.begin);
        if (writable(extent.begin, size)) {
            // Compaction only ever moves bytes towards the buffer's start.
            auto* const write = compact_sequence_whitespace(
                extent.begin, extent.end, const_cast<char*>(extent.begin)
            );
            result.extents[index] = {extent.begin, write};
        } else {
            auto const begin = copied.size();
            // Sized before the copy so the write needs no bounds check per byte.
            copied.resize(begin + size);
            auto* const write =
                compact_sequence_whitespace(extent.begin, extent.end, copied.data() + begin);
            copied.resize(static_cast<size_t>(write - copied.data()));
            offsets.emplace_back(index, std::pair{begin, copied.size()});
        }
    }
    if (!copied.empty()) {
        result.decompressed = std::move(copied);
        for (auto const& [index, span] : offsets) {
            result.extents[index] = {
                result.decompressed.data() + span.first, result.decompressed.data() + span.second
            };
        }
    }
}
 
inline Result<void>
parse_fastx_sequence_file(fastx_sequence_file& result, std::string const& path) {
    auto const data = std::exchange(result.input, {});
    auto& sequence = result.sequence;
    auto& extents = result.extents;
    auto const first = data.find_first_not_of("\r\n");
    if (first != std::string_view::npos && data[first] == '@') {
        fastx_fastq_reader reader{data, path};
        std::vector<std::pair<size_t, size_t>> records;
        while (true) {
            auto const begin = sequence.size();
            auto next = reader.append_next_record(sequence);
            if (!next) return Err(Error::invalid_argument(next.error().message()));
            if (!*next) break;
            if (sequence.size() > begin) records.emplace_back(begin, sequence.size());
        }
        if (sequence.empty()) return Ok();
        for (auto const& [begin, end] : records) {
            extents.push_back({sequence.data() + begin, sequence.data() + end});
        }
        compact_fastx_sequence_extents(result);
        return Ok();
    }
    // Mapped plain bytes are read-only; decompressed bytes can be compacted in place.
    char* out = result.decompressed_target;
    if (out == nullptr) {
        if (data.data() != result.decompressed.data()) result.decompressed.resize(data.size());
        out = result.decompressed.data();
    }
    extents = fastx_fasta_compact_extents(data, out);
    if (extents.empty() && !data.empty()) {
        return Err(Error::invalid_argument("FASTX parse error near: " + path));
    }
    return Ok();
}
 
inline Result<void> pack_fasta_chunks(
    fastx_sequence_file& result,
    std::string_view data,
    char* out,
    uint32_t k,
    std::string const& path
) {
    size_t header = data.find('>');
    while (header != std::string_view::npos && header != 0 && data[header - 1] != '\n' &&
           data[header - 1] != '\r') {
        header = data.find('>', header + 1);
    }
    if (header == std::string_view::npos) {
        if (!data.empty()) return Err(Error::invalid_argument("FASTX parse error near: " + path));
        return Ok();
    }
    std::string scratch, carry;
    fastx_sequence_file chunk;
    uint64_t written = 0;
    bool has_sequence = false;
    auto const* const last = data.data() + data.size();
    while (true) {
        auto const end = fastx_line_end(data, header);
        if (end + 1 >= data.size()) break;
        auto const* at = data.data() + end + 1;
        carry.clear();
        while (at != last) {
            auto const offset = static_cast<size_t>(at - data.data());
            auto const boundary = offset + std::min<size_t>(16384, data.size() - offset);
            auto const limit = boundary == data.size()
                                   ? boundary
                                   : std::min(data.size(), fastx_line_end(data, boundary) + 1);
            auto const* const chunk_end = data.data() + limit;
            scratch.resize(carry.size() + static_cast<size_t>(chunk_end - at));
            std::memcpy(scratch.data(), carry.data(), carry.size());
            char const* stop = nullptr;
            auto* compact_end =
                fasta_sequence_compact(at, chunk_end, scratch.data() + carry.size(), stop);
            has_sequence |= stop != at;
            auto const bases = static_cast<size_t>(compact_end - scratch.data());
            carry.assign(compact_end - std::min<size_t>(k - 1, bases), compact_end);
            chunk.extents.push_back({scratch.data(), compact_end});
            pack_fastx_sequence_file(chunk, k);
            uint64_t words = 0;
            for (auto record : chunk.packed) {
                words = record.word + (record.bases + 15) / 16;
                record.word += written;
                result.packed.push_back(record);
            }
            std::memcpy(out + written * sizeof(uint32_t), scratch.data(), words * sizeof(uint32_t));
            written += words;
            chunk.packed.clear();
            at = stop;
            if (stop != chunk_end || stop == last) break;
        }
        if (at == last) break;
        header = static_cast<size_t>(at - data.data());
    }
    if (!has_sequence && !data.empty()) {
        return Err(Error::invalid_argument("FASTX parse error near: " + path));
    }
    result.words = reinterpret_cast<uint32_t const*>(out);
    result.input = {};
    return Ok();
}
 
inline Result<std::unique_ptr<fastx_sequence_file>> finish_fastx_sequence_file(
    std::unique_ptr<fastx_sequence_file> result,
    std::string_view data,
    std::string const& path,
    bool defer_fasta,
    uint32_t pack_k
) {
    result->input = data;
    if (defer_fasta) {
        auto const first = data.find_first_not_of("\r\n");
        if (first != std::string_view::npos && data[first] == '>') {
            auto const end = fastx_line_end(data, first);
            // Empty first records and malformed inputs stay with the CPU parser.
            if (end + 1 < data.size() && data[end + 1] != '>') return result;
        }
    }
    if (pack_k >= 16) {
        auto const first = data.find_first_not_of("\r\n");
        if (first == std::string_view::npos || data[first] != '@') {
            char* out = result->decompressed_target;
            if (out == nullptr) {
                if (data.data() != result->decompressed.data())
                    result->decompressed.resize(data.size());
                out = result->decompressed.data();
            }
            CUDDL_TRY(pack_fasta_chunks(*result, data, out, pack_k, path));
            return result;
        }
    }
    CUDDL_TRY(parse_fastx_sequence_file(*result, path));
    return result;
}
 
[[nodiscard]] inline Result<std::optional<std::string_view>> inflate_and_pack(
    fastx_sequence_file& result,
    std::string_view input,
    std::string const& path,
    decompression_source source,
    uint32_t k
) {
    thread_local std::unique_ptr<char[]> scratch;
    thread_local size_t capacity = 0;
    // A single-member file's ISIZE is its inflated length. ISIZE sizes nothing for a BGZF stream,
    // which ends in an empty EOF block that reads back as zero, so a hint below the compressed
    // size falls back to an over-estimate; a short buffer doubles and inflates again.
    auto const hint = gzip_size_hint(input);
    size_t wanted = hint >= input.size() ? hint : input.size() * 8;
    size_t produced = 0;
    while (true) {
        if (capacity < wanted) {
            capacity = std::bit_ceil(wanted);
            scratch = std::make_unique_for_overwrite<char[]>(capacity);
        }
        bool insufficient = false;
        produced = CUDDL_TRY(
            gunzip_members_direct(input, {scratch.get(), capacity, {}}, insufficient, path)
        );
        if (!insufficient) break;
        wanted = capacity * 2;
    }
    std::string_view const data{scratch.get(), produced};
    auto const first = data.find_first_not_of("\r\n");
    bool const fastq = first != std::string_view::npos && data[first] == '@';
    if (!fastq) CUDDL_TRY(pack_fasta_chunks(result, data, scratch.get(), k, path));
    uint64_t bytes = produced;
    if (!fastq) {
        uint64_t words = 0;
        for (auto const& record : result.packed) {
            words = std::max(words, record.word + (record.bases + 15) / 16);
        }
        bytes = words * sizeof(uint32_t);
    }
    auto target = source.request(bytes);
    char* out = target.data;
    if (out == nullptr) {
        result.decompressed.resize(bytes);
        out = result.decompressed.data();
    } else {
        result.decompressed_target = out;
        result.decompressed_size = bytes;
        result.storage_owner = std::move(target.owner);
    }
    std::memcpy(out, scratch.get(), bytes);
    if (fastq) return std::optional{std::string_view{out, bytes}};
    result.words = reinterpret_cast<uint32_t const*>(out);
    return std::optional<std::string_view>{};
}
 
inline Result<std::unique_ptr<fastx_sequence_file>> decode_fastx_sequence_file(
    std::unique_ptr<fastx_mapped_file> file,
    std::string const& path,
    decompression_source source,
    bool defer_fasta,
    uint32_t pack_k = 0
) {
    auto result = std::make_unique<fastx_sequence_file>();
    result->file = std::move(file);
    auto data = result->file->data();
    auto& decompressed = result->decompressed;
    if (data.size() >= 2 && static_cast<unsigned char>(data[0]) == 0x1f &&
        static_cast<unsigned char>(data[1]) == 0x8b) {
        if (pack_k >= 16 && !defer_fasta && source.acquire != nullptr) {
            auto const fastq = CUDDL_TRY(inflate_and_pack(*result, data, path, source, pack_k));
            if (!fastq) return result;
            return finish_fastx_sequence_file(std::move(result), *fastq, path, false, pack_k);
        }
        // A single-member file's ISIZE is its inflated length. Asking for that rather than an
        // over-estimate keeps loader buffers small, which matters most when they are
        // page-locked, since the driver pins one allocation at a time. ISIZE sizes
        // nothing for a BGZF stream, which ends in an empty EOF block that reads back as zero, so
        // a hint below the compressed size falls back to the over-estimate the heap path assumes.
        // A concatenated stream advertises only its last member; a short target grows the
        // request, which keeps a multi-member file off the heap path.
        auto const hint = gzip_size_hint(data);
        auto request = hint >= data.size() ? hint : data.size() * 8;
        auto target = source.request(request);
        bool insufficient = false;
        bool inflated = false;
        size_t produced = 0;
        while (target.data != nullptr && target.capacity > 0) {
            auto const written = gunzip_members_direct(data, target, insufficient, path);
            if (!written) return Err(written.error());
            produced = *written;
            if (!insufficient) {
                inflated = true;
                break;
            }
            request = std::max(request * 2, data.size());
            target = source.request(request);
        }
        if (inflated) {
            result->decompressed_target = target.data;
            result->decompressed_size = produced;
            result->storage_owner = std::move(target.owner);
            data = std::string_view{target.data, produced};
        }
        if (result->decompressed_target == nullptr) {
            auto inflated = gunzip_members_into(data, decompressed, path);
            if (!inflated) return Err(inflated.error());
            data = decompressed;
        }
    } else if (source.acquire != nullptr && !defer_fasta) {
        // Plain bytes go to the caller's buffer, as inflated bytes do: the parser compacts them
        // in place and the stager reads them without another copy.
        auto target = source.request(data.size());
        if (target.data != nullptr) {
            std::memcpy(target.data, data.data(), data.size());
            result->decompressed_target = target.data;
            result->decompressed_size = data.size();
            result->storage_owner = std::move(target.owner);
            data = std::string_view{target.data, data.size()};
        }
    }
    return finish_fastx_sequence_file(std::move(result), data, path, defer_fasta, pack_k);
}
 
inline Result<std::unique_ptr<fastx_sequence_file>> load_fastx_sequence_file(
    std::string const& path,
    decompression_source source = {},
    bool defer_fasta = false,
    bool read = false,
    uint32_t pack_k = 0
) {
#if defined(__linux__)
    if (source.acquire != nullptr) {
        int const fd = ::open(path.c_str(), O_RDONLY);
        struct stat status{};
        if (fd == -1 || ::fstat(fd, &status) != 0 || status.st_size < 0) {
            if (fd != -1) ::close(fd);
            return Err(Error::invalid_argument("cannot open FASTX file: " + path));
        }
        auto const size = static_cast<size_t>(status.st_size);
        unsigned char magic[2] = {};
        bool const gzip = ::pread(fd, magic, sizeof(magic), 0) == sizeof(magic) &&
                          magic[0] == 0x1f && magic[1] == 0x8b;
        // Plain bytes go straight into the caller's buffer, as inflated bytes do: the parser
        // compacts them in place and the stager reads them without another copy.
        auto target = gzip ? decompression_target{} : source.request(size);
        if (target.data == nullptr) {
            auto file = fastx_mapped_file::adopt(fd, size, path, read, source.direct_reads);
            if (!file) return Err(Error::invalid_argument("cannot open FASTX file: " + path));
            return decode_fastx_sequence_file(std::move(*file), path, source, defer_fasta, pack_k);
        }
        if (!read_whole_file(fd, target.data, size, source.direct_reads, target.capacity)) {
            return Err(Error::invalid_argument("cannot open FASTX file: " + path));
        }
        auto result = std::make_unique<fastx_sequence_file>();
        result->decompressed_target = target.data;
        result->decompressed_size = size;
        result->storage_owner = std::move(target.owner);
        return finish_fastx_sequence_file(
            std::move(result), {target.data, size}, path, defer_fasta, pack_k
        );
    }
#endif
    auto file = fastx_mapped_file::load(path, read, source.direct_reads);
    if (!file) return Err(Error::invalid_argument("cannot open FASTX file: " + path));
    return decode_fastx_sequence_file(std::move(*file), path, source, defer_fasta, pack_k);
}
 
inline constexpr size_t parallel_read_depth = 64;
 
// Fixed worker pool loading FASTX files ahead of the GPU loop. A pthread spawn and join
// per genome costs tens of microseconds and dominates at high file counts. Readers pull
// ids in order with a bounded number of outstanding files. Consumers may collect by ID
// or by completion, retaining each input's ID. Loader exceptions are
// rethrown by take, matching std::async propagation into the build error handlers.
//
// Readers read whole files ahead of the loaders, which only decode. A parallel file system
// gives its bandwidth only to many reads in flight, and a loader that read its own files left
// storage idle while it inflated: on MOGON's Lustre, 16 loaders sketched 4,000 genomes no node
// had cached in 10.3-11.6 s reading their own files, and in 4.0-4.5 s with readers.
class fastx_load_pool {
   public:
    fastx_load_pool(
        std::span<std::filesystem::path const> paths,
        size_t workers,
        decompression_source source = {},
        size_t window = 0,
        bool defer_fasta = false,
        bool completion_order = false,
        uint32_t pack_k = 0
    )
        : paths_(paths),
          workers_(std::max(size_t{1}, workers)),
          window_(std::max({size_t{1}, window, workers_})),
          source_(source),
          defer_fasta_(defer_fasta),
          pack_k_(pack_k),
          results_(paths.size()),
          errors_(paths.size()),
          ready_(completion_order ? window_ : 0),
          retired_(defer_fasta ? window_ : 0) {
        threads_.reserve(workers_);
        if (defer_fasta_) retirer_ = std::thread([this] { retire_files(); });
        // A reader spends most of its time waiting on storage, so there are as many as the reads
        // a parallel file system needs in flight, however few the loaders.
        auto const readers =
            ablated("read-depth") ? workers_
                                  : std::min(window_, std::max(workers_, parallel_read_depth));
        readers_.reserve(readers);
        for (size_t i = 0; i < readers; ++i) {
            readers_.emplace_back([this] { read_files(); });
        }
        for (size_t i = 0; i < workers_; ++i) threads_.emplace_back([this] { work(); });
    }
    ~fastx_load_pool() {
        {
            std::lock_guard lock(mutex_);
            stop_ = true;
        }
        assign_.notify_all();
        mapped_ready_.notify_all();
        for (auto& thread : readers_) thread.join();
        for (auto& thread : threads_) thread.join();
        // Loaders may be waiting for retirement space, so stop the reaper only after joining them.
        if (retirer_.joinable()) {
            {
                std::lock_guard lock(retired_mutex_);
                retire_stop_ = true;
            }
            retired_ready_.notify_one();
            retirer_.join();
        }
    }
    fastx_load_pool(fastx_load_pool const&) = delete;
    fastx_load_pool& operator=(fastx_load_pool const&) = delete;
 
    [[nodiscard]] Result<std::unique_ptr<fastx_sequence_file>> take(size_t id) {
        std::unique_lock lock(mutex_);
        filled_.wait(lock, [&] { return results_[id].has_value() || errors_[id] != nullptr; });
        ++taken_;
        // One freed window slot needs one loader, not all of them: waking the rest is a syscall
        // and a context switch each, once per genome, at 100,000 genomes.
        assign_.notify_one();
        lock.unlock();
        if (errors_[id] != nullptr) std::rethrow_exception(errors_[id]);
        return std::move(*results_[id]);
    }
 
    [[nodiscard]] Result<std::pair<size_t, std::unique_ptr<fastx_sequence_file>>> take_ready() {
        if (ready_.empty())
            return Err(Error::invalid_argument("completion-order loading is disabled"));
        std::unique_lock lock(mutex_);
        filled_.wait(lock, [&] { return ready_count_ != 0; });
        auto const id = ready_[ready_head_];
        ready_head_ = (ready_head_ + 1) % ready_.size();
        --ready_count_;
        ++taken_;
        assign_.notify_one();
        lock.unlock();
        if (errors_[id] != nullptr) std::rethrow_exception(errors_[id]);
        return std::pair{id, CUDDL_TRY(std::move(*results_[id]))};
    }
 
    [[nodiscard]] Result<std::optional<std::pair<size_t, std::unique_ptr<fastx_sequence_file>>>>
    take_ready_until_limit() {
        {
            std::unique_lock lock(mutex_);
            filled_.wait(lock, [&] { return ready_count_ != 0 || taken_ >= limit_; });
            if (ready_count_ == 0) return std::nullopt;
        }
        return CUDDL_TRY(take_ready());
    }
 
    [[nodiscard]] size_t ready() {
        std::lock_guard lock(mutex_);
        return ready_count_;
    }
 
    size_t cap(size_t limit) {
        {
            std::lock_guard lock(mutex_);
            limit_ = std::clamp(limit, claimed_end_, limit_);
            limit = limit_;
        }
        assign_.notify_all();
        filled_.notify_all();
        mapped_ready_.notify_all();
        return limit;
    }
 
   private:
    struct mapped_input {
        size_t id = 0;
        std::optional<Result<std::unique_ptr<fastx_mapped_file>>> file;
        std::exception_ptr error;
    };
 
    void read_files() {
        while (true) {
            size_t id;
            {
                std::unique_lock lock(mutex_);
                assign_.wait(lock, [&] {
                    return stop_ || mapped_next_ >= limit_ || mapped_next_ - taken_ < window_;
                });
                if (stop_ || mapped_next_ >= limit_) return;
                id = mapped_next_++;
            }
            mapped_input input;
            input.id = id;
            try {
                input.file =
                    fastx_mapped_file::load(paths_[id].string(), true, source_.direct_reads);
            } catch (...) {
                input.error = std::current_exception();
            }
            {
                std::lock_guard lock(mutex_);
                mapped_.push_back(std::move(input));
            }
            mapped_ready_.notify_one();
        }
    }
 
    void retire(std::unique_ptr<fastx_mapped_file> file) {
        std::unique_lock lock(retired_mutex_);
        retired_space_.wait(lock, [&] { return retired_count_ < retired_.size(); });
        retired_[(retired_head_ + retired_count_) % retired_.size()] = std::move(file);
        ++retired_count_;
        lock.unlock();
        retired_ready_.notify_one();
    }
 
    void retire_files() {
        while (true) {
            std::unique_ptr<fastx_mapped_file> file;
            {
                std::unique_lock lock(retired_mutex_);
                retired_ready_.wait(lock, [&] { return retire_stop_ || retired_count_ != 0; });
                if (retired_count_ == 0) return;
                file = std::move(retired_[retired_head_]);
                retired_head_ = (retired_head_ + 1) % retired_.size();
                --retired_count_;
            }
            retired_space_.notify_one();
        }
    }
 
    void work() {
        while (true) {
            size_t id;
            mapped_input input;
            {
                std::unique_lock lock(mutex_);
                // Readers own the window; loaders take files as they are read, so one slow read
                // holds up only itself. Files read past a cap belong to another consumer.
                mapped_ready_.wait(lock, [&] {
                    return stop_ || claimed_ >= limit_ || !mapped_.empty();
                });
                if (stop_ || claimed_ >= limit_) return;
                input = std::move(mapped_.front());
                mapped_.pop_front();
                id = input.id;
                if (id >= limit_) continue;
                ++claimed_;
                claimed_end_ = std::max(claimed_end_, id + 1);
                // Another file may have been read with no loader woken for it.
                if (!mapped_.empty()) mapped_ready_.notify_one();
            }
            try {
                auto const path = paths_[id].string();
                auto loaded = [&]() -> Result<std::unique_ptr<fastx_sequence_file>> {
                    if (input.error) std::rethrow_exception(input.error);
                    if (!*input.file)
                        return Err(Error::invalid_argument("cannot open FASTX file: " + path));
                    return decode_fastx_sequence_file(
                        std::move(**input.file), path, source_, defer_fasta_, pack_k_
                    );
                }();
                if (pack_k_ != 0 && loaded) pack_fastx_sequence_file(**loaded, pack_k_);
                // Inflated or copied bytes no longer reference the file's bytes.
                if (defer_fasta_ && loaded && (*loaded)->file &&
                    ((*loaded)->decompressed_target != nullptr ||
                     !(*loaded)->decompressed.empty())) {
                    retire(std::move((*loaded)->file));
                }
                std::lock_guard lock(mutex_);
                results_[id] = std::move(loaded);
                if (!ready_.empty()) {
                    ready_[(ready_head_ + ready_count_) % ready_.size()] = id;
                    ++ready_count_;
                }
            } catch (...) {
                std::lock_guard lock(mutex_);
                errors_[id] = std::current_exception();
                if (!ready_.empty()) {
                    ready_[(ready_head_ + ready_count_) % ready_.size()] = id;
                    ++ready_count_;
                }
            }
            filled_.notify_all();
        }
    }
    std::span<std::filesystem::path const> paths_;
    size_t workers_;
    size_t window_;
    decompression_source source_;
    bool defer_fasta_;
    uint32_t pack_k_;
    std::vector<std::optional<Result<std::unique_ptr<fastx_sequence_file>>>> results_;
    std::vector<std::exception_ptr> errors_;
    std::vector<size_t> ready_;
    std::vector<std::unique_ptr<fastx_mapped_file>> retired_;
    std::deque<mapped_input> mapped_;
    std::vector<std::thread> threads_;
    std::thread retirer_;
    std::vector<std::thread> readers_;
    std::mutex mutex_;
    std::mutex retired_mutex_;
    std::condition_variable assign_, filled_;
    std::condition_variable retired_ready_, retired_space_;
    std::condition_variable mapped_ready_;
    // Files loaders took, and one past the highest id among them.
    size_t claimed_ = 0, claimed_end_ = 0, taken_ = 0, limit_ = paths_.size();
    size_t ready_head_ = 0, ready_count_ = 0;
    size_t retired_head_ = 0, retired_count_ = 0;
    size_t mapped_next_ = 0;
    bool stop_ = false;
    bool retire_stop_ = false;
};
 
}  // namespace cuddl::detail
 
namespace cuddl {
 
struct reference_build_statistics {
    size_t direct_bytes = 0;
    size_t staged_bytes = 0;
    size_t direct_chunks = 0;
    size_t staged_chunks = 0;
    size_t pinned_buffers = 0;
    unsigned workers = 0;   // loaders the build ran, after its own defaulting
    bool in_place = false;  // the kernels read the loader's own buffers, so nothing was copied
    size_t staging_bytes = 0;
    size_t batches = 0;
    size_t transfers = 0;  // runs copied in one piece, plus records split across pieces
    // Files the gpu backend took without a host loader: inflated by nvCOMP or the
    // decompression engine, or plain FASTA read straight into the GPU's buffers.
    size_t device_files = 0;
    bool measure_resident = false;
    double resident_compute_ms = 0;
};
}  // namespace cuddl