From 0722fc560fed4045c822e02473e2bfe2b70cc87e Mon Sep 17 00:00:00 2001 From: Sebastian Visan Date: Wed, 22 Jul 2026 12:33:55 +0300 Subject: [PATCH 01/12] =?UTF-8?q?=E2=99=BB=EF=B8=8F=20=20implicit=20upper?= =?UTF-8?q?=20bound=20=E2=9A=97=EF=B8=8F=20=20segmentation=20statistics?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- src/mincard.cpp | 38 +++++++++++++++++++++++++++++++++++++- 1 file changed, 37 insertions(+), 1 deletion(-) diff --git a/src/mincard.cpp b/src/mincard.cpp index dd04b85..3aeb65d 100644 --- a/src/mincard.cpp +++ b/src/mincard.cpp @@ -47,9 +47,15 @@ int main(int argc, char* argv[]) { seg_index U; CLI::Option *Uopt = app.add_option("-U,--max-segment-length", U, "Maximum segment length") - ->default_val(31) + ->default_val(0) ->expected(1, numeric_limits::max()); + bool min_size = false; + app.add_flag("--min-size", min_size, "Minimize the size of the segmentation instead of the cardinality"); + + bool stats = false; + app.add_flag("--stats", stats, "Calculate segmentation statistics"); + bool allow_perfect_segments = false; app.add_flag("-p,--perfect-segments", allow_perfect_segments, "In normal mode, additionally consider perfect segments of any length (recommended). With --trivial-vertical and --trivial-horizontal, use the maximal perfect segments and the trivial strategy in-between."); @@ -83,6 +89,22 @@ int main(int argc, char* argv[]) { } catch (const CLI::ParseError &e) { return app.exit(e); } + // No U value specified by the user + if (U == 0) { + if (min_size) { + if (gaps_as_symbols) { + // Implicit upper bound that keeps the segmentation optimal + U = L * 2 - 1; + } + else { + // Arbitrary upper bound so the algorithm is practical + U = L * 4 - 1; + } + } else { + // Default upper bound for min-card + U = 31; + } + } if (L > U) { cerr << "Upper and lower bounds are not compatible!" << endl; return 1; @@ -224,5 +246,19 @@ int main(int argc, char* argv[]) { cerr << " done: " << card << " cardinality, " << size << " gap-aware size" << ((verbose) ? " (" + to_string(duration.count()) + "ms)" : "") << endl; } + // Display segmentation statistics (min/max/avg segment length) + if (stats and segmentation.size() > 0) { + cerr << "Segmentation statistics: "; + int min_segment = c; + int max_segment = 0; + for (auto& segment: segmentation) { + int segment_size = segment.second - segment.first + 1; + min_segment = min(min_segment, segment_size); + max_segment = max(max_segment, segment_size); + } + double avg_segment = double(c) / double(segmentation.size()); + cerr << min_segment << " minimum length, " << max_segment << " maximum length, " << std::setprecision (2) << std::fixed << avg_segment << " average length" << endl; + } + return 0; } From 0d15af108684bf93581fc01db69259f87445a6dc Mon Sep 17 00:00:00 2001 From: Sebastian Visan Date: Thu, 23 Jul 2026 10:43:55 +0300 Subject: [PATCH 02/12] =?UTF-8?q?=E2=99=BB=EF=B8=8F=20=20nice=20pbwt=20api?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- Makefile | 4 +- src/algo.hpp | 7 +- src/mincard.cpp | 6 +- src/msa_chunker.hpp | 205 +++++++++++++++++++++++------------------ src/pbwt.cpp | 169 ++++++++++++++++++++++++++++++++++ src/pbwt.h | 74 +++++++++++++++ src/pbwt.hpp | 217 -------------------------------------------- 7 files changed, 373 insertions(+), 309 deletions(-) create mode 100644 src/pbwt.cpp create mode 100644 src/pbwt.h delete mode 100644 src/pbwt.hpp diff --git a/Makefile b/Makefile index 387af9f..ee9d5ab 100644 --- a/Makefile +++ b/Makefile @@ -9,8 +9,8 @@ HTSLIB_FLAGS=-L $(HTSLIB_LIB) -lhts -Wl,-rpath $(HTSLIB_LIB) OTHER_INCLUDE=ext/ SDSL_INCLUDE=ext/sdsl-lite/include/ -mincard: src/mincard.cpp src/segment.hpp src/RMaxQTree.h src/RMaxQTree.cpp src/msa_chunker.hpp src/trie.hpp src/pbwt.hpp src/algo.hpp - ${CXX} $(CXX_FLAGS) -DVERSION="\"$(VERSION)\"" -o mincard src/mincard.cpp src/RMaxQTree.cpp -I $(HTSLIB_INCLUDE) -I $(OTHER_INCLUDE) -I $(SDSL_INCLUDE) $(HTSLIB_FLAGS) +mincard: src/mincard.cpp src/segment.hpp src/RMaxQTree.h src/RMaxQTree.cpp src/msa_chunker.hpp src/trie.hpp src/pbwt.h src/pbwt.cpp src/algo.hpp + ${CXX} $(CXX_FLAGS) -DVERSION="\"$(VERSION)\"" -o mincard src/mincard.cpp src/RMaxQTree.cpp src/pbwt.cpp -I $(HTSLIB_INCLUDE) -I $(OTHER_INCLUDE) -I $(SDSL_INCLUDE) $(HTSLIB_FLAGS) clean: rm -f mincard diff --git a/src/algo.hpp b/src/algo.hpp index 38a967f..b530404 100644 --- a/src/algo.hpp +++ b/src/algo.hpp @@ -10,7 +10,7 @@ #include "trie.hpp" #include "RMaxQTree.h" // i_type #include "msa_chunker.hpp" -#include "pbwt.hpp" +#include "pbwt.h" using std::vector; using std::pair; @@ -34,9 +34,10 @@ namespace algo { const bool gaps_as_symbols, const bool use_pbwt ) { + static pbwt pb = pbwt(r); if (use_pbwt) { - pbwt::pbwt& PBWT = pbwt::pbwt::instance(r); - return PBWT.compute_meaningful_extensions(idx, r, c, L, U, y); + pb.next(idx.get_column(y - 1)); + return pb.meaningful_extensions(L, U); } vector> L_y; // 1-based indexing diff --git a/src/mincard.cpp b/src/mincard.cpp index 3aeb65d..fddea14 100644 --- a/src/mincard.cpp +++ b/src/mincard.cpp @@ -116,10 +116,14 @@ int main(int argc, char* argv[]) { std::unique_ptr storage; if (column_major) - storage = std::make_unique(inputfile, U); + storage = std::make_unique(inputfile, U, verbose); else storage = std::make_unique(inputfile, U, verbose); msa_chunker::msa_chunker& idx = *storage; + idx.set_row_major(true); + if (use_pbwt) { + idx.set_column_major(true); + } const int r = idx.get_rows(); const int c = idx.get_cols(); diff --git a/src/msa_chunker.hpp b/src/msa_chunker.hpp index 61c6823..278634f 100644 --- a/src/msa_chunker.hpp +++ b/src/msa_chunker.hpp @@ -34,12 +34,68 @@ namespace msa_chunker { * msa can be read into chunks either through a fasta file or a txt file (column major msa) */ class msa_chunker { + protected: + constexpr static msa_pos_t MIN_CHUNK_COLS = 131072; + msa_pos_t rows, cols, max_chunk_cols; + msa_pos_t chunk_start = std::numeric_limits::max(), chunk_cols = -1; + vector msa_rows; + vector msa_cols; + // Which chunks to store in memory + bool row_major = true; + bool column_major = true; + + virtual void load_chunk(msa_pos_t startcol, msa_pos_t length) = 0; + virtual string msa_substr_long(msa_pos_t row, msa_pos_t col, msa_pos_t length) = 0; public: virtual ~msa_chunker() = default; - virtual msa_pos_t get_rows() = 0; - virtual msa_pos_t get_cols() = 0; - virtual string msa_substr(msa_pos_t row, msa_pos_t col, msa_pos_t length) = 0; - virtual char msa_at(msa_pos_t row, msa_pos_t col) = 0; + msa_pos_t get_rows() { + return rows; + } + msa_pos_t get_cols() { + return cols; + } + void set_row_major(bool rm) { + row_major = rm; + msa_rows.clear(); + chunk_start = std::numeric_limits::max(); + chunk_cols = -1; + } + void set_column_major(bool cm) { + column_major = cm; + msa_cols.clear(); + chunk_start = std::numeric_limits::max(); + chunk_cols = -1; + } + // Return the sequence at MSA[row, col..col+length] + string msa_substr(msa_pos_t row, msa_pos_t col, msa_pos_t length) { + assert(0 <= row and row < rows and 0 <= col and col < cols and col + length <= cols); + assert (row_major); + if (length < max_chunk_cols) { + load_chunk(col, length); + return msa_rows[row].substr(col - chunk_start, length); + } else { + // Not in chunk memory + return msa_substr_long(row, col, length); + } + } + // Return the character at MSA[row, col] + char msa_at(msa_pos_t row, msa_pos_t col) { + assert(0 <= row and row < rows and 0 <= col and col < cols); + assert(column_major or row_major); + load_chunk(col, 1); + if(column_major){ + return msa_cols[col - chunk_start][row]; + } else { + return msa_rows[row][col - chunk_start]; + } + } + // Return a column from the MSA + string& get_column(msa_pos_t col) { + assert (0 <= col and col < cols); + assert (column_major); + load_chunk(col, 1); + return msa_cols[col - chunk_start]; + } }; /* @@ -49,12 +105,7 @@ namespace msa_chunker { */ class fasta_chunker : public msa_chunker { private: - constexpr static msa_pos_t MIN_CHUNK_COLS = 131072; - faidx_t *idx = NULL; - msa_pos_t rows, cols, max_chunk_cols; - msa_pos_t chunk_start = std::numeric_limits::max(), chunk_cols = -1; - vector msa_chunk; /* * explicitly load chunk [startcol..startcol+length) into memory @@ -68,22 +119,48 @@ namespace msa_chunker { return; } - msa_chunk.clear(); + msa_rows.clear(); chunk_start = startcol; chunk_cols = min(max_chunk_cols, cols - chunk_start); #ifdef MSA_CHUNKER_DEBUG cerr << "DEBUG: loading chunk [" << chunk_start << ".." << chunk_start + chunk_cols - 1 << "] (0-based) (query was [" << startcol << ".." << startcol + length - 1 << "])" << endl; #endif - for (msa_pos_t r = 0; r < rows; ++r) { - msa_pos_t out_len; - char *s_str = faidx_fetch_seq64(idx, faidx_iseq(idx, r), chunk_start, chunk_start + chunk_cols - 1, &out_len); - assert(out_len == chunk_cols); - msa_chunk.push_back(string(s_str)); - free(s_str); + // Calculate row-major chunk + if (row_major) { + for (msa_pos_t r = 0; r < rows; ++r) { + msa_pos_t out_len; + char *s_str = faidx_fetch_seq64(idx, faidx_iseq(idx, r), chunk_start, chunk_start + chunk_cols - 1, &out_len); + assert(out_len == chunk_cols); + msa_rows.push_back(string(s_str)); + free(s_str); + } + } + + // Calculate column-major chunk + if (column_major) { + // Relies on msa_rows + assert(row_major); + msa_cols.resize(chunk_cols); + for (msa_pos_t j = 0; j < chunk_cols; j++) { + string column = ""; + for (msa_pos_t i = 0; i < rows; i++) { + column.push_back(msa_rows[i][j]); + } + msa_cols[j] = column; + } } } + string msa_substr_long(msa_pos_t row, msa_pos_t col, msa_pos_t length) override { + msa_pos_t out_len; + char *s_str = faidx_fetch_seq64(idx, faidx_iseq(idx, row), col, col + length - 1, &out_len); + assert(out_len == length); + string s(s_str); + free(s_str); + return s; + } + public: fasta_chunker() = delete; @@ -92,6 +169,9 @@ namespace msa_chunker { */ fasta_chunker(const string &fastapath, const msa_pos_t max_qlen, const bool verbose) { max_chunk_cols = max(max_qlen, MIN_CHUNK_COLS); + if(!std::filesystem::is_regular_file(fastapath)){ + throw runtime_error("ERROR: FASTA file could not be found"); + } path fastap(fastapath); path fastaindex(fastapath + ".fai"); path gzip_index(fastapath + ".gzi"); @@ -153,40 +233,6 @@ namespace msa_chunker { cols = c; } - msa_pos_t get_rows() override { - return rows; - } - - msa_pos_t get_cols() override { - return cols; - } - - string msa_substr(msa_pos_t row, msa_pos_t col, msa_pos_t length) override { - assert(0 <= row and row < rows and 0 <= col and col < cols and col + length <= cols); - if (length <= max_chunk_cols) { - load_chunk(col, length); - return msa_chunk[row].substr(col - chunk_start, length); - } else { - msa_chunk.clear(); - chunk_start = std::numeric_limits::max(); - chunk_cols = -1; - - msa_pos_t out_len; - char *s_str = faidx_fetch_seq64(idx, faidx_iseq(idx, row), col, col + length - 1, &out_len); - assert(out_len == length); - string s(s_str); - free(s_str); - return s; - } - } - - // Return the character at MSA[row, col] - char msa_at(msa_pos_t row, msa_pos_t col) override { - assert(0 <= row and row < rows and 0 <= col and col < cols); - load_chunk(col, 1); - return msa_chunk[row][col - chunk_start]; - } - ~fasta_chunker() { fai_destroy(idx); } @@ -197,19 +243,13 @@ namespace msa_chunker { */ class column_chunker : public msa_chunker { private: - constexpr static msa_pos_t MIN_CHUNK_COLS = 131072; - - msa_pos_t rows, cols, max_chunk_cols; - msa_pos_t chunk_start = std::numeric_limits::max(), chunk_cols = -1; - vector msa_transposed_chunk; - vector msa_chunk; // used for faster msa_substr std::streampos matrix_start; ifstream msa_file; /* * explicitly load chunk [startcol..startcol+length) into memory */ - void load_chunk(msa_pos_t startcol, msa_pos_t length, bool transpose) { + void load_chunk(msa_pos_t startcol, msa_pos_t length) { assert(length <= max_chunk_cols); if (chunk_start <= startcol and startcol + length <= chunk_start + chunk_cols) { return; @@ -220,24 +260,40 @@ namespace msa_chunker { chunk_cols = min(max_chunk_cols, cols - chunk_start); std::streampos col_pos = matrix_start + chunk_start * (rows + 1); msa_file.seekg(col_pos); - msa_transposed_chunk.resize(chunk_cols * (rows + 1)); - msa_file.read(msa_transposed_chunk.data(), chunk_cols * (rows + 1)); + vector msa_buffer(chunk_cols * (rows + 1)); + msa_file.read(msa_buffer.data(), chunk_cols * (rows + 1)); + + // Split the column-major chunk + if(column_major) { + msa_cols.resize(chunk_cols); + for (msa_pos_t j = 0; j < chunk_cols; j++) { + auto msa_col = msa_buffer.begin() + j * (rows + 1); + msa_cols[j] = string(msa_col, msa_col + rows); + } + } // Compute original row-major chunk - if(transpose || chunk_start == 0) { - msa_chunk.resize(rows * chunk_cols); + if(row_major) { + msa_rows.resize(rows); for (msa_pos_t i = 0; i < rows; i++) { + string row = ""; for (msa_pos_t j = 0; j < chunk_cols; j++) { - msa_chunk[i * chunk_cols + j] = msa_transposed_chunk[j * (rows + 1) + i]; + row.push_back(msa_buffer[j * (rows + 1) + i]); } + msa_rows[i] = row; } } } + // Not implemented + string msa_substr_long(msa_pos_t row, msa_pos_t col, msa_pos_t length) override { + throw runtime_error("ERROR: substring length exceeds maximum column chunk value"); + } + public: column_chunker() = delete; - column_chunker(const string &msapath, const msa_pos_t max_qlen) { + column_chunker(const string &msapath, const msa_pos_t max_qlen, const bool verbose) { max_chunk_cols = max(MIN_CHUNK_COLS, max_qlen); msa_file = ifstream(msapath, std::ios::binary); if (!msa_file) { @@ -253,29 +309,6 @@ namespace msa_chunker { matrix_start = msa_file.tellg(); } - msa_pos_t get_rows() override { - return rows; - } - - msa_pos_t get_cols() override { - return cols; - } - - // Return the sequence at MSA[row, col..col+length] - string msa_substr(msa_pos_t row, msa_pos_t col, msa_pos_t length) override { - assert(0 <= row and row < rows and 0 <= col and col < cols and col + length <= cols); - load_chunk(col, length, true); - auto msa_row = msa_chunk.begin() + row * chunk_cols; - return string(msa_row + col - chunk_start, msa_row + col - chunk_start + length); - } - - // Return the character at MSA[row, col] - char msa_at(msa_pos_t row, msa_pos_t col) override { - assert(0 <= row and row < rows and 0 <= col and col < cols); - load_chunk(col, 1, false); - return msa_transposed_chunk[(col - chunk_start) * (rows + 1) + row]; - } - ~column_chunker(){ msa_file.close(); } diff --git a/src/pbwt.cpp b/src/pbwt.cpp new file mode 100644 index 0000000..45dec8d --- /dev/null +++ b/src/pbwt.cpp @@ -0,0 +1,169 @@ +#include "pbwt.h" + +#include +#include +#include +#include "rmq.hpp" + +pbwt::pbwt(const long long r): + r(r), k(0), + ak(r + 1, 0), sk(r + 2, 0), tk(r + 2, 0), ek(r + 1, 0), + dy(0), a(r + 2, 0), e(r + 1, 0) +{ + cnt.resize(alphabet_size + 1); + prev.resize(alphabet_size); + // Initial sorted order 1, 2, 3.. + std::iota(ak.begin(), ak.end(), 0); +} + +// If the symbol doesn't exist in the alpabet, add it +unsigned int pbwt::symbol_number(const char symbol) { + auto s = static_cast(symbol); + + if (alphabet[s] == -1) { + alphabet[s] = alphabet_size++; + cnt.push_back(0); + prev.push_back(0); + } + + return alphabet[s]; +} + +long long pbwt::max_ek(long long j, long long i) { + if (j != i) { + ek[j] = std::max(ek[j], max_ek(ak[j], i)); + ak[j] = i + 1; + } + return ek[j]; +} + +void pbwt::next(const string &column) { + // Increment column + k++; + // Symbol frequency + cnt.assign(cnt.size(), 0); + for (long long i = 0; i < r; i++){ + auto c = column[i]; + auto s = symbol_number(c) + 1; + cnt[s]++; + } + int symbols = 0; + for (long long sym = 1; sym <= alphabet_size; sym++) + if(cnt[sym] > 0) + symbols++; + if (k == 1 or symbols > 1) { + // There is at least one new divergence => height change + // Recompute pBWT arrays + prev.assign(prev.size(), 0); + tk.assign(tk.size(), 0); + sk[++dy] = k; + // RMQ for max ek + std::optional> rm; + if constexpr (max_range == RMQ) { + rm.emplace(ek); + } + // Counting sort + for (unsigned int i = 1; i < alphabet_size; i++) + cnt[i] += cnt[i - 1]; + for (long long i = 1; i <= r; i++) { + unsigned int b = symbol_number(column[ak[i] - 1]); + cnt[b]++; + a[cnt[b]] = ak[i]; + if constexpr (max_range == RECURSIVE) { + ak[i] = i + 1; + } + if (prev[b] == 0) + e[cnt[b]] = dy; + else { + // Calculate the range maximum + if constexpr (max_range == NAIVE) { + // O(alphabet) amortized - very fast in practice + e[cnt[b]] = *std::max_element(ek.begin() + prev[b] + 1, ek.begin() + i + 1); + } else if constexpr (max_range == RECURSIVE) { + // O(log alphabet) - solution from the paper + e[cnt[b]] = max_ek(prev[b] + 1, i); + } else { + // RMQ O(1) - best complexity but longer construction + e[cnt[b]] = rm->query(prev[b] + 1, i); + } + } + prev[b] = i; + } + for (long long i = 1; i <= r; i++) { + ak[i] = a[i]; + } + // Calculate frequency array + for (long long i = 1; i <= r; i++) { + tk[e[i]]++; + } + // Shrink arrays sk and tk, array a acts as tmp + long long j = 1; + for (long long i = 1; i <= dy; i++) { + if (tk[i] != 0) { + a[i] = j; + sk[j] = sk[i]; + tk[j] = tk[i]; + j++; + } + } + dy = j - 1; + // Fix ek array + for (long long i = 1; i <= r; i++) { + ek[i] = a[e[i]]; + } + // Calculate heights of extensions + for(long long i = 1; i <= dy; i++){ + tk[dy - i] += tk[dy - i + 1]; + } + } + else { + // Update last element of sk for first row divergence + if (tk[dy] > 1){ + dy++; + } + sk[dy] = k; + ek[1] = dy; + tk[dy] = 1; + } +} + +vector> pbwt::meaningful_extensions( + const long long L, + const long long U +) { + if (k < L) { + // No extension possible + return {}; + } + vector> L_y; + for(long long i = 0; i < dy; i++){ + if (k - sk[dy - i] + 1 < L) { + if (k - sk[dy - i - 1] + 1 > L) { + // Introduce new left extension at k - L + L_y.emplace_back(k - L + 1, tk[dy - i]); + } + } else { + if (k - sk[dy - i] + 1 <= U) { + // Left extension is in range [k - U, k - L] + L_y.emplace_back(sk[dy - i], tk[dy - i]); + } + } + } + // Dummy element, extension at k - U - 1 + L_y.emplace_back(std::max((long long)0, k - U), -1); + return L_y; +} + +long long pbwt::get_column(){ + return k; +} +vector pbwt::get_ak(){ + return ak; +} +vector pbwt::get_dk(){ + vector d(r + 1, 0); + for(int i = 1; i <= r; i++) { + d[i] = sk[ek[i]]; + } + return d; +} diff --git a/src/pbwt.h b/src/pbwt.h new file mode 100644 index 0000000..e7551db --- /dev/null +++ b/src/pbwt.h @@ -0,0 +1,74 @@ +#ifndef PBWT_H +#define PBWT_H + +#include +#include +#include + +using std::vector; +using std::string; +using std::pair; +using std::array; +using std::cerr, std::endl; + +enum MaxRange { + NAIVE, + RECURSIVE, + RMQ +}; + +inline constexpr enum MaxRange max_range = NAIVE; + +class pbwt { + private: + // Number of rows + long long r; + // Current column + long long k = 0; + // Arrays for current column 1-based indexing + vector ak; + vector sk; + vector tk; + vector ek; + // Counting sort arrays + vector cnt; + vector prev; + // Extra space + long long dy; // size of sk + vector a; + vector e; + // Alphabet (may grow with new symbols) + array alphabet = [] { + array alph{}; + alph.fill(-1); + + alph['A'] = 0; + alph['C'] = 1; + alph['G'] = 2; + alph['T'] = 3; + alph['-'] = 4; + + return alph; + }(); + unsigned int alphabet_size = 5; + // Map char -> int + unsigned int symbol_number(const char symbol); + // Recursive function for range max + long long max_ek(long long j, long long i); + public: + // Constructor taking the number of rows + pbwt(const long long r); + ~pbwt() = default; + // 1 indexed algorithm components + long long get_column(); + vector get_ak(); + vector get_dk(); + // Column streaming algorithm + void next(const string &column); + vector> meaningful_extensions( + const long long L, + const long long U + ); +}; + +#endif diff --git a/src/pbwt.hpp b/src/pbwt.hpp deleted file mode 100644 index 6d90a8c..0000000 --- a/src/pbwt.hpp +++ /dev/null @@ -1,217 +0,0 @@ -#ifndef PBWT_HPP -#define PBWT_HPP - -#include -#include -#include -#include -#include -#include -#include -#include "segment.hpp" -#include "msa_chunker.hpp" -#include "rmq.hpp" - -using std::vector; -using std::array; -using std::cerr, std::endl; - -// #define PBWT_DEBUG - -/* - * data structure for computing meaningful left extensions in linear time - */ -namespace pbwt { - typedef msa_chunker::msa_chunker msa_t; - typedef segment::seg_index seg_index; - - enum MaxRange { - NAIVE, - RECURSIVE, - RMQ - }; - - // Range maximum algorithm - inline constexpr enum MaxRange max_range = NAIVE; - - class pbwt { - private: - // Arrays for current column 1-based indexing - vector ak; - vector sk; - vector tk; - vector ek; - // Counting sort arrays - vector cnt; - vector prev; - // Extra space - seg_index dy; // size of sk - vector a; - vector e; - // Alphabet - array alphabet = [] { - array alph{}; - alph.fill(-1); - - alph['A'] = 0; - alph['C'] = 1; - alph['G'] = 2; - alph['T'] = 3; - alph['-'] = 4; - - return alph; - }(); - unsigned int alphabet_size = 5; - // Singleton constructor - pbwt(const seg_index r): - ak(r + 1), sk(r + 2), tk(r + 2), ek(r + 1, 0), - dy(0), a(r + 2), e(r + 1) - { - cnt.resize(alphabet_size + 1); - prev.resize(alphabet_size); - // Initial sorted order 1, 2, 3.. - std::iota(ak.begin(), ak.end(), 0); - } - ~pbwt() = default; - // Add new symbols to alphabet - unsigned int symbol_number(const char symbol) { - auto s = static_cast(symbol); - - if (alphabet[s] == -1) { - alphabet[s] = alphabet_size++; - cnt.push_back(0); - prev.push_back(0); - } - - return alphabet[s]; - } - // Recursive function for range max - seg_index max_ek(seg_index j, seg_index i) { - if (j != i) { - ek[j] = std::max(ek[j], max_ek(ak[j], i)); - ak[j] = i + 1; - } - return ek[j]; - } - public: - // Singleton pattern - static pbwt& instance(const seg_index r) { - static pbwt inst(r); - return inst; - } - // Delete copy and move operations - pbwt(const pbwt&) = delete; - pbwt& operator=(const pbwt&) = delete; - pbwt(pbwt&&) = delete; - pbwt& operator=(pbwt&&) = delete; - // Column streaming algorithm - vector> compute_meaningful_extensions( - msa_t &idx, - const seg_index r, - const seg_index c, - const seg_index L, - const seg_index U, - const seg_index y - ) { - if (y < L) { - return {}; // No extension possible - } - // Symbol frequency - cnt.assign(cnt.size(), 0); - for (seg_index i = 1; i <= r; i++){ - auto c = idx.msa_at(i - 1, y - 1); - auto s = symbol_number(c) + 1; - cnt[s]++; - } - int symbols = 0; - for (seg_index sym = 1; sym <= alphabet_size; sym++) - if(cnt[sym] > 0) - symbols++; - if (y == 1 or symbols > 1) { - // There is at least one new divergence => height change - // Recompute pBWT arrays - prev.assign(prev.size(), 0); - tk.assign(tk.size(), 0); - sk[++dy] = y; - // RMQ for max ek - std::optional> rm; - if constexpr (max_range == RMQ) { - rm.emplace(ek); - } - // Counting sort - for (unsigned int i = 1; i < alphabet_size; i++) - cnt[i] += cnt[i - 1]; - for (seg_index i = 1; i <= r; i++) { - unsigned int b = symbol_number(idx.msa_at(ak[i] - 1, y - 1)); - cnt[b]++; - a[cnt[b]] = ak[i]; - if constexpr (max_range == RECURSIVE) { - ak[i] = i + 1; - } - if (prev[b] == 0) - e[cnt[b]] = dy; - else { - // Calculate the range maximum - if constexpr (max_range == NAIVE) { - // O(alphabet) amortized - very fast in practice - e[cnt[b]] = *std::max_element(ek.begin() + prev[b] + 1, ek.begin() + i + 1); - } else if constexpr (max_range == RECURSIVE) { - // O(log alphabet) - solution from the paper - e[cnt[b]] = max_ek(prev[b] + 1, i); - } else { - // RMQ O(1) - best complexity but longer construction - e[cnt[b]] = rm->query(prev[b] + 1, i); - } - } - prev[b] = i; - } - for (seg_index i = 1; i <= r; i++) { - ak[i] = a[i]; - } - // Calculate frequency array - for (seg_index i = 1; i <= r; i++) { - tk[e[i]]++; - } - // Shrink arrays sk and tk, array a acts as tmp - seg_index j = 1; - for (seg_index i = 1; i <= dy; i++) { - if (tk[i] != 0) { - a[i] = j; - sk[j] = sk[i]; - tk[j] = tk[i]; - j++; - } - } - dy = j - 1; - // Fix ek array - for (seg_index i = 1; i <= r; i++) { - ek[i] = a[e[i]]; - } - // Calculate heights of extensions - for(seg_index i = 1; i <= dy; i++){ - tk[dy - i] += tk[dy - i + 1]; - } - } - else { - // Update last element of sk for first row divergence - if (tk[dy] > 1){ - dy++; - } - sk[dy] = y; - ek[1] = dy; - tk[dy] = 1; - } - // Compute meaningful extensions - vector> L_y; - for(seg_index i = 0; i < dy; i++){ - if(y - sk[dy - i] + 1 >= L and y - sk[dy - i] + 1 <= U) - L_y.emplace_back(sk[dy - i], tk[dy - i]); - } - // Dummy element - L_y.emplace_back(std::max((seg_index)0, y - U), -1); - return L_y; - } - }; -} - -#endif \ No newline at end of file From a762618e13375fb7333fdd5417105a6e938f3603 Mon Sep 17 00:00:00 2001 From: Sebastian Visan Date: Mon, 27 Jul 2026 12:27:36 +0300 Subject: [PATCH 03/12] =?UTF-8?q?=E2=9C=A8=20ring=20buffer=20minsize?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- Makefile | 2 +- src/mincard.cpp | 17 +++++-- src/minsize.hpp | 133 ++++++++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 146 insertions(+), 6 deletions(-) create mode 100644 src/minsize.hpp diff --git a/Makefile b/Makefile index ee9d5ab..c7ae6a5 100644 --- a/Makefile +++ b/Makefile @@ -9,7 +9,7 @@ HTSLIB_FLAGS=-L $(HTSLIB_LIB) -lhts -Wl,-rpath $(HTSLIB_LIB) OTHER_INCLUDE=ext/ SDSL_INCLUDE=ext/sdsl-lite/include/ -mincard: src/mincard.cpp src/segment.hpp src/RMaxQTree.h src/RMaxQTree.cpp src/msa_chunker.hpp src/trie.hpp src/pbwt.h src/pbwt.cpp src/algo.hpp +mincard: src/mincard.cpp src/segment.hpp src/RMaxQTree.h src/RMaxQTree.cpp src/msa_chunker.hpp src/trie.hpp src/pbwt.h src/pbwt.cpp src/algo.hpp src/minsize.hpp ${CXX} $(CXX_FLAGS) -DVERSION="\"$(VERSION)\"" -o mincard src/mincard.cpp src/RMaxQTree.cpp src/pbwt.cpp -I $(HTSLIB_INCLUDE) -I $(OTHER_INCLUDE) -I $(SDSL_INCLUDE) $(HTSLIB_FLAGS) clean: diff --git a/src/mincard.cpp b/src/mincard.cpp index fddea14..58c5d43 100644 --- a/src/mincard.cpp +++ b/src/mincard.cpp @@ -7,6 +7,7 @@ #include #include "algo.hpp" +#include "minsize.hpp" #include "segment.hpp" #include "msa_chunker.hpp" #include "CLI11.hpp" @@ -190,14 +191,20 @@ int main(int argc, char* argv[]) { cerr << " done" << ((verbose) ? " (" + to_string(duration.count()) + "ms)" : "") << endl; } - cerr << "Computing the minimum-cardinality segmentation..." << flush; auto start = high_resolution_clock::now(); - seg_index mincard; - tie(mincard, segmentation) = segment_with_rmq(idx, r, c, L, U, gaps_as_symbols, use_pbwt, L_y, perfect_columns); + seg_index minval; // cardinality or size + if (min_size) { + cerr << "Computing the minimum-size segmentation..." << flush; + minsize::minsize alg(idx, r, c, L, U, gaps_as_symbols, use_pbwt); + tie(minval, segmentation) = alg.segment(); + } else { + cerr << "Computing the minimum-cardinality segmentation..." << flush; + tie(minval, segmentation) = segment_with_rmq(idx, r, c, L, U, gaps_as_symbols, use_pbwt, L_y, perfect_columns); + } auto stop = high_resolution_clock::now(); auto duration = duration_cast(stop - start); - if (mincard != std::numeric_limits::max()) { - cerr << " done: " << segmentation.size() << " segments/ED words, " << mincard << " cardinality" << ((verbose) ? " (" + to_string(duration.count()) + "ms)" : "") << endl; + if (minval != std::numeric_limits::max()) { + cerr << " done: " << segmentation.size() << " segments/ED words, " << minval << (min_size ? " size" : " cardinality") << ((verbose) ? " (" + to_string(duration.count()) + "ms)" : "") << endl; } else { cerr << " done: no valid segmentation found!" << endl; return 1; diff --git a/src/minsize.hpp b/src/minsize.hpp new file mode 100644 index 0000000..30ad6ed --- /dev/null +++ b/src/minsize.hpp @@ -0,0 +1,133 @@ +#ifndef MINSIZE_HPP +#define MINSIZE_HPP + +#include "segment.hpp" +#include "msa_chunker.hpp" +#include "algo.hpp" + +using std::vector, std::tie; + +namespace minsize { + typedef msa_chunker::msa_chunker msa_t; + typedef segment::seg_index seg_index; + + // Class for solving the min-L-size segmentation problem + class minsize { + msa_t &idx; + const seg_index r; + const seg_index c; + const seg_index L; + const seg_index U; + const bool gaps_as_symbols; + const bool use_pbwt; + // Backtracking + vector back; + seg_index n_y; + // Min query + vector> m; + // Returns the index and size for the min size of M[height][left..right] + pair min_query(const seg_index height, const seg_index left, const seg_index right) { + const seg_index l = (max(left - 1, seg_index(0))) % U; + const seg_index r = (right - 1) % U; + seg_index minq = left == 0 ? 0 : numeric_limits::max(); + seg_index x = 0; + // Naive ring buffer + if (l <= r) { + for (seg_index j = l; j <= r; j++) { + if (minq > m[height - 1][j]) { + minq = m[height - 1][j]; + x = left + j - l; + } + } + } else { + for (seg_index j = l; j < U; j++) { + if (minq > m[height - 1][j]) { + minq = m[height - 1][j]; + x = left + j - l; + } + } + for (seg_index j = 0; j <= r; j++) { + if (minq > m[height - 1][j]) { + minq = m[height - 1][j]; + x = left + j + U - l - 1; + } + } + } + return {x, minq}; + } + // Updates the min query data structure for row h, column y + void update(const seg_index h, const seg_index y, const seg_index value) { + // Naive ring buffer + seg_index j = (y - 1) % U; + m[h - 1][j] = value; + } + // Calculate the min-L-size for the next column (y is 1-index) given meaningful left extensions L_yy + void next(const seg_index y, const vector> *L_yy) { + n_y = numeric_limits::max(); + for (size_t j = 0; j + 1 < L_yy->size(); j++) { + seg_index height = (*L_yy)[j].second; + seg_index left = (*L_yy)[j + 1].first; + seg_index right = (*L_yy)[j].first - 1; + seg_index x, s; + tie(x, s) = min_query(height, left, right); + // if (y < 15) + // cerr << "y = " << y << " h = " << height << " l = " << left << " r = " << right << " x = " << x << " s = " << s << " n_y = " << n_y << " small = " << height * y + s << endl; + if (n_y > height * y + s && height * y + s > 0) { + n_y = height * y + s; + back[y] = x; + } + } + for (seg_index height = 1; height <= r; height++) { + update(height, y, n_y - height * y); + } + } + public: + minsize( + msa_t &idx, + const seg_index r, + const seg_index c, + const seg_index L, + const seg_index U, + const bool gaps_as_symbols, + const bool use_pbwt + ): idx(idx), r(r), c(c), L(L), U(U), gaps_as_symbols(gaps_as_symbols), use_pbwt(use_pbwt), + n_y(numeric_limits::max()), back(c + 1, 0), + m(r, vector(U, 0)) {} + /* find the minimum-size segmentation of MSA[1..r][1..c] (indexed by + * idx) respecting lower bound L, upper bound U, gaps_as_symbols strategy, + * optionally given the meaningful left extensions L_y + * returns: pair {size, segmentation}, with size = + * numeric_limits::max() if no segmentation exists */ + pair>> segment( + const vector>> &L_y = {} + ) { + // Algorithm for gaps as symbols + // TODO: Add trie algo for gaps + for (i_type y = 1; y <= c; y++) { + // compute L_y if it was not given in input + const vector> *L_yy; + vector> L_yy_on_the_fly; + if (L_y.size() > 0) { + L_yy = &(*(L_y.begin() + y)); + } else { + L_yy_on_the_fly = algo::compute_meaningful_extensions(idx, r, c, L, U, y, gaps_as_symbols, use_pbwt); + L_yy = &L_yy_on_the_fly; + } + // Algorithm for column y + next(y, L_yy); + } + + // trace back + vector> segments; + for (i_type pos = c; pos > 0; pos = back[pos]) { + segments.emplace_back(back[pos] + 1, pos); + } + reverse(segments.begin(), segments.end()); + + return {n_y, segments}; + } + + }; +} + +#endif From 42c15b4ba69d44157bc015a9f833140d70adaae2 Mon Sep 17 00:00:00 2001 From: Sebastian Visan Date: Tue, 28 Jul 2026 10:51:29 +0300 Subject: [PATCH 04/12] =?UTF-8?q?=E2=9A=97=EF=B8=8F=20=20rmaxqtree=20min?= =?UTF-8?q?=20query?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- src/RMaxQTree.cpp | 18 +++--- src/minsize.hpp | 137 +++++++++++++++++++++++++++++++++++----------- 2 files changed, 114 insertions(+), 41 deletions(-) diff --git a/src/RMaxQTree.cpp b/src/RMaxQTree.cpp index 2aa65d0..78b1a3b 100644 --- a/src/RMaxQTree.cpp +++ b/src/RMaxQTree.cpp @@ -24,10 +24,10 @@ void RMaxQTree::init(i_type node, i_type b, i_type e, i_type *keys) { init(2 * node, b, (b + e) / 2, keys); init(2 * node + 1, (b + e) / 2 + 1, e, keys); // propagate up - if (tree[2 * node].Cj >= tree[2 * node + 1].Cj) - tree[node] = tree[2 * node]; - else + if (tree[2 * node + 1].Cj >= tree[2 * node].Cj) tree[node] = tree[2 * node + 1]; + else + tree[node] = tree[2 * node]; } } @@ -49,10 +49,10 @@ void RMaxQTree::updateTree(i_type key, i_type j, i_type Cj, i_type node, i_type else updateTree(key, j, Cj, 2 * node + 1, mid + 1, e); // And propagate back up - if (tree[2 * node].Cj > tree[2 * node + 1].Cj) // NR: propagate right value - tree[node] = tree[2 * node]; - else + if (tree[2 * node + 1].Cj > tree[2 * node].Cj) // NR: propagate left value tree[node] = tree[2 * node + 1]; + else + tree[node] = tree[2 * node]; } } @@ -73,10 +73,10 @@ std::pair RMaxQTree::queryTree(i_type i, i_type j, i_type node, i return right; if (right.second == negative_infinity) return left; - if (left.second <= right.second) // NR: prefer higher key - return right; - else + if (right.second <= left.second) // NR: prefer lower key return left; + else + return right; } // Empty constructor for creating arrays diff --git a/src/minsize.hpp b/src/minsize.hpp index 30ad6ed..0b68d47 100644 --- a/src/minsize.hpp +++ b/src/minsize.hpp @@ -4,6 +4,9 @@ #include "segment.hpp" #include "msa_chunker.hpp" #include "algo.hpp" +#include "RMaxQTree.h" + +#include using std::vector, std::tie; @@ -11,6 +14,15 @@ namespace minsize { typedef msa_chunker::msa_chunker msa_t; typedef segment::seg_index seg_index; + // Which data structure to use for the sliding window minimum range queries + enum MinRange { + RING, // Naive ring buffer + TREE, // RMaxQTree + QUEUE // RMQueue + }; + + inline constexpr enum MinRange min_range = TREE; + // Class for solving the min-L-size segmentation problem class minsize { msa_t &idx; @@ -20,65 +32,109 @@ namespace minsize { const seg_index U; const bool gaps_as_symbols; const bool use_pbwt; - // Backtracking + // n[y] is min size ending at column y + vector n; + // Traceback vector back; - seg_index n_y; // Min query - vector> m; + vector> ring; + vector rmqtree; + vector keys; // Returns the index and size for the min size of M[height][left..right] pair min_query(const seg_index height, const seg_index left, const seg_index right) { const seg_index l = (max(left - 1, seg_index(0))) % U; - const seg_index r = (right - 1) % U; + const seg_index r = (max(right - 1, seg_index(0))) % U; seg_index minq = left == 0 ? 0 : numeric_limits::max(); seg_index x = 0; // Naive ring buffer - if (l <= r) { - for (seg_index j = l; j <= r; j++) { - if (minq > m[height - 1][j]) { - minq = m[height - 1][j]; - x = left + j - l; + if constexpr (min_range == RING) { + if (l <= r) { + for (seg_index j = l; j <= r; j++) { + if (minq > ring[height - 1][j]) { + minq = ring[height - 1][j]; + x = left + j - l; + } + } + } else { + for (seg_index j = l; j < U; j++) { + if (minq > ring[height - 1][j]) { + minq = ring[height - 1][j]; + x = left + j - l; + } + } + for (seg_index j = 0; j <= r; j++) { + if (minq > ring[height - 1][j]) { + minq = ring[height - 1][j]; + x = left + j + U - l - 1; + } } } - } else { - for (seg_index j = l; j < U; j++) { - if (minq > m[height - 1][j]) { - minq = m[height - 1][j]; - x = left + j - l; + } + if constexpr (min_range == TREE) { + if (l <= r) { + auto [x1, minq1] = rmqtree[height - 1].query(l + 1, r + 1); + minq1 *= -1; + x1--; + if (minq1 < minq) { + minq = minq1; + x = left + x1 - l; + } + } else { + auto [x1, minq1] = rmqtree[height - 1].query(l + 1, U); + minq1 *= -1; + x1--; + auto [x2, minq2] = rmqtree[height - 1].query(1, r + 1); + minq2 *= -1; + x2--; + if (minq1 < minq) { + minq = minq1; + x = left + x1 - l; } - } - for (seg_index j = 0; j <= r; j++) { - if (minq > m[height - 1][j]) { - minq = m[height - 1][j]; - x = left + j + U - l - 1; + if (minq2 < minq) { + minq = minq2; + x = left + x2 + U - l - 1; } } } + if constexpr (min_range == QUEUE) { + cerr << "Not yet implemented"; + } return {x, minq}; } // Updates the min query data structure for row h, column y void update(const seg_index h, const seg_index y, const seg_index value) { // Naive ring buffer - seg_index j = (y - 1) % U; - m[h - 1][j] = value; + if constexpr (min_range == RING) { + seg_index j = (y - 1) % U; + ring[h - 1][j] = value; + } + if constexpr (min_range == TREE) { + seg_index j = (y - 1) % U; + rmqtree[h - 1].update(j + 1, j + 1, -value); + } + if constexpr (min_range == QUEUE) { + cerr << "Not yet implemented"; + } } // Calculate the min-L-size for the next column (y is 1-index) given meaningful left extensions L_yy void next(const seg_index y, const vector> *L_yy) { - n_y = numeric_limits::max(); for (size_t j = 0; j + 1 < L_yy->size(); j++) { seg_index height = (*L_yy)[j].second; seg_index left = (*L_yy)[j + 1].first; seg_index right = (*L_yy)[j].first - 1; - seg_index x, s; - tie(x, s) = min_query(height, left, right); - // if (y < 15) - // cerr << "y = " << y << " h = " << height << " l = " << left << " r = " << right << " x = " << x << " s = " << s << " n_y = " << n_y << " small = " << height * y + s << endl; - if (n_y > height * y + s && height * y + s > 0) { - n_y = height * y + s; + seg_index x, s, x2, s2; + + tie(x, s) = (left == right) + ? pair{ left, n[left] - height * left } + : min_query(height, left, right); + + if (n[y] > height * y + s && height * y + s > 0) { + n[y] = height * y + s; back[y] = x; } } for (seg_index height = 1; height <= r; height++) { - update(height, y, n_y - height * y); + update(height, y, n[y] - height * y); } } public: @@ -91,8 +147,25 @@ namespace minsize { const bool gaps_as_symbols, const bool use_pbwt ): idx(idx), r(r), c(c), L(L), U(U), gaps_as_symbols(gaps_as_symbols), use_pbwt(use_pbwt), - n_y(numeric_limits::max()), back(c + 1, 0), - m(r, vector(U, 0)) {} + n(c + 1, numeric_limits::max()), back(c + 1, 0) { + n[0] = 0; + if constexpr (min_range == RING) { + ring = vector>(r, vector(U, 0)); + } + if constexpr (min_range == TREE) { + rmqtree.resize(r); + keys.resize(U + 1); + for (i_type i = 0; i <= U; ++i) + keys[i] = i; + for(auto &rmq: rmqtree) { + rmq.fillRMaxQTree(keys.data(), U + 1); + rmq.update(0, 0, 0); + } + } + if constexpr (min_range == QUEUE) { + cerr << "Not yet implemented"; + } + } /* find the minimum-size segmentation of MSA[1..r][1..c] (indexed by * idx) respecting lower bound L, upper bound U, gaps_as_symbols strategy, * optionally given the meaningful left extensions L_y @@ -124,7 +197,7 @@ namespace minsize { } reverse(segments.begin(), segments.end()); - return {n_y, segments}; + return {n[c], segments}; } }; From f0ea5cecd1f02c2e2597f4350bf2542b3e1d97fa Mon Sep 17 00:00:00 2001 From: Sebastian Visan Date: Fri, 31 Jul 2026 17:08:54 +0300 Subject: [PATCH 05/12] =?UTF-8?q?=E2=9C=A8=20range=20minimum=20queue?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- .gitignore | 1 + Makefile | 7 +- src/minsize.hpp | 63 ++++--- src/rmqueue.cpp | 382 +++++++++++++++++++++++++++++++++++++++++++ src/rmqueue.h | 67 ++++++++ src/test_rmqueue.cpp | 107 ++++++++++++ 6 files changed, 598 insertions(+), 29 deletions(-) create mode 100644 src/rmqueue.cpp create mode 100644 src/rmqueue.h create mode 100644 src/test_rmqueue.cpp diff --git a/.gitignore b/.gitignore index 305efb3..3006dce 100644 --- a/.gitignore +++ b/.gitignore @@ -1,4 +1,5 @@ mincard +test_rmqueue test/*.eds test/*.gfa test/*.fai diff --git a/Makefile b/Makefile index c7ae6a5..96ca1b8 100644 --- a/Makefile +++ b/Makefile @@ -9,8 +9,11 @@ HTSLIB_FLAGS=-L $(HTSLIB_LIB) -lhts -Wl,-rpath $(HTSLIB_LIB) OTHER_INCLUDE=ext/ SDSL_INCLUDE=ext/sdsl-lite/include/ -mincard: src/mincard.cpp src/segment.hpp src/RMaxQTree.h src/RMaxQTree.cpp src/msa_chunker.hpp src/trie.hpp src/pbwt.h src/pbwt.cpp src/algo.hpp src/minsize.hpp - ${CXX} $(CXX_FLAGS) -DVERSION="\"$(VERSION)\"" -o mincard src/mincard.cpp src/RMaxQTree.cpp src/pbwt.cpp -I $(HTSLIB_INCLUDE) -I $(OTHER_INCLUDE) -I $(SDSL_INCLUDE) $(HTSLIB_FLAGS) +mincard: src/mincard.cpp src/segment.hpp src/rmqueue.h src/rmqueue.cpp src/RMaxQTree.h src/RMaxQTree.cpp src/msa_chunker.hpp src/trie.hpp src/pbwt.h src/pbwt.cpp src/algo.hpp src/minsize.hpp + ${CXX} $(CXX_FLAGS) -DVERSION="\"$(VERSION)\"" -o mincard src/mincard.cpp src/rmqueue.cpp src/RMaxQTree.cpp src/pbwt.cpp -I $(HTSLIB_INCLUDE) -I $(OTHER_INCLUDE) -I $(SDSL_INCLUDE) $(HTSLIB_FLAGS) clean: rm -f mincard + +test_rmqueue: src/test_rmqueue.cpp src/rmqueue.cpp + ${CXX} $(CXX_FLAGS) -o test_rmqueue src/test_rmqueue.cpp src/rmqueue.cpp diff --git a/src/minsize.hpp b/src/minsize.hpp index 0b68d47..a932e45 100644 --- a/src/minsize.hpp +++ b/src/minsize.hpp @@ -5,6 +5,7 @@ #include "msa_chunker.hpp" #include "algo.hpp" #include "RMaxQTree.h" +#include "rmqueue.h" #include @@ -21,7 +22,7 @@ namespace minsize { QUEUE // RMQueue }; - inline constexpr enum MinRange min_range = TREE; + inline constexpr enum MinRange min_range = QUEUE; // Class for solving the min-L-size segmentation problem class minsize { @@ -37,42 +38,43 @@ namespace minsize { // Traceback vector back; // Min query - vector> ring; - vector rmqtree; + vector> rings; + vector rmqtrees; + vector rmqueues; vector keys; // Returns the index and size for the min size of M[height][left..right] - pair min_query(const seg_index height, const seg_index left, const seg_index right) { + pair min_query(const seg_index height, const seg_index y, const seg_index left, const seg_index right) { const seg_index l = (max(left - 1, seg_index(0))) % U; const seg_index r = (max(right - 1, seg_index(0))) % U; seg_index minq = left == 0 ? 0 : numeric_limits::max(); seg_index x = 0; // Naive ring buffer - if constexpr (min_range == RING) { + if constexpr (min_range == RING) { if (l <= r) { for (seg_index j = l; j <= r; j++) { - if (minq > ring[height - 1][j]) { - minq = ring[height - 1][j]; + if (minq > rings[height - 1][j]) { + minq = rings[height - 1][j]; x = left + j - l; } } } else { for (seg_index j = l; j < U; j++) { - if (minq > ring[height - 1][j]) { - minq = ring[height - 1][j]; + if (minq > rings[height - 1][j]) { + minq = rings[height - 1][j]; x = left + j - l; } } for (seg_index j = 0; j <= r; j++) { - if (minq > ring[height - 1][j]) { - minq = ring[height - 1][j]; - x = left + j + U - l - 1; + if (minq > rings[height - 1][j]) { + minq = rings[height - 1][j]; + x = left + j + U - l; } } } } if constexpr (min_range == TREE) { if (l <= r) { - auto [x1, minq1] = rmqtree[height - 1].query(l + 1, r + 1); + auto [x1, minq1] = rmqtrees[height - 1].query(l + 1, r + 1); minq1 *= -1; x1--; if (minq1 < minq) { @@ -80,10 +82,10 @@ namespace minsize { x = left + x1 - l; } } else { - auto [x1, minq1] = rmqtree[height - 1].query(l + 1, U); + auto [x1, minq1] = rmqtrees[height - 1].query(l + 1, U); minq1 *= -1; x1--; - auto [x2, minq2] = rmqtree[height - 1].query(1, r + 1); + auto [x2, minq2] = rmqtrees[height - 1].query(1, r + 1); minq2 *= -1; x2--; if (minq1 < minq) { @@ -92,12 +94,15 @@ namespace minsize { } if (minq2 < minq) { minq = minq2; - x = left + x2 + U - l - 1; + x = left + x2 + U - l; } } } if constexpr (min_range == QUEUE) { - cerr << "Not yet implemented"; + long long queue_col = std::max(0LL, y - U); + long long x_q = rmqueues[height - 1].query(left - queue_col, right - queue_col); + minq = rmqueues[height - 1].get(x_q); + x = x_q + queue_col; } return {x, minq}; } @@ -106,14 +111,16 @@ namespace minsize { // Naive ring buffer if constexpr (min_range == RING) { seg_index j = (y - 1) % U; - ring[h - 1][j] = value; + rings[h - 1][j] = value; } if constexpr (min_range == TREE) { seg_index j = (y - 1) % U; - rmqtree[h - 1].update(j + 1, j + 1, -value); + rmqtrees[h - 1].update(j + 1, j + 1, -value); } if constexpr (min_range == QUEUE) { - cerr << "Not yet implemented"; + rmqueues[h - 1].push(value); + if (y >= U) + rmqueues[h - 1].pop(); } } // Calculate the min-L-size for the next column (y is 1-index) given meaningful left extensions L_yy @@ -122,12 +129,11 @@ namespace minsize { seg_index height = (*L_yy)[j].second; seg_index left = (*L_yy)[j + 1].first; seg_index right = (*L_yy)[j].first - 1; - seg_index x, s, x2, s2; + seg_index x, s; tie(x, s) = (left == right) ? pair{ left, n[left] - height * left } - : min_query(height, left, right); - + : min_query(height, y, left, right); if (n[y] > height * y + s && height * y + s > 0) { n[y] = height * y + s; back[y] = x; @@ -150,20 +156,23 @@ namespace minsize { n(c + 1, numeric_limits::max()), back(c + 1, 0) { n[0] = 0; if constexpr (min_range == RING) { - ring = vector>(r, vector(U, 0)); + rings = vector>(r, vector(U, 0)); } if constexpr (min_range == TREE) { - rmqtree.resize(r); + rmqtrees.resize(r); keys.resize(U + 1); for (i_type i = 0; i <= U; ++i) keys[i] = i; - for(auto &rmq: rmqtree) { + for(auto &rmq: rmqtrees) { rmq.fillRMaxQTree(keys.data(), U + 1); rmq.update(0, 0, 0); } } if constexpr (min_range == QUEUE) { - cerr << "Not yet implemented"; + rmqueues = vector(r, RMQueue(U + 1)); + for(auto &rmq: rmqueues) { + rmq.push(0); + } } } /* find the minimum-size segmentation of MSA[1..r][1..c] (indexed by diff --git a/src/rmqueue.cpp b/src/rmqueue.cpp new file mode 100644 index 0000000..112f1f0 --- /dev/null +++ b/src/rmqueue.cpp @@ -0,0 +1,382 @@ +#include "rmqueue.h" + +#include +#include +#include + +using std::cout; + +// Simple function for floor of base 2 logarithm +long long log2(long long x) { + int y = 0; + while (x >>= 1) { + y++; + } + return y; +} + +// Ceil division operation +long long ceil_div(long long x, long long y) { + if (x % y == 0) { + return x / y; + } else { + return x / y + 1; + } +} + +// API + +RMQueue::RMQueue(long long n): n(n), logn(log2(n)) { + b = ceil_div(logn, 4); + blocks = ceil_div(n, b) + 1; + q_size = b * blocks; + Q = vector(q_size, std::numeric_limits::max()); + // Initialize logarithm array + lsize = 1 << ceil_div(logn, 2); // O(sqrt(n)) + L.resize(lsize + 1); + L[0] = 0; + long long i = 1, lg = 0, count = 1; + while(i <= lsize) { + L[i++] = lg; + count--; + if(count == 0){ + lg++; + count = 1 << lg; + } + } + // Block aligned queries + p = vector>(logn + 1, vector(blocks + 1, -1)); + // Precompute tables + compute_tj(); + compute_lca(); + // In block queries + cart = vector(blocks, 0); + spine = vector(b + 1, 0); +} + +void RMQueue::pop() { + start++; + if (start >= b) { + start = 0; + begin = plus(begin, 1); + } +} + +void RMQueue::push(long long x) { + setq(end * b + stop, x); + // Update C array + update_cartesian_tree(x); + stop++; + if (stop >= b) { + stop = 0; + // Update arrays with new block + update_pk(end); + // Next block + end = plus(end, 1); + // Clear cartesian tree at that block + cart[end] = 0; + spine_i = 0; + } +} + +// Get element at index i in Q +long long RMQueue::get(long long i) { + if (i < 0 || i > n) { + return -1; + } + return getq(i + start + begin * b); +} + +long long RMQueue::query(long long left, long long right) { + if (left < 0 || right > n || left > right) { + return -1; + } + long long i = left + start; + long long j = right + start; + // Inclusive with first/last block element + long long b1 = i / b; // block of i + long long b2 = j / b; // block of j + long long offset = begin * b + start; + long long x = -1; + if (b1 == b2) { + x = in_block_query(b1, i % b, j % b); + } else { + long long q1 = in_block_query(b1, i % b, b - 1); + long long q2 = block_aligned_query(b1 + 1, b2 - 1); + long long q3 = in_block_query(b2, 0, j % b); + x = minx(minx(q1, q2), q3); + } + return (x + q_size - offset) % q_size; +} + +// Util methods + +long long RMQueue::block_aligned_query(long long b1, long long b2) { + // Check if query is empty + if (b2 < b1) + return -1; + // Block query [b1, b2] + long long k = log(b2 - b1); + long long x1 = p[k][plus(b2, begin)]; + long long x2 = p[k][plus(b1 - 1 + (1 << k), begin)]; + return minx(x2, x1); +} + +long long RMQueue::in_block_query(long long block, long long i1, long long i2) { + if (i1 == i2) + return plus(block, begin) * b + i1; + if (plus(block, begin) != end) { + // Complete block + return plus(block, begin) * b + lca[cart[plus(block, begin)]][1 + i1][1 + i2] - 1; + } + // Last block is partial => cartesian tree with #stop number of vertices + // Get a valid cartesian tree that starts the same + long long t = (cart[plus(block, begin)] << ((b - stop) * 2)) + (1 << (b - stop)) - 1; + return plus(block, begin) * b + lca[t][1 + i1][1 + i2] - 1; +} + +long long RMQueue::getq(long long i) { + return Q[i % q_size]; +} + +void RMQueue::setq(long long i, long long x) { + Q[i % q_size] = x; +} + +long long RMQueue::minx(long long x1, long long x2) { + if (x1 == -1) + return x2; + if (x2 == -1) + return x1; + if(getq(x1) <= getq(x2)) { + return x1; + } + return x2; +} + +long long RMQueue::plus(long long x, long long y) { + return (x + y + blocks * 10) % blocks; +} + +long long RMQueue::log(long long x) { + if(x <= lsize) + return L[x]; + return L[x / lsize] + ceil_div(logn, 2); +} + +void RMQueue::update_pk(long long block) { + long long min_b = std::numeric_limits::max(); + long long l = -1; // index + for (int i = 0; i < b; i++) { + long long val = getq(block * b + i); + if (min_b > val) { + min_b = val; + l = block * b + i; + } + p[0][block] = l; + for (int k = 1; k <= logn; k++) { + long long l1 = p[k - 1][block]; + long long l2 = p[k - 1][plus(block, - (1 << (k - 1)))]; + p[k][block] = minx(l2, l1); + } + } +} + +void RMQueue::compute_tj() { + tj.resize(b); // [1,b-1] + tj[0] = vector>(1, vector(1, 1)); + vector right_spine(b * 2, -1); + right_spine[0] = 0; + for (int j = 1; j < b; j++) { + tj[j].resize(1 << (2 * j)); + // All cartesian trees x with j vertices + // Balanced parentheses 0 = (, 1 = ) + for (int x = 0; x < (1 << (2 * j)); x++) { + long long i = 0; // right_spine index + long long open = 0; + long long v = 0; // tree node + bool ok = true; + for (int bit = 2 * j - 1; bit >= 0; bit--) { + if (open == -1) { + ok = false; + break; + } + if (open == 0) { + // Start new spine + i = 0; + } + if ((x >> bit) & 1) { + open--; + } else { + open++; + // Add v to spine + right_spine[++i] = ++v; + } + } + if (ok && open == 0) { + // Integer x describes a cartesian tree + tj[j][x] = vector(j + 1, -1LL); + long long new_x = (x << 2) + 3; // add 2 ) at the end + for (int i2 = 0; i2 <= i; i2 ++) { + // set ( at position i2 from the end + tj[j][x][right_spine[i2]] = new_x & ~(1 << (i2 + 1)); + } + } + } + } +} + +int calculate_level(vector& parent, int node) { + int level = 0; + while(parent[node] != 0) { + level++; + node = parent[node]; + } + return level; +} + +void RMQueue::compute_lca() { + cout << "Computing LCA b = " << b << "\n"; + lca.resize(1 << (2 * b)); + vector parent(b + 1, -1); + vector level(b + 1, -1); + vector stack_trace(b + 1, -1); + // All cartesian trees with vertices + for (int x = 0; x < (1 << (2 * b)); x++) { + long long open = 0; + long long i = -1; // index for stack-trace + long long v = 0; // current number of ( + long long last = -1; // last closed node + bool ok = true; + for (int bit = 2 * b - 1; bit >= 0; bit--) { + if ((x >> bit) & 1) { + open--; + } else { + open++; + } + if (open < 0 || open > b) { + ok = false; + break; + } + if ((x >> bit) & 1) { + // Close last node from stack trace + if (last != -1) { + parent[last] = stack_trace[i]; + } + last = stack_trace[i]; + i--; + } else { + // Add next v to trace + stack_trace[++i] = ++v; + if (last != -1) { + parent[last] = v; + last = -1; + } + } + } + if (ok && open == 0) { + // Integer x describes a cartesian tree + lca[x] = vector(b + 1, vector(b + 1, -1LL)); + // Simply calculate lca of each 2 nodes + parent[last] = 0; // last node closed is root + for (int i = 1; i <= b; i++) { + level[i] = calculate_level(parent, i); + } + for (int p = 1; p <= b; p++) { + for (int q = 1; q <= b; q++) { + // LCA for p & q + int parent_p = p; + int parent_q = q; + while (parent_p != parent_q) { + if (level[parent_p] > level[parent_q]) + parent_p = parent[parent_p]; + else + parent_q = parent[parent_q]; + } + lca[x][p][q] = parent_p; + } + } + } + } +} + +void RMQueue::update_cartesian_tree(long long x) { + // Find where to insert x on the spine + while(spine_i > 0 && getq(end * b + spine[spine_i] - 1) > x) { + spine_i--; + } + // Update the cartesian tree by inserting the new node + cart[end] = tj[stop][cart[end]][spine[spine_i]]; + // Place the new node on the right spine + spine_i++; + spine[spine_i] = stop + 1; +} + +// Debugging + +void print_n_bits(long long x, int n) { + for (int i = n - 1; i >= 0; i--) { + cout << ((x >> i) & 1LL); + } +} + +void RMQueue::debug_tables() { + cout << "Debugging Tj table:\n"; + for (int j = 1; j < b; j++) { + cout << "Tj[" << j << "]:\n"; + for (int x = 0; x < (1 << (2 * j)); x++) { + if(!tj[j][x].empty()) { + cout << "x = "; + print_n_bits(x, j * 2); + cout << "\n"; + for (int h = 0; h <= j; h++) { + if (tj[j][x][h] != -1) { + cout << " " << h << " -> "; + print_n_bits(tj[j][x][h], j * 2 + 2); + cout << "\n"; + } + } + } + } + } + cout << "Debugging LCA table:\n"; + for (int x = 0; x < (1 << (2 * b)); x++) { + if(!lca[x].empty()) { + cout << "LCA for x = "; + print_n_bits(x, b * 2); + cout << "\n"; + for (int p = 1; p <= b; p++) { + for (int q = p; q <= b; q++) { + cout << p << " " << q << " -> " << lca[x][p][q] << "\n"; + } + } + } + } +} + +void RMQueue::debug_trees() { + cout << "Debugging cartesian trees:\n"; + if (end >= begin) { + // No ring buffer wrapping + for (int i = begin; i < end; i++) { + cout << "c[" << i << "] = "; + print_n_bits(cart[i], 2 * b); + cout << "\n"; + } + } else { + for (int i = begin; i < blocks; i++) { + cout << "c[" << i << "] = "; + print_n_bits(cart[i], 2 * b); + cout << "\n"; + } + for (int i = 0; i < end; i++) { + cout << "c[" << i << "] = "; + print_n_bits(cart[i], 2 * b); + cout << "\n"; + } + } + // Current block + cout << "current cartesian tree "; + print_n_bits(cart[end], 2 * stop); + cout << "\n"; +} \ No newline at end of file diff --git a/src/rmqueue.h b/src/rmqueue.h new file mode 100644 index 0000000..0b9825c --- /dev/null +++ b/src/rmqueue.h @@ -0,0 +1,67 @@ +#ifndef RMQUEUE_H +#define RMQUEUE_H + +#include + +using std::vector; + +/* Range Minimum Queue + * mantains a queue of at most N integers and supports operations: + * pop in O(1) + * - deletes the first element of the queue + * push (element) in amortized O(1) + * - inserts the element at the back of the queue + * query [left, right] in O(1) + * - returns the minimum in the specified interval + */ +class RMQueue { + private: + // Ring buffer for the actual queue elements + vector Q; + long long getq(long long i); + void setq(long long i, long long x); + long long minx(long long x1, long long x2); // index of minimum from 2 ind + // Invariants + long long n; + long long logn; + long long b; // block length + long long blocks; + long long q_size; + // Current element + long long begin = 0; // block + long long start = 0; // shift + long long end = 0; // block + long long stop = 0; // in-block shift + // Modulo addition + long long plus(long long x, long long y); + // Precomputed logarithm array + long long lsize; + vector L; + long long log(long long x); + // Queries + long long block_aligned_query(long long b1, long long b2); + long long in_block_query(long long block, long long i1, long long i2); + // Block-aligned + // p[k][h] is the index of the minimum element in Q[(h-2^k)b+1, hb] + vector> p; + void update_pk(long long block); + // In-block + vector>> tj; + void compute_tj(); + vector>> lca; + void compute_lca(); + vector cart; + vector spine; + long long spine_i = 0; + void update_cartesian_tree(long long x); + public: + RMQueue(long long n); + void pop(); + void push(long long x); + long long query(long long left, long long right); + long long get(long long i); + void debug_tables(); + void debug_trees(); +}; + +#endif \ No newline at end of file diff --git a/src/test_rmqueue.cpp b/src/test_rmqueue.cpp new file mode 100644 index 0000000..3a74607 --- /dev/null +++ b/src/test_rmqueue.cpp @@ -0,0 +1,107 @@ +#include "rmqueue.h" +#include +#include + +using std::cout; + +void test_b_1() { + // Small queue + RMQueue rmq(7); + rmq.push(100); + rmq.push(12); + rmq.push(50); + rmq.push(75); + rmq.push(7); + // rmq can be indexed + assert(rmq.get(0) == 100); + assert(rmq.get(1) == 12); + assert(rmq.get(4) == 7); + // test queries + assert(rmq.query(0, 4) == 4); + assert(rmq.query(0, 3) == 1); + assert(rmq.query(2, 2) == 2); + assert(rmq.query(3, 4) == 4); + // remove element 100 + rmq.pop(); + assert(rmq.get(0) == 12); + assert(rmq.query(0, 3) == 3); + assert(rmq.query(0, 2) == 0); + assert(rmq.query(1, 2) == 1); + // add more elements + rmq.push(21); + rmq.push(5); + rmq.pop(); + rmq.push(2014); + rmq.pop(); + rmq.pop(); + rmq.push(-3); + rmq.push(9); + rmq.pop(); + rmq.push(90); + rmq.push(6); + // queue 21 5 2014 -3 9 90 6 + assert(rmq.get(0) == 21); + assert(rmq.get(1) == 5); + assert(rmq.get(6) == 6); + assert(rmq.query(0, 6) == 3); + assert(rmq.query(0, 2) == 1); + assert(rmq.query(4, 5) == 4); +} + +void test_b_3() { + RMQueue rmq(1<<12); + // rmq.debug_tables(); + + for (int i = 0; i < 600; i++) + rmq.push(i); + + assert(rmq.query(1, 32) == 1); + assert(rmq.query(20, 400) == 20); + assert(rmq.query(3, 4) == 3); + + for (int i = 0; i < 300; i++) + rmq.pop(); + + assert(rmq.get(0) == 300); + assert(rmq.get(rmq.query(2, 5)) == 302); + + for (int i = 600; i < (1<<12); i++) + rmq.push(i); + + for (int i = 0; i < 200; i++) + rmq.pop(); + + for (int i = 1<<12; i < (1<<12) + 250; i++) + rmq.push(i); + + assert(rmq.query(0, 12) == 0); + assert(rmq.query(25, 26) == 25); + assert(rmq.get(rmq.query(25, 26)) == 525); + + for (int i = 0; i < 3500 ; i++) + rmq.pop(); + + assert(rmq.query(14, 122) == 14); + assert(rmq.get(rmq.query(14, 122)) == 4014); +} + +void test_equal_mins() { + RMQueue rmq(10); + rmq.push(622); + rmq.push(620); + rmq.push(620); + rmq.push(620); + rmq.push(621); + rmq.push(622); + rmq.push(620); + assert(rmq.query(0, 0) == 0); + assert(rmq.query(0, 6) == 1); // first minimum index +} + +int main() { + cout << "Basic RMQueue Testing\n"; + test_b_1(); + test_b_3(); + test_equal_mins(); + cout << "All tests passed!\n"; +} \ No newline at end of file From 4e235f64980f375d64a92d8cd9c4897a4cfc5faf Mon Sep 17 00:00:00 2001 From: Sebastian Visan Date: Mon, 3 Aug 2026 11:03:49 +0300 Subject: [PATCH 06/12] =?UTF-8?q?=E2=9A=97=EF=B8=8F=20=20minsize=20experim?= =?UTF-8?q?ent?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- .gitignore | 7 ++++ experiments/minsize/README.md | 18 ++++++++ experiments/minsize/input/.gitkeep | 0 experiments/minsize/run_experiment.sh | 60 +++++++++++++++++++++++++++ src/rmqueue.cpp | 1 - 5 files changed, 85 insertions(+), 1 deletion(-) create mode 100644 experiments/minsize/README.md create mode 100644 experiments/minsize/input/.gitkeep create mode 100755 experiments/minsize/run_experiment.sh diff --git a/.gitignore b/.gitignore index 3006dce..d0b7ae7 100644 --- a/.gitignore +++ b/.gitignore @@ -26,4 +26,11 @@ experiments/pbwt/mincard-naive experiments/pbwt/mincard-recursive experiments/pbwt/mincard-rmq +experiments/minsize/input/* +!experiments/minsize/input/.gitkeep +experiments/minsize/output +experiments/minsize/mincard-ring +experiments/minsize/mincard-tree +experiments/minsize/mincard-queue + experiments/transpose/output \ No newline at end of file diff --git a/experiments/minsize/README.md b/experiments/minsize/README.md new file mode 100644 index 0000000..a56e429 --- /dev/null +++ b/experiments/minsize/README.md @@ -0,0 +1,18 @@ +# min range experiment + +Simple minsize experiment to compare the different ways of calculating the range minimum query for the M matrix rows. + +Get covid dataset as described [here](../covid/README.md) + +Compile mincard with each algorithm (ring buffer, rmaxtree, rmqueue) by changing [this line](../../src/minsize.hpp#L25) and then running + +```sh + make -C ../../ + mv ../../mincard ./mincard-ring +``` + +Do this with the names [ring, tree, queue]. Finally run the experiment script (running times will be collected in output/results.csv) + +```sh + ./run_experiment.sh +``` diff --git a/experiments/minsize/input/.gitkeep b/experiments/minsize/input/.gitkeep new file mode 100644 index 0000000..e69de29 diff --git a/experiments/minsize/run_experiment.sh b/experiments/minsize/run_experiment.sh new file mode 100755 index 0000000..486f398 --- /dev/null +++ b/experiments/minsize/run_experiment.sh @@ -0,0 +1,60 @@ +#!/bin/bash +set -euo pipefail +thisfolder=$( cd -- "$( dirname -- "${BASH_SOURCE[0]}" )" &> /dev/null && pwd ) # https://stackoverflow.com/questions/59895/how-do-i-get-the-directory-where-a-bash-script-is-located-from-within-the-script +cd "$thisfolder" + +mincard_ring="$thisfolder/mincard-ring" +mincard_tree="$thisfolder/mincard-tree" +mincard_queue="$thisfolder/mincard-queue" +inputmsa="$thisfolder/input/covid_100000.fa.gz" +usrbintimeformat="%e" + +outdir="output" +mkdir -p "$outdir" +rm -f "$outdir"/* +cd "$outdir" +ln -s "$inputmsa" msa.fa + +L_values=(1 2 4 8 16 32 64 128 256 512) +printf "Running minsize with pBWT and different range min strategies (ring, rmaxtree, rmqueue) on different L values %s \n" "${L_values[*]}" + +tmpfile=$(mktemp) +trap 'rm -f "$tmpfile"' EXIT + +RED='\033[0;31m' +GREEN='\033[0;32m' +BLUE='\033[0;34m' +PURPLE='\033[0;35m' +NC='\033[0m' + +# Write results to table +echo -e "L value,Ring Buffer,RMaxTree,RMQueue" > results.csv + +# minsize with pbwt (gaps as symbols strategy) +for L in "${L_values[@]}" +do + printf "\n${BLUE}Comparison for L = $L${NC}\n" + # ring buffer + /usr/bin/time -f"$usrbintimeformat" -o "$tmpfile" "$mincard_ring" msa.fa -v --gaps-as-symbols --pbwt -L $L --min-size -o minsize_ring_L${L}.eds + t1=$(<"$tmpfile") + # rmaxtree + /usr/bin/time -f"$usrbintimeformat" -o "$tmpfile" "$mincard_tree" msa.fa -v --gaps-as-symbols --pbwt -L $L --min-size -o minsize_tree_L${L}.eds + t2=$(<"$tmpfile") + # rmqueue + /usr/bin/time -f"$usrbintimeformat" -o "$tmpfile" "$mincard_queue" msa.fa -v --gaps-as-symbols --pbwt -L $L --min-size -o minsize_queue_L${L}.eds + t3=$(<"$tmpfile") + + if cmp -s minsize_ring_L${L}.eds minsize_tree_L${L}.eds && + cmp -s minsize_ring_L${L}.eds minsize_queue_L${L}.eds + then + printf "${GREEN}Outputs are identical for L = ${L}!${NC}\n" + printf "${PURPLE}ring buffer %.3fs vs. rmaxtree suffix trie %.3fs vs. rmqueue %.3fs${NC}\n" "$t1" "$t2" "$t3" + else + printf "${RED}Outputs are different for ${name} on U = ${U}!${NC}\n" + exit 1 + fi + + echo -e "${L},${t1},${t2},${t3}" >> results.csv +done + +printf "\n${GREEN}Finished running the experiment!\n${BLUE}Wrote algorithm times to results.csv${NC}\n" diff --git a/src/rmqueue.cpp b/src/rmqueue.cpp index 112f1f0..1eb545d 100644 --- a/src/rmqueue.cpp +++ b/src/rmqueue.cpp @@ -236,7 +236,6 @@ int calculate_level(vector& parent, int node) { } void RMQueue::compute_lca() { - cout << "Computing LCA b = " << b << "\n"; lca.resize(1 << (2 * b)); vector parent(b + 1, -1); vector level(b + 1, -1); From 18d89c039dce9728e9160d2ef9e386f31e35fbb3 Mon Sep 17 00:00:00 2001 From: Sebastian Visan Date: Mon, 3 Aug 2026 12:00:56 +0300 Subject: [PATCH 07/12] =?UTF-8?q?=E2=9A=A1=EF=B8=8F=20semi-dynamic=20rmq?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- src/algo.hpp | 40 +++++++++++++++++++++++----------------- 1 file changed, 23 insertions(+), 17 deletions(-) diff --git a/src/algo.hpp b/src/algo.hpp index b530404..7cb0497 100644 --- a/src/algo.hpp +++ b/src/algo.hpp @@ -8,7 +8,7 @@ #include "segment.hpp" #include "trie.hpp" -#include "RMaxQTree.h" // i_type +#include "rmqueue.h" #include "msa_chunker.hpp" #include "pbwt.h" @@ -253,20 +253,16 @@ namespace algo { seg_index perfect_first = -1, perfect_m = numeric_limits::max(); const bool allow_perfect_segments = (perfect_columns.size() > 0); - RMaxQTree rmq; // TODO use range min query data structure - vector keys(c + 1); // keys is passed to rmq - for (i_type i = 0; i <= c; ++i) - keys[i] = i; - rmq.fillRMaxQTree(keys.data(), c + 1); - + RMQueue rmq(U + 1); m[0] = 0; - rmq.update(0, 0, 0); + rmq.push(0); + if (allow_perfect_segments) { perfect_m = numeric_limits::max(); perfect_first = -1; } - for (i_type y = 1; y <= c; ++y) { + for (seg_index y = 1; y <= c; ++y) { m[y] = numeric_limits::max(); // compute L_y if it was not given in input @@ -284,14 +280,22 @@ namespace algo { // optimal solution using L_y for (size_t j = 0; j + 1 < L_yy->size(); ++j) { - i_type l = (*L_yy)[j + 1].first; - i_type r = (*L_yy)[j].first - 1; + seg_index l = (*L_yy)[j + 1].first; + seg_index r = (*L_yy)[j].first - 1; if (l > r) continue; - auto [x, neg_mx] = rmq.query(l, r); - if (x == -1) continue; - if (neg_mx == -numeric_limits::max()) continue; - i_type candidate = (*L_yy)[j].second + m[x]; + seg_index x = l; + if (L > 1) { + // m_y is not monotone non-decreasing + seg_index queue_col = std::max(0LL, y - U); + x = rmq.query(l - queue_col, r - queue_col); + auto mx = rmq.get(x); + + if (x == -1) continue; + if (mx == numeric_limits::max()) continue; + x += queue_col; + } + seg_index candidate = (*L_yy)[j].second + m[x]; assert(candidate >= 0); if (candidate < m[y]) { @@ -308,7 +312,9 @@ namespace algo { } } - rmq.update(y, y, -m[y]); + rmq.push(m[y]); + if (y >= U) + rmq.pop(); // update perfect-segment run if (allow_perfect_segments) { @@ -324,7 +330,7 @@ namespace algo { // trace back vector> segments; - for (i_type pos = c; pos > 0; pos = back[pos]) { + for (seg_index pos = c; pos > 0; pos = back[pos]) { segments.emplace_back(back[pos] + 1, pos); } reverse(segments.begin(), segments.end()); From 0e84ea4082ad3e379ee0aaff85b445f9dd18abb8 Mon Sep 17 00:00:00 2001 From: Sebastian Visan Date: Thu, 6 Aug 2026 12:22:53 +0300 Subject: [PATCH 08/12] =?UTF-8?q?=E2=99=BB=EF=B8=8F=20=20quiet=20flag?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- src/mincard.cpp | 118 ++++++++++++++++++++++++++++++++++---------- src/minsize.hpp | 2 +- src/msa_chunker.hpp | 41 +++++++++------ 3 files changed, 118 insertions(+), 43 deletions(-) diff --git a/src/mincard.cpp b/src/mincard.cpp index 58c5d43..09f9a90 100644 --- a/src/mincard.cpp +++ b/src/mincard.cpp @@ -23,8 +23,6 @@ using algo::seg_index, algo::compute_perfect_columns, algo::compute_all_meaningf typedef segment::seg_size_t seg_size_t; using segment::segment_stream_gfa, segment::segment_stream_eds, segment::segment_stream_no_output; -bool verbose = false; - int main(int argc, char* argv[]) { CLI::App app{"mincard version " + string(VERSION) + " — build Elastic Degenerate Strings (EDSes) from multiple sequence alignments in FASTA format"}; argv = app.ensure_utf8(argv); @@ -79,6 +77,9 @@ int main(int argc, char* argv[]) { bool verbose = false; app.add_flag("-v,--verbose", verbose, "Print running times"); + bool quiet = false; + app.add_flag("-q,--quiet", quiet, "Print only cardinality and size"); + bool use_pbwt = false; app.add_flag("--pbwt", use_pbwt, "Compute meaningful extensions using positional Burrows-Wheeler Transform"); @@ -114,12 +115,22 @@ int main(int argc, char* argv[]) { cerr << "pBWT only works with the gaps as symbols strategy! Add flag --gaps-as-symbols" << endl; return 1; } + if (verbose and quiet) { + cerr << "Quitet flag is ignored if verbose is set" << endl; + } + int verbosity = 1; + if (quiet) { + verbosity = 0; + } + if (verbose) { + verbosity = 2; + } std::unique_ptr storage; if (column_major) - storage = std::make_unique(inputfile, U, verbose); + storage = std::make_unique(inputfile, U, verbosity); else - storage = std::make_unique(inputfile, U, verbose); + storage = std::make_unique(inputfile, U, verbosity); msa_chunker::msa_chunker& idx = *storage; idx.set_row_major(true); if (use_pbwt) { @@ -128,7 +139,9 @@ int main(int argc, char* argv[]) { const int r = idx.get_rows(); const int c = idx.get_cols(); - cerr << "Processing MSA[1.." << r << ",1.." << c << "] (\"" << inputfile << "\")" << endl; + if (verbosity > 0) { + cerr << "Processing MSA[1.." << r << ",1.." << c << "] (\"" << inputfile << "\")" << endl; + } vector perfect_columns = {}; if (allow_perfect_segments) { @@ -137,12 +150,17 @@ int main(int argc, char* argv[]) { tie(p, perfect_columns) = compute_perfect_columns(idx, r, c); auto stop = high_resolution_clock::now(); auto duration = duration_cast(stop - start); - cerr << "MSA contains " << p << "/" << c << " (" << setprecision(4) << (double) 100 * p / c << "%) perfect columns" <<((verbose) ? " (" + to_string(duration.count()) + "ms)" : "") << endl; + if (verbosity > 0) { + cerr << "MSA contains " << p << "/" << c << " (" << setprecision(4) << (double) 100 * p / c << "%) perfect columns" + << ((verbosity > 1) ? " (" + to_string(duration.count()) + "ms)" : "") << endl; + } } vector> segmentation; // 1-based segments [x..y] if (trivial_segmentation) { - cerr << "Computing the S^¦¦¦ segmentation..." << flush; + if (verbosity > 0) { + cerr << "Computing the S^¦¦¦ segmentation..." << flush; + } auto start = high_resolution_clock::now(); segmentation.reserve(c); for (seg_index i = 1; i <= c; ++i) { @@ -157,9 +175,13 @@ int main(int argc, char* argv[]) { } auto stop = high_resolution_clock::now(); auto duration = duration_cast(stop - start); - cerr << " done: " << segmentation.size() << " segments/ED words" << ((verbose) ? " (" + to_string(duration.count()) + "ms)" : "") << endl; + if (verbosity > 0) { + cerr << " done: " << segmentation.size() << " segments/ED words" << ((verbosity > 1) ? " (" + to_string(duration.count()) + "ms)" : "") << endl; + } } else if (no_segmentation) { - cerr << "Computing the S^≡ segmentation..." << flush; + if (verbosity > 0) { + cerr << "Computing the S^≡ segmentation..." << flush; + } auto start = high_resolution_clock::now(); if (!allow_perfect_segments) { segmentation.push_back({ 1, c }); @@ -177,34 +199,48 @@ int main(int argc, char* argv[]) { } auto stop = high_resolution_clock::now(); auto duration = duration_cast(stop - start); - cerr << " done: " << segmentation.size() << " segments/ED words" << ((verbose) ? " (" + to_string(duration.count()) + "ms)" : "") << endl; + if (verbosity > 0) { + cerr << " done: " << segmentation.size() << " segments/ED words" << ((verbosity > 1) ? " (" + to_string(duration.count()) + "ms)" : "") << endl; + } } else { - cerr << "The allowed segments are" << ((allow_perfect_segments) ? " perfect segments and those" : "") << " of length [" << L << ".." << U << "]" << endl; - + if (verbosity > 0) { + cerr << "The allowed segments are" << ((allow_perfect_segments) ? " perfect segments and those" : "") << " of length [" << L << ".." << U << "]" << endl; + } vector>> L_y; if (preprocess) { - cerr << "Computing the meaningful extensions..." << flush; + if (verbosity > 0) { + cerr << "Computing the meaningful extensions..." << flush; + } auto start = high_resolution_clock::now(); L_y = compute_all_meaningful_extensions(idx, r, c, L, U, gaps_as_symbols, use_pbwt); auto stop = high_resolution_clock::now(); auto duration = duration_cast(stop - start); - cerr << " done" << ((verbose) ? " (" + to_string(duration.count()) + "ms)" : "") << endl; + if (verbosity > 0) { + cerr << " done" << ((verbosity > 1) ? " (" + to_string(duration.count()) + "ms)" : "") << endl; + } } auto start = high_resolution_clock::now(); seg_index minval; // cardinality or size if (min_size) { - cerr << "Computing the minimum-size segmentation..." << flush; + if (verbosity > 0) { + cerr << "Computing the minimum-size segmentation..." << flush; + } minsize::minsize alg(idx, r, c, L, U, gaps_as_symbols, use_pbwt); tie(minval, segmentation) = alg.segment(); } else { - cerr << "Computing the minimum-cardinality segmentation..." << flush; + if (verbosity > 0) { + cerr << "Computing the minimum-cardinality segmentation..." << flush; + } tie(minval, segmentation) = segment_with_rmq(idx, r, c, L, U, gaps_as_symbols, use_pbwt, L_y, perfect_columns); } auto stop = high_resolution_clock::now(); auto duration = duration_cast(stop - start); if (minval != std::numeric_limits::max()) { - cerr << " done: " << segmentation.size() << " segments/ED words, " << minval << (min_size ? " size" : " cardinality") << ((verbose) ? " (" + to_string(duration.count()) + "ms)" : "") << endl; + if (verbosity > 0) { + cerr << " done: " << segmentation.size() << " segments/ED words, " << minval << (min_size ? " size" : " cardinality") + << ((verbosity > 1) ? " (" + to_string(duration.count()) + "ms)" : "") << endl; + } } else { cerr << " done: no valid segmentation found!" << endl; return 1; @@ -216,10 +252,14 @@ int main(int argc, char* argv[]) { ostream *out; ofstream outfile; if (outputgfafile == "-") { - cerr << "Streaming the block graph to stdout..." << flush; + if (verbosity > 0) { + cerr << "Streaming the block graph to stdout..." << flush; + } out = &cout; } else { - cerr << "Streaming the block graph to \"" << outputgfafile << "\"..." << flush; + if (verbosity > 0) { + cerr << "Streaming the block graph to \"" << outputgfafile << "\"..." << flush; + } outfile = ofstream(outputgfafile); out = &outfile; } @@ -228,16 +268,24 @@ int main(int argc, char* argv[]) { auto stop = high_resolution_clock::now(); auto duration = duration_cast(stop - start); outfile.close(); - cerr << " done: " << card << " cardinality, " << size << " gap-aware size" << ((verbose) ? " (" + to_string(duration.count()) + "ms)" : "") << endl; + if (verbosity > 0) { + cerr << " done: " << card << " cardinality, " << size << " gap-aware size" << ((verbosity > 1) ? " (" + to_string(duration.count()) + "ms)" : "") << endl; + } else { + cerr << card << " " << size << endl; + } } if (outputedsfile != "") { ostream *out; ofstream outfile; if (outputedsfile == "-") { - cerr << "Streaming the EDS to stdout..." << flush; + if (verbosity > 0) { + cerr << "Streaming the EDS to stdout..." << flush; + } out = &cout; } else { - cerr << "Streaming the EDS to \"" << outputedsfile << "\"..." << flush; + if (verbosity > 0) { + cerr << "Streaming the EDS to \"" << outputedsfile << "\"..." << flush; + } outfile = ofstream(outputedsfile); out = &outfile; } @@ -246,20 +294,32 @@ int main(int argc, char* argv[]) { auto stop = high_resolution_clock::now(); auto duration = duration_cast(stop - start); outfile.close(); - cerr << " done: " << card << " cardinality, " << size << " gap-aware size" << ((verbose) ? " (" + to_string(duration.count()) + "ms)" : "") << endl; + if (verbosity > 0) { + cerr << " done: " << card << " cardinality, " << size << " gap-aware size" << ((verbosity > 1) ? " (" + to_string(duration.count()) + "ms)" : "") << endl; + } else { + cerr << card << " " << size << endl; + } } if (outputgfafile == "" and outputedsfile == "") { - cerr << "Computing the EDS stats (no output selected)..." << flush; + if (verbosity > 0) { + cerr << "Computing the EDS stats (no output selected)..." << flush; + } auto start = high_resolution_clock::now(); tie(card, size) = segment_stream_no_output(idx, r, c, segmentation); auto stop = high_resolution_clock::now(); auto duration = duration_cast(stop - start); - cerr << " done: " << card << " cardinality, " << size << " gap-aware size" << ((verbose) ? " (" + to_string(duration.count()) + "ms)" : "") << endl; + if (verbosity > 0) { + cerr << " done: " << card << " cardinality, " << size << " gap-aware size" << ((verbosity > 1) ? " (" + to_string(duration.count()) + "ms)" : "") << endl; + } else { + cerr << card << " " << size << endl; + } } // Display segmentation statistics (min/max/avg segment length) if (stats and segmentation.size() > 0) { - cerr << "Segmentation statistics: "; + if (verbosity > 0) { + cerr << "Segmentation statistics: "; + } int min_segment = c; int max_segment = 0; for (auto& segment: segmentation) { @@ -268,7 +328,11 @@ int main(int argc, char* argv[]) { max_segment = max(max_segment, segment_size); } double avg_segment = double(c) / double(segmentation.size()); - cerr << min_segment << " minimum length, " << max_segment << " maximum length, " << std::setprecision (2) << std::fixed << avg_segment << " average length" << endl; + if (verbosity > 0) { + cerr << min_segment << " minimum length, " << max_segment << " maximum length, " << std::setprecision (2) << std::fixed << avg_segment << " average length" << endl; + } else { + cerr << min_segment << " " << max_segment << " " << std::setprecision (2) << std::fixed << avg_segment << endl; + } } return 0; diff --git a/src/minsize.hpp b/src/minsize.hpp index a932e45..53d7164 100644 --- a/src/minsize.hpp +++ b/src/minsize.hpp @@ -22,7 +22,7 @@ namespace minsize { QUEUE // RMQueue }; - inline constexpr enum MinRange min_range = QUEUE; + inline constexpr enum MinRange min_range = RING; // Class for solving the min-L-size segmentation problem class minsize { diff --git a/src/msa_chunker.hpp b/src/msa_chunker.hpp index 278634f..e0cab9c 100644 --- a/src/msa_chunker.hpp +++ b/src/msa_chunker.hpp @@ -167,7 +167,7 @@ namespace msa_chunker { /* * index a given (gzipped) FASTA file */ - fasta_chunker(const string &fastapath, const msa_pos_t max_qlen, const bool verbose) { + fasta_chunker(const string &fastapath, const msa_pos_t max_qlen, const int verbosity) { max_chunk_cols = max(max_qlen, MIN_CHUNK_COLS); if(!std::filesystem::is_regular_file(fastapath)){ throw runtime_error("ERROR: FASTA file could not be found"); @@ -186,36 +186,45 @@ namespace msa_chunker { if (!exists(fastaindex) or (compressed and !exists(gzip_index))) { auto start = high_resolution_clock::now(); - cerr << "Index" << ((compressed) ? "es " : " ") << fastaindex << ((compressed) ? " or \"" + gzip_index.string() + "\"" : "") << " not found, generating the index" << ((compressed) ? "es" : "") << "..." << flush; + if (verbosity > 0) { + cerr << "Index" << ((compressed) ? "es " : " ") << fastaindex << ((compressed) ? " or \"" + gzip_index.string() + "\"" : "") + << " not found, generating the index" << ((compressed) ? "es" : "") << "..." << flush; + } if (fai_build3(fastap.c_str(), fastaindex.c_str(), gzip_index.c_str()) == -1) { - cerr << "\nERROR: failed to create index" << endl; - exit(1); + throw runtime_error("ERROR:failed to create index"); } auto stop = high_resolution_clock::now(); auto duration = duration_cast(stop - start); - cerr << " done" << ((verbose) ? " (" + to_string(duration.count()) + "ms)" : "") << endl; + if (verbosity > 0) { + cerr << " done" << ((verbosity > 1) ? " (" + to_string(duration.count()) + "ms)" : "") << endl; + } } else if (last_write_time(fastaindex) < last_write_time(fastap) or (compressed and last_write_time(gzip_index) < last_write_time(fastap))) { auto start = high_resolution_clock::now(); - cerr << "Index" << ((compressed) ? "es " : " ") << fastaindex << ((compressed) ? " or \"" + gzip_index.string() + "\" are" : " is") << " older than MSA, regenerating the index" << ((compressed) ? "es" : "") << "..." << flush; + if (verbosity > 0) { + cerr << "Index" << ((compressed) ? "es " : " ") << fastaindex << ((compressed) ? " or \"" + gzip_index.string() + "\" are" : " is") + << " older than MSA, regenerating the index" << ((compressed) ? "es" : "") << "..." << flush; + } remove(fastaindex); if (compressed and exists(gzip_index)) remove(gzip_index); if (fai_build3(fastap.c_str(), fastaindex.c_str(), gzip_index.c_str()) == -1) { - cerr << "\nERROR: failed to create index" << endl; - exit(1); + throw runtime_error("ERROR:failed to create index"); } auto stop = high_resolution_clock::now(); auto duration = duration_cast(stop - start); - cerr << " done" << ((verbose) ? " (" + to_string(duration.count()) + "ms)" : "") << endl; + if (verbosity > 0) { + cerr << " done" << ((verbosity > 1) ? " (" + to_string(duration.count()) + "ms)" : "") << endl; + } } else { - cerr << "Index" << ((compressed) ? "es " : " ") << fastaindex << ((compressed) ? " and \"" + gzip_index.string() + "\"" : "") << " found" << endl; + if (verbosity > 0) { + cerr << "Index" << ((compressed) ? "es " : " ") << fastaindex << ((compressed) ? " and \"" + gzip_index.string() + "\"" : "") << " found" << endl; + } } assert(exists(fastaindex) and (!compressed or exists(gzip_index))); if (!(idx = fai_load3(fastap.c_str(), fastaindex.c_str(), gzip_index.c_str(), FAI_NONE))) { - cerr << "\nERROR: failed to create index" << endl; - exit(1); + throw runtime_error("ERROR:failed to create index"); } rows = faidx_nseq(idx); @@ -226,8 +235,7 @@ namespace msa_chunker { if (c == -1) { c = seq_len; } else if (seq_len != c) { - cerr << "ERROR: MSA has rows of different length! (" << string(seq_name) << ")" << endl; - exit(1); + throw runtime_error("ERROR: MSA has rows of different length! (" + string(seq_name) + ")"); } } cols = c; @@ -293,7 +301,7 @@ namespace msa_chunker { public: column_chunker() = delete; - column_chunker(const string &msapath, const msa_pos_t max_qlen, const bool verbose) { + column_chunker(const string &msapath, const msa_pos_t max_qlen, const int verbosity) { max_chunk_cols = max(MIN_CHUNK_COLS, max_qlen); msa_file = ifstream(msapath, std::ios::binary); if (!msa_file) { @@ -306,6 +314,9 @@ namespace msa_chunker { if (!(iss >> cols >> rows) || (iss >> extra)) { throw runtime_error("ERROR: first line must contain exactly two integers: columns and rows"); } + if (verbosity > 0) { + cerr << "Found column major matrix file" << endl; + } matrix_start = msa_file.tellg(); } From 54f11fa47394b075dd40a85a6d372d3b2b479a33 Mon Sep 17 00:00:00 2001 From: Sebastian Visan Date: Mon, 10 Aug 2026 13:13:06 +0300 Subject: [PATCH 09/12] =?UTF-8?q?=E2=9C=A8=20minsize=20gaps=20algorithm?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- src/minsize.hpp | 90 ++++++++++++++++++++++++++++++++++++++++++------- 1 file changed, 77 insertions(+), 13 deletions(-) diff --git a/src/minsize.hpp b/src/minsize.hpp index 53d7164..074a24b 100644 --- a/src/minsize.hpp +++ b/src/minsize.hpp @@ -143,6 +143,64 @@ namespace minsize { update(height, y, n[y] - height * y); } } + // Calculate the min-L-size for the next column (y is 1-index) using susffix tree algo + void next_gaps(const seg_index y) { + set reverse_unique_chunk; + for (seg_index i = 0; i < r; i++) { + string s = idx.msa_substr(i, y - min(U, y), min(U, y)); + if (!gaps_as_symbols) + s.erase(remove(s.begin(), s.end(), '-'), s.end()); + reverse(s.begin(), s.end()); + reverse_unique_chunk.insert(move(s)); + } + + trie::trie T(reverse_unique_chunk); + reverse_unique_chunk.clear(); + set reverse_chunk; + for (seg_index i = 0; i < r; i++) { + string s = idx.msa_substr(i, y - min(U, y), min(U, y)); + reverse(s.begin(), s.end()); + reverse_chunk.insert(s); + } + + vector counts(T.nodes() + 1, 0); + vector node_length(T.nodes() + 1, 0); // length from the root + counts[T.nodes()] = reverse_chunk.size(); // current count of the root + vector active_node(reverse_chunk.size(), T.nodes()); + for (seg_index len = 1; len <= min(U, y); len++) { + auto it = reverse_chunk.begin(); + for (seg_index i = 0; i < reverse_chunk.size(); ++i, ++it) { + if ((*it)[len - 1] != '-') { + counts[active_node[i]] -= 1; + auto active_length = node_length[active_node[i]]; + if (active_node[i] == T.nodes()) { + active_node[i] = T.child(T.root(), (*it)[len - 1]); + } else { + active_node[i] = T.child(active_node[i], (*it)[len - 1]); + } + assert(active_node[i] != trie::trie::null); + counts[active_node[i]] += 1; + node_length[active_node[i]] = active_length + 1; + } + } + assert(it == reverse_chunk.end()); + // Calculate the total length of the active nodes (segmentation size) + // Update min-size arrays n and back + if (len >= L) { + long long segment_size = 0; + for (seg_index i = 0; i <= T.nodes(); i++) { + if (counts[i] > 0) { + // empty string has size 1 + segment_size += max(1ULL, node_length[i]); + } + } + if (n[y] > n[y - len] + segment_size && n[y - len] + segment_size > 0) { + n[y] = n[y - len] + segment_size; + back[y] = y - len; + } + } + } + } public: minsize( msa_t &idx, @@ -183,20 +241,26 @@ namespace minsize { pair>> segment( const vector>> &L_y = {} ) { - // Algorithm for gaps as symbols - // TODO: Add trie algo for gaps - for (i_type y = 1; y <= c; y++) { - // compute L_y if it was not given in input - const vector> *L_yy; - vector> L_yy_on_the_fly; - if (L_y.size() > 0) { - L_yy = &(*(L_y.begin() + y)); - } else { - L_yy_on_the_fly = algo::compute_meaningful_extensions(idx, r, c, L, U, y, gaps_as_symbols, use_pbwt); - L_yy = &L_yy_on_the_fly; + if (gaps_as_symbols) { + // algorithm for gaps as symbols from the paper + for (i_type y = 1; y <= c; y++) { + // compute L_y if it was not given in input + const vector> *L_yy; + vector> L_yy_on_the_fly; + if (L_y.size() > 0) { + L_yy = &(*(L_y.begin() + y)); + } else { + L_yy_on_the_fly = algo::compute_meaningful_extensions(idx, r, c, L, U, y, gaps_as_symbols, use_pbwt); + L_yy = &L_yy_on_the_fly; + } + // Algorithm for column y + next(y, L_yy); + } + } else { + // suffix trie algorithm for gaps + for (i_type y = 1; y <= c; y++) { + next_gaps(y); } - // Algorithm for column y - next(y, L_yy); } // trace back From cbcf5273ccefbf178084ba56870cc68656a3fbcd Mon Sep 17 00:00:00 2001 From: Sebastian Visan Date: Mon, 10 Aug 2026 18:02:23 +0300 Subject: [PATCH 10/12] =?UTF-8?q?=E2=9A=97=EF=B8=8F=20=20cardinality=20vs?= =?UTF-8?q?=20size=20experiment?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- .gitignore | 3 +- experiments/card_vs_size/README.md | 7 ++ experiments/card_vs_size/run_experiment.sh | 104 +++++++++++++++++++++ 3 files changed, 113 insertions(+), 1 deletion(-) create mode 100644 experiments/card_vs_size/README.md create mode 100755 experiments/card_vs_size/run_experiment.sh diff --git a/.gitignore b/.gitignore index d0b7ae7..2d9c68a 100644 --- a/.gitignore +++ b/.gitignore @@ -33,4 +33,5 @@ experiments/minsize/mincard-ring experiments/minsize/mincard-tree experiments/minsize/mincard-queue -experiments/transpose/output \ No newline at end of file +experiments/transpose/output +experiments/card_vs_size/output diff --git a/experiments/card_vs_size/README.md b/experiments/card_vs_size/README.md new file mode 100644 index 0000000..6a11c8e --- /dev/null +++ b/experiments/card_vs_size/README.md @@ -0,0 +1,7 @@ +# cardinality vs size experiment + +This experiment compares the resulting eds stats and running times for all 4 algorithms (mincard/minsize x gaps/gapless) with various bounds (L and U) on all 3 datasets. + +Get the ecoli, covid, and chr datasets in their respective experiments. Then, simply run the experiment. + +This experiment does not stream the eds-es, though it could with a minor change in the script. The results are collected automatically in output/results.csv diff --git a/experiments/card_vs_size/run_experiment.sh b/experiments/card_vs_size/run_experiment.sh new file mode 100755 index 0000000..7962c76 --- /dev/null +++ b/experiments/card_vs_size/run_experiment.sh @@ -0,0 +1,104 @@ +#!/bin/bash +set -euo pipefail +thisfolder=$( cd -- "$( dirname -- "${BASH_SOURCE[0]}" )" &> /dev/null && pwd ) +cd "$thisfolder" +export LC_NUMERIC="en_US.UTF-8" + +mincard="$thisfolder/../../mincard" + +# MSA fasta files from the 3 datasets +input_ecoli_msa="$thisfolder/../e_coli_sim/input/msa.fa" +input_covid_msa="$thisfolder/../covid/input/covid_100000.fa.gz" +input_chr19_msa="$thisfolder/../chr19/input/chr19_100.aligned.uppercase.fa" + +# Create links to msa files +outdir="output" +mkdir -p "$outdir" +rm -f "$outdir"/* +cd "$outdir" + +ln -s "$input_ecoli_msa" ecoli_msa.fa +ln -s "$input_covid_msa" covid_msa.fa +ln -s "$input_chr19_msa" chr19_msa.fa + +# Keep them in an array for easy iteration +declare -A msa + +datasets=( + ecoli + covid + chr19 +) + +msa[ecoli]=ecoli_msa.fa +msa[covid]=covid_msa.fa +msa[chr19]=chr19_msa.fa + +# Store algorithm flags +declare -A flags + +algos=( + mincard_gaps + mincard_gapless + minsize_gaps + minsize_gapless +) + +flags[mincard_gaps]="" +flags[mincard_gapless]="--gaps-as-symbols --pbwt" +flags[minsize_gaps]="--min-size" +flags[minsize_gapless]="--min-size --gaps-as-symbols --pbwt" + +# Create temporary files for collecting results +outfile=$(mktemp) +timefile=$(mktemp) +trap 'rm -f "$outfile"' EXIT +trap 'rm -f "$timefile"' EXIT + +usrbintimeformat="%e" + +RED='\033[0;31m' +GREEN='\033[0;32m' +BLUE='\033[0;34m' +PURPLE='\033[0;35m' +NC='\033[0m' + +# Iterate through different [L, U] bounds + +values=("1 16" "16 32" "32 64") +printf "Running mincard/minsize with gaps/gapless strategies on MSAs with LU bounds" +for pair in "${values[@]}"; do + printf " [%s]" "${pair/ /, }" +done +printf "\n" + + +# Write results to csv table +echo -e "Dataset,Algorithm,L,U,Cardinality,Size,Avg. Segment,Time" > results.csv + +for dataset in "${datasets[@]}"; do + input_msa="${msa[$dataset]}" + + for pair in "${values[@]}"; do + read lower upper <<< "$pair" + printf "\nDataset: ${BLUE}${dataset}${NC} | L = ${BLUE}${lower}${NC} U = ${BLUE}${upper}${NC}\n" + + for algo in "${algos[@]}"; do + /usr/bin/time -f"$usrbintimeformat" -o "$timefile" \ + "$mincard" "$input_msa" \ + -L "$lower" -U "$upper" \ + -q --stats \ + ${flags[$algo]} > "$outfile" 2>&1; + + time=$(cat "$timefile") + read card size <<< "$(head -n1 "$outfile")" + read min max avg <<< "$(tail -n1 "$outfile")" + + printf "Algorithm: ${GREEN}${algo}${NC} | cardinality = ${PURPLE}${card}${NC} | gap-aware size = ${PURPLE}${size}${NC} | segment sizes: min = ${PURPLE}${min}${NC}, max = ${PURPLE}${max}${NC}, avg = ${PURPLE}${avg}${NC} | ${BLUE}${time}s${NC}\n" + echo -e "${dataset},${algo},${lower},${upper},${card},${size},${avg},${time}" >> results.csv + done + done +done + +printf "\n${GREEN}Finished running the experiment!\n" +printf "${BLUE}Wrote algorithm times and stats to results.csv${NC}\n" From af78afc066c9bccf7c6695fa1256d429dacc88e2 Mon Sep 17 00:00:00 2001 From: Sebastian Visan Date: Thu, 27 Aug 2026 11:54:23 +0300 Subject: [PATCH 11/12] =?UTF-8?q?=E2=99=BB=EF=B8=8F=20=20separate=20minsiz?= =?UTF-8?q?e=20tool?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- .gitignore | 1 + Makefile | 5 ++- experiments/card_vs_size/README.md | 6 ++++ experiments/card_vs_size/run_experiment.sh | 13 ++++++-- experiments/minsize/README.md | 6 ++-- experiments/minsize/run_experiment.sh | 12 +++---- src/mincard.cpp | 38 +++++++++++++--------- 7 files changed, 52 insertions(+), 29 deletions(-) diff --git a/.gitignore b/.gitignore index 2d9c68a..0cd479d 100644 --- a/.gitignore +++ b/.gitignore @@ -1,4 +1,5 @@ mincard +minsize test_rmqueue test/*.eds test/*.gfa diff --git a/Makefile b/Makefile index 96ca1b8..1988fb0 100644 --- a/Makefile +++ b/Makefile @@ -12,8 +12,11 @@ SDSL_INCLUDE=ext/sdsl-lite/include/ mincard: src/mincard.cpp src/segment.hpp src/rmqueue.h src/rmqueue.cpp src/RMaxQTree.h src/RMaxQTree.cpp src/msa_chunker.hpp src/trie.hpp src/pbwt.h src/pbwt.cpp src/algo.hpp src/minsize.hpp ${CXX} $(CXX_FLAGS) -DVERSION="\"$(VERSION)\"" -o mincard src/mincard.cpp src/rmqueue.cpp src/RMaxQTree.cpp src/pbwt.cpp -I $(HTSLIB_INCLUDE) -I $(OTHER_INCLUDE) -I $(SDSL_INCLUDE) $(HTSLIB_FLAGS) +minsize: src/mincard.cpp src/segment.hpp src/rmqueue.h src/rmqueue.cpp src/RMaxQTree.h src/RMaxQTree.cpp src/msa_chunker.hpp src/trie.hpp src/pbwt.h src/pbwt.cpp src/algo.hpp src/minsize.hpp + ${CXX} $(CXX_FLAGS) -DMETRIC=SIZE -DVERSION="\"$(VERSION)\"" -o minsize src/mincard.cpp src/rmqueue.cpp src/RMaxQTree.cpp src/pbwt.cpp -I $(HTSLIB_INCLUDE) -I $(OTHER_INCLUDE) -I $(SDSL_INCLUDE) $(HTSLIB_FLAGS) + clean: - rm -f mincard + rm -f mincard && rm -f minsize && rm -f test_rmqueue test_rmqueue: src/test_rmqueue.cpp src/rmqueue.cpp ${CXX} $(CXX_FLAGS) -o test_rmqueue src/test_rmqueue.cpp src/rmqueue.cpp diff --git a/experiments/card_vs_size/README.md b/experiments/card_vs_size/README.md index 6a11c8e..db936da 100644 --- a/experiments/card_vs_size/README.md +++ b/experiments/card_vs_size/README.md @@ -2,6 +2,12 @@ This experiment compares the resulting eds stats and running times for all 4 algorithms (mincard/minsize x gaps/gapless) with various bounds (L and U) on all 3 datasets. +Compile `mincard` and `minsize`: +``` +make -C ../../ +make minsize -C ../../ +``` + Get the ecoli, covid, and chr datasets in their respective experiments. Then, simply run the experiment. This experiment does not stream the eds-es, though it could with a minor change in the script. The results are collected automatically in output/results.csv diff --git a/experiments/card_vs_size/run_experiment.sh b/experiments/card_vs_size/run_experiment.sh index 7962c76..27643ce 100755 --- a/experiments/card_vs_size/run_experiment.sh +++ b/experiments/card_vs_size/run_experiment.sh @@ -5,6 +5,7 @@ cd "$thisfolder" export LC_NUMERIC="en_US.UTF-8" mincard="$thisfolder/../../mincard" +minsize="$thisfolder/../../minsize" # MSA fasta files from the 3 datasets input_ecoli_msa="$thisfolder/../e_coli_sim/input/msa.fa" @@ -36,6 +37,7 @@ msa[chr19]=chr19_msa.fa # Store algorithm flags declare -A flags +declare -A alg algos=( mincard_gaps @@ -45,9 +47,14 @@ algos=( ) flags[mincard_gaps]="" +flags[minsize_gaps]="" flags[mincard_gapless]="--gaps-as-symbols --pbwt" -flags[minsize_gaps]="--min-size" -flags[minsize_gapless]="--min-size --gaps-as-symbols --pbwt" +flags[minsize_gapless]="--gaps-as-symbols --pbwt" + +alg[mincard_gaps]="$mincard" +alg[mincard_gapless]="$mincard" +alg[minsize_gaps]="$minsize" +alg[minsize_gapless]="$minsize" # Create temporary files for collecting results outfile=$(mktemp) @@ -85,7 +92,7 @@ for dataset in "${datasets[@]}"; do for algo in "${algos[@]}"; do /usr/bin/time -f"$usrbintimeformat" -o "$timefile" \ - "$mincard" "$input_msa" \ + ${alg[$algo]} "$input_msa" \ -L "$lower" -U "$upper" \ -q --stats \ ${flags[$algo]} > "$outfile" 2>&1; diff --git a/experiments/minsize/README.md b/experiments/minsize/README.md index a56e429..f38f0e2 100644 --- a/experiments/minsize/README.md +++ b/experiments/minsize/README.md @@ -4,11 +4,11 @@ Simple minsize experiment to compare the different ways of calculating the range Get covid dataset as described [here](../covid/README.md) -Compile mincard with each algorithm (ring buffer, rmaxtree, rmqueue) by changing [this line](../../src/minsize.hpp#L25) and then running +Compile `minsize` with each algorithm (ring buffer, rmaxtree, rmqueue) by changing [this line](../../src/minsize.hpp#L25) and then running ```sh - make -C ../../ - mv ../../mincard ./mincard-ring + make minsize -C ../../ + mv ../../minsize ./minsize-ring ``` Do this with the names [ring, tree, queue]. Finally run the experiment script (running times will be collected in output/results.csv) diff --git a/experiments/minsize/run_experiment.sh b/experiments/minsize/run_experiment.sh index 486f398..0ed1178 100755 --- a/experiments/minsize/run_experiment.sh +++ b/experiments/minsize/run_experiment.sh @@ -3,9 +3,9 @@ set -euo pipefail thisfolder=$( cd -- "$( dirname -- "${BASH_SOURCE[0]}" )" &> /dev/null && pwd ) # https://stackoverflow.com/questions/59895/how-do-i-get-the-directory-where-a-bash-script-is-located-from-within-the-script cd "$thisfolder" -mincard_ring="$thisfolder/mincard-ring" -mincard_tree="$thisfolder/mincard-tree" -mincard_queue="$thisfolder/mincard-queue" +minsize_ring="$thisfolder/minsize-ring" +minsize_tree="$thisfolder/minsize-tree" +minsize_queue="$thisfolder/minsize-queue" inputmsa="$thisfolder/input/covid_100000.fa.gz" usrbintimeformat="%e" @@ -35,13 +35,13 @@ for L in "${L_values[@]}" do printf "\n${BLUE}Comparison for L = $L${NC}\n" # ring buffer - /usr/bin/time -f"$usrbintimeformat" -o "$tmpfile" "$mincard_ring" msa.fa -v --gaps-as-symbols --pbwt -L $L --min-size -o minsize_ring_L${L}.eds + /usr/bin/time -f"$usrbintimeformat" -o "$tmpfile" "$minsize_ring" msa.fa -v --gaps-as-symbols --pbwt -L $L --min-size -o minsize_ring_L${L}.eds t1=$(<"$tmpfile") # rmaxtree - /usr/bin/time -f"$usrbintimeformat" -o "$tmpfile" "$mincard_tree" msa.fa -v --gaps-as-symbols --pbwt -L $L --min-size -o minsize_tree_L${L}.eds + /usr/bin/time -f"$usrbintimeformat" -o "$tmpfile" "$minsize_tree" msa.fa -v --gaps-as-symbols --pbwt -L $L --min-size -o minsize_tree_L${L}.eds t2=$(<"$tmpfile") # rmqueue - /usr/bin/time -f"$usrbintimeformat" -o "$tmpfile" "$mincard_queue" msa.fa -v --gaps-as-symbols --pbwt -L $L --min-size -o minsize_queue_L${L}.eds + /usr/bin/time -f"$usrbintimeformat" -o "$tmpfile" "$minsize_queue" msa.fa -v --gaps-as-symbols --pbwt -L $L --min-size -o minsize_queue_L${L}.eds t3=$(<"$tmpfile") if cmp -s minsize_ring_L${L}.eds minsize_tree_L${L}.eds && diff --git a/src/mincard.cpp b/src/mincard.cpp index 09f9a90..3f74078 100644 --- a/src/mincard.cpp +++ b/src/mincard.cpp @@ -23,6 +23,13 @@ using algo::seg_index, algo::compute_perfect_columns, algo::compute_all_meaningf typedef segment::seg_size_t seg_size_t; using segment::segment_stream_gfa, segment::segment_stream_eds, segment::segment_stream_no_output; +// The current tool optimizes 2 metrics: cardinality and size +#ifndef METRIC +#define METRIC CARDINALITY +#endif +#define CARDINALITY 0 +#define SIZE 1 + int main(int argc, char* argv[]) { CLI::App app{"mincard version " + string(VERSION) + " — build Elastic Degenerate Strings (EDSes) from multiple sequence alignments in FASTA format"}; argv = app.ensure_utf8(argv); @@ -49,9 +56,6 @@ int main(int argc, char* argv[]) { ->default_val(0) ->expected(1, numeric_limits::max()); - bool min_size = false; - app.add_flag("--min-size", min_size, "Minimize the size of the segmentation instead of the cardinality"); - bool stats = false; app.add_flag("--stats", stats, "Calculate segmentation statistics"); @@ -93,7 +97,11 @@ int main(int argc, char* argv[]) { } // No U value specified by the user if (U == 0) { - if (min_size) { + #if METRIC == CARDINALITY + // Default upper bound for min-card + U = 31; + #endif + #if METRIC == SIZE if (gaps_as_symbols) { // Implicit upper bound that keeps the segmentation optimal U = L * 2 - 1; @@ -102,10 +110,7 @@ int main(int argc, char* argv[]) { // Arbitrary upper bound so the algorithm is practical U = L * 4 - 1; } - } else { - // Default upper bound for min-card - U = 31; - } + #endif } if (L > U) { cerr << "Upper and lower bounds are not compatible!" << endl; @@ -222,23 +227,24 @@ int main(int argc, char* argv[]) { auto start = high_resolution_clock::now(); seg_index minval; // cardinality or size - if (min_size) { + #if METRIC == CARDINALITY + if (verbosity > 0) { + cerr << "Computing the minimum-cardinality segmentation..." << flush; + } + tie(minval, segmentation) = segment_with_rmq(idx, r, c, L, U, gaps_as_symbols, use_pbwt, L_y, perfect_columns); + #endif + #if METRIC == SIZE if (verbosity > 0) { cerr << "Computing the minimum-size segmentation..." << flush; } minsize::minsize alg(idx, r, c, L, U, gaps_as_symbols, use_pbwt); tie(minval, segmentation) = alg.segment(); - } else { - if (verbosity > 0) { - cerr << "Computing the minimum-cardinality segmentation..." << flush; - } - tie(minval, segmentation) = segment_with_rmq(idx, r, c, L, U, gaps_as_symbols, use_pbwt, L_y, perfect_columns); - } + #endif auto stop = high_resolution_clock::now(); auto duration = duration_cast(stop - start); if (minval != std::numeric_limits::max()) { if (verbosity > 0) { - cerr << " done: " << segmentation.size() << " segments/ED words, " << minval << (min_size ? " size" : " cardinality") + cerr << " done: " << segmentation.size() << " segments/ED words, " << minval << (METRIC == CARDINALITY ? " cardinality" : " size") << ((verbosity > 1) ? " (" + to_string(duration.count()) + "ms)" : "") << endl; } } else { From 3b1bc7e0437dec64a06c075e6378b2ef041f8564 Mon Sep 17 00:00:00 2001 From: Sebastian Visan Date: Thu, 27 Aug 2026 13:30:57 +0300 Subject: [PATCH 12/12] =?UTF-8?q?=E2=99=BB=EF=B8=8F=20=20set=20gaps-as-sym?= =?UTF-8?q?bols=20and=20pbwt=20as=20default?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- experiments/card_vs_size/run_experiment.sh | 8 ++++---- experiments/chr19/bench_pbwt.sh | 4 ++-- experiments/chr19/run_experiment.sh | 8 ++++---- experiments/covid/bench_pbwt.sh | 4 ++-- experiments/covid/run_experiment.sh | 12 ++++++------ experiments/e_coli_sim/bench_pbwt.sh | 4 ++-- experiments/e_coli_sim/run_experiment.sh | 8 ++++---- experiments/minsize/run_experiment.sh | 6 +++--- experiments/pbwt/run_experiment.sh | 6 +++--- experiments/transpose/run_experiment.sh | 8 ++++---- src/mincard.cpp | 18 ++++++++---------- 11 files changed, 42 insertions(+), 44 deletions(-) diff --git a/experiments/card_vs_size/run_experiment.sh b/experiments/card_vs_size/run_experiment.sh index 27643ce..20b561d 100755 --- a/experiments/card_vs_size/run_experiment.sh +++ b/experiments/card_vs_size/run_experiment.sh @@ -46,10 +46,10 @@ algos=( minsize_gapless ) -flags[mincard_gaps]="" -flags[minsize_gaps]="" -flags[mincard_gapless]="--gaps-as-symbols --pbwt" -flags[minsize_gapless]="--gaps-as-symbols --pbwt" +flags[mincard_gaps]="--gaps" +flags[minsize_gaps]="--gaps" +flags[mincard_gapless]="" +flags[minsize_gapless]="" alg[mincard_gaps]="$mincard" alg[mincard_gapless]="$mincard" diff --git a/experiments/chr19/bench_pbwt.sh b/experiments/chr19/bench_pbwt.sh index 15f7fc0..0a35be1 100755 --- a/experiments/chr19/bench_pbwt.sh +++ b/experiments/chr19/bench_pbwt.sh @@ -36,11 +36,11 @@ for U in "${U_values[@]}" do printf "\n${BLUE}Comparison for U = $U${NC}\n" # trie run - /usr/bin/time -f"$usrbintimeformat" -o "$tmpfile" "$mincard" msa.fa -v --gaps-as-symbols -U $U -o mincard_U${U}.eds + /usr/bin/time -f"$usrbintimeformat" -o "$tmpfile" "$mincard" msa.fa -v --trie -U $U -o mincard_U${U}.eds t1=$(<"$tmpfile") trie_times+=("$t1") # pbwt run - /usr/bin/time -f"$usrbintimeformat" -o "$tmpfile" "$mincard" msa.fa --gaps-as-symbols --pbwt -U $U -o mincard_U${U}_pbwt.eds + /usr/bin/time -f"$usrbintimeformat" -o "$tmpfile" "$mincard" msa.fa -v -U $U -o mincard_U${U}_pbwt.eds t2=$(<"$tmpfile") pbwt_times+=("$t2") diff --git a/experiments/chr19/run_experiment.sh b/experiments/chr19/run_experiment.sh index c49f4db..a19e4d6 100755 --- a/experiments/chr19/run_experiment.sh +++ b/experiments/chr19/run_experiment.sh @@ -17,20 +17,20 @@ ln -s $inputmsa msa.fa # mincard for U in 4 8 16 do - /usr/bin/time -f"$usrbintimeformat" $mincard msa.fa -v -U $U -o mincard_U${U}.eds + /usr/bin/time -f"$usrbintimeformat" $mincard msa.fa -v --gaps -U $U -o mincard_U${U}.eds done # mincard with perfect segments for U in 4 8 16 do - /usr/bin/time -f"$usrbintimeformat" $mincard msa.fa -v -U $U --perfect-segments -o mincard_U${U}_p.eds + /usr/bin/time -f"$usrbintimeformat" $mincard msa.fa -v --gaps -U $U --perfect-segments -o mincard_U${U}_p.eds done # mincard trivial S^||| -/usr/bin/time -f"$usrbintimeformat" $mincard msa.fa -v --trivial-vertical -o mincard_t.eds +/usr/bin/time -f"$usrbintimeformat" $mincard msa.fa -v --gaps --trivial-vertical -o mincard_t.eds # mincard trivial S^≡ with perfect segments -/usr/bin/time -f"$usrbintimeformat" $mincard msa.fa -v --trivial-horizontal --perfect-segments -o mincard_np.eds +/usr/bin/time -f"$usrbintimeformat" $mincard msa.fa -v --gaps --trivial-horizontal --perfect-segments -o mincard_np.eds exit # msatoeds heuristics, they require >= 100GB RAM diff --git a/experiments/covid/bench_pbwt.sh b/experiments/covid/bench_pbwt.sh index e2fb9d7..66b3d0f 100755 --- a/experiments/covid/bench_pbwt.sh +++ b/experiments/covid/bench_pbwt.sh @@ -35,11 +35,11 @@ for U in "${U_values[@]}" do printf "\n${BLUE}Comparison for U = $U${NC}\n" # trie run - /usr/bin/time -f"$usrbintimeformat" -o "$tmpfile" "$mincard" msa.fa --gaps-as-symbols -U $U -o mincard_U${U}.eds + /usr/bin/time -f"$usrbintimeformat" -o "$tmpfile" "$mincard" msa.fa --trie -U $U -o mincard_U${U}.eds t1=$(<"$tmpfile") trie_times+=("$t1") # pbwt run - /usr/bin/time -f"$usrbintimeformat" -o "$tmpfile" "$mincard" msa.fa --gaps-as-symbols --pbwt -U $U -o mincard_U${U}_pbwt.eds + /usr/bin/time -f"$usrbintimeformat" -o "$tmpfile" "$mincard" msa.fa -U $U -o mincard_U${U}_pbwt.eds t2=$(<"$tmpfile") pbwt_times+=("$t2") diff --git a/experiments/covid/run_experiment.sh b/experiments/covid/run_experiment.sh index 8724fc3..395e498 100755 --- a/experiments/covid/run_experiment.sh +++ b/experiments/covid/run_experiment.sh @@ -20,21 +20,21 @@ do # mincard for U in 4 8 16 32 64 128 256 512 do - /usr/bin/time -f"$usrbintimeformat" $mincard msa.fa.gz -v -U $U -o ${base}_mincard_U${U}.eds # plain + /usr/bin/time -f"$usrbintimeformat" $mincard msa.fa.gz -v --gaps -U $U -o ${base}_mincard_U${U}.eds # plain rm msa.fa.gz.fai msa.fa.gz.gzi - /usr/bin/time -f"$usrbintimeformat" $mincard msa.fa.gz -v -U $U --perfect-segments -o ${base}_mincard_U${U}_p.eds # perfect segments + /usr/bin/time -f"$usrbintimeformat" $mincard msa.fa.gz -v --gaps -U $U --perfect-segments -o ${base}_mincard_U${U}_p.eds # perfect segments rm msa.fa.gz.fai msa.fa.gz.gzi - /usr/bin/time -f"$usrbintimeformat" $mincard msa.fa.gz --preprocess -v -U $U -o ${base}_mincard_U${U}.eds # preprocess + /usr/bin/time -f"$usrbintimeformat" $mincard msa.fa.gz --gaps --preprocess -v -U $U -o ${base}_mincard_U${U}.eds # preprocess rm msa.fa.gz.fai msa.fa.gz.gzi - /usr/bin/time -f"$usrbintimeformat" $mincard msa.fa.gz -v -U $U --perfect-segments --preprocess -o ${base}_mincard_U${U}_p.eds + /usr/bin/time -f"$usrbintimeformat" $mincard msa.fa.gz -v --gaps -U $U --perfect-segments --preprocess -o ${base}_mincard_U${U}_p.eds rm msa.fa.gz.fai msa.fa.gz.gzi done # mincard trivial S^||| - /usr/bin/time -f"$usrbintimeformat" $mincard msa.fa.gz -v --trivial-vertical -o ${base}_mincard_t.eds + /usr/bin/time -f"$usrbintimeformat" $mincard msa.fa.gz -v --gaps --trivial-vertical -o ${base}_mincard_t.eds # mincard trivial S^≡ with perfect segments - /usr/bin/time -f"$usrbintimeformat" $mincard msa.fa.gz -v --trivial-horizontal --perfect-segments -o ${base}_mincard_np.eds + /usr/bin/time -f"$usrbintimeformat" $mincard msa.fa.gz -v --gaps --trivial-horizontal --perfect-segments -o ${base}_mincard_np.eds # msatoeds heuristics for strat in trivial greedy double-greedy diff --git a/experiments/e_coli_sim/bench_pbwt.sh b/experiments/e_coli_sim/bench_pbwt.sh index d9a8a81..1f48841 100755 --- a/experiments/e_coli_sim/bench_pbwt.sh +++ b/experiments/e_coli_sim/bench_pbwt.sh @@ -35,11 +35,11 @@ for U in "${U_values[@]}" do printf "\n${BLUE}Comparison for U = $U${NC}\n" # trie run - /usr/bin/time -f"$usrbintimeformat" -o "$tmpfile" "$mincard" msa.fa --gaps-as-symbols -U $U -o mincard_U${U}.eds + /usr/bin/time -f"$usrbintimeformat" -o "$tmpfile" "$mincard" msa.fa --trie -U $U -o mincard_U${U}.eds t1=$(<"$tmpfile") trie_times+=("$t1") # pbwt run - /usr/bin/time -f"$usrbintimeformat" -o "$tmpfile" "$mincard" msa.fa --gaps-as-symbols --pbwt -U $U -o mincard_U${U}_pbwt.eds + /usr/bin/time -f"$usrbintimeformat" -o "$tmpfile" "$mincard" msa.fa -U $U -o mincard_U${U}_pbwt.eds t2=$(<"$tmpfile") pbwt_times+=("$t2") diff --git a/experiments/e_coli_sim/run_experiment.sh b/experiments/e_coli_sim/run_experiment.sh index c67ce6d..d6538ca 100755 --- a/experiments/e_coli_sim/run_experiment.sh +++ b/experiments/e_coli_sim/run_experiment.sh @@ -16,20 +16,20 @@ ln -s $inputmsa msa.fa # mincard for U in 4 8 16 32 64 do - /usr/bin/time -f"$usrbintimeformat" $mincard msa.fa -v -U $U -o mincard_U${U}.eds + /usr/bin/time -f"$usrbintimeformat" $mincard msa.fa -v --gaps -U $U -o mincard_U${U}.eds done # mincard with perfect segments for U in 4 8 16 32 64 do - /usr/bin/time -f"$usrbintimeformat" $mincard msa.fa -v -U $U --perfect-segments -o mincard_U${U}_p.eds + /usr/bin/time -f"$usrbintimeformat" $mincard msa.fa -v --gaps -U $U --perfect-segments -o mincard_U${U}_p.eds done # mincard trivial S^||| -/usr/bin/time -f"$usrbintimeformat" $mincard msa.fa -v --trivial-vertical -o mincard_t.eds +/usr/bin/time -f"$usrbintimeformat" $mincard msa.fa -v --gaps --trivial-vertical -o mincard_t.eds # mincard trivial S^≡ with perfect segments -/usr/bin/time -f"$usrbintimeformat" $mincard msa.fa -v --trivial-horizontal --perfect-segments -o mincard_np.eds +/usr/bin/time -f"$usrbintimeformat" $mincard msa.fa -v --gaps --trivial-horizontal --perfect-segments -o mincard_np.eds # msatoeds heuristics for strat in trivial greedy double-greedy diff --git a/experiments/minsize/run_experiment.sh b/experiments/minsize/run_experiment.sh index 0ed1178..132a1fb 100755 --- a/experiments/minsize/run_experiment.sh +++ b/experiments/minsize/run_experiment.sh @@ -35,13 +35,13 @@ for L in "${L_values[@]}" do printf "\n${BLUE}Comparison for L = $L${NC}\n" # ring buffer - /usr/bin/time -f"$usrbintimeformat" -o "$tmpfile" "$minsize_ring" msa.fa -v --gaps-as-symbols --pbwt -L $L --min-size -o minsize_ring_L${L}.eds + /usr/bin/time -f"$usrbintimeformat" -o "$tmpfile" "$minsize_ring" msa.fa -v -L $L --min-size -o minsize_ring_L${L}.eds t1=$(<"$tmpfile") # rmaxtree - /usr/bin/time -f"$usrbintimeformat" -o "$tmpfile" "$minsize_tree" msa.fa -v --gaps-as-symbols --pbwt -L $L --min-size -o minsize_tree_L${L}.eds + /usr/bin/time -f"$usrbintimeformat" -o "$tmpfile" "$minsize_tree" msa.fa -v -L $L --min-size -o minsize_tree_L${L}.eds t2=$(<"$tmpfile") # rmqueue - /usr/bin/time -f"$usrbintimeformat" -o "$tmpfile" "$minsize_queue" msa.fa -v --gaps-as-symbols --pbwt -L $L --min-size -o minsize_queue_L${L}.eds + /usr/bin/time -f"$usrbintimeformat" -o "$tmpfile" "$minsize_queue" msa.fa -v -L $L --min-size -o minsize_queue_L${L}.eds t3=$(<"$tmpfile") if cmp -s minsize_ring_L${L}.eds minsize_tree_L${L}.eds && diff --git a/experiments/pbwt/run_experiment.sh b/experiments/pbwt/run_experiment.sh index 78f67d2..1dbb971 100755 --- a/experiments/pbwt/run_experiment.sh +++ b/experiments/pbwt/run_experiment.sh @@ -36,15 +36,15 @@ for U in "${U_values[@]}" do printf "\n${BLUE}Comparison for U = $U${NC}\n" # naive - /usr/bin/time -f"$usrbintimeformat" -o "$tmpfile" "$mincard_naive" msa.fa --gaps-as-symbols --pbwt -U $U -o mincard_naive_U${U}.eds + /usr/bin/time -f"$usrbintimeformat" -o "$tmpfile" "$mincard_naive" msa.fa -U $U -o mincard_naive_U${U}.eds t1=$(<"$tmpfile") naive_times+=("$t1") # recursive - /usr/bin/time -f"$usrbintimeformat" -o "$tmpfile" "$mincard_recursive" msa.fa --gaps-as-symbols --pbwt -U $U -o mincard_recursive_U${U}.eds + /usr/bin/time -f"$usrbintimeformat" -o "$tmpfile" "$mincard_recursive" msa.fa -U $U -o mincard_recursive_U${U}.eds t2=$(<"$tmpfile") recursive_times+=("$t2") # rmq - /usr/bin/time -f"$usrbintimeformat" -o "$tmpfile" "$mincard_rmq" msa.fa --gaps-as-symbols --pbwt -U $U -o mincard_rmq_U${U}.eds + /usr/bin/time -f"$usrbintimeformat" -o "$tmpfile" "$mincard_rmq" msa.fa -U $U -o mincard_rmq_U${U}.eds t3=$(<"$tmpfile") rmq_times+=("$t3") diff --git a/experiments/transpose/run_experiment.sh b/experiments/transpose/run_experiment.sh index 54b61bd..f7d63ef 100755 --- a/experiments/transpose/run_experiment.sh +++ b/experiments/transpose/run_experiment.sh @@ -72,14 +72,14 @@ for name in "${datasets[@]}"; do printf "\n${BLUE}Dataset: ${name} | U = ${U}${NC}\n" # suffix trie - /usr/bin/time -f"$usrbintimeformat" -o "$tmpfile" "$mincard" "$input_msa_row" -v --gaps-as-symbols -U $U -o mincard_trie_U${U}.eds + /usr/bin/time -f"$usrbintimeformat" -o "$tmpfile" "$mincard" "$input_msa_row" -v --trie -U $U -o mincard_trie_U${U}.eds t1=$(<"$tmpfile") - /usr/bin/time -f"$usrbintimeformat" -o "$tmpfile" "$mincard" "$input_msa_col" -v --gaps-as-symbols --column-major -U $U -o mincard_cm_trie_U${U}.eds + /usr/bin/time -f"$usrbintimeformat" -o "$tmpfile" "$mincard" "$input_msa_col" -v --trie --column-major -U $U -o mincard_cm_trie_U${U}.eds t2=$(<"$tmpfile") # pbwt - /usr/bin/time -f"$usrbintimeformat" -o "$tmpfile" "$mincard" "$input_msa_row" -v --gaps-as-symbols --pbwt -U $U -o mincard_pbwt_U${U}.eds + /usr/bin/time -f"$usrbintimeformat" -o "$tmpfile" "$mincard" "$input_msa_row" -v -U $U -o mincard_pbwt_U${U}.eds t3=$(<"$tmpfile") - /usr/bin/time -f"$usrbintimeformat" -o "$tmpfile" "$mincard" "$input_msa_col" -v --gaps-as-symbols --pbwt --column-major -U $U -o mincard_cm_pbwt_U${U}.eds + /usr/bin/time -f"$usrbintimeformat" -o "$tmpfile" "$mincard" "$input_msa_col" -v --column-major -U $U -o mincard_cm_pbwt_U${U}.eds t4=$(<"$tmpfile") if cmp -s mincard_trie_U${U}.eds mincard_cm_trie_U${U}.eds && diff --git a/src/mincard.cpp b/src/mincard.cpp index 3f74078..a4757e1 100644 --- a/src/mincard.cpp +++ b/src/mincard.cpp @@ -70,10 +70,13 @@ int main(int argc, char* argv[]) { CLI::Option *nsopt = app.add_flag("-n,--trivial-horizontal", no_segmentation, "Use trivial S^≡ segmentation (no segmentation)") ->excludes(Lopt)->excludes(Uopt)->excludes(tsopt); - bool gaps_as_symbols = false; - app.add_flag("--gaps-as-symbols", gaps_as_symbols, "In preprocessing the MSA, consider gaps '-' as normal symbols") - ->excludes(tsopt)->excludes(nsopt); + bool gaps = false; + app.add_flag("--gaps", gaps, "In preprocessing the MSA, remove gap symbols '-'") + ->excludes(tsopt)->excludes(nsopt); + bool naive_trie = false; + app.add_flag("--trie", naive_trie, "Compute meaningful extensions using naive Suffix Trie instead of positional Burrows-Wheeler Transform"); + bool preprocess = false; app.add_flag("--preprocess", preprocess, "Compute all meaningful extensions before segmenting") ->excludes(tsopt)->excludes(nsopt); @@ -84,9 +87,6 @@ int main(int argc, char* argv[]) { bool quiet = false; app.add_flag("-q,--quiet", quiet, "Print only cardinality and size"); - bool use_pbwt = false; - app.add_flag("--pbwt", use_pbwt, "Compute meaningful extensions using positional Burrows-Wheeler Transform"); - bool column_major = false; app.add_flag("--column-major", column_major, "Read msa in column-major format for faster column streaming"); @@ -95,6 +95,8 @@ int main(int argc, char* argv[]) { } catch (const CLI::ParseError &e) { return app.exit(e); } + bool gaps_as_symbols = !gaps; + bool use_pbwt = gaps_as_symbols && !naive_trie; // No U value specified by the user if (U == 0) { #if METRIC == CARDINALITY @@ -116,10 +118,6 @@ int main(int argc, char* argv[]) { cerr << "Upper and lower bounds are not compatible!" << endl; return 1; } - if(use_pbwt and !gaps_as_symbols){ - cerr << "pBWT only works with the gaps as symbols strategy! Add flag --gaps-as-symbols" << endl; - return 1; - } if (verbose and quiet) { cerr << "Quitet flag is ignored if verbose is set" << endl; }