Streamline CLI, data structures, and filtering logic

cpp/include/seqwin/filter.hpp

@@ -12,7 +12,7 @@ namespace seqwin {
 using Subgraphs = std::vector<std::vector<std::uint64_t>>;
 
 /**
- * @brief Graph-filtering configurations.
+ * @brief Part of Seqwin configurations.
  */
 struct FilterConfig {
     std::optional<double> penalty_th;
@@ -28,13 +28,17 @@ struct FilterConfig {
 };
 
 /**
- * @brief Graph-filtering results.
+ * @brief Includes filtered graph arrays, low-penalty subgraphs and calculated values.
  */
 struct FilterResult {
     NoInitArray<Kmer> kmers;
     NoInitArray<Node> nodes;
     std::vector<Edge> edges;
     Subgraphs subgraphs;
+    std::size_t total_tar;
+    std::size_t total_neg;
+    double e_absence_tar;
+    double e_presence_neg;
     double penalty_th;
     double edge_weight_th;
     std::size_t min_nodes;

cpp/include/seqwin/graph.hpp

@@ -2,13 +2,90 @@
 
 #include <cstddef>
 #include <cstdint>
+#include <memory>
 #include <string>
+#include <utility>
 #include <vector>
 
-#include "seqwin/no_init_array.hpp"
-
 namespace seqwin {
 
+/**
+ * @brief Fixed-size owning array that avoids value-initializing elements.
+ *
+ * Unlike `std::vector<T>(n)` or `std::make_unique<T[]>(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 <typename T>
+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* 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]; }
+
+    void swap(NoInitArray& other) noexcept
+    {
+        std::swap(size_, other.size_);
+        std::swap(data_, other.data_);
+    }
+
+    friend void swap(NoInitArray& a, NoInitArray& b) noexcept { a.swap(b); }
+
+    void reset() noexcept
+    {
+        data_.reset();
+        size_ = 0;
+    }
+
+private:
+    std::size_t size_ = 0;
+    std::unique_ptr<T[]> data_;
+};
+
 /**
  * @brief Location metadata for a minimizer.
  */
@@ -32,11 +109,20 @@ struct Node {
     std::size_t start;
     /** End of the half-open range for this node's minimizer entries. */
     std::size_t stop;
-    /** Node scoring placeholder. */
+    /**
+     * Number of target assemblies containing this node's minimizer.
+     * Initialized to 0; populated during the filtering phase.
+     */
     std::uint32_t n_tar = 0;
-    /** Node scoring placeholder. */
+    /**
+     * Number of non-target assemblies containing this node's minimizer.
+     * Initialized to 0; populated during the filtering phase.
+     */
     std::uint32_t n_neg = 0;
-    /** Node scoring placeholder. */
+    /**
+     * Node penalty score.
+     * Initialized to 0; calculated from `n_tar` and `n_neg`.
+     */
     double penalty = 0.0;
 };
 

cpp/src/bindings/python_bindings.cpp

@@ -174,6 +174,8 @@ PYBIND11_MODULE(_core, m) {
                 array_to_numpy(std::move(result.nodes)),
                 array_to_numpy(std::move(result.edges)),
                 std::move(result.subgraphs),
+                result.total_tar,
+                result.total_neg,
                 result.penalty_th,
                 result.edge_weight_th,
                 result.min_nodes,

cpp/src/seqwin/filter.cpp

@@ -4,7 +4,6 @@
 #include <cmath>
 #include <cstddef>
 #include <cstdint>
-#include <limits>
 #include <sstream>
 #include <stdexcept>
 #include <string>
@@ -69,25 +68,21 @@ double expected_presence(
     return sum / static_cast<double>(count);
 }
 
-FilterResult calculate_thresholds(
-    const Node* nodes,
-    std::size_t n_nodes,
+/** @brief Calculate thresholds and add them to `result`. */
+void calculate_thresholds(
     const bool* is_targets,
     std::size_t n_assemblies,
     const double* jaccard,
     std::size_t jaccard_rows,
     std::size_t jaccard_cols,
-    const FilterConfig& config
+    const FilterConfig& config,
+    FilterResult& result
 ) {
-    const std::size_t total_tar = std::count(is_targets, is_targets + n_assemblies, true);
-    const auto total_neg = n_assemblies - total_tar;
-
     double penalty_th;
     if (config.penalty_th) {
         penalty_th = *config.penalty_th;
         internal::log_python("Penalty threshold is provided (--penalty-th), skip auto estimation", "warning");
     } else {
-        internal::log_python(" - Calculating penalty threshold...");
         // Consider k-mers in target assemblies:
         double e_absence_tar; // their expected absence in target assemblies
         double e_presence_neg; // their expected presence in non-target assemblies
@@ -98,23 +93,9 @@ FilterResult calculate_thresholds(
             e_absence_tar = 1.0 - expected_presence(jaccard, n_assemblies, is_targets, true);
             e_presence_neg = expected_presence(jaccard, n_assemblies, is_targets, false);
         } else {
-            // Calculate expected presence from minimizer sketches
-            // For all k-mers in targets, calculate their average presence in targets or non-targets
-            double sum_n_tar = 0.0; // Number of k-mers in all targets
-            double sum_presence_tar = 0.0;
-            double sum_presence_neg = 0.0;
-            for (std::size_t i = 0; i < n_nodes; ++i) {
-                const double node_n_tar = nodes[i].n_tar;
-                const double node_n_neg = nodes[i].n_neg;
-                sum_n_tar += node_n_tar;
-                sum_presence_tar += (node_n_tar / total_tar) * node_n_tar;
-                sum_presence_neg += (node_n_neg / total_neg) * node_n_tar;
-            }
-            if (sum_n_tar == 0.0) {
-                throw std::invalid_argument("No target minimizers are available for threshold estimation");
-            }
-            e_absence_tar = 1.0 - sum_presence_tar / sum_n_tar;
-            e_presence_neg = sum_presence_neg / sum_n_tar;
+            // Use values calculated by `get_penalty()`
+            e_absence_tar = result.e_absence_tar;
+            e_presence_neg = result.e_presence_neg;
         }
         internal::log_python(" - Expected k-mer absence in targets: " + format_value(e_absence_tar, 5));
         internal::log_python(" - Expected k-mer presence in non-targets: " + format_value(e_presence_neg, 5));
@@ -134,7 +115,7 @@ FilterResult calculate_thresholds(
     // Consider N as the number of assemblies that include a certain k-mer. Since we want k-mers with
     // penalty lower than penalty_th, based on the definition of penalty, N ≥ (1 - penalty_th) * total_tar.
     // So edge weight threshold is calculated based on the lower bound of N, times a multiplier < 1.
-    const double edge_weight_th = config.edge_w_th_mul * (1.0 - penalty_th) * total_tar;
+    const double edge_weight_th = config.edge_w_th_mul * (1.0 - penalty_th) * result.total_tar;
 
     // Calculate size range of subgraphs
     const std::size_t gap_len = (config.windowsize + 1) / 2;
@@ -153,12 +134,10 @@ FilterResult calculate_thresholds(
         );
     }
 
-    FilterResult result;
     result.penalty_th = penalty_th;
     result.edge_weight_th = edge_weight_th;
     result.min_nodes = min_nodes;
     result.max_nodes = max_nodes;
-    return result;
 }
 
 } // namespace
@@ -179,11 +158,11 @@ FilterResult filter(
     const FilterConfig& config
 ) {
     internal::log_python(" - Calculating node penalty scores...");
-    internal::get_penalty(
+    auto result = internal::get_penalty(
         kmers, nodes, n_nodes, record_offsets, n_record_offsets, is_targets, n_assemblies, config.n_cpu
     );
-    auto result = calculate_thresholds(
-        nodes, n_nodes, is_targets, n_assemblies, jaccard, jaccard_rows, jaccard_cols, config
+    calculate_thresholds(
+        is_targets, n_assemblies, jaccard, jaccard_rows, jaccard_cols, config, result
     );
 
     internal::log_python(" - Filtering graph edges and nodes...");
@@ -200,8 +179,7 @@ FilterResult filter(
     );
 
     auto [subgraphs, used_nodes] = internal::get_subgraphs(
-        pruned.nodes, pruned.edges, result.penalty_th, result.min_nodes,
-        result.max_nodes.value_or(std::numeric_limits<std::size_t>::max())
+        pruned.nodes, pruned.edges, result.penalty_th, result.min_nodes, result.max_nodes
     );
     if (subgraphs.empty()) {
         throw std::runtime_error("No low-penalty subgraph was found. Try decrease --stringency, or increase --penalty-th");