cuddl::a48 Namespace Reference

Namespaces
namespace	detail

Classes
struct	database
	Fully decoded A48 reference database. More...
struct	database_metadata
	Per-database construction metadata decoded from the A48 file-level header. More...
struct	oracle_result
	Independent host oracle for one query over a decoded reference database. More...
struct	record
	One decoded DDL score row plus its stable record identity and construction metadata. More...
struct	record_metadata
	Metadata attached to one decoded A48 record. More...
struct	row_summary
	Exact bucket-classification counts for one decoded row pair. More...

Functions
std::optional< uint64_t >	decode_a48_token (std::string_view token)
	Decodes one A48 (base-64, ASCII-offset-48) token to its integer value.
std::string	encode_a48_token (uint64_t value)
	Encodes an integer as an A48 (base-64, ASCII-offset-48) token.
uint16_t	promote_score (uint16_t loaded, int32_t offset, uint8_t exponent_bits)
	Reconstructs one absolute 16-bit bucket maxima from an A48-decoded value.
std::vector< std::string_view >	split_fields (std::string_view text, char separator, bool keep_trailing_empty=false)
	Splits text on separator.
bool	is_record_header (std::string_view key)
	True when key is a record-level header (part of a record block) rather than a file-level header that precedes every record.
Result< record_metadata >	parse_record_headers (std::string_view block)
	Parses a record header block (one "#key\tvalue" per line) into record_metadata.
Result< std::vector< uint16_t > >	decode_a48_row (std::string_view row, record_metadata const &metadata, uint8_t exponent_bits)
	Decodes one A48 data row into absolute 16-bit bucket scores.
Result< database >	decode_a48_tsv (std::string_view input)
	Decodes an uncompressed A48-encoded DDL TSV into score rows and metadata.
Result< database >	decode_a48_tsv_parallel (std::string_view input, unsigned requested_threads=std::thread::hardware_concurrency())
	Multithreaded variant of decode_a48_tsv.
oracle_result	exhaustive_oracle (std::vector< uint16_t > const &rows, std::vector< uint16_t > const &query, size_t bucket_count, uint32_t minimum_matches)
	Exact match counts and pair summaries for one flat query row.
oracle_result	exhaustive_oracle (std::vector< std::vector< uint16_t > > const &rows, std::vector< uint16_t > const &query, uint32_t minimum_matches)
	Convenience overload over decoded records instead of one flat score array.

Function Documentation

decode_a48_row()

Result< std::vector< uint16_t > > cuddl::a48::decode_a48_row ( std::string_view  row, record_metadata const &  metadata, uint8_t  exponent_bits  )

inline

Decodes one A48 data row into absolute 16-bit bucket scores.

Definition at line 264 of file a48.hpp.

                                                                                           {
    auto const fields = split_fields(row, '\t', true);
    std::vector<uint16_t> scores;
    scores.reserve(fields.size());
    for (auto const field : fields) {
        auto const loaded = decode_a48_token(field);
        if (!loaded) {
            return Err(
                Error::invalid_argument(
                    "A48 data row contains a character outside the A48 alphabet"
                )
            );
        }
        scores.push_back(
            promote_score(static_cast<uint16_t>(*loaded), metadata.offset, exponent_bits)
        );
    }
    return Result<std::vector<uint16_t>>::ok(std::move(scores));
}

Here is the call graph for this function:

Here is the caller graph for this function:

decode_a48_token()

std::optional< uint64_t > cuddl::a48::decode_a48_token ( std::string_view  token )

inline

Decodes one A48 (base-64, ASCII-offset-48) token to its integer value.

Returns std::nullopt when a character is outside the A48 alphabet (ASCII 48..111), i.e. the token is malformed.

Definition at line 84 of file a48.hpp.

                                                                                  {
    if (token.empty()) {
        return uint64_t{0};
    }
    uint64_t value = 0;
    for (char const character : token) {
        auto const code = static_cast<uint8_t>(character);
        if (code < 48U || code > 111U) {
            return std::nullopt;
        }
        value = (value << 6U) | static_cast<uint64_t>(code - 48U);
    }
    return value;
}

Here is the caller graph for this function:

decode_a48_tsv()

Result< database > cuddl::a48::decode_a48_tsv ( std::string_view  input )

inline

Decodes an uncompressed A48-encoded DDL TSV into score rows and metadata.

The decoder is intentionally strict about the structural invariants that would otherwise yield silent garbage: A48 alphabet violations, an omitted #len disagreeing with the data row, a bucket count that changes between records, an absent #exponent on a relative-encoded file, and records without any bucket values are all reported as malformed input. It does not, however, depend on the large official RefSeq asset, so a small deterministic fixture can exercise the code paths in unit tests.

Definition at line 292 of file a48.hpp.

                                                                           {
    database result;
    uint8_t exponent_bits = 6U;
    bool has_exponent = false;
    size_t expected_bucket_count = 0;
    bool bucket_count_set = false;
 
    std::string header_block;
    size_t ordinal = 0;
    size_t begin = 0;
    while (begin <= input.size()) {
        auto const end = input.find('\n', begin);
        auto line =
            input.substr(begin, end == std::string_view::npos ? input.size() - begin : end - begin);
        if (!line.empty() && line.back() == '\r') {
            line.remove_suffix(1);
        }
 
        if (line.empty()) {
            if (end == std::string_view::npos) {
                break;
            }
            begin = end + 1U;
            continue;
        }
 
        if (line.front() == '#') {
            auto const fields = split_fields(line, '\t');
            if (fields.empty()) {
                break;
            }
            auto key = fields.front();
            key.remove_prefix(1);  // drop '#'
            auto const value = fields.size() > 1U ? fields[1] : std::string_view{};
 
            if (key == "k") {
                auto const parsed = detail::parse_field<uint32_t>(
                    [](std::string const& s) { return std::stoul(s); }, value, "k"
                );
                if (!parsed) {
                    return Err(parsed.error());
                }
                result.metadata.kmer_length = *parsed;
                result.metadata.has_kmer_length = true;
            } else if (key == "seed") {
                auto const parsed = detail::parse_field<uint64_t>(
                    [](std::string const& s) { return std::stoull(s); }, value, "seed"
                );
                if (!parsed) {
                    return Err(parsed.error());
                }
                result.metadata.seed = *parsed;
                result.metadata.has_seed = true;
            } else if (key == "exponent") {
                auto const parsed = detail::parse_field<uint32_t>(
                    [](std::string const& s) { return std::stoul(s); }, value, "exponent"
                );
                if (!parsed) {
                    return Err(parsed.error());
                }
                if (*parsed < 1U || *parsed > 15U) {
                    return Err(Error::invalid_argument("A48 exponent must be in [1, 15]"));
                }
                exponent_bits = *parsed;
                result.metadata.exponent_bits = *parsed;
                result.metadata.has_exponent = true;
                has_exponent = true;
            } else if (key == "blacklist") {
                result.metadata.blacklist = std::string(value);
            }
 
            if (is_record_header(key)) {
                if (!header_block.empty()) {
                    header_block.push_back('\n');
                }
                header_block.append(line);
            }
 
            if (end == std::string_view::npos) {
                break;
            }
            begin = end + 1U;
            continue;
        }
 
        // A data row terminates the currently accumulated record header block.
        if (header_block.empty()) {
            return Err(Error::invalid_argument("A48 data row appears before any record header"));
        }
        auto meta = CUDDL_TRY(parse_record_headers(header_block));
        header_block.clear();
 
        if (meta.offset >= 0 && !has_exponent) {
            return Err(
                Error::invalid_argument("relative-encoded A48 record requires a #exponent header")
            );
        }
        auto scores = CUDDL_TRY(decode_a48_row(line, meta, exponent_bits));
        if (scores.empty()) {
            return Err(Error::invalid_argument("A48 record has no bucket values"));
        }
        if (!bucket_count_set) {
            expected_bucket_count = scores.size();
            bucket_count_set = true;
        } else if (scores.size() != expected_bucket_count) {
            return Err(Error::invalid_argument("A48 data row width changed between records"));
        }
        if (meta.len != 0U && meta.len != scores.size()) {
            return Err(Error::invalid_argument("#len disagrees with the A48 data row width"));
        }
 
        record decoded;
        decoded.ordinal = static_cast<uint32_t>(ordinal);
        decoded.metadata = meta;
        decoded.scores = std::move(scores);
        result.records.push_back(std::move(decoded));
        ++ordinal;
 
        if (end == std::string_view::npos) {
            break;
        }
        begin = end + 1U;
    }
 
    if (!header_block.empty()) {
        return Err(Error::invalid_argument("A48 file ends inside a record header block"));
    }
    if (result.records.empty()) {
        return Err(Error::invalid_argument("A48 file contains no records"));
    }
 
    result.metadata.record_count = result.records.size();
    return Result<database>::ok(std::move(result));
}

Here is the call graph for this function:

decode_a48_tsv_parallel()

Result< database > cuddl::a48::decode_a48_tsv_parallel ( std::string_view  input, unsigned  requested_threads = std::thread::hardware_concurrency()  )

inline

Multithreaded variant of decode_a48_tsv.

One pass discovers file-level metadata, record header blocks, and data-row spans. The first record is decoded synchronously to establish the bucket width; the remaining independent records are decoded by a fixed pool and written directly to their ordinal slots, so order and strict validation are identical to the single-threaded decoder.

requested_threads defaults to std::thread::hardware_concurrency().

One undecoded record: its header block, data row, and file ordinal.

Definition at line 435 of file a48.hpp.

  {
    struct record_job {
        std::string header;
        std::string_view data;
        uint32_t ordinal{};
    };
 
    database result;
    std::vector<record_job> jobs;
    uint8_t exponent_bits = 6U;
    bool has_exponent = false;
    size_t expected_bucket_count = 0;
    bool bucket_count_set = false;
 
    std::string header_block;
    size_t ordinal = 0;
    size_t begin = 0;
    while (begin <= input.size()) {
        auto const end = input.find('\n', begin);
        auto line =
            input.substr(begin, end == std::string_view::npos ? input.size() - begin : end - begin);
        if (!line.empty() && line.back() == '\r') {
            line.remove_suffix(1);
        }
        if (line.empty()) {
            if (end == std::string_view::npos) {
                break;
            }
            begin = end + 1U;
            continue;
        }
 
        if (line.front() == '#') {
            auto const fields = split_fields(line, '\t');
            if (fields.empty()) {
                break;
            }
            auto key = fields.front();
            key.remove_prefix(1);
            auto const value = fields.size() > 1U ? fields[1] : std::string_view{};
            if (key == "k") {
                auto const parsed = detail::parse_field<uint32_t>(
                    [](std::string const& s) { return std::stoul(s); }, value, "k"
                );
                if (!parsed) {
                    return Err(parsed.error());
                }
                result.metadata.kmer_length = *parsed;
                result.metadata.has_kmer_length = true;
            } else if (key == "seed") {
                auto const parsed = detail::parse_field<uint64_t>(
                    [](std::string const& s) { return std::stoull(s); }, value, "seed"
                );
                if (!parsed) {
                    return Err(parsed.error());
                }
                result.metadata.seed = *parsed;
                result.metadata.has_seed = true;
            } else if (key == "exponent") {
                auto const parsed = detail::parse_field<uint32_t>(
                    [](std::string const& s) { return std::stoul(s); }, value, "exponent"
                );
                if (!parsed) {
                    return Err(parsed.error());
                }
                if (*parsed < 1U || *parsed > 15U) {
                    return Err(Error::invalid_argument("A48 exponent must be in [1, 15]"));
                }
                exponent_bits = *parsed;
                result.metadata.exponent_bits = *parsed;
                result.metadata.has_exponent = true;
                has_exponent = true;
            } else if (key == "blacklist") {
                result.metadata.blacklist = std::string(value);
            }
 
            if (is_record_header(key)) {
                if (!header_block.empty()) {
                    header_block.push_back('\n');
                }
                header_block.append(line);
            }
            if (end == std::string_view::npos) {
                break;
            }
            begin = end + 1U;
            continue;
        }
 
        if (header_block.empty()) {
            return Err(Error::invalid_argument("A48 data row appears before any record header"));
        }
        auto const record_ordinal = static_cast<uint32_t>(ordinal);
        if (record_ordinal == 0U) {
            auto meta = CUDDL_TRY(parse_record_headers(header_block));
            header_block.clear();
            if (meta.offset >= 0 && !has_exponent) {
                return Err(
                    Error::invalid_argument(
                        "relative-encoded A48 record requires a #exponent header"
                    )
                );
            }
            auto scores = CUDDL_TRY(decode_a48_row(line, meta, exponent_bits));
            if (scores.empty()) {
                return Err(Error::invalid_argument("A48 record has no bucket values"));
            }
            expected_bucket_count = scores.size();
            bucket_count_set = true;
            if (meta.len != 0U && meta.len != scores.size()) {
                return Err(Error::invalid_argument("#len disagrees with the A48 data row width"));
            }
            record decoded;
            decoded.ordinal = record_ordinal;
            decoded.metadata = std::move(meta);
            decoded.scores = std::move(scores);
            result.records.push_back(std::move(decoded));
        } else {
            jobs.push_back(
                record_job{std::move(header_block), line, static_cast<uint32_t>(ordinal)}
            );
            header_block.clear();
        }
        ++ordinal;
 
        if (end == std::string_view::npos) {
            break;
        }
        begin = end + 1U;
    }
 
    if (!header_block.empty()) {
        return Err(Error::invalid_argument("A48 file ends inside a record header block"));
    }
    if (result.records.empty()) {
        return Err(Error::invalid_argument("A48 file contains no records"));
    }
 
    result.records.resize(result.records.size() + jobs.size());
    auto const requested = requested_threads;
    auto const worker_count = static_cast<size_t>(std::max(1U, requested));
    std::vector<std::thread> workers;
    std::vector<std::exception_ptr> errors(worker_count);
    workers.reserve(worker_count);
 
    for (size_t worker = 0; worker < worker_count; ++worker) {
        auto const chunk_begin = jobs.size() * worker / worker_count;
        auto const chunk_end = jobs.size() * (worker + 1U) / worker_count;
        if (chunk_begin == chunk_end) {
            continue;
        }
        workers.emplace_back([&, chunk_begin, chunk_end, worker]() {
            try {
                for (size_t i = chunk_begin; i < chunk_end; ++i) {
                    auto const& job = jobs[i];
                    auto meta = parse_record_headers(job.header);
                    if (!meta) {
                        throw std::runtime_error(meta.error().message());
                    }
                    if (meta->offset >= 0 && !has_exponent) {
                        throw std::runtime_error(
                            "relative-encoded A48 record requires a #exponent header"
                        );
                    }
                    auto scores = decode_a48_row(job.data, *meta, exponent_bits);
                    if (!scores) {
                        throw std::runtime_error(scores.error().message());
                    }
                    if (scores->empty()) {
                        throw std::runtime_error("A48 record has no bucket values");
                    }
                    if (!bucket_count_set || scores->size() != expected_bucket_count) {
                        throw std::runtime_error("A48 data row width changed between records");
                    }
                    if (meta->len != 0U && meta->len != scores->size()) {
                        throw std::runtime_error("#len disagrees with the A48 data row width");
                    }
                    record decoded;
                    decoded.ordinal = job.ordinal;
                    decoded.metadata = std::move(*meta);
                    decoded.scores = std::move(*scores);
                    result.records[job.ordinal] = std::move(decoded);
                }
            } catch (...) {
                errors[worker] = std::current_exception();
            }
        });
    }
    for (auto& worker : workers) {
        worker.join();
    }
    for (auto const& error : errors) {
        if (error) {
            try {
                std::rethrow_exception(error);
            } catch (std::exception const& exception) {
                return Err(Error::invalid_argument(exception.what()));
            }
        }
    }
 
    result.metadata.record_count = result.records.size();
    return Result<database>::ok(std::move(result));
}

Here is the call graph for this function:

encode_a48_token()

std::string cuddl::a48::encode_a48_token ( uint64_t  value )

inline

Encodes an integer as an A48 (base-64, ASCII-offset-48) token.

Definition at line 100 of file a48.hpp.

                                                                {
    if (value == 0U) {
        return std::string(1, '0');
    }
    std::string reversed;
    while (value != 0U) {
        reversed.push_back(static_cast<char>(48U + (value & 0x3fU)));
        value >>= 6U;
    }
    return {reversed.rbegin(), reversed.rend()};
}

exhaustive_oracle() [1/2]

oracle_result cuddl::a48::exhaustive_oracle ( std::vector< std::vector< uint16_t > > const &  rows, std::vector< uint16_t > const &  query, uint32_t  minimum_matches  )

inline

Convenience overload over decoded records instead of one flat score array.

Definition at line 718 of file a48.hpp.

  {
    std::vector<uint16_t> flat;
    auto const bucket_count = rows.empty() ? size_t{0} : rows.front().size();
    flat.reserve(rows.size() * bucket_count);
    for (auto const& row : rows) {
        flat.insert(flat.end(), row.begin(), row.end());
    }
    return exhaustive_oracle(flat, query, bucket_count, minimum_matches);
}

Here is the call graph for this function:

exhaustive_oracle() [2/2]

oracle_result cuddl::a48::exhaustive_oracle ( std::vector< uint16_t > const &  rows, std::vector< uint16_t > const &  query, size_t  bucket_count, uint32_t  minimum_matches  )

inline

Exact match counts and pair summaries for one flat query row.

Definition at line 680 of file a48.hpp.

  {
    auto const reference_count = rows.size() / bucket_count;
    oracle_result result;
    result.match_counts.assign(reference_count, 0U);
    result.summaries.assign(reference_count, row_summary{});
    for (size_t reference = 0; reference < reference_count; ++reference) {
        auto summary = row_summary{};
        auto const* reference_row = rows.data() + reference * bucket_count;
        auto const* query_row = query.data();
        for (size_t bucket = 0; bucket < bucket_count; ++bucket) {
            auto const query_score = query_row[bucket];
            auto const reference_score = reference_row[bucket];
            if (query_score == reference_score) {
                if (query_score == 0U) {
                    ++summary.both_empty;
                } else {
                    ++summary.equal;
                    ++result.match_counts[reference];
                }
            } else if (query_score < reference_score) {
                ++summary.lower;
            } else {
                ++summary.higher;
            }
        }
        if (result.match_counts[reference] >= minimum_matches) {
            result.summaries[reference] = summary;
        }
    }
    return result;
}

Here is the caller graph for this function:

is_record_header()

bool cuddl::a48::is_record_header ( std::string_view  key )

inline

True when key is a record-level header (part of a record block) rather than a file-level header that precedes every record.

Definition at line 170 of file a48.hpp.

                                                               {
    return key == "id" || key == "tid" || key == "name" || key == "file" || key == "bases" ||
           key == "contigs" || key == "gc" || key == "origin" || key == "lineage" || key == "len" ||
           key == "offset";
}

Here is the caller graph for this function:

parse_record_headers()

Result< record_metadata > cuddl::a48::parse_record_headers ( std::string_view  block )

inline

Parses a record header block (one "#key\tvalue" per line) into record_metadata.

Definition at line 177 of file a48.hpp.

                                                                                        {
    record_metadata meta;
    for (auto const line : split_fields(block, '\n')) {
        if (line.empty() || line.front() != '#') {
            continue;
        }
        auto const fields = split_fields(line, '\t');
        if (fields.empty()) {
            continue;
        }
        auto key = fields.front();
        key.remove_prefix(1);  // drop '#'
        if (fields.size() < 2U) {
            continue;
        }
        auto const value = fields[1];
        if (key == "id") {
            auto const parsed = detail::parse_field<int64_t>(
                [](std::string const& s) { return std::stoll(s); }, value, "id"
            );
            if (!parsed) {
                return Err(parsed.error());
            }
            meta.id = *parsed;
        } else if (key == "tid") {
            auto const parsed = detail::parse_field<int64_t>(
                [](std::string const& s) { return std::stoll(s); }, value, "tid"
            );
            if (!parsed) {
                return Err(parsed.error());
            }
            meta.tax_id = *parsed;
        } else if (key == "name") {
            meta.name = std::string(value);
        } else if (key == "file") {
            meta.filename = std::string(value);
        } else if (key == "bases") {
            auto const parsed = detail::parse_field<int64_t>(
                [](std::string const& s) { return std::stoll(s); }, value, "bases"
            );
            if (!parsed) {
                return Err(parsed.error());
            }
            meta.bases = *parsed;
        } else if (key == "contigs") {
            auto const parsed = detail::parse_field<int32_t>(
                [](std::string const& s) { return std::stol(s); }, value, "contigs"
            );
            if (!parsed) {
                return Err(parsed.error());
            }
            meta.contigs = *parsed;
        } else if (key == "gc") {
            auto const parsed = detail::parse_field<float>(
                [](std::string const& s) { return std::stof(s); }, value, "gc"
            );
            if (!parsed) {
                return Err(parsed.error());
            }
            meta.gc = *parsed;
        } else if (key == "origin") {
            meta.origin = std::string(value);
        } else if (key == "lineage") {
            meta.lineage = std::string(value);
        } else if (key == "len") {
            auto const parsed = detail::parse_field<uint32_t>(
                [](std::string const& s) { return std::stoul(s); }, value, "len"
            );
            if (!parsed) {
                return Err(parsed.error());
            }
            meta.len = *parsed;
        } else if (key == "offset") {
            auto const parsed = detail::parse_field<int32_t>(
                [](std::string const& s) { return std::stol(s); }, value, "offset"
            );
            if (!parsed) {
                return Err(parsed.error());
            }
            meta.offset = *parsed;
        }
    }
    return Result<record_metadata>::ok(std::move(meta));
}

Here is the call graph for this function:

Here is the caller graph for this function:

promote_score()

uint16_t cuddl::a48::promote_score ( uint16_t  loaded, int32_t  offset, uint8_t  exponent_bits  )

inline

Reconstructs one absolute 16-bit bucket maxima from an A48-decoded value.

Matches BBTools' DynamicDemiLog.fromArray: when offset is >= 0 the stored value is relative to the record's global NLZ and must be promoted back into the full 16-bit space; when offset is -1 the stored value is already an absolute score. exponent_bits selects the split between the exponent (NLZ tier) and mantissa fields.

Definition at line 119 of file a48.hpp.

                                                                      {
    if (loaded == 0U || offset < 0) {
        return loaded;
    }
    auto const mantissa_bits = 16U - exponent_bits;
    auto const mask = (1U << mantissa_bits) - 1U;
    auto const rel_nlz = static_cast<uint32_t>(loaded >> mantissa_bits);
    auto const abs_nlz = rel_nlz + static_cast<uint32_t>(offset);
    return static_cast<uint16_t>((abs_nlz << mantissa_bits) | (loaded & mask));
}

Here is the caller graph for this function:

split_fields()

std::vector< std::string_view > cuddl::a48::split_fields ( std::string_view  text, char  separator, bool  keep_trailing_empty = false  )

inline

Splits text on separator.

With keep_trailing_empty the split is exactly BBTools' LineParser1 term layout: a line ending in a separator produces one final empty field, which parseLongA48 decodes as zero. Header parsing drops that field instead (a trailing empty value is not a distinct header).

Definition at line 136 of file a48.hpp.

                                                                                    {
    std::vector<std::string_view> fields;
    size_t begin = 0;
    while (true) {
        auto const end = text.find(separator, begin);
        if (end == std::string_view::npos) {
            fields.push_back(text.substr(begin));
            break;
        }
        fields.push_back(text.substr(begin, end - begin));
        begin = end + 1U;
    }
    if (!keep_trailing_empty && !fields.empty() && fields.back().empty()) {
        fields.pop_back();
    }
    return fields;
}

Here is the caller graph for this function: