fastx_sequence_scan.hpp

Go to the documentation of this file.

#pragma once
 
#include <cstdint>
#include <string>
#include <string_view>
#include <utility>
#include <vector>
 
#include <cusbf/Alphabet.cuh>
#include <cusbf/detail/fastx_buffer_reader.hpp>
 
namespace cusbf::detail {
 
struct fastx_sequence_extent {
    const char* begin;
    const char* end;
};
 
[[nodiscard]] inline constexpr bool fastx_is_sequence_whitespace(char ch) noexcept {
    return ch == '\n' || ch == '\r' || ch == ' ' || ch == '\t';
}
 
[[nodiscard]] inline std::vector<fastx_sequence_extent> fastx_fasta_extents(std::string_view data) {
    std::vector<fastx_sequence_extent> extents;
    size_t sequence = std::string_view::npos;
    size_t header = data.find('>');
    while (header != std::string_view::npos) {
        if (header != 0 && data[header - 1] != '\n' && data[header - 1] != '\r') {
            header = data.find('>', header + 1);
            continue;
        }
        if (sequence != std::string_view::npos && header > sequence) {
            extents.push_back({data.data() + sequence, data.data() + header});
        }
        auto const end = fastx_line_end(data, header);
        sequence = end < data.size() ? end + 1 : data.size();
        header = data.find('>', sequence);
    }
    if (sequence != std::string_view::npos && sequence < data.size()) {
        extents.push_back({data.data() + sequence, data.data() + data.size()});
    }
    return extents;
}
 
struct fastx_sequence_span {
    std::string prefix;
    std::vector<std::string_view> segments;
    uint64_t total_bytes = 0;
};
 
namespace {
 
inline bool scan_step_back(
    std::vector<fastx_sequence_extent> const& extents, size_t& extent, const char*& position
) {
    if (position > extents[extent].begin) {
        --position;
        return true;
    }
    if (extent > 0) {
        --extent;
        position = extents[extent].end - 1;
        return true;
    }
    return false;
}
 
inline bool scan_step_forward(
    std::vector<fastx_sequence_extent> const& extents, size_t& extent, const char*& position
) {
    if (position + 1 < extents[extent].end) {
        ++position;
        return true;
    }
    if (extent + 1 < extents.size()) {
        ++extent;
        position = extents[extent].begin;
        return true;
    }
    return false;
}
 
}  // namespace
 
[[nodiscard]] inline std::vector<fastx_sequence_span> fastx_split_sequence_spans(
    std::vector<fastx_sequence_extent> const& extents, uint32_t seed_bases, uint32_t count
) {
    std::vector<fastx_sequence_span> spans;
    if (extents.empty() || count == 0U) {
        return spans;
    }
    size_t total = 0;
    for (auto const& extent : extents) {
        total += static_cast<size_t>(extent.end - extent.begin);
    }
    if (total == 0) {
        return spans;
    }
    if (static_cast<size_t>(count) > total) {
        count = static_cast<uint32_t>(total);
    }
    auto const per = total / count + (total % count != 0U);
    spans.reserve(count);
 
    for (uint32_t t = 0; t < count; ++t) {
        auto const span_begin = static_cast<size_t>(t) * per;
        if (span_begin >= total) {
            break;
        }
        auto const span_end = span_begin + per < total ? span_begin + per : total;
 
        fastx_sequence_span span;
        span.total_bytes = span_end - span_begin;
 
        // Locate the extent and offset of the span start.
        size_t extent = 0;
        size_t offset = span_begin;
        while (extent < extents.size() &&
               offset >= static_cast<size_t>(extents[extent].end - extents[extent].begin)) {
            offset -= static_cast<size_t>(extents[extent].end - extents[extent].begin);
            ++extent;
        }
 
        // Seed: walk back up to seed_bases valid bases (or until an invalid base / EOF).
        size_t seed_extent = extent;
        const char* seed_start = extents[extent].begin + offset;
        {
            size_t e = extent;
            const char* q = seed_start;
            uint32_t collected = 0;
            while (collected < seed_bases && scan_step_back(extents, e, q)) {
                auto const ch = *q;
                if (fastx_is_sequence_whitespace(ch)) {
                    continue;
                }
                if (DnaAlphabet::encode(&ch) == DnaAlphabet::invalidSymbol) {
                    seed_extent = e;
                    seed_start = q;
                    scan_step_forward(extents, seed_extent, seed_start);
                    break;
                }
                ++collected;
                seed_extent = e;
                seed_start = q;
            }
        }
 
        // Materialise the prefix from the seed start to the span start.
        {
            size_t e = seed_extent;
            const char* q = seed_start;
            const char* target = extents[extent].begin + offset;
            while (e < extent || q < target) {
                if (q >= extents[e].end) {
                    ++e;
                    q = extents[e].begin;
                    continue;
                }
                span.prefix.push_back(*q);
                ++q;
            }
        }
 
        // Collect the span's segments.
        {
            size_t remaining = span_end - span_begin;
            const char* q = extents[extent].begin + offset;
            size_t e = extent;
            while (remaining > 0 && e < extents.size()) {
                auto const limit = static_cast<size_t>(extents[e].end - q);
                auto const step = limit < remaining ? limit : remaining;
                span.segments.emplace_back(q, step);
                remaining -= step;
                if (remaining > 0) {
                    ++e;
                    if (e < extents.size()) {
                        q = extents[e].begin;
                    }
                }
            }
        }
 
        spans.push_back(std::move(span));
    }
    return spans;
}
 
}  // namespace cusbf::detail