diff --git a/cpp/CMakeLists.txt b/cpp/CMakeLists.txt index b7d21e7..943358d 100644 --- a/cpp/CMakeLists.txt +++ b/cpp/CMakeLists.txt @@ -41,14 +41,16 @@ endif() Python_add_library(_core MODULE WITH_SOABI src/bindings/python_bindings.cpp - src/seqwin/graph.cpp - src/seqwin/helpers.cpp + src/seqwin/build.cpp + src/seqwin/build_internals.cpp + src/seqwin/filter.cpp + src/utils/fasta_reader.cpp + src/utils/logging.cpp vendor/btllib/minimizer.cpp - src/seqwin/fasta_reader.cpp vendor/btllib/status.cpp ) -target_include_directories(_core PRIVATE include vendor "${PYBIND11_INCLUDE}") +target_include_directories(_core PRIVATE include src vendor "${PYBIND11_INCLUDE}") target_link_libraries(_core PRIVATE Python::Module ZLIB::ZLIBSTATIC) # GCC 8's libstdc++ keeps std::filesystem symbols in a separate library. diff --git a/cpp/include/seqwin/build.hpp b/cpp/include/seqwin/build.hpp new file mode 100644 index 0000000..f3218cd --- /dev/null +++ b/cpp/include/seqwin/build.hpp @@ -0,0 +1,34 @@ +#pragma once + +#include +#include +#include + +#include "seqwin/graph.hpp" + +namespace seqwin { + +/** + * @brief Build a minimizer graph from assembly FASTA files. + * + * `assembly_paths` and `is_targets` are parallel lists. + * + * @param assembly_paths Paths to input assemblies in FASTA format (plain or gzipped). + * @param kmerlen K-mer length for minimizer sketch. + * @param windowsize Window size for minimizer sketch. + * @param is_targets Whether each assembly is a target assembly. + * @param n_cpu Number of worker threads to use. + * @param low_memory Recompute minimizers in a second pass to reduce peak memory. + * @return Minimizer graph. + * @throws `std::runtime_error` If input sizes are inconsistent or counts exceed supported ranges. + */ +Graph build( + const std::vector& assembly_paths, + std::size_t kmerlen, + std::size_t windowsize, + const std::vector& is_targets, + std::size_t n_cpu = 1, + bool low_memory = false +); + +} // namespace seqwin diff --git a/cpp/include/seqwin/filter.hpp b/cpp/include/seqwin/filter.hpp new file mode 100644 index 0000000..b2a4405 --- /dev/null +++ b/cpp/include/seqwin/filter.hpp @@ -0,0 +1,18 @@ +#pragma once + +#include +#include +#include + +#include "seqwin/graph.hpp" + +namespace seqwin { + +Graph filter_kmers( + const Kmer* kmers, + const Node* nodes, + std::size_t n_nodes, + std::vector used_hashes +); + +} // namespace seqwin diff --git a/cpp/include/seqwin/graph.hpp b/cpp/include/seqwin/graph.hpp index 422b204..c0ec30f 100644 --- a/cpp/include/seqwin/graph.hpp +++ b/cpp/include/seqwin/graph.hpp @@ -68,27 +68,4 @@ struct Graph { std::vector> ids_by_assembly; }; -/** - * @brief Build a minimizer graph from assembly FASTA files. - * - * `assembly_paths` and `is_targets` are parallel lists. - * - * @param assembly_paths Paths to input assemblies in FASTA format (plain or gzipped). - * @param kmerlen K-mer length for minimizer sketch. - * @param windowsize Window size for minimizer sketch. - * @param is_targets Whether each assembly is a target assembly. - * @param n_cpu Number of worker threads to use. - * @param low_memory Recompute minimizers in a second pass to reduce peak memory. - * @return Minimizer graph. - * @throws `std::runtime_error` If input sizes are inconsistent or counts exceed supported ranges. - */ -Graph build( - const std::vector& assembly_paths, - std::size_t kmerlen, - std::size_t windowsize, - const std::vector& is_targets, - std::size_t n_cpu = 1, - bool low_memory = false -); - } // namespace seqwin diff --git a/cpp/src/bindings/python_bindings.cpp b/cpp/src/bindings/python_bindings.cpp index 6d45028..db443e6 100644 --- a/cpp/src/bindings/python_bindings.cpp +++ b/cpp/src/bindings/python_bindings.cpp @@ -9,8 +9,9 @@ #include #include +#include "seqwin/build.hpp" +#include "seqwin/filter.hpp" #include "seqwin/graph.hpp" -#include "seqwin/helpers.hpp" namespace py = pybind11; diff --git a/cpp/src/seqwin/graph.cpp b/cpp/src/seqwin/build.cpp similarity index 82% rename from cpp/src/seqwin/graph.cpp rename to cpp/src/seqwin/build.cpp index e0f2b40..622d7c5 100644 --- a/cpp/src/seqwin/graph.cpp +++ b/cpp/src/seqwin/build.cpp @@ -1,25 +1,61 @@ -#include "seqwin/graph.hpp" +#include "seqwin/build.hpp" #include #include -#include +#include #include #include #include +#include #include #include #include #include "btllib/minimizer.hpp" -#include "seqwin/fasta_reader.hpp" -#include "seqwin/helpers.hpp" -#include "seqwin/thread_pool.hpp" + +#include "seqwin/build_internals.hpp" +#include "utils/fasta_reader.hpp" +#include "utils/logging.hpp" +#include "utils/thread_pool.hpp" namespace seqwin { +namespace internal { namespace { constexpr std::size_t map_reserve_divisor = 100; +constexpr std::size_t plain_fasta_seq_len_per_byte = 1; +constexpr std::size_t gz_fasta_seq_len_per_byte = 4; + +bool ends_with(std::string_view text, std::string_view suffix) +{ + return text.size() >= suffix.size() && + text.substr(text.size() - suffix.size()) == suffix; +} + +std::size_t estimate_minimizer_count( + const std::vector& assembly_paths, + std::size_t start_assembly, + std::size_t end_assembly, + std::size_t windowsize +) { + // Reserve-only heuristic, not correctness-critical + std::size_t est_total_seq_len = 0; + for (std::size_t assembly_i = start_assembly; assembly_i < end_assembly; ++assembly_i) { + const auto& assembly_path = assembly_paths[assembly_i]; + const std::size_t seq_len_per_byte = ends_with(assembly_path, ".gz") + ? gz_fasta_seq_len_per_byte + : plain_fasta_seq_len_per_byte; + std::error_code ec; + const auto file_bytes = std::filesystem::file_size(assembly_path, ec); + if (ec) { + continue; + } + est_total_seq_len += static_cast(file_bytes * seq_len_per_byte); + } + return (2 * est_total_seq_len) / (windowsize + 1); +} + struct RawKmer { std::uint64_t hash; Kmer kmer; @@ -27,7 +63,6 @@ struct RawKmer { struct NodeState { std::size_t count = 0; - std::size_t start = 0; std::size_t cursor = 0; std::uint32_t n_tar = 0; std::uint32_t n_neg = 0; @@ -67,16 +102,15 @@ ThreadGraph build_worker( bool low_memory ) { // Estimate total minimizer count in all assemblies - const auto n_kmers_est = seqwin::est_kmer_number( - std::vector( - assembly_paths.begin() + static_cast(start_assembly), - assembly_paths.begin() + static_cast(end_assembly) - ), + const auto n_minimizers_est = estimate_minimizer_count( + assembly_paths, + start_assembly, + end_assembly, windowsize ); std::vector raw_kmers; if (!low_memory) { - raw_kmers.reserve(n_kmers_est); + raw_kmers.reserve(n_minimizers_est); } ThreadGraph graph; @@ -89,14 +123,14 @@ ThreadGraph build_worker( // Reserving for unordered_map will actually allocate physical memory std::unordered_map node_map; std::unordered_map edge_map; - const auto n_map_entries_est = n_kmers_est / map_reserve_divisor; + const auto n_map_entries_est = n_minimizers_est / map_reserve_divisor; node_map.reserve(n_map_entries_est); edge_map.reserve(n_map_entries_est); for (std::size_t assembly_i = start_assembly; assembly_i < end_assembly; ++assembly_i) { const bool is_target = is_targets[assembly_i]; - auto records = seqwin::read_fasta(assembly_paths[assembly_i]); + auto records = read_fasta(assembly_paths[assembly_i]); std::vector record_ids; record_ids.reserve(records.size()); @@ -176,7 +210,6 @@ ThreadGraph build_worker( std::size_t cursor = 0; for (auto& [hash, state] : node_map) { (void)hash; - state.start = cursor; state.cursor = cursor; cursor += state.count; } @@ -191,9 +224,11 @@ ThreadGraph build_worker( graph.nodes = NoInitArray(node_map.size()); std::size_t node_i = 0; for (const auto& [hash, state] : node_map) { + const std::size_t start = low_memory ? 0 : state.cursor - state.count; + graph.nodes[node_i++] = ThreadNode{ hash, - state.start, + start, state.count, state.n_tar, state.n_neg, @@ -233,7 +268,7 @@ NoInitArray recompute_kmers( for (std::size_t assembly_i = graph.start_assembly; assembly_i < graph.end_assembly; ++assembly_i) { - auto records = seqwin::read_fasta(assembly_paths[assembly_i]); + auto records = read_fasta(assembly_paths[assembly_i]); for (std::size_t record_i = 0; record_i < records.size(); ++record_i) { const auto record_idx = record_offsets[assembly_i] + record_i; @@ -272,6 +307,7 @@ NoInitArray recompute_kmers( } } // namespace +} // namespace internal Graph build( const std::vector& assembly_paths, @@ -294,8 +330,8 @@ Graph build( n_workers = std::min(n_workers, n_assemblies); } - ThreadPool pool(n_workers); // Avoid spawning threads every time - std::vector graphs(n_workers); + internal::ThreadPool pool(n_workers); // Avoid spawning threads every time + std::vector graphs(n_workers); const std::size_t base = n_assemblies / n_workers; const std::size_t rem = n_assemblies % n_workers; @@ -304,7 +340,7 @@ Graph build( for (std::size_t thread_id = start; thread_id < end; ++thread_id) { std::size_t start_assembly = thread_id * base + std::min(thread_id, rem); std::size_t end_assembly = start_assembly + base + (thread_id < rem ? 1 : 0); - graphs[thread_id] = build_worker( + graphs[thread_id] = internal::build_worker( assembly_paths, kmerlen, windowsize, @@ -317,15 +353,15 @@ Graph build( } }); - auto [graph, kmer_maps] = merge_thread_graphs( + auto [graph, kmer_maps] = internal::merge_thread_graphs( graphs, n_assemblies, pool, low_memory ); if (low_memory) { - log_python(" - Recomputing minimizers for low-memory mode..."); - graph.kmers = recompute_kmers( + internal::log_python(" - Recomputing minimizers for low-memory mode..."); + graph.kmers = internal::recompute_kmers( assembly_paths, kmerlen, windowsize, diff --git a/cpp/src/seqwin/helpers.cpp b/cpp/src/seqwin/build_internals.cpp similarity index 81% rename from cpp/src/seqwin/helpers.cpp rename to cpp/src/seqwin/build_internals.cpp index cc0a28c..27cf0e7 100644 --- a/cpp/src/seqwin/helpers.cpp +++ b/cpp/src/seqwin/build_internals.cpp @@ -1,5 +1,4 @@ -#include "seqwin/helpers.hpp" -#include "seqwin/thread_pool.hpp" +#include "seqwin/build_internals.hpp" #include #include @@ -9,11 +8,10 @@ #include #include -#include +#include "utils/logging.hpp" +#include "utils/thread_pool.hpp" -namespace py = pybind11; - -namespace seqwin { +namespace seqwin::internal { namespace { /** @@ -298,28 +296,6 @@ static void merge_edges(NoInitArray& edges, ThreadPool& pool) } // namespace -void log_python(const std::string& message, const std::string& level) -{ - py::gil_scoped_acquire acquire; - - py::object logging = py::module_::import("logging"); - py::object logger = logging.attr("getLogger")(); - - if (level == "debug") { - logger.attr("debug")(message); - } else if (level == "info") { - logger.attr("info")(message); - } else if (level == "warning" || level == "warn") { - logger.attr("warning")(message); - } else if (level == "error") { - logger.attr("error")(message); - } else if (level == "critical") { - logger.attr("critical")(message); - } else { - logger.attr("info")(message); - } -} - std::pair merge_thread_graphs( std::vector& graphs, std::size_t n_assemblies, @@ -423,68 +399,4 @@ std::pair merge_thread_graphs( }; } -Graph filter_kmers( - const Kmer* kmers, - const Node* nodes, - std::size_t n_nodes, - std::vector used_hashes -) { - std::sort(used_hashes.begin(), used_hashes.end()); - - Graph graph; - - std::vector used_node_indices; - used_node_indices.reserve(used_hashes.size()); - - std::size_t n_kmers = 0; - std::size_t node_i = 0; - std::size_t used_i = 0; - while (node_i < n_nodes && used_i < used_hashes.size()) { - const auto node_hash = nodes[node_i].hash; - const auto used_hash = used_hashes[used_i]; - - if (node_hash < used_hash) { - ++node_i; - continue; - } - if (used_hash < node_hash) { - ++used_i; - continue; - } - - used_node_indices.push_back(node_i); - n_kmers += nodes[node_i].stop - nodes[node_i].start; - ++node_i; - ++used_i; - } - - graph.nodes = NoInitArray(used_node_indices.size()); - graph.kmers = NoInitArray(n_kmers); - - std::size_t new_start = 0; - for (std::size_t out_node_i = 0; out_node_i < used_node_indices.size(); ++out_node_i) { - const auto in_node_i = used_node_indices[out_node_i]; - const Node& old_node = nodes[in_node_i]; - - const auto old_start = old_node.start; - const auto old_stop = old_node.stop; - const auto size = old_stop - old_start; - - Node new_node = old_node; - new_node.start = new_start; - new_node.stop = new_start + size; - graph.nodes[out_node_i] = new_node; - - for (std::size_t k = 0; k < size; ++k) { - const auto out_i = new_start + k; - const auto in_i = old_start + k; - graph.kmers[out_i] = kmers[in_i]; - } - - new_start += size; - } - - return graph; -} - -} // namespace seqwin +} // namespace seqwin::internal diff --git a/cpp/include/seqwin/helpers.hpp b/cpp/src/seqwin/build_internals.hpp similarity index 85% rename from cpp/include/seqwin/helpers.hpp rename to cpp/src/seqwin/build_internals.hpp index ac5f895..aea68ce 100644 --- a/cpp/include/seqwin/helpers.hpp +++ b/cpp/src/seqwin/build_internals.hpp @@ -9,7 +9,7 @@ #include "seqwin/graph.hpp" -namespace seqwin { +namespace seqwin::internal { class ThreadPool; @@ -65,17 +65,6 @@ struct ThreadGraph { */ using KmerMaps = std::vector>; -/** - * @brief Emit a message through Python's logging module. - * - * @param message Message to log. - * @param level Logging level: `debug`, `info`, `warning`, `error`, or `critical`. - */ -void log_python( - const std::string& message, - const std::string& level = "info" -); - /** * @brief Merge thread-local minimizer graphs into a single graph. * @@ -92,11 +81,4 @@ std::pair merge_thread_graphs( bool low_memory ); -Graph filter_kmers( - const Kmer* kmers, - const Node* nodes, - std::size_t n_nodes, - std::vector used_hashes -); - -} // namespace seqwin +} // namespace seqwin::internal diff --git a/cpp/src/seqwin/filter.cpp b/cpp/src/seqwin/filter.cpp new file mode 100644 index 0000000..fbbc0fa --- /dev/null +++ b/cpp/src/seqwin/filter.cpp @@ -0,0 +1,74 @@ +#include "seqwin/filter.hpp" + +#include +#include +#include +#include + +namespace seqwin { + +Graph filter_kmers( + const Kmer* kmers, + const Node* nodes, + std::size_t n_nodes, + std::vector used_hashes +) { + std::sort(used_hashes.begin(), used_hashes.end()); + + Graph graph; + + std::vector used_node_indices; + used_node_indices.reserve(used_hashes.size()); + + std::size_t n_kmers = 0; + std::size_t node_i = 0; + std::size_t used_i = 0; + while (node_i < n_nodes && used_i < used_hashes.size()) { + const auto node_hash = nodes[node_i].hash; + const auto used_hash = used_hashes[used_i]; + + if (node_hash < used_hash) { + ++node_i; + continue; + } + if (used_hash < node_hash) { + ++used_i; + continue; + } + + used_node_indices.push_back(node_i); + n_kmers += nodes[node_i].stop - nodes[node_i].start; + ++node_i; + ++used_i; + } + + graph.nodes = NoInitArray(used_node_indices.size()); + graph.kmers = NoInitArray(n_kmers); + + std::size_t new_start = 0; + for (std::size_t out_node_i = 0; out_node_i < used_node_indices.size(); ++out_node_i) { + const auto in_node_i = used_node_indices[out_node_i]; + const Node& old_node = nodes[in_node_i]; + + const auto old_start = old_node.start; + const auto old_stop = old_node.stop; + const auto size = old_stop - old_start; + + Node new_node = old_node; + new_node.start = new_start; + new_node.stop = new_start + size; + graph.nodes[out_node_i] = new_node; + + for (std::size_t k = 0; k < size; ++k) { + const auto out_i = new_start + k; + const auto in_i = old_start + k; + graph.kmers[out_i] = kmers[in_i]; + } + + new_start += size; + } + + return graph; +} + +} // namespace seqwin diff --git a/cpp/src/seqwin/fasta_reader.cpp b/cpp/src/utils/fasta_reader.cpp similarity index 86% rename from cpp/src/seqwin/fasta_reader.cpp rename to cpp/src/utils/fasta_reader.cpp index 720dc32..90b4ee0 100644 --- a/cpp/src/seqwin/fasta_reader.cpp +++ b/cpp/src/utils/fasta_reader.cpp @@ -1,9 +1,8 @@ -#include "seqwin/fasta_reader.hpp" +#include "utils/fasta_reader.hpp" -#include #include +#include #include -#include #include #include #include @@ -13,11 +12,9 @@ #include -namespace seqwin { +namespace seqwin::internal { namespace { -constexpr std::size_t plain_fasta_seq_len_per_byte = 1; -constexpr std::size_t gz_fasta_seq_len_per_byte = 4; constexpr std::size_t gz_read_buf_size = 1U << 16; bool ends_with(std::string_view text, std::string_view suffix) @@ -215,24 +212,4 @@ std::vector read_fasta(const std::string& assembly_path) return read_plain_fasta(assembly_path); } -std::size_t est_kmer_number( - const std::vector& assembly_paths, - std::size_t windowsize -) { - // Reserve-only heuristic, not correctness-critical. - std::size_t est_total_seq_len = 0; - for (const auto& assembly_path : assembly_paths) { - const std::size_t seq_len_per_byte = ends_with(assembly_path, ".gz") - ? gz_fasta_seq_len_per_byte - : plain_fasta_seq_len_per_byte; - std::error_code ec; - const auto file_bytes = std::filesystem::file_size(assembly_path, ec); - if (ec) { - continue; - } - est_total_seq_len += static_cast(file_bytes * seq_len_per_byte); - } - return (2 * est_total_seq_len) / (windowsize + 1); -} - -} // namespace seqwin +} // namespace seqwin::internal diff --git a/cpp/include/seqwin/fasta_reader.hpp b/cpp/src/utils/fasta_reader.hpp similarity index 51% rename from cpp/include/seqwin/fasta_reader.hpp rename to cpp/src/utils/fasta_reader.hpp index 9db3b3c..ef4934e 100644 --- a/cpp/include/seqwin/fasta_reader.hpp +++ b/cpp/src/utils/fasta_reader.hpp @@ -4,7 +4,7 @@ #include #include -namespace seqwin { +namespace seqwin::internal { struct FastaRecord { std::string id; @@ -20,16 +20,4 @@ struct FastaRecord { */ std::vector read_fasta(const std::string& assembly_path); -/** - * @brief Estimate the number of minimizers based on the size of assembly files. - * - * @param assembly_paths Paths to assembly FASTA files (plain or gzipped). - * @param windowsize Minimizer window size used for the estimate. - * @return Estimated total number of minimizers for all assembly files. - */ -std::size_t est_kmer_number( - const std::vector& assembly_paths, - std::size_t windowsize -); - -} // namespace seqwin +} // namespace seqwin::internal diff --git a/cpp/src/utils/logging.cpp b/cpp/src/utils/logging.cpp new file mode 100644 index 0000000..62cba8b --- /dev/null +++ b/cpp/src/utils/logging.cpp @@ -0,0 +1,31 @@ +#include "utils/logging.hpp" + +#include + +namespace py = pybind11; + +namespace seqwin::internal { + +void log_python(const std::string& message, const std::string& level) +{ + py::gil_scoped_acquire acquire; + + py::object logging = py::module_::import("logging"); + py::object logger = logging.attr("getLogger")(); + + if (level == "debug") { + logger.attr("debug")(message); + } else if (level == "info") { + logger.attr("info")(message); + } else if (level == "warning" || level == "warn") { + logger.attr("warning")(message); + } else if (level == "error") { + logger.attr("error")(message); + } else if (level == "critical") { + logger.attr("critical")(message); + } else { + logger.attr("info")(message); + } +} + +} // namespace seqwin::internal diff --git a/cpp/src/utils/logging.hpp b/cpp/src/utils/logging.hpp new file mode 100644 index 0000000..2421f8b --- /dev/null +++ b/cpp/src/utils/logging.hpp @@ -0,0 +1,18 @@ +#pragma once + +#include + +namespace seqwin::internal { + +/** + * @brief Emit a message through Python's logging module. + * + * @param message Message to log. + * @param level Logging level: `debug`, `info`, `warning`, `error`, or `critical`. + */ +void log_python( + const std::string& message, + const std::string& level = "info" +); + +} // namespace seqwin::internal diff --git a/cpp/include/seqwin/thread_pool.hpp b/cpp/src/utils/thread_pool.hpp similarity index 98% rename from cpp/include/seqwin/thread_pool.hpp rename to cpp/src/utils/thread_pool.hpp index dba4672..9a0d3aa 100644 --- a/cpp/include/seqwin/thread_pool.hpp +++ b/cpp/src/utils/thread_pool.hpp @@ -11,7 +11,7 @@ #include #include -namespace seqwin { +namespace seqwin::internal { /** * @brief Simple fixed-size worker pool for parallel range processing. @@ -155,4 +155,4 @@ class ThreadPool { bool stopping_ = false; }; -} // namespace seqwin +} // namespace seqwin::internal