Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
30 changes: 22 additions & 8 deletions cpp/include/seqwin/fasta_reader.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -6,16 +6,30 @@

namespace seqwin {

struct FastaRecord
{
std::string id;
std::string sequence;
struct FastaRecord {
std::string id;
std::string sequence;
};

std::vector<FastaRecord>
read_fasta(const std::string& assembly_path);
/**
* @brief Read records from a plain-text or gzip-compressed FASTA file.
*
* @param assembly_path Path to a FASTA or `.gz` FASTA file.
* @return Vector of parsed FASTA records.
* @throws `std::runtime_error` If the file cannot be opened or the FASTA format is invalid.
*/
std::vector<FastaRecord> read_fasta(const std::string& assembly_path);

std::size_t
est_kmer_number(const std::vector<std::string>& assembly_paths, std::size_t windowsize);
/**
* @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<std::string>& assembly_paths,
std::size_t windowsize
);

} // namespace seqwin
52 changes: 51 additions & 1 deletion cpp/include/seqwin/graph.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -9,35 +9,85 @@

namespace seqwin {

/**
* @brief Location metadata for a minimizer.
*/
struct Kmer {
/** 0-based position of the minimizer within its FASTA record. */
std::uint32_t pos;
/** 0-based index of the FASTA record within the assembly. */
std::uint16_t record_idx;
/** Assembly index assigned to the source assembly. */
std::uint16_t assembly_idx;
};

/**
* @brief Minimizer graph node for one unique minimizer hash.
*
* The `[start, stop)` range identifies minimizers with this hash.
* Before merging, it selects entries in `ThreadGraph.idx`, whose values are indices into `ThreadGraph.kmers`.
* After merging, it indexes directly into the parallel `Graph.kmers` and `Graph.idx` arrays.
*/
struct Node {
/** Hash value of the minimizer represented by this node. */
std::uint64_t hash;
/** Start of the half-open range for this node's minimizer entries. */
std::uint64_t start;
/** End of the half-open range for this node's minimizer entries. */
std::uint64_t stop;
/** Number of target assemblies containing this minimizer hash. */
std::uint32_t n_tar;
/** Number of non-target assemblies containing this minimizer hash. */
std::uint32_t n_neg;
double penalty; // Used as thread_id before merging from different threads
/**
* Node penalty score used for downstream graph filtering (set to 0.0).
* Temporarily stores thread ID before thread-local graphs are merged.
*/
double penalty;
};

/**
* @brief Undirected weighted edge between two minimizers.
*/
struct Edge {
/** Smaller endpoint hash of the undirected edge. */
std::uint64_t first;
/** Larger endpoint hash of the undirected edge. */
std::uint64_t second;
/** Number of assemblies where the endpoints are adjacent. */
std::uint64_t weight;
};

/**
* @brief Container for the minimizer graph returned by `build()`.
*/
struct Graph {
/** Minimizer occurrences in all assemblies, grouped and sorted by hash. */
NoInitArray<Kmer> kmers;
/** Parallel to `kmers` and stores each minimizer's original generation index, ordered by genomic position. */
NoInitArray<std::uint64_t> idx;
/** Sorted by hash. */
NoInitArray<Node> nodes;
/** Sorted by hash. */
NoInitArray<Edge> edges;
/** FASTA record IDs of each assembly. */
std::vector<std::vector<std::string>> ids_by_assembly;
};

/**
* @brief Build a minimizer graph from assembly FASTA files.
*
* `assembly_paths`, `assembly_indices`, 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 assembly_indices Assembly indices, expected to be `0..N-1`.
* @param is_targets Whether each assembly is a target assembly.
* @param n_cpu Number of worker threads to use.
* @return Minimizer graph.
* @throws `std::runtime_error` If input sizes are inconsistent or indices exceed supported ranges.
*/
Graph build(
const std::vector<std::string>& assembly_paths,
std::size_t kmerlen,
Expand Down
38 changes: 32 additions & 6 deletions cpp/include/seqwin/helpers.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -9,28 +9,54 @@

namespace seqwin {

/**
* ThreadPool.parallel_for(n_items, fn) splits [0, n_items) into contiguous
* chunks and calls fn(start, end, worker_id) once per non-empty chunk.
* Callers should let parallel_for() choose chunk boundaries.
*/
class ThreadPool;

/**
* @brief Partial minimizer graph built by one worker thread.
*/
struct ThreadGraph {
std::vector<Kmer> kmers; // Needs push_back() support (exact kmer count is not known)
/**
* Minimizer occurrences generated by this worker, stored in genomic position order.
* Uses `std::vector` for `push_back()` support (exact kmer count is not known).
*/
std::vector<Kmer> kmers;
/**
* Indices into `kmers`, grouped by minimizer hash.
* The order is stable within each hash group, preserving the original genomic
* position order of minimizers with the same hash.
*/
NoInitArray<std::uint64_t> idx;
/** Unsorted. */
NoInitArray<Node> nodes;
/** Unsorted. */
NoInitArray<Edge> edges;
/** FASTA record IDs of each assembly. */
std::vector<std::vector<std::string>> ids_by_assembly;
/** Total number of minimizers generated by this worker. */
std::size_t n_kmers = 0;
/** Global index of the first assembly assigned to this worker. */
std::size_t start_assembly = 0;
};

/**
* @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.
*
* @param graphs Thread-local graphs produced during parallel graph construction.
* @param n_assemblies Total number of input assemblies.
* @param pool Thread pool used for parallel merge steps.
* @return Fully merged minimizer graph.
*/
Graph merge_thread_graphs(
std::vector<ThreadGraph>& graphs,
std::size_t n_assemblies,
Expand Down
2 changes: 1 addition & 1 deletion cpp/include/seqwin/no_init_array.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,7 @@
namespace seqwin {

/**
* A fixed-size owning array.
* @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,
Expand Down
13 changes: 13 additions & 0 deletions cpp/include/seqwin/thread_pool.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -13,6 +13,9 @@

namespace seqwin {

/**
* @brief Simple fixed-size worker pool for parallel range processing.
*/
class ThreadPool {
public:
explicit ThreadPool(std::size_t n_workers)
Expand Down Expand Up @@ -45,6 +48,16 @@ class ThreadPool {

std::size_t size() const noexcept { return n_workers_; }

/**
* @brief Run a function over contiguous chunks of `[0, n_items)`.
*
* The callable receives `(start, end, worker_id)`. Any exception thrown by a
* worker is captured and rethrown on the calling thread.
*
* @tparam Fn Callable type.
* @param n_items Number of items in the range.
* @param fn Function invoked once for each non-empty chunk.
*/
template <typename Fn>
void parallel_for(std::size_t n_items, Fn&& fn)
{
Expand Down
14 changes: 7 additions & 7 deletions cpp/src/bindings/python_bindings.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -71,22 +71,22 @@ PYBIND11_MODULE(_core, m) {
owner->edges.data(),
capsule
);
py::list all_idx_to_id;
py::list ids_by_assembly;
for (const auto& ids : owner->ids_by_assembly) {
py::tuple idx_to_id(ids.size());
py::tuple ids_tuple(ids.size());
for (std::size_t i = 0; i < ids.size(); ++i) {
idx_to_id[i] = ids[i];
ids_tuple[i] = ids[i];
}
all_idx_to_id.append(idx_to_id);
ids_by_assembly.append(ids_tuple);
}

return py::make_tuple(kmers, idx, nodes, edges, all_idx_to_id);
return py::make_tuple(kmers, idx, nodes, edges, ids_by_assembly);
},
py::arg("assembly_paths"),
py::arg("kmerlen"),
py::arg("windowsize"),
py::arg("assembly_idx"),
py::arg("is_target"),
py::arg("assembly_indices"),
py::arg("is_targets"),
py::arg("n_cpu") = 1
);

Expand Down
Loading
Loading