From d637b3f8f4282cc98e2569be0e38376108a8bfbc Mon Sep 17 00:00:00 2001 From: mxwang66 <54554363+mxwang66@users.noreply.github.com> Date: Mon, 25 May 2026 01:57:52 -0500 Subject: [PATCH 1/2] Optimize `kmers` and `idx` allocation with `NoInitArray` Avoid value-initializing large arrays. --- cpp/include/seqwin/graph.hpp | 8 +-- cpp/include/seqwin/helpers.hpp | 12 ++++- cpp/include/seqwin/no_init_array.hpp | 78 ++++++++++++++++++++++++++++ cpp/src/seqwin/graph.cpp | 8 +-- cpp/src/seqwin/helpers.cpp | 57 ++++++++++++-------- 5 files changed, 133 insertions(+), 30 deletions(-) create mode 100644 cpp/include/seqwin/no_init_array.hpp diff --git a/cpp/include/seqwin/graph.hpp b/cpp/include/seqwin/graph.hpp index 6706f42..06a88ae 100644 --- a/cpp/include/seqwin/graph.hpp +++ b/cpp/include/seqwin/graph.hpp @@ -5,6 +5,8 @@ #include #include +#include "seqwin/no_init_array.hpp" + namespace seqwin { struct Kmer { @@ -29,13 +31,11 @@ struct Edge { }; struct Graph { - std::vector kmers; - std::vector idx; + NoInitArray kmers; + NoInitArray idx; std::vector nodes; std::vector edges; std::vector> ids_by_assembly; - std::size_t n_kmers = 0; - std::size_t start_assembly = 0; }; Graph build( diff --git a/cpp/include/seqwin/helpers.hpp b/cpp/include/seqwin/helpers.hpp index 9b0fe8c..125f9ad 100644 --- a/cpp/include/seqwin/helpers.hpp +++ b/cpp/include/seqwin/helpers.hpp @@ -16,13 +16,23 @@ namespace seqwin { */ class ThreadPool; +struct ThreadGraph { + std::vector kmers; + NoInitArray idx; + std::vector nodes; + std::vector edges; + std::vector> ids_by_assembly; + std::size_t n_kmers = 0; + std::size_t start_assembly = 0; +}; + void log_python( const std::string& message, const std::string& level = "info" ); Graph merge_thread_graphs( - std::vector& graphs, + std::vector& graphs, std::size_t n_assemblies, ThreadPool& pool ); diff --git a/cpp/include/seqwin/no_init_array.hpp b/cpp/include/seqwin/no_init_array.hpp new file mode 100644 index 0000000..026a38b --- /dev/null +++ b/cpp/include/seqwin/no_init_array.hpp @@ -0,0 +1,78 @@ +#pragma once + +#include +#include +#include + +namespace seqwin { + +/** + * A fixed-size owning array. + * + * Unlike `std::vector(n) or std::make_unique(n)`, this class allocates + * with `new T[n]`. For scalar and trivially default-initialized element types, + * this avoids value-initializing every element, which can be expensive for very + * large arrays. + * + * Important: + * - Every element must be assigned before it is read. + * - This is not a full `std::vector` replacement. + * - It intentionally provides no `resize()`, `reserve()`, `push_back()`, or `capacity()`. + */ +template +class NoInitArray { +public: + NoInitArray() noexcept = default; + + explicit NoInitArray(std::size_t size) + : size_(size) + , data_(size == 0 ? nullptr : new T[size]) + {} + + NoInitArray(const NoInitArray&) = delete; + NoInitArray& operator=(const NoInitArray&) = delete; + + NoInitArray(NoInitArray&& other) noexcept + : size_(other.size_) + , data_(std::move(other.data_)) + { + other.size_ = 0; + } + + NoInitArray& operator=(NoInitArray&& other) noexcept + { + if (this != &other) { + data_ = std::move(other.data_); + size_ = other.size_; + other.size_ = 0; + } + return *this; + } + + std::size_t size() const noexcept { return size_; } + bool empty() const noexcept { return size_ == 0; } + + T* data() noexcept { return data_.get(); } + const T* data() const noexcept { return data_.get(); } + + T& operator[](std::size_t i) noexcept { return data_[i]; } + const T& operator[](std::size_t i) const noexcept { return data_[i]; } + + void swap(NoInitArray& other) noexcept + { + std::swap(size_, other.size_); + std::swap(data_, other.data_); + } + + void reset() noexcept + { + data_.reset(); + size_ = 0; + } + +private: + std::size_t size_ = 0; + std::unique_ptr data_; +}; + +} // namespace seqwin diff --git a/cpp/src/seqwin/graph.cpp b/cpp/src/seqwin/graph.cpp index 19b1a35..430d799 100644 --- a/cpp/src/seqwin/graph.cpp +++ b/cpp/src/seqwin/graph.cpp @@ -43,7 +43,7 @@ struct EdgeState { std::size_t last_seen_assembly = std::numeric_limits::max(); }; -Graph build_worker( +ThreadGraph build_worker( const std::vector& assembly_paths, std::size_t kmerlen, std::size_t windowsize, @@ -61,7 +61,7 @@ Graph build_worker( windowsize ); - Graph graph; + ThreadGraph graph; graph.kmers.reserve(n_kmers_est); graph.ids_by_assembly.resize(end_assembly - start_assembly); graph.start_assembly = start_assembly; @@ -151,7 +151,7 @@ Graph build_worker( } std::unordered_map().swap(edge_map); - graph.idx.resize(graph.n_kmers); + graph.idx = NoInitArray(graph.n_kmers); std::uint64_t cursor = 0; for (auto& [hash, state] : node_map) { (void)hash; @@ -210,7 +210,7 @@ Graph build( } ThreadPool pool(n_workers); - std::vector graphs(n_workers); + std::vector graphs(n_workers); const std::size_t base = n_assemblies / n_workers; const std::size_t rem = n_assemblies % n_workers; diff --git a/cpp/src/seqwin/helpers.cpp b/cpp/src/seqwin/helpers.cpp index fc1959a..4bb5348 100644 --- a/cpp/src/seqwin/helpers.cpp +++ b/cpp/src/seqwin/helpers.cpp @@ -21,7 +21,7 @@ struct IdxSegment { }; template -std::vector concat(std::vector& graphs, MemberPtr member, ThreadPool& pool) +std::vector concat(std::vector& graphs, MemberPtr member, ThreadPool& pool) { std::vector out; if (graphs.empty()) { @@ -54,14 +54,14 @@ std::vector concat(std::vector& graphs, MemberPtr member, ThreadPool& return out; } -std::vector concat_nodes(std::vector& graphs, ThreadPool& pool) +std::vector concat_nodes(std::vector& graphs, ThreadPool& pool) { - return concat(graphs, &Graph::nodes, pool); + return concat(graphs, &ThreadGraph::nodes, pool); } -std::vector concat_edges(std::vector& graphs, ThreadPool& pool) +std::vector concat_edges(std::vector& graphs, ThreadPool& pool) { - return concat(graphs, &Graph::edges, pool); + return concat(graphs, &ThreadGraph::edges, pool); } // 1. Parallel LSD radix sort based on hash (stable) @@ -154,13 +154,13 @@ static std::vector merge_nodes( return idx_segments; } -static std::vector merge_kmers( +static NoInitArray merge_kmers( const std::vector& idx_segments, - std::vector& graphs, + std::vector& graphs, std::size_t total_kmers, ThreadPool& pool ) { - std::vector kmers(total_kmers); + NoInitArray kmers(total_kmers); pool.parallel_for(idx_segments.size(), [&](std::size_t start, std::size_t end, std::size_t) { for (std::size_t s = start; s < end; ++s) { @@ -177,14 +177,14 @@ static std::vector merge_kmers( return kmers; } -static std::vector merge_idx( +static NoInitArray merge_idx( const std::vector& idx_segments, const std::vector& kmer_offsets, - const std::vector& graphs, + const std::vector& graphs, std::size_t total_kmers, ThreadPool& pool ) { - std::vector idx(total_kmers); + NoInitArray idx(total_kmers); pool.parallel_for(idx_segments.size(), [&](std::size_t start, std::size_t end, std::size_t) { for (std::size_t s = start; s < end; ++s) { @@ -290,7 +290,7 @@ void log_python(const std::string& message, const std::string& level) } Graph merge_thread_graphs( - std::vector& graphs, + std::vector& graphs, std::size_t n_assemblies, ThreadPool& pool ) { @@ -316,7 +316,7 @@ Graph merge_thread_graphs( graph.n_kmers, pool ); - std::vector().swap(graph.idx); + graph.idx.reset(); return { std::move(kmers), @@ -350,7 +350,7 @@ Graph merge_thread_graphs( auto idx = merge_idx(idx_segments, kmer_offsets, graphs, total_kmers, pool); for (auto& graph : graphs) { - std::vector().swap(graph.idx); + graph.idx.reset(); } std::vector> ids_by_assembly(n_assemblies); @@ -360,7 +360,7 @@ Graph merge_thread_graphs( } } - std::vector().swap(graphs); + std::vector().swap(graphs); return { std::move(kmers), @@ -383,10 +383,12 @@ Graph filter_kmers( Graph graph; graph.nodes.reserve(used_hashes.size()); + 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; - std::uint64_t new_start = 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]; @@ -400,6 +402,17 @@ Graph filter_kmers( continue; } + used_node_indices.push_back(node_i); + n_kmers += static_cast(nodes[node_i].stop - nodes[node_i].start); + ++node_i; + ++used_i; + } + + graph.kmers = NoInitArray(n_kmers); + graph.idx = NoInitArray(n_kmers); + + std::uint64_t new_start = 0; + for (const auto node_i : used_node_indices) { const Node& old_node = nodes[node_i]; const auto old_start = old_node.start; const auto old_stop = old_node.stop; @@ -410,12 +423,14 @@ Graph filter_kmers( new_node.stop = new_start + size; graph.nodes.push_back(new_node); - graph.kmers.insert(graph.kmers.end(), kmers + old_start, kmers + old_stop); - graph.idx.insert(graph.idx.end(), idx + old_start, idx + old_stop); + for (std::uint64_t k = 0; k < size; ++k) { + const auto out_i = static_cast(new_start + k); + const auto in_i = static_cast(old_start + k); + graph.kmers[out_i] = kmers[in_i]; + graph.idx[out_i] = idx[in_i]; + } new_start += size; - ++node_i; - ++used_i; } return graph; From 8dcb6ded86abd799a19fb0004871b5b7cc04366e Mon Sep 17 00:00:00 2001 From: mxwang66 <54554363+mxwang66@users.noreply.github.com> Date: Mon, 25 May 2026 16:49:38 -0500 Subject: [PATCH 2/2] Use `NoInitArray` for `nodes` and `edges` --- cpp/include/seqwin/graph.hpp | 4 +- cpp/include/seqwin/helpers.hpp | 6 +- cpp/include/seqwin/no_init_array.hpp | 8 ++ cpp/src/seqwin/graph.cpp | 4 +- cpp/src/seqwin/helpers.cpp | 122 +++++++++++++++++---------- 5 files changed, 91 insertions(+), 53 deletions(-) diff --git a/cpp/include/seqwin/graph.hpp b/cpp/include/seqwin/graph.hpp index 06a88ae..50c870d 100644 --- a/cpp/include/seqwin/graph.hpp +++ b/cpp/include/seqwin/graph.hpp @@ -33,8 +33,8 @@ struct Edge { struct Graph { NoInitArray kmers; NoInitArray idx; - std::vector nodes; - std::vector edges; + NoInitArray nodes; + NoInitArray edges; std::vector> ids_by_assembly; }; diff --git a/cpp/include/seqwin/helpers.hpp b/cpp/include/seqwin/helpers.hpp index 125f9ad..b881c3f 100644 --- a/cpp/include/seqwin/helpers.hpp +++ b/cpp/include/seqwin/helpers.hpp @@ -17,10 +17,10 @@ namespace seqwin { class ThreadPool; struct ThreadGraph { - std::vector kmers; + std::vector kmers; // Needs push_back() support (exact kmer count is not known) NoInitArray idx; - std::vector nodes; - std::vector edges; + NoInitArray nodes; + NoInitArray edges; std::vector> ids_by_assembly; std::size_t n_kmers = 0; std::size_t start_assembly = 0; diff --git a/cpp/include/seqwin/no_init_array.hpp b/cpp/include/seqwin/no_init_array.hpp index 026a38b..5130705 100644 --- a/cpp/include/seqwin/no_init_array.hpp +++ b/cpp/include/seqwin/no_init_array.hpp @@ -54,6 +54,12 @@ class NoInitArray { T* data() noexcept { return data_.get(); } const T* data() const noexcept { return data_.get(); } + T* begin() noexcept { return data_.get(); } + T* end() noexcept { return data_.get() + size_; } + const T* begin() const noexcept { return data_.get(); } + const T* end() const noexcept { return data_.get() + size_; } + const T* cbegin() const noexcept { return data_.get(); } + const T* cend() const noexcept { return data_.get() + size_; } T& operator[](std::size_t i) noexcept { return data_[i]; } const T& operator[](std::size_t i) const noexcept { return data_[i]; } @@ -64,6 +70,8 @@ class NoInitArray { std::swap(data_, other.data_); } + friend void swap(NoInitArray& a, NoInitArray& b) noexcept { a.swap(b); } + void reset() noexcept { data_.reset(); diff --git a/cpp/src/seqwin/graph.cpp b/cpp/src/seqwin/graph.cpp index 430d799..881a54b 100644 --- a/cpp/src/seqwin/graph.cpp +++ b/cpp/src/seqwin/graph.cpp @@ -144,7 +144,7 @@ ThreadGraph build_worker( } // Create edges first to reduce peak memory - graph.edges.resize(edge_map.size()); + graph.edges = NoInitArray(edge_map.size()); std::size_t edge_i = 0; for (const auto& [key, state] : edge_map) { graph.edges[edge_i++] = Edge{key.first, key.second, state.weight}; @@ -165,7 +165,7 @@ ThreadGraph build_worker( graph.idx[node_it->second.cursor++] = static_cast(i); } - graph.nodes.resize(node_map.size()); + graph.nodes = NoInitArray(node_map.size()); std::size_t node_i = 0; for (const auto& [hash, state] : node_map) { const auto stop = state.start + state.count; diff --git a/cpp/src/seqwin/helpers.cpp b/cpp/src/seqwin/helpers.cpp index 4bb5348..48799ff 100644 --- a/cpp/src/seqwin/helpers.cpp +++ b/cpp/src/seqwin/helpers.cpp @@ -21,11 +21,10 @@ struct IdxSegment { }; template -std::vector concat(std::vector& graphs, MemberPtr member, ThreadPool& pool) +NoInitArray concat(std::vector& graphs, MemberPtr member, ThreadPool& pool) { - std::vector out; if (graphs.empty()) { - return out; + return {}; } std::vector offsets(graphs.size(), 0); @@ -34,7 +33,7 @@ std::vector concat(std::vector& graphs, MemberPtr member, Thread offsets[i] = cursor; cursor += (graphs[i].*member).size(); } - out.resize(cursor); + NoInitArray out(cursor); pool.parallel_for(graphs.size(), [&](std::size_t start, std::size_t end, std::size_t) { for (std::size_t i = start; i < end; ++i) { @@ -47,28 +46,30 @@ std::vector concat(std::vector& graphs, MemberPtr member, Thread out.begin() + static_cast(offsets[i]) ); } - std::vector().swap(local); + local.reset(); } }); return out; } -std::vector concat_nodes(std::vector& graphs, ThreadPool& pool) +NoInitArray concat_nodes(std::vector& graphs, ThreadPool& pool) { return concat(graphs, &ThreadGraph::nodes, pool); } -std::vector concat_edges(std::vector& graphs, ThreadPool& pool) +NoInitArray concat_edges(std::vector& graphs, ThreadPool& pool) { return concat(graphs, &ThreadGraph::edges, pool); } -// 1. Parallel LSD radix sort based on hash (stable) -// 2. Merge nodes with the same hash -// 3. Build segment metadata for final materialization +/** + * 1. Parallel LSD radix sort based on hash (stable) + * 2. Merge nodes with the same hash + * 3. Build segment metadata for final materialization + */ static std::vector merge_nodes( - std::vector& nodes, + NoInitArray& nodes, ThreadPool& pool ) { const std::size_t n_nodes = nodes.size(); @@ -76,9 +77,9 @@ static std::vector merge_nodes( return {}; } - std::vector buf(nodes.size()); - auto* src = &nodes; - auto* dst = &buf; + NoInitArray buf(nodes.size()); + auto* src = nodes.data(); + auto* dst = buf.data(); std::vector counts(pool.size() * 65536); for (std::size_t shift = 0; shift < 64; shift += 16) { @@ -87,7 +88,7 @@ static std::vector merge_nodes( pool.parallel_for(n_nodes, [&](std::size_t start, std::size_t end, std::size_t t) { auto* local_counts = counts.data() + t * 65536; for (std::size_t i = start; i < end; ++i) { - const auto bucket = ((*src)[i].hash >> shift) & 0xFFFFULL; + const auto bucket = (src[i].hash >> shift) & 0xFFFFULL; ++local_counts[static_cast(bucket)]; } }); @@ -105,39 +106,52 @@ static std::vector merge_nodes( pool.parallel_for(n_nodes, [&](std::size_t start, std::size_t end, std::size_t t) { auto* local_offsets = counts.data() + t * 65536; for (std::size_t i = start; i < end; ++i) { - const auto bucket = ((*src)[i].hash >> shift) & 0xFFFFULL; + const auto bucket = (src[i].hash >> shift) & 0xFFFFULL; const auto pos = local_offsets[static_cast(bucket)]++; - (*dst)[pos] = (*src)[i]; + dst[pos] = src[i]; } }); std::swap(src, dst); } + // Determine final node count + std::size_t unique_count = 0; + std::size_t i = 0; + while (i < n_nodes) { + const auto hash = src[i].hash; + while (i < n_nodes && src[i].hash == hash) { + ++i; + } + ++unique_count; + } + buf.reset(); + buf = NoInitArray(unique_count); + // Aggregate nodes and keep idx ranges std::vector idx_segments; idx_segments.reserve(n_nodes); std::uint64_t n_kmers = 0; std::size_t write_i = 0; - std::size_t i = 0; + i = 0; while (i < n_nodes) { - const auto hash = (*src)[i].hash; + const auto hash = src[i].hash; std::uint32_t n_tar = 0; std::uint32_t n_neg = 0; const auto start = n_kmers; - while (i < n_nodes && (*src)[i].hash == hash) { - n_tar += (*src)[i].n_tar; - n_neg += (*src)[i].n_neg; + while (i < n_nodes && src[i].hash == hash) { + n_tar += src[i].n_tar; + n_neg += src[i].n_neg; - const auto local_start = (*src)[i].start; - const auto local_stop = (*src)[i].stop; + const auto local_start = src[i].start; + const auto local_stop = src[i].stop; const auto length = local_stop - local_start; if (length != 0) { idx_segments.push_back(IdxSegment{ - static_cast((*src)[i].penalty), // thread_id + static_cast(src[i].penalty), // thread_id static_cast(local_start), static_cast(n_kmers), static_cast(length) @@ -148,9 +162,9 @@ static std::vector merge_nodes( } const auto stop = n_kmers; - nodes[write_i++] = Node{hash, start, stop, n_tar, n_neg, 0.0}; + buf[write_i++] = Node{hash, start, stop, n_tar, n_neg, 0.0}; } - nodes.resize(write_i); + nodes = std::move(buf); return idx_segments; } @@ -200,17 +214,17 @@ static NoInitArray merge_idx( return idx; } -static void merge_edges(std::vector& edges, ThreadPool& pool) +static void merge_edges(NoInitArray& edges, ThreadPool& pool) { const std::size_t n_edges = edges.size(); if (n_edges == 0) { - edges.clear(); + edges.reset(); return; } - std::vector buf(edges.size()); - auto* src = &edges; - auto* dst = &buf; + NoInitArray buf(edges.size()); + auto* src = edges.data(); + auto* dst = buf.data(); std::vector counts(pool.size() * 65536); for (auto key : {&Edge::second, &Edge::first}) { @@ -220,7 +234,7 @@ static void merge_edges(std::vector& edges, ThreadPool& pool) pool.parallel_for(n_edges, [&](std::size_t start, std::size_t end, std::size_t t) { auto* local_counts = counts.data() + t * 65536; for (std::size_t i = start; i < end; ++i) { - const auto bucket = (((*src)[i].*key) >> shift) & 0xFFFFULL; + const auto bucket = ((src[i].*key) >> shift) & 0xFFFFULL; ++local_counts[static_cast(bucket)]; } }); @@ -238,9 +252,9 @@ static void merge_edges(std::vector& edges, ThreadPool& pool) pool.parallel_for(n_edges, [&](std::size_t start, std::size_t end, std::size_t t) { auto* local_offsets = counts.data() + t * 65536; for (std::size_t i = start; i < end; ++i) { - const auto bucket = (((*src)[i].*key) >> shift) & 0xFFFFULL; + const auto bucket = ((src[i].*key) >> shift) & 0xFFFFULL; const auto pos = local_offsets[static_cast(bucket)]++; - (*dst)[pos] = (*src)[i]; + dst[pos] = src[i]; } }); @@ -248,21 +262,35 @@ static void merge_edges(std::vector& edges, ThreadPool& pool) } } - std::size_t write_i = 0; + // Determine final edge count + std::size_t unique_count = 0; std::size_t i = 0; while (i < n_edges) { - const auto first = (*src)[i].first; - const auto second = (*src)[i].second; + const auto first = src[i].first; + const auto second = src[i].second; + while (i < n_edges && src[i].first == first && src[i].second == second) { + ++i; + } + ++unique_count; + } + buf.reset(); + buf = NoInitArray(unique_count); + + std::size_t write_i = 0; + i = 0; + while (i < n_edges) { + const auto first = src[i].first; + const auto second = src[i].second; std::uint64_t weight = 0; - while (i < n_edges && (*src)[i].first == first && (*src)[i].second == second) { - weight += (*src)[i].weight; + while (i < n_edges && src[i].first == first && src[i].second == second) { + weight += src[i].weight; ++i; } - edges[write_i++] = Edge{first, second, weight}; + buf[write_i++] = Edge{first, second, weight}; } - edges.resize(write_i); + edges = std::move(buf); } } // namespace @@ -381,7 +409,6 @@ Graph filter_kmers( std::sort(used_hashes.begin(), used_hashes.end()); Graph graph; - graph.nodes.reserve(used_hashes.size()); std::vector used_node_indices; used_node_indices.reserve(used_hashes.size()); @@ -408,12 +435,15 @@ Graph filter_kmers( ++used_i; } + graph.nodes = NoInitArray(used_node_indices.size()); graph.kmers = NoInitArray(n_kmers); graph.idx = NoInitArray(n_kmers); std::uint64_t new_start = 0; - for (const auto node_i : used_node_indices) { - const Node& old_node = nodes[node_i]; + 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; @@ -421,7 +451,7 @@ Graph filter_kmers( Node new_node = old_node; new_node.start = new_start; new_node.stop = new_start + size; - graph.nodes.push_back(new_node); + graph.nodes[out_node_i] = new_node; for (std::uint64_t k = 0; k < size; ++k) { const auto out_i = static_cast(new_start + k);