From aad39bd847ea71c086f6bec29147e8bf1d978d4c Mon Sep 17 00:00:00 2001
From: Mikhail Karasikov <mike-karas@ya.ru>
Date: Sat, 28 Nov 2020 23:39:04 +0100
Subject: [PATCH 01/23] row-diff transform for relation counts

---
 .../src/annotation/annotation_converters.cpp  |  27 +++-
 .../src/annotation/annotation_converters.hpp  |   3 +-
 .../annotate_column_compressed.hpp            |   3 +-
 metagraph/src/annotation/row_diff_builder.cpp | 147 ++++++++++++++++--
 metagraph/src/annotation/row_diff_builder.hpp |   3 +-
 metagraph/src/cli/config/config.cpp           |   1 +
 metagraph/src/cli/transform_annotation.cpp    |   3 +-
 7 files changed, 163 insertions(+), 24 deletions(-)

diff --git a/metagraph/src/annotation/annotation_converters.cpp b/metagraph/src/annotation/annotation_converters.cpp
index c9d65ac9da..27e89446a6 100644
--- a/metagraph/src/annotation/annotation_converters.cpp
+++ b/metagraph/src/annotation/annotation_converters.cpp
@@ -1122,7 +1122,8 @@ void convert_to_row_diff(const std::vector<std::string> &files,
                          size_t mem_bytes,
                          uint32_t max_path_length,
                          std::filesystem::path dest_dir,
-                         bool optimize) {
+                         bool optimize,
+                         bool with_counts) {
     if (!files.size())
         return;
 
@@ -1142,9 +1143,24 @@ void convert_to_row_diff(const std::vector<std::string> &files,
         size_t mem_bytes_left = mem_bytes;
         std::vector<std::string> file_batch;
         for ( ; i < files.size(); ++i) {
-            // also include two buffers (fwd and back) for each column transformed
-            uint64_t file_size = std::filesystem::file_size(files[i])
-                                + ROW_DIFF_BUFFER_SIZE * sizeof(uint64_t) * 2;
+            // take annotation size
+            uint64_t file_size = std::filesystem::file_size(files[i]);
+            if (!with_counts) {
+                // also include two buffers (fwd and back) for each column transformed
+                file_size += ROW_DIFF_BUFFER_SIZE * sizeof(uint64_t) * 2;
+            } else {
+                // each buffer will store pairs (idx, count)
+                file_size += ROW_DIFF_BUFFER_SIZE * (sizeof(uint64_t) * 2) * 2;
+                // also add k-mer counts
+                const auto &counts_fname = files[i] + ".counts";
+                file_size += std::filesystem::file_size(counts_fname);
+                if (!std::filesystem::exists(counts_fname)) {
+                    logger->warn("Could not find counts for annotation {}, skipped",
+                                 counts_fname);
+                    continue;
+                }
+            }
+
             if (file_size > mem_bytes) {
                 logger->warn(
                         "Not enough memory to process {}, requires {} MB, skipped",
@@ -1180,7 +1196,8 @@ void convert_to_row_diff(const std::vector<std::string> &files,
 
         convert_batch_to_row_diff(
                 graph_fname, graph_fname + kRowDiffAnchorExt + (optimize ? "" : ".unopt"),
-                file_batch, dest_dir, row_reduction_fname, ROW_DIFF_BUFFER_SIZE);
+                file_batch, dest_dir, row_reduction_fname, ROW_DIFF_BUFFER_SIZE,
+                with_counts);
 
         logger->trace("Batch transformed in {} sec", timer.elapsed());
     }
diff --git a/metagraph/src/annotation/annotation_converters.hpp b/metagraph/src/annotation/annotation_converters.hpp
index d4ed5294a6..ee0d828c3f 100644
--- a/metagraph/src/annotation/annotation_converters.hpp
+++ b/metagraph/src/annotation/annotation_converters.hpp
@@ -108,7 +108,8 @@ void convert_to_row_diff(const std::vector<std::string> &files,
                          size_t mem_bytes,
                          uint32_t max_path_length,
                          std::filesystem::path dest_dir,
-                         bool optimize = false);
+                         bool optimize = false,
+                         bool with_counts = false);
 
 void convert_row_diff_to_col_compressed(const std::vector<std::string> &files,
                                         const std::string &outfbase);
diff --git a/metagraph/src/annotation/representation/column_compressed/annotate_column_compressed.hpp b/metagraph/src/annotation/representation/column_compressed/annotate_column_compressed.hpp
index db832a88a7..3b089d17ca 100644
--- a/metagraph/src/annotation/representation/column_compressed/annotate_column_compressed.hpp
+++ b/metagraph/src/annotation/representation/column_compressed/annotate_column_compressed.hpp
@@ -25,6 +25,7 @@ class ColumnCompressed : public MultiLabelEncoded<Label> {
   public:
     using Index = typename MultiLabelEncoded<Label>::Index;
     using VLabels = typename MultiLabelEncoded<Label>::VLabels;
+    using counts_container = sdsl::int_vector<>;
 
     ColumnCompressed(uint64_t num_rows = 0,
                      size_t num_columns_cached = 1);
@@ -118,7 +119,7 @@ class ColumnCompressed : public MultiLabelEncoded<Label> {
                               bitmap_builder*,
                               caches::LRUCachePolicy<size_t>> cached_columns_;
 
-    std::vector<sdsl::int_vector<>> relation_counts_;
+    std::vector<counts_container> relation_counts_;
 
     using MultiLabelEncoded<Label>::label_encoder_;
 };
diff --git a/metagraph/src/annotation/row_diff_builder.cpp b/metagraph/src/annotation/row_diff_builder.cpp
index 914cc0fd48..23d4845881 100644
--- a/metagraph/src/annotation/row_diff_builder.cpp
+++ b/metagraph/src/annotation/row_diff_builder.cpp
@@ -11,10 +11,14 @@
 #include "common/vectors/bit_vector_sd.hpp"
 #include "graph/annotated_dbg.hpp"
 
+
 constexpr uint64_t BLOCK_SIZE = 1 << 25;
 constexpr uint64_t ROW_REDUCTION_WIDTH = 32;
 constexpr uint32_t MAX_NUM_FILES_OPEN = 2000;
 
+const uint8_t kCountBits = 9;
+
+
 namespace mtg {
 namespace annot {
 
@@ -263,7 +267,42 @@ void traverse_anno_chunked(
     assert(pred_boundary_it == pred_boundary.end());
 }
 
+template <typename T = uint64_t>
+void convert_batch_to_row_diff(const std::string &pred_succ_fprefix,
+                               const std::string &anchors_fname,
+                               const std::vector<std::string> &source_files,
+                               const std::filesystem::path &dest_dir,
+                               const std::string &row_reduction_fname,
+                               uint64_t buf_size);
 
+void convert_batch_to_row_diff(const std::string &pred_succ_fprefix,
+                               const std::string &anchors_fname,
+                               const std::vector<std::string> &source_files,
+                               const std::filesystem::path &dest_dir,
+                               const std::string &row_reduction_fname,
+                               uint64_t buf_size,
+                               bool with_counts) {
+    if (with_counts) {
+        convert_batch_to_row_diff<std::pair<uint64_t, uint64_t>>(
+                pred_succ_fprefix,
+                anchors_fname,
+                source_files,
+                dest_dir,
+                row_reduction_fname,
+                buf_size);
+    } else {
+        convert_batch_to_row_diff<uint64_t>(
+                pred_succ_fprefix,
+                anchors_fname,
+                source_files,
+                dest_dir,
+                row_reduction_fname,
+                buf_size);
+    }
+}
+
+// 'T' is either a row index, or a pair (row index, relation count)
+template <typename T>
 void convert_batch_to_row_diff(const std::string &pred_succ_fprefix,
                                const std::string &anchors_fname,
                                const std::vector<std::string> &source_files,
@@ -273,6 +312,8 @@ void convert_batch_to_row_diff(const std::string &pred_succ_fprefix,
     if (source_files.empty())
         return;
 
+    constexpr bool with_counts = utils::is_pair_v<T>;
+
     uint32_t num_threads = get_num_threads();
 
     anchor_bv_type anchor;
@@ -282,6 +323,7 @@ void convert_batch_to_row_diff(const std::string &pred_succ_fprefix,
     }
 
     std::vector<annot::ColumnCompressed<>> sources(source_files.size());
+    std::vector<std::vector<annot::ColumnCompressed<>::counts_container>> counts(source_files.size());
 
     #pragma omp parallel for num_threads(num_threads)
     for (size_t i = 0; i < source_files.size(); ++i) {
@@ -293,6 +335,26 @@ void convert_batch_to_row_diff(const std::string &pred_succ_fprefix,
                           anchors_fname, source_files[i], anchor.size(), sources[i].num_objects());
             std::exit(1);
         }
+
+        if (!with_counts)
+            continue;
+
+        const auto &counts_fname = source_files[i] + ".counts";
+        std::ifstream instream(counts_fname, std::ios::binary);
+        if (!instream.good()) {
+            logger->error("Could not open file with relation counts {}", counts_fname);
+            exit(1);
+        }
+
+        counts[i].resize(sources[i].num_labels());
+        for (size_t j = 0; j < counts[i].size(); ++j) {
+            try {
+                counts[i][j].load(instream);
+            } catch (...) {
+                logger->error("Couldn't read relation counts from {}", counts_fname);
+                exit(1);
+            }
+        }
     }
     logger->trace("Done loading {} annotations", sources.size());
 
@@ -302,19 +364,19 @@ void convert_batch_to_row_diff(const std::string &pred_succ_fprefix,
 
     // stores the row indices that were set because of differences to incoming/outgoing
     // edges, for each of the sources, per chunk. set_rows_fwd is already sorted
-    std::vector<std::vector<std::vector<uint64_t>>> set_rows_bwd(sources.size());
-    std::vector<std::vector<std::vector<uint64_t>>> set_rows_fwd(sources.size());
+    std::vector<std::vector<std::vector<T>>> set_rows_bwd(sources.size());
+    std::vector<std::vector<std::vector<T>>> set_rows_fwd(sources.size());
     std::vector<std::vector<uint64_t>> row_diff_bits(sources.size());
     std::vector<std::vector<uint64_t>> num_chunks(sources.size());
 
     auto tmp_file = [&](size_t s, size_t j, size_t chunk) {
         return tmp_path/fmt::format("{}/col_{}_{}/chunk_{}", s / 100, s, j, chunk);
     };
-    auto dump_chunk_to_disk = [&](const std::vector<uint64_t> &v,
+    auto dump_chunk_to_disk = [&](const std::vector<T> &v,
                                   size_t s, size_t j, size_t chunk) {
         assert(std::is_sorted(v.begin(), v.end()) && "all bits in chunks must be sorted");
         std::ofstream f(tmp_file(s, j, chunk), std::ios::binary | std::ios::app);
-        f.write(reinterpret_cast<const char *>(v.data()), v.size() * sizeof(uint64_t));
+        f.write(reinterpret_cast<const char *>(v.data()), v.size() * sizeof(T));
         row_diff_bits[s][j] += v.size();
     };
 
@@ -359,6 +421,20 @@ void convert_batch_to_row_diff(const std::string &pred_succ_fprefix,
     // total number of set bits in the original rows
     std::vector<uint32_t> row_nbits_block;
 
+    // get bit at position |i| or its relation count
+    auto get_value = [&](const bit_vector &col,
+                         size_t s, size_t j, uint64_t i) {
+        if constexpr(!with_counts) {
+            return col[i];
+        } else {
+            if (uint64_t rk = col.conditional_rank1(i)) {
+                return (uint64_t)counts[s][j][rk - 1];
+            } else {
+                return (uint64_t)0;
+            }
+        }
+    };
+
     traverse_anno_chunked(
             "Compute diffs", anchor.size(), pred_succ_fprefix, sources,
             [&](uint64_t chunk_size) {
@@ -371,11 +447,24 @@ void convert_batch_to_row_diff(const std::string &pred_succ_fprefix,
 
                 __atomic_add_fetch(&row_nbits_block[chunk_idx], 1, __ATOMIC_RELAXED);
 
+                // get bits for these positions (or relation counts, hence uint64_t)
+                uint64_t curr_value;
+                if constexpr(with_counts) {
+                    curr_value = counts[source_idx][j][source_col.rank1(row_idx) - 1];
+                } else {
+                    curr_value = 1;
+                }
                 // add current bit if this node is an anchor
-                // or if the successor has zero bit
-                if (anchor[row_idx] || !source_col[succ]) {
+                // or if the successor has zero diff
+                bool is_anchor = anchor[row_idx];
+                uint64_t succ_value = is_anchor ? 0 : get_value(source_col, source_idx, j, succ);
+                if (is_anchor || succ_value != curr_value) {
                     auto &v = set_rows_fwd[source_idx][j];
-                    v.push_back(row_idx);
+                    if constexpr(with_counts) {
+                        v.emplace_back(row_idx, curr_value - succ_value);
+                    } else {
+                        v.push_back(row_idx);
+                    }
 
                     if (v.size() == v.capacity()) {
                         // dump chunk to disk
@@ -386,9 +475,17 @@ void convert_batch_to_row_diff(const std::string &pred_succ_fprefix,
 
                 // check non-anchor predecessor nodes and add them if they are zero
                 for (const uint64_t *pred_p = pred_begin; pred_p < pred_end; ++pred_p) {
-                    if (!source_col[*pred_p] && !anchor[*pred_p]) {
+                    if (anchor[*pred_p])
+                        continue;
+
+                    uint64_t pred_value = get_value(source_col, source_idx, j, *pred_p);
+                    if (!pred_value && curr_value) {
                         auto &v = set_rows_bwd[source_idx][j];
-                        v.push_back(*pred_p);
+                        if constexpr(with_counts) {
+                            v.emplace_back(*pred_p, pred_value - curr_value);
+                        } else {
+                            v.push_back(*pred_p);
+                        }
 
                         if (v.size() == v.capacity()) {
                             std::sort(v.begin(), v.end());
@@ -441,6 +538,8 @@ void convert_batch_to_row_diff(const std::string &pred_succ_fprefix,
 
     // free memory occupied by original columns
     sources.clear();
+    counts.clear();
+    counts.resize(label_encoders.size());
 
     std::vector<std::unique_ptr<RowDiffColumnAnnotator>> row_diff(label_encoders.size());
 
@@ -459,6 +558,10 @@ void convert_batch_to_row_diff(const std::string &pred_succ_fprefix,
         std::vector<std::unique_ptr<bit_vector>> columns(label_encoders[l_idx].size());
 
         for (size_t j = 0; j < label_encoders[l_idx].size(); ++j) {
+            if constexpr(with_counts) {
+                counts[l_idx].emplace_back(row_diff_bits[l_idx][j], 0, kCountBits);
+            }
+
             auto call_ones = [&](const std::function<void(uint64_t)> &call) {
                 std::vector<std::string> filenames;
                 for (uint32_t chunk = 0; chunk < num_chunks[l_idx][j]; ++chunk) {
@@ -490,20 +593,20 @@ void convert_batch_to_row_diff(const std::string &pred_succ_fprefix,
                         new_chunks.push_back(to_merge.at(0) + "_");
                         std::ofstream f(new_chunks.back(), std::ios::binary);
 
-                        std::vector<uint64_t> buf;
+                        std::vector<T> buf;
                         buf.reserve(1024 * 1024);
 
-                        common::merge_files<uint64_t>(to_merge, [&](uint64_t i) {
+                        common::merge_files<T>(to_merge, [&](T i) {
                             buf.push_back(i);
                             if (buf.size() == buf.capacity()) {
                                 f.write(reinterpret_cast<const char *>(buf.data()),
-                                        buf.size() * sizeof(decltype(buf.front())));
+                                        buf.size() * sizeof(T));
                                 buf.resize(0);
                             }
                         });
                         if (buf.size()) {
                             f.write(reinterpret_cast<const char *>(buf.data()),
-                                    buf.size() * sizeof(decltype(buf.front())));
+                                    buf.size() * sizeof(T));
                         }
 
                         // remove merged chunks
@@ -515,7 +618,14 @@ void convert_batch_to_row_diff(const std::string &pred_succ_fprefix,
                 }
 
                 //TODO: benchmark using Elias-Fano encoders + merger
-                common::merge_files<uint64_t>(filenames, call);
+                uint64_t rk = 0;
+                common::merge_files<T>(filenames, [&](T v) {
+                    call(utils::get_first(v));
+                    if constexpr(with_counts) {
+                        counts[l_idx][j][rk] = v.second;
+                    }
+                    rk++;
+                });
             };
             columns[j] = std::make_unique<bit_vector_sd>(call_ones, anchor.size(),
                                                          row_diff_bits[l_idx][j]);
@@ -528,10 +638,17 @@ void convert_batch_to_row_diff(const std::string &pred_succ_fprefix,
         auto fpath = dest_dir/std::filesystem::path(source_files[l_idx])
                                 .filename()
                                 .replace_extension()
-                                .replace_extension(RowDiffColumnAnnotator::kExtension);
+                                .replace_extension(row_diff[l_idx]->kExtension);
 
         row_diff[l_idx]->serialize(fpath);
 
+        if constexpr(with_counts) {
+            std::ofstream outstream(std::string(fpath) + ".counts", std::ios::binary);
+            for (size_t j = 0; j < label_encoders[l_idx].size(); ++j) {
+                counts[l_idx][j].serialize(outstream);
+            }
+        }
+
         logger->trace("Serialized {}", fpath);
     }
     logger->trace("Removing temp directory: {}", tmp_path);
diff --git a/metagraph/src/annotation/row_diff_builder.hpp b/metagraph/src/annotation/row_diff_builder.hpp
index 50fd4bd93e..2f2fa7671d 100644
--- a/metagraph/src/annotation/row_diff_builder.hpp
+++ b/metagraph/src/annotation/row_diff_builder.hpp
@@ -22,7 +22,8 @@ void convert_batch_to_row_diff(const std::string &pred_succ_fprefix,
                                const std::vector<std::string> &source_files,
                                const std::filesystem::path &dest_dir,
                                const std::string &row_reduction_fname,
-                               uint64_t buf_size);
+                               uint64_t buf_size,
+                               bool with_counts = false);
 
 void optimize_anchors_in_row_diff(const std::string &graph_fname,
                                   const std::filesystem::path &dest_dir,
diff --git a/metagraph/src/cli/config/config.cpp b/metagraph/src/cli/config/config.cpp
index c048be6719..b2a8409fc7 100644
--- a/metagraph/src/cli/config/config.cpp
+++ b/metagraph/src/cli/config/config.cpp
@@ -1033,6 +1033,7 @@ void Config::print_usage(const std::string &prog_name, IdentityType identity) {
             fprintf(stderr, "\t   --parallel-nodes [INT] \tnumber of nodes processed in parallel in brwt tree [n_threads]\n");
             fprintf(stderr, "\t   --max-path-length [INT] \tmaximum path length in row_diff annotation [50]\n");
             fprintf(stderr, "\t   --optimize \t\t\toptimize anchors in row_diff annotation [off]\n");
+            fprintf(stderr, "\t   --count-kmers \t\tadd k-mer counts to the row_diff annotation [off]\n");
         } break;
         case RELAX_BRWT: {
             fprintf(stderr, "Usage: %s relax_brwt -o <annotation-basename> [options] ANNOTATOR\n\n", prog_name.c_str());
diff --git a/metagraph/src/cli/transform_annotation.cpp b/metagraph/src/cli/transform_annotation.cpp
index 9524a353df..a0e372d3f2 100644
--- a/metagraph/src/cli/transform_annotation.cpp
+++ b/metagraph/src/cli/transform_annotation.cpp
@@ -390,7 +390,8 @@ int transform_annotation(Config *config) {
             case Config::RowDiff: {
                 auto out_dir = std::filesystem::path(config->outfbase).remove_filename();
                 convert_to_row_diff(files, config->infbase, config->memory_available * 1e9,
-                                    config->max_path_length, out_dir, config->optimize);
+                                    config->max_path_length, out_dir, config->optimize,
+                                    config->count_kmers);
                 break;
             }
             case Config::RowCompressed: {

From 7e432e4a53a1d5d55a57fd0065425e1b5ea01837 Mon Sep 17 00:00:00 2001
From: Mikhail Karasikov <mike-karas@ya.ru>
Date: Fri, 16 Apr 2021 18:48:27 +0200
Subject: [PATCH 02/23] row-diff transform for counts

---
 .../src/annotation/annotation_converters.cpp  |  16 +-
 .../src/annotation/annotation_converters.hpp  |   4 +-
 metagraph/src/annotation/row_diff_builder.cpp | 158 +++++++++++++++---
 metagraph/src/annotation/row_diff_builder.hpp |   3 +-
 metagraph/src/cli/config/config.cpp           |   1 +
 metagraph/src/cli/transform_annotation.cpp    |   2 +-
 .../src/common/elias_fano/elias_fano.cpp      |  39 ++++-
 .../src/common/elias_fano/elias_fano.hpp      |   4 +
 8 files changed, 197 insertions(+), 30 deletions(-)

diff --git a/metagraph/src/annotation/annotation_converters.cpp b/metagraph/src/annotation/annotation_converters.cpp
index 2426bbfe2c..8bd86d2edf 100644
--- a/metagraph/src/annotation/annotation_converters.cpp
+++ b/metagraph/src/annotation/annotation_converters.cpp
@@ -1146,7 +1146,8 @@ void convert_to_row_diff(const std::vector<std::string> &files,
                          fs::path out_dir,
                          fs::path swap_dir,
                          RowDiffStage construction_stage,
-                         fs::path count_vector_fname) {
+                         fs::path count_vector_fname,
+                         bool with_counts) {
     if (out_dir.empty())
         out_dir = "./";
 
@@ -1184,6 +1185,16 @@ void convert_to_row_diff(const std::vector<std::string> &files,
         for ( ; i < files.size(); ++i) {
             // also add some space for buffers for each column
             uint64_t file_size = fs::file_size(files[i]) + ROW_DIFF_BUFFER_BYTES;
+            if (with_counts) {
+                // also add k-mer counts
+                const auto &counts_fname = files[i] + ".counts";
+                file_size += fs::file_size(counts_fname);
+                if (!fs::exists(counts_fname)) {
+                    logger->warn("Could not find counts for annotation {}, skipped",
+                                 counts_fname);
+                    continue;
+                }
+            }
             if (file_size > mem_bytes) {
                 logger->warn("Not enough memory to process {}, requires {} MB, skipped",
                              files[i], file_size / 1e6);
@@ -1214,7 +1225,8 @@ void convert_to_row_diff(const std::vector<std::string> &files,
         } else {
             convert_batch_to_row_diff(graph_fname,
                     file_batch, out_dir, swap_dir, count_vector_fname, ROW_DIFF_BUFFER_BYTES,
-                    construction_stage == RowDiffStage::COMPUTE_REDUCTION);
+                    construction_stage == RowDiffStage::COMPUTE_REDUCTION,
+                    with_counts);
         }
 
         logger->trace("Batch processed in {} sec", timer.elapsed());
diff --git a/metagraph/src/annotation/annotation_converters.hpp b/metagraph/src/annotation/annotation_converters.hpp
index 41221a9e42..b670cce74f 100644
--- a/metagraph/src/annotation/annotation_converters.hpp
+++ b/metagraph/src/annotation/annotation_converters.hpp
@@ -104,6 +104,7 @@ void convert_to_row_annotator(const ColumnCompressed<Label> &annotator,
  * @param out_dir directory where the transformed columns will be dumped. Filenames are
  * kept, extension is changed from 'column.annodbg' to 'row_diff.annodbg'
  * @param swap_dir directory for temporary files
+ * @param with_counts row-diff transform for k-mer counts
  */
 enum class RowDiffStage { COUNT_LABELS = 0, COMPUTE_REDUCTION, CONVERT };
 void convert_to_row_diff(const std::vector<std::string> &files,
@@ -113,7 +114,8 @@ void convert_to_row_diff(const std::vector<std::string> &files,
                          std::filesystem::path out_dir,
                          std::filesystem::path swap_dir,
                          RowDiffStage construction_stage,
-                         std::filesystem::path count_vector_fname = "");
+                         std::filesystem::path count_vector_fname = "",
+                         bool with_counts = false);
 
 void convert_row_diff_to_col_compressed(const std::vector<std::string> &files,
                                         const std::string &outfbase);
diff --git a/metagraph/src/annotation/row_diff_builder.cpp b/metagraph/src/annotation/row_diff_builder.cpp
index d1ca53cf83..95fb0dd25d 100644
--- a/metagraph/src/annotation/row_diff_builder.cpp
+++ b/metagraph/src/annotation/row_diff_builder.cpp
@@ -628,6 +628,46 @@ void traverse_anno_chunked(
 }
 
 
+template <typename T = uint64_t>
+void convert_batch_to_row_diff(const std::string &pred_succ_fprefix,
+                               const std::vector<std::string> &source_files,
+                               const fs::path &col_out_dir,
+                               const fs::path &swap_dir,
+                               const std::string &row_reduction_fname,
+                               uint64_t buf_size_bytes,
+                               bool compute_row_reduction);
+
+void convert_batch_to_row_diff(const std::string &pred_succ_fprefix,
+                               const std::vector<std::string> &source_files,
+                               const fs::path &col_out_dir,
+                               const fs::path &swap_dir,
+                               const std::string &row_reduction_fname,
+                               uint64_t buf_size_bytes,
+                               bool compute_row_reduction,
+                               bool with_counts) {
+    if (with_counts) {
+        convert_batch_to_row_diff<std::pair<uint64_t, uint64_t>>(
+                pred_succ_fprefix,
+                source_files,
+                col_out_dir,
+                swap_dir,
+                row_reduction_fname,
+                buf_size_bytes,
+                compute_row_reduction);
+    } else {
+        convert_batch_to_row_diff<uint64_t>(
+                pred_succ_fprefix,
+                source_files,
+                col_out_dir,
+                swap_dir,
+                row_reduction_fname,
+                buf_size_bytes,
+                compute_row_reduction);
+    }
+}
+
+// 'T' is either a row index, or a pair (row index, relation count)
+template <typename T>
 void convert_batch_to_row_diff(const std::string &pred_succ_fprefix,
                                const std::vector<std::string> &source_files,
                                const fs::path &col_out_dir,
@@ -647,6 +687,33 @@ void convert_batch_to_row_diff(const std::string &pred_succ_fprefix,
         return;
     }
 
+    constexpr bool with_counts = utils::is_pair_v<T>;
+    std::vector<std::vector<sdsl::int_vector<>>> counts;
+    if (with_counts) {
+        #pragma omp parallel for num_threads(get_num_threads())
+        for (size_t i = 0; i < source_files.size(); ++i) {
+            const auto &counts_fname = utils::remove_suffix(source_files[i],
+                                                            ColumnCompressed<>::kExtension)
+                                            + ColumnCompressed<>::kCountExtension;
+            std::ifstream instream(counts_fname, std::ios::binary);
+            if (!instream.good()) {
+                logger->error("Could not open file with relation counts {}", counts_fname);
+                exit(1);
+            }
+
+            counts[i].resize(sources[i].num_labels());
+            for (size_t j = 0; j < counts[i].size(); ++j) {
+                try {
+                    counts[i][j].load(instream);
+                } catch (...) {
+                    logger->error("Couldn't read relation counts from {}", counts_fname);
+                    exit(1);
+                }
+            }
+        }
+        logger->trace("Done loading relation counts", counts.size());
+    }
+
     anchor_bv_type anchor;
     if (!compute_row_reduction) {
         const std::string anchors_fname = pred_succ_fprefix + kRowDiffAnchorExt;
@@ -667,26 +734,26 @@ void convert_batch_to_row_diff(const std::string &pred_succ_fprefix,
 
     // stores the row indices that were set because of differences to incoming/outgoing
     // edges, for each of the sources, per chunk. set_rows_fwd is already sorted
-    std::vector<std::vector<std::vector<uint64_t>>> set_rows_bwd(sources.size());
-    std::vector<std::vector<std::vector<uint64_t>>> set_rows_fwd(sources.size());
+    std::vector<std::vector<std::vector<T>>> set_rows_bwd(sources.size());
+    std::vector<std::vector<std::vector<T>>> set_rows_fwd(sources.size());
     std::vector<std::vector<uint64_t>> row_diff_bits(sources.size());
     std::vector<std::vector<uint64_t>> num_chunks(sources.size());
 
     auto tmp_file = [&](size_t s, size_t j, size_t chunk) {
         return tmp_path/fmt::format("{}/col_{}_{}/chunk_{}", s / 100, s, j, chunk);
     };
-    auto dump_chunk_to_disk = [&](const std::vector<uint64_t> &v,
+    auto dump_chunk_to_disk = [&](const std::vector<T> &v,
                                   size_t s, size_t j, size_t chunk) {
         assert(std::is_sorted(v.begin(), v.end()) && "all bits in chunks must be sorted");
-        Encoder<uint64_t>::append_block(v, tmp_file(s, j, chunk));
+        Encoder<T>::append_block(v, tmp_file(s, j, chunk));
         row_diff_bits[s][j] += v.size();
     };
 
     // In the first stage, only one buffer is created per column (`bwd`).
     // In the last stage, two buffers (`fwd` and `bwd`) are created per column.
     const uint64_t buf_size = compute_row_reduction
-                                ? buf_size_bytes / sizeof(uint64_t)
-                                : buf_size_bytes / sizeof(uint64_t) / 2;
+                                ? buf_size_bytes / sizeof(T)
+                                : buf_size_bytes / sizeof(T) / 2;
 
     #pragma omp parallel for num_threads(num_threads)
     for (size_t s = 0; s < sources.size(); ++s) {
@@ -740,6 +807,20 @@ void convert_batch_to_row_diff(const std::string &pred_succ_fprefix,
     // buffer for writing previous block while populating next row_nbits_block
     std::vector<uint32_t> row_nbits_block_other;
 
+    // get bit at position |i| or its relation count
+    auto get_value = [&](const bit_vector &col,
+                         size_t s, size_t j, uint64_t i) -> uint64_t {
+        if constexpr(!with_counts) {
+            return col[i];
+        } else {
+            if (uint64_t rk = col.conditional_rank1(i)) {
+                return counts[s][j][rk - 1];
+            } else {
+                return 0;
+            }
+        }
+    };
+
     traverse_anno_chunked(
             num_rows, pred_succ_fprefix, sources,
             [&](uint64_t chunk_size) {
@@ -751,18 +832,32 @@ void convert_batch_to_row_diff(const std::string &pred_succ_fprefix,
                     const uint64_t *succ,
                     const uint64_t *pred_begin, const uint64_t *pred_end) {
 
+                // get bits for these positions (or relation counts, hence uint64_t)
+                uint64_t curr_value;
+                if constexpr(with_counts) {
+                    curr_value = counts[source_idx][j][source_col.rank1(row_idx) - 1];
+                } else {
+                    curr_value = 1;
+                }
+
                 if (compute_row_reduction) {
-                    if (succ && source_col[*succ]) {
-                        // reduction (no bit in row-diff)
+                    if (succ && curr_value == get_value(source_col, source_idx, j, *succ)) {
+                        // reduction (zero diff)
                         __atomic_add_fetch(&row_nbits_block[chunk_idx], 1, __ATOMIC_RELAXED);
                     }
                 } else {
                     // add current bit if this node is an anchor
-                    // or if the successor has zero bit
-                    if (anchor[row_idx] || !source_col[*succ]) {
+                    // or if the successor has zero diff
+                    bool is_anchor = anchor[row_idx];
+                    uint64_t succ_value = is_anchor ? 0 : get_value(source_col, source_idx, j, *succ);
+                    if (is_anchor || succ_value != curr_value) {
                         // no reduction, we must keep the bit
                         auto &v = set_rows_fwd[source_idx][j];
-                        v.push_back(row_idx);
+                        if constexpr(with_counts) {
+                            v.emplace_back(row_idx, curr_value - succ_value);
+                        } else {
+                            v.push_back(row_idx);
+                        }
 
                         if (v.size() == v.capacity()) {
                             // dump chunk to disk
@@ -774,9 +869,13 @@ void convert_batch_to_row_diff(const std::string &pred_succ_fprefix,
 
                 // check non-anchor predecessor nodes and add them if they are zero
                 for (const uint64_t *pred_p = pred_begin; pred_p < pred_end; ++pred_p) {
-                    if (!source_col[*pred_p] && (compute_row_reduction || !anchor[*pred_p])) {
+                    if (curr_value && !source_col[*pred_p] && (compute_row_reduction || !anchor[*pred_p])) {
                         auto &v = set_rows_bwd[source_idx][j];
-                        v.push_back(*pred_p);
+                        if constexpr(with_counts) {
+                            v.emplace_back(*pred_p, -curr_value);
+                        } else {
+                            v.push_back(*pred_p);
+                        }
 
                         if (v.size() == v.capacity()) {
                             std::sort(v.begin(), v.end());
@@ -853,6 +952,12 @@ void convert_batch_to_row_diff(const std::string &pred_succ_fprefix,
         std::vector<std::unique_ptr<bit_vector>> columns(label_encoders[l_idx].size());
 
         for (size_t j = 0; j < label_encoders[l_idx].size(); ++j) {
+            if constexpr(with_counts) {
+                // diff counts may be negative, hence we need wider integers
+                counts[l_idx][j] = sdsl::int_vector<>(row_diff_bits[l_idx][j], 0,
+                                                      counts[l_idx][j].width() * 2);
+            }
+
             auto call_ones = [&](const std::function<void(uint64_t)> &call) {
                 std::vector<std::string> filenames;
                 // skip chunk with fwd bits which have already been counted if stage 1
@@ -884,24 +989,31 @@ void convert_batch_to_row_diff(const std::string &pred_succ_fprefix,
                         assert(to_merge.size() < files_open_per_thread);
 
                         new_chunks.push_back(to_merge.at(0) + "_");
-                        std::vector<uint64_t> buf;
+                        std::vector<T> buf;
                         buf.reserve(BUFFER_SIZE);
 
-                        elias_fano::merge_files<uint64_t>(to_merge, [&](uint64_t i) {
+                        elias_fano::merge_files<T>(to_merge, [&](T i) {
                             buf.push_back(i);
                             if (buf.size() == buf.capacity()) {
-                                Encoder<uint64_t>::append_block(buf, new_chunks.back());
+                                Encoder<T>::append_block(buf, new_chunks.back());
                                 buf.resize(0);
                             }
                         });
                         if (buf.size()) {
-                            Encoder<uint64_t>::append_block(buf, new_chunks.back());
+                            Encoder<T>::append_block(buf, new_chunks.back());
                         }
                     }
                     filenames.swap(new_chunks);
                 }
 
-                elias_fano::merge_files<uint64_t>(filenames, call);
+                uint64_t rk = 0;
+                elias_fano::merge_files<T>(filenames, [&](T v) {
+                    call(utils::get_first(v));
+                    if constexpr(with_counts) {
+                        counts[l_idx][j][rk] = v.second;
+                    }
+                    rk++;
+                });
             };
             columns[j] = std::make_unique<bit_vector_sd>(call_ones, num_rows,
                                                          row_diff_bits[l_idx][j]);
@@ -915,10 +1027,18 @@ void convert_batch_to_row_diff(const std::string &pred_succ_fprefix,
             auto fpath = col_out_dir/fs::path(source_files[l_idx])
                                     .filename()
                                     .replace_extension()
-                                    .replace_extension(RowDiffColumnAnnotator::kExtension);
+                                    .replace_extension(row_diff[l_idx]->kExtension);
 
             row_diff[l_idx]->serialize(fpath);
 
+            // TODO: IntRowDiff::kExtension
+            if constexpr(with_counts) {
+                std::ofstream outstream(std::string(fpath) + ".counts", std::ios::binary);
+                for (size_t j = 0; j < label_encoders[l_idx].size(); ++j) {
+                    counts[l_idx][j].serialize(outstream);
+                }
+            }
+
             logger->trace("Serialized {}", fpath);
         }
     }
diff --git a/metagraph/src/annotation/row_diff_builder.hpp b/metagraph/src/annotation/row_diff_builder.hpp
index 2e6aeb887b..fcbf046f38 100644
--- a/metagraph/src/annotation/row_diff_builder.hpp
+++ b/metagraph/src/annotation/row_diff_builder.hpp
@@ -34,7 +34,8 @@ void convert_batch_to_row_diff(const std::string &pred_succ_fprefix,
                                const std::filesystem::path &swap_dir,
                                const std::string &row_reduction_fname,
                                uint64_t buf_size_bytes,
-                               bool compute_row_reduction = true);
+                               bool compute_row_reduction = true,
+                               bool with_counts = false);
 
 } // namespace annot
 } // namespace mtg
diff --git a/metagraph/src/cli/config/config.cpp b/metagraph/src/cli/config/config.cpp
index b60978cfd1..2ba62ebcf0 100644
--- a/metagraph/src/cli/config/config.cpp
+++ b/metagraph/src/cli/config/config.cpp
@@ -1113,6 +1113,7 @@ void Config::print_usage(const std::string &prog_name, IdentityType identity) {
             fprintf(stderr, "\t   --row-diff-stage [0|1|2] \tstage of the row_diff construction [0]\n");
             fprintf(stderr, "\t   --max-path-length [INT] \tmaximum path length in row_diff annotation [100]\n");
             fprintf(stderr, "\t-i --infile-base [STR] \t\tgraph for generating succ/pred/anchors (for row_diff types) []\n");
+            fprintf(stderr, "\t   --count-kmers \t\tadd k-mer counts to the row_diff annotation [off]\n");
             fprintf(stderr, "\n");
             fprintf(stderr, "\t   --parallel-nodes [INT] \tnumber of nodes processed in parallel in brwt tree [n_threads]\n");
             fprintf(stderr, "\n");
diff --git a/metagraph/src/cli/transform_annotation.cpp b/metagraph/src/cli/transform_annotation.cpp
index 617774a948..8142455100 100644
--- a/metagraph/src/cli/transform_annotation.cpp
+++ b/metagraph/src/cli/transform_annotation.cpp
@@ -616,7 +616,7 @@ int transform_annotation(Config *config) {
                 convert_to_row_diff(files, config->infbase, config->memory_available * 1e9,
                                     config->max_path_length, out_dir, config->tmp_dir,
                                     static_cast<RowDiffStage>(config->row_diff_stage),
-                                    config->outfbase);
+                                    config->outfbase, config->count_kmers);
                 break;
             }
             case Config::RowCompressed: {
diff --git a/metagraph/src/common/elias_fano/elias_fano.cpp b/metagraph/src/common/elias_fano/elias_fano.cpp
index f9f2663fc3..794f34760a 100644
--- a/metagraph/src/common/elias_fano/elias_fano.cpp
+++ b/metagraph/src/common/elias_fano/elias_fano.cpp
@@ -406,11 +406,11 @@ void EliasFanoEncoder<T>::init(size_t size, T max_value) {
         upper_.resize(num_upper_bytes_ + sizeof(T) - 1, 0);
     }
 
-    sink_->write(reinterpret_cast<char *>(&size), sizeof(size_t));
-    sink_->write(reinterpret_cast<char *>(&offset_), sizeof(T));
-    sink_->write(reinterpret_cast<char *>(&num_lower_bits_), 1);
-    sink_->write(reinterpret_cast<char *>(&num_lower_bytes_), sizeof(size_t));
-    sink_->write(reinterpret_cast<char *>(&num_upper_bytes_), sizeof(size_t));
+    sink_->write(reinterpret_cast<const char *>(&size), sizeof(size_t));
+    sink_->write(reinterpret_cast<const char *>(&offset_), sizeof(T));
+    sink_->write(reinterpret_cast<const char *>(&num_lower_bits_), 1);
+    sink_->write(reinterpret_cast<const char *>(&num_lower_bytes_), sizeof(size_t));
+    sink_->write(reinterpret_cast<const char *>(&num_upper_bytes_), sizeof(size_t));
 }
 
 template <typename T>
@@ -668,6 +668,31 @@ EliasFanoEncoderBuffered<std::pair<T, C>>::~EliasFanoEncoderBuffered() {
     assert(!sink_second_.is_open());
 }
 
+/** Append sorted array #data to EF-coded #out_fname */
+template <typename T, typename C>
+size_t EliasFanoEncoderBuffered<std::pair<T, C>>
+::append_block(const std::vector<std::pair<T, C>> &data,
+               const std::string &file_name) {
+    std::vector<T> ts;
+    std::vector<C> cs;
+    ts.reserve(data.size());
+    cs.reserve(data.size());
+    for (const auto &[v, c] : data) {
+        ts.push_back(v);
+        cs.push_back(c);
+    }
+    size_t total_size = EliasFanoEncoderBuffered<T>::append_block(ts, file_name);
+
+    std::ofstream sink_second(file_name + ".count", std::ios::binary|std::ios::app);
+    if (!sink_second) {
+        logger->error("Unable to open {} for writing", file_name + ".count");
+        std::exit(EXIT_FAILURE);
+    }
+    sink_second.write(reinterpret_cast<const char *>(cs.data()), cs.size() * sizeof(C));
+
+    return total_size + cs.size() * sizeof(C);
+}
+
 template <typename T, typename C>
 size_t EliasFanoEncoderBuffered<std::pair<T, C>>::finish() {
     encode_chunk();
@@ -677,7 +702,7 @@ size_t EliasFanoEncoderBuffered<std::pair<T, C>>::finish() {
 
 template <typename T, typename C>
 void EliasFanoEncoderBuffered<std::pair<T, C>>::encode_chunk() {
-    sink_second_.write(reinterpret_cast<char *>(buffer_second_.data()),
+    sink_second_.write(reinterpret_cast<const char *>(buffer_second_.data()),
                        buffer_second_.size() * sizeof(C));
     buffer_second_.resize(0);
 }
@@ -689,6 +714,7 @@ template class EliasFanoDecoder<sdsl::uint256_t>;
 template class EliasFanoDecoder<std::pair<uint64_t, uint8_t>>;
 template class EliasFanoDecoder<std::pair<uint64_t, uint16_t>>;
 template class EliasFanoDecoder<std::pair<uint64_t, uint32_t>>;
+template class EliasFanoDecoder<std::pair<uint64_t, uint64_t>>;
 template class EliasFanoDecoder<std::pair<sdsl::uint128_t, uint8_t>>;
 template class EliasFanoDecoder<std::pair<sdsl::uint128_t, uint16_t>>;
 template class EliasFanoDecoder<std::pair<sdsl::uint128_t, uint32_t>>;
@@ -702,6 +728,7 @@ template class EliasFanoEncoderBuffered<sdsl::uint256_t>;
 template class EliasFanoEncoderBuffered<std::pair<uint64_t, uint8_t>>;
 template class EliasFanoEncoderBuffered<std::pair<uint64_t, uint16_t>>;
 template class EliasFanoEncoderBuffered<std::pair<uint64_t, uint32_t>>;
+template class EliasFanoEncoderBuffered<std::pair<uint64_t, uint64_t>>;
 template class EliasFanoEncoderBuffered<std::pair<sdsl::uint128_t, uint8_t>>;
 template class EliasFanoEncoderBuffered<std::pair<sdsl::uint128_t, uint16_t>>;
 template class EliasFanoEncoderBuffered<std::pair<sdsl::uint128_t, uint32_t>>;
diff --git a/metagraph/src/common/elias_fano/elias_fano.hpp b/metagraph/src/common/elias_fano/elias_fano.hpp
index e504958859..77a1249da6 100644
--- a/metagraph/src/common/elias_fano/elias_fano.hpp
+++ b/metagraph/src/common/elias_fano/elias_fano.hpp
@@ -235,6 +235,10 @@ class EliasFanoEncoderBuffered<std::pair<T, C>> {
 
     size_t finish();
 
+    /** Append sorted array #data to EF-coded #out_fname */
+    static size_t append_block(const std::vector<std::pair<T, C>> &data,
+                               const std::string &file_name);
+
   private:
     void encode_chunk();
 

From 9b9d551ff892a9b2a6bff8d07efb776798ad0cd4 Mon Sep 17 00:00:00 2001
From: Mikhail Karasikov <mike-karas@ya.ru>
Date: Sat, 17 Apr 2021 13:23:17 +0200
Subject: [PATCH 03/23] minor: removed constexpr

---
 metagraph/src/annotation/row_diff_builder.cpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/metagraph/src/annotation/row_diff_builder.cpp b/metagraph/src/annotation/row_diff_builder.cpp
index 95fb0dd25d..d5627ef943 100644
--- a/metagraph/src/annotation/row_diff_builder.cpp
+++ b/metagraph/src/annotation/row_diff_builder.cpp
@@ -810,7 +810,7 @@ void convert_batch_to_row_diff(const std::string &pred_succ_fprefix,
     // get bit at position |i| or its relation count
     auto get_value = [&](const bit_vector &col,
                          size_t s, size_t j, uint64_t i) -> uint64_t {
-        if constexpr(!with_counts) {
+        if (!with_counts) {
             return col[i];
         } else {
             if (uint64_t rk = col.conditional_rank1(i)) {

From 832f7d51ed0676678b20e1cff77cf66fe5e2019b Mon Sep 17 00:00:00 2001
From: Mikhail Karasikov <mike-karas@ya.ru>
Date: Sun, 18 Apr 2021 13:56:36 +0200
Subject: [PATCH 04/23] minor

---
 .../src/annotation/annotation_converters.cpp  |  4 +++-
 metagraph/src/annotation/row_diff_builder.cpp | 20 +++++--------------
 2 files changed, 8 insertions(+), 16 deletions(-)

diff --git a/metagraph/src/annotation/annotation_converters.cpp b/metagraph/src/annotation/annotation_converters.cpp
index 8bd86d2edf..b082d25f26 100644
--- a/metagraph/src/annotation/annotation_converters.cpp
+++ b/metagraph/src/annotation/annotation_converters.cpp
@@ -1187,7 +1187,9 @@ void convert_to_row_diff(const std::vector<std::string> &files,
             uint64_t file_size = fs::file_size(files[i]) + ROW_DIFF_BUFFER_BYTES;
             if (with_counts) {
                 // also add k-mer counts
-                const auto &counts_fname = files[i] + ".counts";
+                const auto &counts_fname
+                    = utils::remove_suffix(files[i], ColumnCompressed<>::kExtension)
+                                                + ColumnCompressed<>::kCountExtension;
                 file_size += fs::file_size(counts_fname);
                 if (!fs::exists(counts_fname)) {
                     logger->warn("Could not find counts for annotation {}, skipped",
diff --git a/metagraph/src/annotation/row_diff_builder.cpp b/metagraph/src/annotation/row_diff_builder.cpp
index d5627ef943..05cd62e097 100644
--- a/metagraph/src/annotation/row_diff_builder.cpp
+++ b/metagraph/src/annotation/row_diff_builder.cpp
@@ -647,22 +647,12 @@ void convert_batch_to_row_diff(const std::string &pred_succ_fprefix,
                                bool with_counts) {
     if (with_counts) {
         convert_batch_to_row_diff<std::pair<uint64_t, uint64_t>>(
-                pred_succ_fprefix,
-                source_files,
-                col_out_dir,
-                swap_dir,
-                row_reduction_fname,
-                buf_size_bytes,
-                compute_row_reduction);
+                pred_succ_fprefix, source_files, col_out_dir, swap_dir,
+                row_reduction_fname, buf_size_bytes, compute_row_reduction);
     } else {
         convert_batch_to_row_diff<uint64_t>(
-                pred_succ_fprefix,
-                source_files,
-                col_out_dir,
-                swap_dir,
-                row_reduction_fname,
-                buf_size_bytes,
-                compute_row_reduction);
+                pred_succ_fprefix, source_files, col_out_dir, swap_dir,
+                row_reduction_fname, buf_size_bytes, compute_row_reduction);
     }
 }
 
@@ -955,7 +945,7 @@ void convert_batch_to_row_diff(const std::string &pred_succ_fprefix,
             if constexpr(with_counts) {
                 // diff counts may be negative, hence we need wider integers
                 counts[l_idx][j] = sdsl::int_vector<>(row_diff_bits[l_idx][j], 0,
-                                                      counts[l_idx][j].width() * 2);
+                                                      std::min(counts[l_idx][j].width() * 2, 64));
             }
 
             auto call_ones = [&](const std::function<void(uint64_t)> &call) {

From 938e709ebbeaf6fa5c02bd85e09ce3d4e8d23add Mon Sep 17 00:00:00 2001
From: Mikhail Karasikov <mike-karas@ya.ru>
Date: Tue, 20 Apr 2021 10:27:06 +0200
Subject: [PATCH 05/23] cleanup

---
 .../annotation/binary_matrix/row_diff/row_diff.cpp | 14 --------------
 .../annotation/binary_matrix/row_diff/row_diff.hpp | 14 +++++---------
 2 files changed, 5 insertions(+), 23 deletions(-)

diff --git a/metagraph/src/annotation/binary_matrix/row_diff/row_diff.cpp b/metagraph/src/annotation/binary_matrix/row_diff/row_diff.cpp
index fed992595a..4dda7cab78 100644
--- a/metagraph/src/annotation/binary_matrix/row_diff/row_diff.cpp
+++ b/metagraph/src/annotation/binary_matrix/row_diff/row_diff.cpp
@@ -11,20 +11,6 @@ namespace mtg {
 namespace annot {
 namespace binmat {
 
-template <class BaseMatrix>
-void RowDiff<BaseMatrix>::serialize(const std::string &filename) const {
-    std::ofstream f(filename, ios::binary);
-    serialize(f);
-    f.close();
-}
-
-template <class BaseMatrix>
-bool RowDiff<BaseMatrix>::load(const std::string &filename) {
-    std::ifstream f(filename, ios::binary);
-    bool result = load(f);
-    return result;
-}
-
 template <class BaseMatrix>
 void RowDiff<BaseMatrix>::load_anchor(const std::string &filename) {
     if (!std::filesystem::exists(filename)) {
diff --git a/metagraph/src/annotation/binary_matrix/row_diff/row_diff.hpp b/metagraph/src/annotation/binary_matrix/row_diff/row_diff.hpp
index 49d4b4a3c9..7dce0c4456 100644
--- a/metagraph/src/annotation/binary_matrix/row_diff/row_diff.hpp
+++ b/metagraph/src/annotation/binary_matrix/row_diff/row_diff.hpp
@@ -27,7 +27,7 @@ const std::string kRowDiffForkSuccExt = ".rd_succ";
 const size_t RD_PATH_RESERVE_SIZE = 2;
 
 
-class IRowDiff : public BinaryMatrix {
+class IRowDiff {
   public:
     virtual ~IRowDiff() {}
 };
@@ -54,7 +54,7 @@ class IRowDiff : public BinaryMatrix {
 // does not instantiate the virtual methods in this class, so I had to move definitions
 // to the header (gcc works fine)
 template <class BaseMatrix>
-class RowDiff : public IRowDiff {
+class RowDiff : public IRowDiff, public BinaryMatrix {
   public:
     using anchor_bv_type = bit_vector_small;
     using fork_succ_bv_type = bit_vector_small;
@@ -90,9 +90,6 @@ class RowDiff : public IRowDiff {
     bool load(std::istream &f) override;
     void serialize(std::ostream &f) const override;
 
-    void serialize(const std::string &filename) const;
-    bool load(const std::string &filename);
-
     void load_fork_succ(const std::string &filename);
     void load_anchor(const std::string &filename);
     const anchor_bv_type& anchor() const { return anchor_; }
@@ -100,8 +97,6 @@ class RowDiff : public IRowDiff {
     const BaseMatrix& diffs() const { return diffs_; }
     BaseMatrix& diffs() { return diffs_; }
 
-    Vector<uint64_t> get_diff(uint64_t node_id) const { return diffs_.get_row(node_id); }
-
   private:
     static void add_diff(const Vector<uint64_t> &diff, Vector<uint64_t> *row);
 
@@ -131,6 +126,7 @@ std::vector<BinaryMatrix::Row> RowDiff<BaseMatrix>::get_column(Column column) co
     assert(anchor_.size() == diffs_.num_rows() && "anchors must be loaded");
     assert(!fork_succ_.size() || fork_succ_.size() == graph_->get_boss().get_last().size());
 
+    // TODO: implement a more efficient algorithm
     std::vector<Row> result;
     for (Row row = 0; row < num_rows(); ++row) {
         if (get(row, column))
@@ -144,7 +140,7 @@ BinaryMatrix::SetBitPositions RowDiff<BaseMatrix>::get_row(Row row) const {
     assert(anchor_.size() == diffs_.num_rows() && "anchors must be loaded");
     assert(!fork_succ_.size() || fork_succ_.size() == graph_->get_boss().get_last().size());
 
-    Vector<uint64_t> result = get_diff(row);
+    Vector<uint64_t> result = diffs_.get_row(row);
     std::sort(result.begin(), result.end());
 
     uint64_t boss_edge = graph_->kmer_to_boss_index(
@@ -158,7 +154,7 @@ BinaryMatrix::SetBitPositions RowDiff<BaseMatrix>::get_row(Row row) const {
         row = graph::AnnotatedSequenceGraph::graph_to_anno_index(
                 graph_->boss_to_kmer_index(boss_edge));
 
-        auto diff_row = get_diff(row);
+        auto diff_row = diffs_.get_row(row);
         std::sort(diff_row.begin(), diff_row.end());
         add_diff(diff_row, &result);
     }

From 2b7256dbc588c1e7f9716c31add1ed1b2e5ff281 Mon Sep 17 00:00:00 2001
From: Mikhail Karasikov <mike-karas@ya.ru>
Date: Tue, 20 Apr 2021 10:47:11 +0200
Subject: [PATCH 06/23] added IntRowDiff -- a representation for rd-transformed
 integer matrices

---
 .../int_matrix/row_diff/int_row_diff.hpp      | 304 ++++++++++++++++++
 .../static_annotators_def.hpp                 |   5 +
 2 files changed, 309 insertions(+)
 create mode 100644 metagraph/src/annotation/int_matrix/row_diff/int_row_diff.hpp

diff --git a/metagraph/src/annotation/int_matrix/row_diff/int_row_diff.hpp b/metagraph/src/annotation/int_matrix/row_diff/int_row_diff.hpp
new file mode 100644
index 0000000000..8c84c3789c
--- /dev/null
+++ b/metagraph/src/annotation/int_matrix/row_diff/int_row_diff.hpp
@@ -0,0 +1,304 @@
+#ifndef __INT_ROW_DIFF_HPP__
+#define __INT_ROW_DIFF_HPP__
+
+#include <algorithm>
+#include <iostream>
+#include <string>
+#include <vector>
+
+#include "common/vectors/bit_vector_adaptive.hpp"
+#include "common/vector_map.hpp"
+#include "common/vector.hpp"
+#include "common/logger.hpp"
+#include "common/utils/template_utils.hpp"
+#include "graph/annotated_dbg.hpp"
+#include "graph/representation/succinct/boss.hpp"
+#include "graph/representation/succinct/dbg_succinct.hpp"
+#include "annotation/binary_matrix/row_diff/row_diff.hpp"
+#include "annotation/int_matrix/base/int_matrix.hpp"
+
+
+namespace mtg {
+namespace annot {
+namespace matrix {
+
+const size_t RD_PATH_RESERVE_SIZE = 2;
+
+
+template <class BaseMatrix>
+class IntRowDiff : public binmat::IRowDiff, public IntMatrix {
+  public:
+    using anchor_bv_type = bit_vector_small;
+    using fork_succ_bv_type = bit_vector_small;
+    static_assert(std::is_convertible<IntRowDiff*, IntMatrix*>::value);
+
+    IntRowDiff() {}
+
+    IntRowDiff(const graph::DBGSuccinct *graph, BaseMatrix&& diff)
+        : graph_(graph), diffs_(std::move(diff)) {}
+
+    bool get(Row i, Column j) const override;
+    std::vector<Row> get_column(Column j) const override;
+    SetBitPositions get_row(Row i) const override;
+    std::vector<SetBitPositions> get_rows(const std::vector<Row> &rows) const override;
+    // query integer values
+    RowValues get_row_values(Row i) const override;
+    std::vector<RowValues> get_row_values(const std::vector<Row> &rows) const override;
+
+    uint64_t num_columns() const override { return diffs_.num_columns(); }
+    uint64_t num_relations() const override { return diffs_.num_relations(); }
+    uint64_t num_rows() const override { return diffs_.num_rows(); }
+
+    bool load(std::istream &in) override;
+    void serialize(std::ostream &out) const override;
+
+    void load_fork_succ(const std::string &filename);
+    void load_anchor(const std::string &filename);
+    void set_graph(const graph::DBGSuccinct *graph) { graph_ = graph; }
+
+    const anchor_bv_type& anchor() const { return anchor_; }
+    const graph::DBGSuccinct& graph() const { return *graph_; }
+    const BaseMatrix& diffs() const { return diffs_; }
+    BaseMatrix& diffs() { return diffs_; }
+
+  private:
+    static void add_diff(const RowValues &diff, RowValues *row);
+
+    const graph::DBGSuccinct *graph_ = nullptr;
+
+    BaseMatrix diffs_;
+    anchor_bv_type anchor_;
+    fork_succ_bv_type fork_succ_;
+};
+
+
+template <class BaseMatrix>
+bool IntRowDiff<BaseMatrix>::get(Row i, Column j) const {
+    SetBitPositions set_bits = get_row(i);
+    auto v = std::lower_bound(set_bits.begin(), set_bits.end(), j);
+    return v != set_bits.end() && *v == j;
+}
+
+template <class BaseMatrix>
+std::vector<IntMatrix::Row> IntRowDiff<BaseMatrix>::get_column(Column j) const {
+    assert(anchor_.size() == diffs_.num_rows() && "anchors must be loaded");
+    assert(!fork_succ_.size() || fork_succ_.size() == graph_->get_boss().get_last().size());
+
+    // TODO: implement a more efficient algorithm
+    std::vector<Row> result;
+    for (Row i = 0; i < num_rows(); ++i) {
+        if (get(i, j))
+            result.push_back(i);
+    }
+    return result;
+}
+
+template <class BaseMatrix>
+IntMatrix::SetBitPositions IntRowDiff<BaseMatrix>::get_row(Row i) const {
+    RowValues row = get_row_values(i);
+    SetBitPositions result(row.size());
+    for (size_t k = 0; k < row.size(); ++k) {
+        result[k] = row[k].first;
+    }
+    return result;
+}
+
+template <class BaseMatrix>
+IntMatrix::RowValues IntRowDiff<BaseMatrix>::get_row_values(Row row) const {
+    assert(anchor_.size() == diffs_.num_rows() && "anchors must be loaded");
+    assert(!fork_succ_.size() || fork_succ_.size() == graph_->get_boss().get_last().size());
+
+    RowValues result = diffs_.get_row_values(row);
+    std::sort(result.begin(), result.end());
+
+    uint64_t boss_edge = graph_->kmer_to_boss_index(
+            graph::AnnotatedSequenceGraph::anno_to_graph_index(row));
+    const graph::boss::BOSS &boss = graph_->get_boss();
+    const bit_vector &rd_succ = fork_succ_.size() ? fork_succ_ : boss.get_last();
+
+    while (!anchor_[row]) {
+        boss_edge = boss.row_diff_successor(boss_edge, rd_succ);
+
+        row = graph::AnnotatedSequenceGraph::graph_to_anno_index(
+                graph_->boss_to_kmer_index(boss_edge));
+
+        RowValues diff_row = diffs_.get_row_values(row);
+        std::sort(diff_row.begin(), diff_row.end());
+        add_diff(diff_row, &result);
+    }
+
+    return result;
+}
+
+template <class BaseMatrix>
+std::vector<IntMatrix::SetBitPositions>
+IntRowDiff<BaseMatrix>::get_rows(const std::vector<Row> &row_ids) const {
+    std::vector<SetBitPositions> result;
+    result.reserve(row_ids.size());
+
+    for (auto&& row : get_row_values(row_ids)) {
+        result.emplace_back(row.size());
+        for (size_t k = 0; k < row.size(); ++k) {
+            result.back()[k] = row[k].first;
+        }
+        row = RowValues();
+    }
+
+    return result;
+}
+
+template <class BaseMatrix>
+std::vector<IntMatrix::RowValues>
+IntRowDiff<BaseMatrix>::get_row_values(const std::vector<Row> &row_ids) const {
+    assert(anchor_.size() == diffs_.num_rows() && "anchors must be loaded");
+    assert(!fork_succ_.size() || fork_succ_.size() == graph_->get_boss().get_last().size());
+
+    // diff rows annotating nodes along the row-diff paths
+    std::vector<Row> rd_ids;
+    rd_ids.reserve(row_ids.size() * RD_PATH_RESERVE_SIZE);
+
+    // map row index to its index in |rd_rows|
+    VectorMap<Row, size_t> node_to_rd;
+    node_to_rd.reserve(row_ids.size() * RD_PATH_RESERVE_SIZE);
+
+    // Truncated row-diff paths, indexes to |rd_rows|.
+    // The last index in each path points to an anchor or to a row which had
+    // been reached before, and thus, will be reconstructed before this one.
+    std::vector<std::vector<size_t>> rd_paths_trunc(row_ids.size());
+
+    const graph::boss::BOSS &boss = graph_->get_boss();
+    const bit_vector &rd_succ = fork_succ_.size() ? fork_succ_ : boss.get_last();
+
+    for (size_t i = 0; i < row_ids.size(); ++i) {
+        Row row = row_ids[i];
+
+        graph::boss::BOSS::edge_index boss_edge = graph_->kmer_to_boss_index(
+                graph::AnnotatedSequenceGraph::anno_to_graph_index(row));
+
+        while (true) {
+            row = graph::AnnotatedSequenceGraph::graph_to_anno_index(
+                    graph_->boss_to_kmer_index(boss_edge));
+
+            auto [it, is_new] = node_to_rd.try_emplace(row, rd_ids.size());
+            rd_paths_trunc[i].push_back(it.value());
+
+            // If a node had been reached before, we interrupt the diff path.
+            // The annotation for that node will have been reconstructed earlier
+            // than for other nodes in this path as well. Thus, we will start
+            // reconstruction from that node and don't need its successors.
+            if (!is_new)
+                break;
+
+            rd_ids.push_back(row);
+
+            if (anchor_[row])
+                break;
+
+            boss_edge = boss.row_diff_successor(boss_edge, rd_succ);
+        }
+    }
+
+    node_to_rd = VectorMap<Row, size_t>();
+
+    std::vector<RowValues> rd_rows = diffs_.get_row_values(rd_ids);
+
+    rd_ids = std::vector<Row>();
+
+    // reconstruct annotation rows from row-diff
+    std::vector<RowValues> rows(row_ids.size());
+
+    for (size_t i = 0; i < row_ids.size(); ++i) {
+        RowValues &result = rows[i];
+        // propagate back and reconstruct full annotations for predecessors
+        for (auto it = rd_paths_trunc[i].rbegin(); it != rd_paths_trunc[i].rend(); ++it) {
+            std::sort(rd_rows[*it].begin(), rd_rows[*it].end());
+            add_diff(rd_rows[*it], &result);
+            // replace diff row with full reconstructed annotation
+            rd_rows[*it] = result;
+        }
+    }
+
+    return rows;
+}
+
+template <class BaseMatrix>
+bool IntRowDiff<BaseMatrix>::load(std::istream &in) {
+    std::string version(4, '\0');
+    in.read(version.data(), 4);
+    return anchor_.load(in) && fork_succ_.load(in) && diffs_.load(in);
+}
+
+template <class BaseMatrix>
+void IntRowDiff<BaseMatrix>::serialize(std::ostream &out) const {
+    out.write("v2.0", 4);
+    anchor_.serialize(out);
+    fork_succ_.serialize(out);
+    diffs_.serialize(out);
+}
+
+template <class BaseMatrix>
+void IntRowDiff<BaseMatrix>::add_diff(const RowValues &diff, RowValues *row) {
+    assert(std::is_sorted(row->begin(), row->end()));
+    assert(std::is_sorted(diff.begin(), diff.end()));
+
+    if (diff.empty())
+        return;
+
+    RowValues result;
+    result.reserve(row->size() + diff.size());
+
+    auto it = row->begin();
+    auto it2 = diff.begin();
+    while (it != row->end() && it2 != diff.end()) {
+        if (it->first < it2->first) {
+            result.push_back(*it);
+            ++it;
+        } else if (it->first > it2->first) {
+            result.push_back(*it2);
+            ++it2;
+        } else {
+            result.push_back(*it);
+            result.back().second += it2->second;
+            ++it;
+            ++it2;
+        }
+    }
+    std::copy(it, row->end(), std::back_inserter(result));
+    std::copy(it2, diff.end(), std::back_inserter(result));
+
+    row->swap(result);
+}
+
+template <class BaseMatrix>
+void IntRowDiff<BaseMatrix>::load_anchor(const std::string &filename) {
+    if (!std::filesystem::exists(filename)) {
+        common::logger->error("Can't read anchor file: {}", filename);
+        std::exit(1);
+    }
+    std::ifstream in(filename, ios::binary);
+    if (!in.good()) {
+        common::logger->error("Could not open anchor file {}", filename);
+        std::exit(1);
+    }
+    anchor_.load(in);
+}
+
+template <class BaseMatrix>
+void IntRowDiff<BaseMatrix>::load_fork_succ(const std::string &filename) {
+    if (!std::filesystem::exists(filename)) {
+        common::logger->error("Can't read fork successor file: {}", filename);
+        std::exit(1);
+    }
+    std::ifstream in(filename, ios::binary);
+    if (!in.good()) {
+        common::logger->error("Could not open fork successor file {}", filename);
+        std::exit(1);
+    }
+    fork_succ_.load(in);
+}
+
+} // namespace matrix
+} // namespace annot
+} // namespace mtg
+
+#endif // __INT_ROW_DIFF_HPP__
diff --git a/metagraph/src/annotation/representation/annotation_matrix/static_annotators_def.hpp b/metagraph/src/annotation/representation/annotation_matrix/static_annotators_def.hpp
index fc762d2e84..ed676b6d59 100644
--- a/metagraph/src/annotation/representation/annotation_matrix/static_annotators_def.hpp
+++ b/metagraph/src/annotation/representation/annotation_matrix/static_annotators_def.hpp
@@ -14,6 +14,7 @@
 #include "annotation/binary_matrix/row_sparse/row_sparse.hpp"
 #include "annotation/binary_matrix/row_vector/unique_row_binmat.hpp"
 #include "annotation/int_matrix/extended/csc_matrix.hpp"
+#include "annotation/int_matrix/row_diff/int_row_diff.hpp"
 
 
 namespace mtg {
@@ -43,6 +44,8 @@ typedef StaticBinRelAnnotator<binmat::RowDiff<binmat::RowSparse>, std::string> R
 
 typedef StaticBinRelAnnotator<matrix::CSCMatrix<binmat::BRWT, sdsl::dac_vector_dp<>>, std::string> IntMultiBRWTAnnotator;
 
+typedef StaticBinRelAnnotator<matrix::IntRowDiff<matrix::CSCMatrix<binmat::BRWT, sdsl::dac_vector_dp<>>>, std::string> IntRowDiffBRWTAnnotator;
+
 
 template <>
 inline const std::string RowFlatAnnotator::kExtension = ".flat.annodbg";
@@ -68,6 +71,8 @@ template <>
 inline const std::string RowDiffRowSparseAnnotator::kExtension = ".row_diff_sparse.annodbg";
 template <>
 inline const std::string IntMultiBRWTAnnotator::kExtension = ".int_brwt.annodbg";
+template <>
+inline const std::string IntRowDiffBRWTAnnotator::kExtension = ".rd_int_brwt.annodbg";
 
 } // namespace annot
 } // namespace mtg

From f4a89762248c88582199b255968379d25127c5d1 Mon Sep 17 00:00:00 2001
From: Mikhail Karasikov <mike-karas@ya.ru>
Date: Tue, 20 Apr 2021 16:16:05 +0200
Subject: [PATCH 07/23] transform to IntRowDiff<Int-Multi-BRWT>

---
 metagraph/integration_tests/test_query.py     |  1 +
 .../annotation_matrix/annotation_matrix.cpp   |  2 +
 .../static_annotators_def.hpp                 |  2 +-
 metagraph/src/cli/annotate.cpp                | 14 ++++---
 metagraph/src/cli/config/config.cpp           |  7 +++-
 metagraph/src/cli/config/config.hpp           |  1 +
 metagraph/src/cli/load/load_annotation.cpp    |  7 ++++
 metagraph/src/cli/transform_annotation.cpp    | 41 ++++++++++++++++++-
 8 files changed, 66 insertions(+), 9 deletions(-)

diff --git a/metagraph/integration_tests/test_query.py b/metagraph/integration_tests/test_query.py
index 1143ce164c..0732ee4b1a 100644
--- a/metagraph/integration_tests/test_query.py
+++ b/metagraph/integration_tests/test_query.py
@@ -30,6 +30,7 @@
                        'row_diff_sparse': '.row_diff_sparse.annodbg',
                        'rb_brwt': '.rb_brwt.annodbg',
                        'brwt': '.brwt.annodbg',
+                       'int_brwt': '.int_brwt.annodbg',
                        'rbfish': '.rbfish.annodbg',
                        'flat': '.flat.annodbg'}
 
diff --git a/metagraph/src/annotation/representation/annotation_matrix/annotation_matrix.cpp b/metagraph/src/annotation/representation/annotation_matrix/annotation_matrix.cpp
index fa96dde45c..9b855056f8 100644
--- a/metagraph/src/annotation/representation/annotation_matrix/annotation_matrix.cpp
+++ b/metagraph/src/annotation/representation/annotation_matrix/annotation_matrix.cpp
@@ -226,5 +226,7 @@ template class StaticBinRelAnnotator<binmat::RowDiff<binmat::RowSparse>, std::st
 
 template class StaticBinRelAnnotator<matrix::CSCMatrix<binmat::BRWT, sdsl::dac_vector_dp<>>, std::string>;
 
+template class StaticBinRelAnnotator<matrix::IntRowDiff<matrix::CSCMatrix<binmat::BRWT, sdsl::dac_vector_dp<>>>, std::string>;
+
 } // namespace annot
 } // namespace mtg
diff --git a/metagraph/src/annotation/representation/annotation_matrix/static_annotators_def.hpp b/metagraph/src/annotation/representation/annotation_matrix/static_annotators_def.hpp
index ed676b6d59..c6c3b7fd62 100644
--- a/metagraph/src/annotation/representation/annotation_matrix/static_annotators_def.hpp
+++ b/metagraph/src/annotation/representation/annotation_matrix/static_annotators_def.hpp
@@ -72,7 +72,7 @@ inline const std::string RowDiffRowSparseAnnotator::kExtension = ".row_diff_spar
 template <>
 inline const std::string IntMultiBRWTAnnotator::kExtension = ".int_brwt.annodbg";
 template <>
-inline const std::string IntRowDiffBRWTAnnotator::kExtension = ".rd_int_brwt.annodbg";
+inline const std::string IntRowDiffBRWTAnnotator::kExtension = ".row_diff_int_brwt.annodbg";
 
 } // namespace annot
 } // namespace mtg
diff --git a/metagraph/src/cli/annotate.cpp b/metagraph/src/cli/annotate.cpp
index 7526835268..65237c6072 100644
--- a/metagraph/src/cli/annotate.cpp
+++ b/metagraph/src/cli/annotate.cpp
@@ -154,15 +154,17 @@ void add_kmer_counts(const std::string &file,
     // remember the number of base labels to remove those unique to each sequence quickly
     const size_t num_base_labels = labels.size();
 
-    mtg::common::logger->trace("Parsing k-mer counts from '{}'",
-        utils::remove_suffix(file, ".gz", ".fasta") + ".kmer_counts.gz"
-    );
+    const std::string &counts_fname
+            = utils::remove_suffix(file, ".gz", ".fasta") + ".kmer_counts.gz";
+
+    logger->trace("Parsing k-mer counts from '{}'", counts_fname);
+
     read_extended_fasta_file_critical<uint32_t>(file, "kmer_counts",
         [&](size_t k, const kseq_t *read_stream, const uint32_t *kmer_counts) {
             if (k != graph.get_k()) {
-                mtg::common::logger->error("File '{}' contains counts for k-mers of "
-                                          "length {} but graph is constructed with k={}",
-                                          file, k, graph.get_k());
+                logger->error("File '{}' contains counts for k-mers of "
+                              "length {} but graph is constructed with k={}",
+                              file, k, graph.get_k());
                 exit(1);
             }
             assert(read_stream->seq.l >= k && "sequences can't be shorter than k-mers");
diff --git a/metagraph/src/cli/config/config.cpp b/metagraph/src/cli/config/config.cpp
index 2ba62ebcf0..08f8c9d9d3 100644
--- a/metagraph/src/cli/config/config.cpp
+++ b/metagraph/src/cli/config/config.cpp
@@ -561,6 +561,7 @@ Config::Config(int argc, char *argv[]) {
     if (identity == TRANSFORM_ANNOTATION) {
         const bool to_row_diff = anno_type == RowDiff
                                     || anno_type == RowDiffBRWT
+                                    || anno_type == IntRowDiffBRWT
                                     || anno_type == RowDiffRowSparse;
         if (to_row_diff && !infbase.size()) {
             std::cerr << "Path to graph must be passed with '-i <GRAPH>'" << std::endl;
@@ -678,6 +679,8 @@ std::string Config::annotype_to_string(AnnotationType state) {
             return "row_sparse";
         case IntBRWT:
             return "int_brwt";
+        case IntRowDiffBRWT:
+            return "row_diff_int_brwt";
     }
     throw std::runtime_error("Never happens");
 }
@@ -709,6 +712,8 @@ Config::AnnotationType Config::string_to_annotype(const std::string &string) {
         return AnnotationType::RowSparse;
     } else if (string == "int_brwt") {
         return AnnotationType::IntBRWT;
+    } else if (string == "row_diff_int_brwt") {
+        return AnnotationType::IntRowDiffBRWT;
     } else {
         std::cerr << "Error: unknown annotation representation" << std::endl;
         exit(1);
@@ -771,7 +776,7 @@ DeBruijnGraph::Mode Config::string_to_graphmode(const std::string &string) {
 
 void Config::print_usage(const std::string &prog_name, IdentityType identity) {
     const char annotation_list[] = "\t\t( column, brwt, rb_brwt, int_brwt,\n"
-                                   "\t\t  row_diff, row_diff_brwt, row_diff_sparse,\n"
+                                   "\t\t  row_diff, row_diff_brwt, row_diff_sparse, row_diff_int_brwt,\n"
                                    "\t\t  row, flat, row_sparse, rbfish, bin_rel_wt, bin_rel_wt_sdsl )";
 
     switch (identity) {
diff --git a/metagraph/src/cli/config/config.hpp b/metagraph/src/cli/config/config.hpp
index 7d9318bea5..a1fd1dc19e 100644
--- a/metagraph/src/cli/config/config.hpp
+++ b/metagraph/src/cli/config/config.hpp
@@ -192,6 +192,7 @@ class Config {
         RBFish,
         RbBRWT,
         IntBRWT,
+        IntRowDiffBRWT,
     };
 
     enum GraphType {
diff --git a/metagraph/src/cli/load/load_annotation.cpp b/metagraph/src/cli/load/load_annotation.cpp
index 401eb6453f..9482b47fd8 100644
--- a/metagraph/src/cli/load/load_annotation.cpp
+++ b/metagraph/src/cli/load/load_annotation.cpp
@@ -57,6 +57,9 @@ Config::AnnotationType parse_annotation_type(const std::string &filename) {
     } else if (utils::ends_with(filename, annot::IntMultiBRWTAnnotator::kExtension)) {
         return Config::AnnotationType::IntBRWT;
 
+    } else if (utils::ends_with(filename, annot::IntRowDiffBRWTAnnotator::kExtension)) {
+        return Config::AnnotationType::IntRowDiffBRWT;
+
     } else {
         logger->error("Unknown annotation format in '{}'", filename);
         exit(1);
@@ -130,6 +133,10 @@ initialize_annotation(Config::AnnotationType anno_type,
             annotation.reset(new annot::IntMultiBRWTAnnotator());
             break;
         }
+        case Config::IntRowDiffBRWT: {
+            annotation.reset(new annot::IntRowDiffBRWTAnnotator());
+            break;
+        }
     }
 
     return annotation;
diff --git a/metagraph/src/cli/transform_annotation.cpp b/metagraph/src/cli/transform_annotation.cpp
index 8142455100..0a0c2b508a 100644
--- a/metagraph/src/cli/transform_annotation.cpp
+++ b/metagraph/src/cli/transform_annotation.cpp
@@ -685,6 +685,40 @@ int transform_annotation(Config *config) {
                 rb_brwt_annotator->serialize(config->outfbase);
                 break;
             }
+            case Config::IntRowDiffBRWT: {
+                assert(config->infbase.size());
+                const std::string anchors_file = config->infbase + annot::binmat::kRowDiffAnchorExt;
+                if (!std::filesystem::exists(anchors_file)) {
+                    logger->error("Anchor bitmap {} does not exist. Run the row_diff"
+                                  " transform followed by anchor optimization.", anchors_file);
+                    std::exit(1);
+                }
+                const std::string fork_succ_file = config->infbase + annot::binmat::kRowDiffForkSuccExt;
+                if (!std::filesystem::exists(fork_succ_file)) {
+                    logger->error("Fork successor bitmap {} does not exist", fork_succ_file);
+                    std::exit(1);
+                }
+                auto int_annotation = convert_to_IntMultiBRWT(files, *config, timer);
+                logger->trace("Annotation converted in {} sec", timer.elapsed());
+
+                using CSCMatrix = matrix::CSCMatrix<binmat::BRWT, sdsl::dac_vector_dp<>>;
+
+                IntRowDiffBRWTAnnotator annotation(
+                        std::make_unique<matrix::IntRowDiff<CSCMatrix>>(nullptr,
+                                std::move(*int_annotation.release_matrix())),
+                        int_annotation.get_label_encoder());
+
+                const_cast<matrix::IntRowDiff<CSCMatrix> &>(annotation.get_matrix())
+                        .load_anchor(anchors_file);
+                const_cast<matrix::IntRowDiff<CSCMatrix> &>(annotation.get_matrix())
+                        .load_fork_succ(fork_succ_file);
+
+                logger->trace("RowDiff support bitmaps loaded");
+
+                annotation.serialize(config->outfbase);
+                logger->trace("Serialized to {}", config->outfbase);
+                break;
+            }
         }
 
     } else if (input_anno_type == Config::RowDiff) {
@@ -859,6 +893,9 @@ int relax_multi_brwt(Config *config) {
         case Config::IntBRWT:
             annotator = std::make_unique<IntMultiBRWTAnnotator>();
             break;
+        case Config::IntRowDiffBRWT:
+            annotator = std::make_unique<IntRowDiffBRWTAnnotator>();
+            break;
         default:
             logger->error("Relaxation only supported for BRWT and RowDiffBRWT");
             exit(1);
@@ -878,7 +915,9 @@ int relax_multi_brwt(Config *config) {
             ? dynamic_cast<MultiBRWTAnnotator &>(*annotator).get_matrix()
             : (anno_type == Config::IntBRWT
                 ? dynamic_cast<IntMultiBRWTAnnotator &>(*annotator).get_matrix().get_binary_matrix()
-                : dynamic_cast<RowDiffBRWTAnnotator &>(*annotator).get_matrix().diffs());
+                : (anno_type == Config::IntRowDiffBRWT
+                    ? dynamic_cast<IntRowDiffBRWTAnnotator &>(*annotator).get_matrix().diffs().get_binary_matrix()
+                    : dynamic_cast<RowDiffBRWTAnnotator &>(*annotator).get_matrix().diffs()));
     relax_BRWT(const_cast<binmat::BRWT *>(&matrix), config->relax_arity_brwt,
                get_num_threads());
 

From 6f6260b5738b0309c9a5b4ca4fb00d64a98675ef Mon Sep 17 00:00:00 2001
From: Mikhail Karasikov <mike-karas@ya.ru>
Date: Tue, 20 Apr 2021 18:21:59 +0200
Subject: [PATCH 08/23] increment k-mer counts from any input source

---
 metagraph/src/cli/annotate.cpp | 77 ++++++++++++++++++++++++----------
 1 file changed, 55 insertions(+), 22 deletions(-)

diff --git a/metagraph/src/cli/annotate.cpp b/metagraph/src/cli/annotate.cpp
index 65237c6072..778d1f495b 100644
--- a/metagraph/src/cli/annotate.cpp
+++ b/metagraph/src/cli/annotate.cpp
@@ -246,28 +246,61 @@ void annotate_data(std::shared_ptr<graph::DeBruijnGraph> graph,
     if (config.count_kmers) {
         // add k-mer counts
         for (const auto &file : files) {
-            add_kmer_counts(
-                file,
-                anno_graph->get_graph(),
-                forward_and_reverse,
-                config.filename_anno,
-                config.annotate_sequence_headers,
-                config.fasta_anno_comment_delim,
-                config.fasta_header_delimiter,
-                config.anno_labels,
-                [&](std::string sequence,
-                            std::vector<std::string> labels,
-                            std::vector<uint64_t> kmer_counts) {
-                    thread_pool.enqueue(
-                        [&](std::string &sequence,
-                                std::vector<std::string> &labels,
-                                std::vector<uint64_t> &kmer_counts) {
-                            anno_graph->add_kmer_counts(sequence, labels, std::move(kmer_counts));
-                        },
-                        std::move(sequence), std::move(labels), std::move(kmer_counts)
-                    );
-                }
-            );
+            const std::string &counts_fname
+                    = utils::remove_suffix(file, ".gz", ".fasta") + ".kmer_counts.gz";
+
+            if (fs::exists(counts_fname)) {
+                add_kmer_counts(
+                    file,
+                    anno_graph->get_graph(),
+                    forward_and_reverse,
+                    config.filename_anno,
+                    config.annotate_sequence_headers,
+                    config.fasta_anno_comment_delim,
+                    config.fasta_header_delimiter,
+                    config.anno_labels,
+                    [&](std::string sequence,
+                                std::vector<std::string> labels,
+                                std::vector<uint64_t> kmer_counts) {
+                        thread_pool.enqueue(
+                            [&](std::string &sequence,
+                                    std::vector<std::string> &labels,
+                                    std::vector<uint64_t> &kmer_counts) {
+                                anno_graph->add_kmer_counts(sequence, labels, std::move(kmer_counts));
+                            },
+                            std::move(sequence), std::move(labels), std::move(kmer_counts)
+                        );
+                    }
+                );
+            } else {
+                logger->warn("No k-mer counts found at '{}'. Every input k-mer"
+                             " will have count 1.", counts_fname);
+                call_annotations(
+                    file,
+                    config.refpath,
+                    anno_graph->get_graph(),
+                    forward_and_reverse,
+                    config.min_count,
+                    config.max_count,
+                    config.filename_anno,
+                    config.annotate_sequence_headers,
+                    config.fasta_anno_comment_delim,
+                    config.fasta_header_delimiter,
+                    config.anno_labels,
+                    [&](std::string sequence, auto labels) {
+                        thread_pool.enqueue(
+                            [&anno_graph](const std::string &sequence, const auto &labels) {
+                                const size_t k = anno_graph->get_graph().get_k();
+                                if (sequence.size() >= k) {
+                                    anno_graph->add_kmer_counts(sequence, labels,
+                                        std::vector<uint64_t>(sequence.size() - k + 1, 1));
+                                }
+                            },
+                            std::move(sequence), std::move(labels)
+                        );
+                    }
+                );
+            }
         }
     }
 

From c753da94941764db1ec627c27a3aef6059db88aa Mon Sep 17 00:00:00 2001
From: Mikhail Karasikov <mike-karas@ya.ru>
Date: Wed, 21 Apr 2021 11:05:21 +0200
Subject: [PATCH 09/23] pack int-row-diff into column-compressed

---
 .../annotate_column_compressed.cpp            | 14 +++++++
 .../annotate_column_compressed.hpp            |  7 ++++
 metagraph/src/annotation/row_diff_builder.cpp | 41 +++++++++++--------
 3 files changed, 44 insertions(+), 18 deletions(-)

diff --git a/metagraph/src/annotation/representation/column_compressed/annotate_column_compressed.cpp b/metagraph/src/annotation/representation/column_compressed/annotate_column_compressed.cpp
index 18c9cffece..de34b9d61d 100644
--- a/metagraph/src/annotation/representation/column_compressed/annotate_column_compressed.cpp
+++ b/metagraph/src/annotation/representation/column_compressed/annotate_column_compressed.cpp
@@ -55,6 +55,20 @@ ColumnCompressed<Label>::ColumnCompressed(sdsl::bit_vector&& column,
     flushed_ = false;
 }
 
+template <typename Label>
+ColumnCompressed<Label>::ColumnCompressed(std::vector<std::unique_ptr<bit_vector>>&& columns,
+                                          const LabelEncoder<Label> &label_encoder,
+                                          size_t num_columns_cached,
+                                          const std::string &swap_dir,
+                                          uint64_t buffer_size_bytes,
+                                          uint8_t count_width)
+      : ColumnCompressed(columns.at(0)->size(),
+                         num_columns_cached, swap_dir, buffer_size_bytes, count_width) {
+    bitmatrix_ = std::move(columns);
+    label_encoder_ = label_encoder;
+    flushed_ = true;
+}
+
 template <typename Label>
 ColumnCompressed<Label>::~ColumnCompressed() {
     // Note: this is needed to make sure that everything is flushed to bitmatrix_
diff --git a/metagraph/src/annotation/representation/column_compressed/annotate_column_compressed.hpp b/metagraph/src/annotation/representation/column_compressed/annotate_column_compressed.hpp
index c4c03310d6..f052fd99e4 100644
--- a/metagraph/src/annotation/representation/column_compressed/annotate_column_compressed.hpp
+++ b/metagraph/src/annotation/representation/column_compressed/annotate_column_compressed.hpp
@@ -42,6 +42,13 @@ class ColumnCompressed : public MultiLabelEncoded<Label> {
                      uint64_t buffer_size_bytes = 1e9,
                      uint8_t count_width = 8);
 
+    ColumnCompressed(std::vector<std::unique_ptr<bit_vector>>&& columns,
+                     const LabelEncoder<Label> &label_encoder,
+                     size_t num_columns_cached = 1,
+                     const std::string &swap_dir = "",
+                     uint64_t buffer_size_bytes = 1e9,
+                     uint8_t count_width = 8);
+
     ColumnCompressed(const ColumnCompressed&) = delete;
     ColumnCompressed& operator=(const ColumnCompressed&) = delete;
 
diff --git a/metagraph/src/annotation/row_diff_builder.cpp b/metagraph/src/annotation/row_diff_builder.cpp
index 05cd62e097..2c979ed577 100644
--- a/metagraph/src/annotation/row_diff_builder.cpp
+++ b/metagraph/src/annotation/row_diff_builder.cpp
@@ -680,6 +680,7 @@ void convert_batch_to_row_diff(const std::string &pred_succ_fprefix,
     constexpr bool with_counts = utils::is_pair_v<T>;
     std::vector<std::vector<sdsl::int_vector<>>> counts;
     if (with_counts) {
+        counts.resize(source_files.size());
         #pragma omp parallel for num_threads(get_num_threads())
         for (size_t i = 0; i < source_files.size(); ++i) {
             const auto &counts_fname = utils::remove_suffix(source_files[i],
@@ -925,7 +926,7 @@ void convert_batch_to_row_diff(const std::string &pred_succ_fprefix,
     // free memory occupied by original columns
     sources.clear();
 
-    std::vector<std::unique_ptr<RowDiffColumnAnnotator>> row_diff(label_encoders.size());
+    std::vector<std::vector<std::unique_ptr<bit_vector>>> diff_columns(label_encoders.size());
 
     const uint32_t files_open_per_thread
             = MAX_NUM_FILES_OPEN / std::max((uint32_t)1, num_threads);
@@ -1009,27 +1010,31 @@ void convert_batch_to_row_diff(const std::string &pred_succ_fprefix,
                                                          row_diff_bits[l_idx][j]);
         }
 
-        row_diff[l_idx] = std::make_unique<RowDiffColumnAnnotator>(
-                std::make_unique<RowDiff<ColumnMajor>>(nullptr, ColumnMajor(std::move(columns))),
-                std::move(label_encoders[l_idx]));
+        if (compute_row_reduction) {
+            diff_columns[l_idx] = std::move(columns);
 
-        if (!compute_row_reduction) {
-            auto fpath = col_out_dir/fs::path(source_files[l_idx])
-                                    .filename()
-                                    .replace_extension()
-                                    .replace_extension(row_diff[l_idx]->kExtension);
-
-            row_diff[l_idx]->serialize(fpath);
-
-            // TODO: IntRowDiff::kExtension
+        } else {
+            auto fpath = col_out_dir/fs::path(source_files[l_idx]).filename();
             if constexpr(with_counts) {
-                std::ofstream outstream(std::string(fpath) + ".counts", std::ios::binary);
+                fpath.replace_extension().replace_extension(ColumnCompressed<>::kExtension);
+                ColumnCompressed<>(std::move(columns), label_encoders[l_idx]).serialize(fpath);
+                logger->trace("Serialized {}", fpath);
+
+                fpath.replace_extension().replace_extension(ColumnCompressed<>::kCountExtension);
+                std::ofstream outstream(fpath, std::ios::binary);
                 for (size_t j = 0; j < label_encoders[l_idx].size(); ++j) {
                     counts[l_idx][j].serialize(outstream);
                 }
-            }
+                logger->trace("Serialized {}", fpath);
+                counts[l_idx].clear();
 
-            logger->trace("Serialized {}", fpath);
+            } else {
+                fpath.replace_extension().replace_extension(RowDiffColumnAnnotator::kExtension);
+                RowDiffColumnAnnotator(
+                        std::make_unique<RowDiff<ColumnMajor>>(nullptr, ColumnMajor(std::move(columns))),
+                        std::move(label_encoders[l_idx])).serialize(fpath);
+                logger->trace("Serialized {}", fpath);
+            }
         }
     }
 
@@ -1045,8 +1050,8 @@ void convert_batch_to_row_diff(const std::string &pred_succ_fprefix,
         row_nbits_block.assign(std::min(BLOCK_SIZE, row_reduction.size() - chunk), 0);
 
         #pragma omp parallel for num_threads(num_threads) schedule(dynamic)
-        for (size_t l_idx = 0; l_idx < row_diff.size(); ++l_idx) {
-            for (const auto &col_ptr : row_diff[l_idx]->get_matrix().diffs().data()) {
+        for (size_t l_idx = 0; l_idx < diff_columns.size(); ++l_idx) {
+            for (const auto &col_ptr : diff_columns[l_idx]) {
                 col_ptr->call_ones_in_range(chunk, chunk + row_nbits_block.size(),
                     [&](uint64_t i) {
                         __atomic_add_fetch(&row_nbits_block[i - chunk], 1, __ATOMIC_RELAXED);

From df8e577a8791d97f0c3028cc7e5feb122d4f6876 Mon Sep 17 00:00:00 2001
From: Mikhail Karasikov <mike-karas@ya.ru>
Date: Wed, 21 Apr 2021 11:05:37 +0200
Subject: [PATCH 10/23] fixed ci

---
 .github/workflows/main.yml | 11 +----------
 1 file changed, 1 insertion(+), 10 deletions(-)

diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml
index 25ef854a78..2792dde3c7 100644
--- a/.github/workflows/main.yml
+++ b/.github/workflows/main.yml
@@ -50,16 +50,7 @@ jobs:
     - name: install dependencies
       run: |
         # If clang, use the default version. Otherwise the compiler install with apt-get.
-        if [[ "${{ matrix.compiler }}" == "g++-8" ]]; then
-            sudo apt-get --allow-downgrades install \
-                g++-8=8-20180414-1ubuntu2 \
-                gcc-8-base=8-20180414-1ubuntu2 \
-                gcc-8=8-20180414-1ubuntu2 \
-                libstdc++-8-dev=8-20180414-1ubuntu2 \
-                cpp-8=8-20180414-1ubuntu2 \
-                libgcc-8-dev=8-20180414-1ubuntu2 \
-                libmpx2=8-20180414-1ubuntu2
-        elif [[ "${{ matrix.compiler }}" != "clang" ]]; then
+        if [[ "${{ matrix.compiler }}" != "clang" ]]; then
             sudo apt-get install ${{ matrix.compiler }}
         fi
 

From 5a062e53e0548bcd8e0c716a9728058b1b917b88 Mon Sep 17 00:00:00 2001
From: Mikhail Karasikov <mike-karas@ya.ru>
Date: Wed, 21 Apr 2021 12:21:52 +0200
Subject: [PATCH 11/23] unified row-diff annotations and their loading

---
 .../binary_matrix/row_diff/row_diff.hpp       | 23 +++++-----
 .../int_matrix/row_diff/int_row_diff.hpp      |  9 ++--
 metagraph/src/annotation/row_diff_builder.cpp |  6 +--
 .../src/cli/load/load_annotated_graph.cpp     | 44 +++++++------------
 4 files changed, 36 insertions(+), 46 deletions(-)

diff --git a/metagraph/src/annotation/binary_matrix/row_diff/row_diff.hpp b/metagraph/src/annotation/binary_matrix/row_diff/row_diff.hpp
index 7dce0c4456..feb66a3ba5 100644
--- a/metagraph/src/annotation/binary_matrix/row_diff/row_diff.hpp
+++ b/metagraph/src/annotation/binary_matrix/row_diff/row_diff.hpp
@@ -30,6 +30,12 @@ const size_t RD_PATH_RESERVE_SIZE = 2;
 class IRowDiff {
   public:
     virtual ~IRowDiff() {}
+
+    const graph::DBGSuccinct* graph() const { return graph_; }
+    void set_graph(const graph::DBGSuccinct *graph) { graph_ = graph; }
+
+  protected:
+    const graph::DBGSuccinct *graph_ = nullptr;
 };
 
 /**
@@ -62,19 +68,14 @@ class RowDiff : public IRowDiff, public BinaryMatrix {
     RowDiff() {}
 
     RowDiff(const graph::DBGSuccinct *graph, BaseMatrix &&diff)
-        : graph_(graph), diffs_(std::move(diff)) {}
-
-    uint64_t num_columns() const override { return diffs_.num_columns(); }
-
-    const graph::DBGSuccinct* graph() const { return graph_; }
+          : diffs_(std::move(diff)) { graph_ = graph; }
 
     /**
-     * Returns the number of set bits in the matrix.
+     * Returns the number of set bits in the row-diff transformed matrix.
      */
     uint64_t num_relations() const override { return diffs_.num_relations(); }
-
+    uint64_t num_columns() const override { return diffs_.num_columns(); }
     uint64_t num_rows() const override { return diffs_.num_rows(); }
-    void set_graph(const graph::DBGSuccinct *graph) { graph_ = graph; }
 
     bool get(Row row, Column column) const override;
 
@@ -100,8 +101,6 @@ class RowDiff : public IRowDiff, public BinaryMatrix {
   private:
     static void add_diff(const Vector<uint64_t> &diff, Vector<uint64_t> *row);
 
-    const graph::DBGSuccinct *graph_ = nullptr;
-
     BaseMatrix diffs_;
     anchor_bv_type anchor_;
     fork_succ_bv_type fork_succ_;
@@ -109,6 +108,7 @@ class RowDiff : public IRowDiff, public BinaryMatrix {
 
 template <class BaseMatrix>
 bool RowDiff<BaseMatrix>::get(Row row, Column column) const {
+    assert(graph_ && "graph must be loaded");
     assert(anchor_.size() == diffs_.num_rows() && "anchors must be loaded");
     assert(!fork_succ_.size() || fork_succ_.size() == graph_->get_boss().get_last().size());
 
@@ -123,6 +123,7 @@ bool RowDiff<BaseMatrix>::get(Row row, Column column) const {
  */
 template <class BaseMatrix>
 std::vector<BinaryMatrix::Row> RowDiff<BaseMatrix>::get_column(Column column) const {
+    assert(graph_ && "graph must be loaded");
     assert(anchor_.size() == diffs_.num_rows() && "anchors must be loaded");
     assert(!fork_succ_.size() || fork_succ_.size() == graph_->get_boss().get_last().size());
 
@@ -137,6 +138,7 @@ std::vector<BinaryMatrix::Row> RowDiff<BaseMatrix>::get_column(Column column) co
 
 template <class BaseMatrix>
 BinaryMatrix::SetBitPositions RowDiff<BaseMatrix>::get_row(Row row) const {
+    assert(graph_ && "graph must be loaded");
     assert(anchor_.size() == diffs_.num_rows() && "anchors must be loaded");
     assert(!fork_succ_.size() || fork_succ_.size() == graph_->get_boss().get_last().size());
 
@@ -165,6 +167,7 @@ BinaryMatrix::SetBitPositions RowDiff<BaseMatrix>::get_row(Row row) const {
 template <class BaseMatrix>
 std::vector<BinaryMatrix::SetBitPositions>
 RowDiff<BaseMatrix>::get_rows(const std::vector<Row> &row_ids) const {
+    assert(graph_ && "graph must be loaded");
     assert(anchor_.size() == diffs_.num_rows() && "anchors must be loaded");
     assert(!fork_succ_.size() || fork_succ_.size() == graph_->get_boss().get_last().size());
 
diff --git a/metagraph/src/annotation/int_matrix/row_diff/int_row_diff.hpp b/metagraph/src/annotation/int_matrix/row_diff/int_row_diff.hpp
index 8c84c3789c..5deb819e3f 100644
--- a/metagraph/src/annotation/int_matrix/row_diff/int_row_diff.hpp
+++ b/metagraph/src/annotation/int_matrix/row_diff/int_row_diff.hpp
@@ -35,7 +35,7 @@ class IntRowDiff : public binmat::IRowDiff, public IntMatrix {
     IntRowDiff() {}
 
     IntRowDiff(const graph::DBGSuccinct *graph, BaseMatrix&& diff)
-        : graph_(graph), diffs_(std::move(diff)) {}
+        : diffs_(std::move(diff)) { graph_ = graph; }
 
     bool get(Row i, Column j) const override;
     std::vector<Row> get_column(Column j) const override;
@@ -54,18 +54,14 @@ class IntRowDiff : public binmat::IRowDiff, public IntMatrix {
 
     void load_fork_succ(const std::string &filename);
     void load_anchor(const std::string &filename);
-    void set_graph(const graph::DBGSuccinct *graph) { graph_ = graph; }
 
     const anchor_bv_type& anchor() const { return anchor_; }
-    const graph::DBGSuccinct& graph() const { return *graph_; }
     const BaseMatrix& diffs() const { return diffs_; }
     BaseMatrix& diffs() { return diffs_; }
 
   private:
     static void add_diff(const RowValues &diff, RowValues *row);
 
-    const graph::DBGSuccinct *graph_ = nullptr;
-
     BaseMatrix diffs_;
     anchor_bv_type anchor_;
     fork_succ_bv_type fork_succ_;
@@ -81,6 +77,7 @@ bool IntRowDiff<BaseMatrix>::get(Row i, Column j) const {
 
 template <class BaseMatrix>
 std::vector<IntMatrix::Row> IntRowDiff<BaseMatrix>::get_column(Column j) const {
+    assert(graph_ && "graph must be loaded");
     assert(anchor_.size() == diffs_.num_rows() && "anchors must be loaded");
     assert(!fork_succ_.size() || fork_succ_.size() == graph_->get_boss().get_last().size());
 
@@ -105,6 +102,7 @@ IntMatrix::SetBitPositions IntRowDiff<BaseMatrix>::get_row(Row i) const {
 
 template <class BaseMatrix>
 IntMatrix::RowValues IntRowDiff<BaseMatrix>::get_row_values(Row row) const {
+    assert(graph_ && "graph must be loaded");
     assert(anchor_.size() == diffs_.num_rows() && "anchors must be loaded");
     assert(!fork_succ_.size() || fork_succ_.size() == graph_->get_boss().get_last().size());
 
@@ -150,6 +148,7 @@ IntRowDiff<BaseMatrix>::get_rows(const std::vector<Row> &row_ids) const {
 template <class BaseMatrix>
 std::vector<IntMatrix::RowValues>
 IntRowDiff<BaseMatrix>::get_row_values(const std::vector<Row> &row_ids) const {
+    assert(graph_ && "graph must be loaded");
     assert(anchor_.size() == diffs_.num_rows() && "anchors must be loaded");
     assert(!fork_succ_.size() || fork_succ_.size() == graph_->get_boss().get_last().size());
 
diff --git a/metagraph/src/annotation/row_diff_builder.cpp b/metagraph/src/annotation/row_diff_builder.cpp
index 2c979ed577..1cba34f208 100644
--- a/metagraph/src/annotation/row_diff_builder.cpp
+++ b/metagraph/src/annotation/row_diff_builder.cpp
@@ -8,7 +8,7 @@
 #include "common/threads/threading.hpp"
 #include "common/elias_fano/elias_fano_merger.hpp"
 #include "common/utils/file_utils.hpp"
-#include "common/vectors/bit_vector_sd.hpp"
+#include "common/vectors/bit_vector_sdsl.hpp"
 #include "graph/annotated_dbg.hpp"
 
 const uint64_t BLOCK_SIZE = 1 << 25;
@@ -1006,8 +1006,8 @@ void convert_batch_to_row_diff(const std::string &pred_succ_fprefix,
                     rk++;
                 });
             };
-            columns[j] = std::make_unique<bit_vector_sd>(call_ones, num_rows,
-                                                         row_diff_bits[l_idx][j]);
+            columns[j] = std::make_unique<bit_vector_smart>(call_ones, num_rows,
+                                                            row_diff_bits[l_idx][j]);
         }
 
         if (compute_row_reduction) {
diff --git a/metagraph/src/cli/load/load_annotated_graph.cpp b/metagraph/src/cli/load/load_annotated_graph.cpp
index 353b8df6d8..b59b11575b 100644
--- a/metagraph/src/cli/load/load_annotated_graph.cpp
+++ b/metagraph/src/cli/load/load_annotated_graph.cpp
@@ -41,39 +41,27 @@ std::unique_ptr<AnnotatedDBG> initialize_annotated_dbg(std::shared_ptr<DeBruijnG
                           config.infbase);
             exit(1);
         }
-        const Config::AnnotationType input_anno_type
-                = parse_annotation_type(config.infbase_annotators.at(0));
         // row_diff annotation is special, as it must know the graph structure
-        if (input_anno_type == Config::AnnotationType::RowDiff
-            || input_anno_type == Config::AnnotationType::RowDiffBRWT
-            || input_anno_type == Config::AnnotationType::RowDiffRowSparse) {
-            auto dbg_graph = dynamic_cast<const DBGSuccinct *>(graph.get());
-
-            if (!dbg_graph) {
-                if (auto *canonical = dynamic_cast<CanonicalDBG *>(graph.get()))
-                    dbg_graph = dynamic_cast<const DBGSuccinct *>(&canonical->get_graph());
+        using namespace annot::binmat;
+        BinaryMatrix &matrix = const_cast<BinaryMatrix &>(annotation_temp->get_matrix());
+        if (IRowDiff *row_diff = dynamic_cast<IRowDiff*>(&matrix)) {
+            const DBGSuccinct *dbg_graph;
+            if (auto *canonical = dynamic_cast<CanonicalDBG *>(graph.get())) {
+                dbg_graph = dynamic_cast<const DBGSuccinct *>(&canonical->get_graph());
+            } else {
+                dbg_graph = dynamic_cast<const DBGSuccinct *>(graph.get());
             }
-
             if (!dbg_graph) {
-                logger->error(
-                        "Only succinct de Bruijn graph representations are "
-                        "supported for row-diff annotations");
+                logger->error("Only succinct de Bruijn graph representations"
+                              " are supported for row-diff annotations");
                 std::exit(1);
             }
-            // this is really ugly, but the alternative is to add a set_graph method to
-            // all annotations
-            using namespace annot::binmat;
-            BinaryMatrix &matrix = const_cast<BinaryMatrix &>(annotation_temp->get_matrix());
-            if (input_anno_type == Config::AnnotationType::RowDiff) {
-                dynamic_cast<RowDiff<ColumnMajor> &>(matrix).set_graph(dbg_graph);
-                std::string anchor_fname = config.infbase + kRowDiffAnchorExt;
-                dynamic_cast<RowDiff<ColumnMajor> &>(matrix).load_anchor(anchor_fname);
-                std::string fork_succ_fname = config.infbase + kRowDiffForkSuccExt;
-                dynamic_cast<RowDiff<ColumnMajor> &>(matrix).load_fork_succ(fork_succ_fname);
-            } else if (input_anno_type == Config::AnnotationType::RowDiffRowSparse) {
-                dynamic_cast<RowDiff<RowSparse> &>(matrix).set_graph(dbg_graph);
-            } else {
-                dynamic_cast<RowDiff<BRWT> &>(matrix).set_graph(dbg_graph);
+
+            row_diff->set_graph(dbg_graph);
+
+            if (auto *row_diff_column = dynamic_cast<RowDiff<ColumnMajor> *>(&matrix)) {
+                row_diff_column->load_anchor(config.infbase + kRowDiffAnchorExt);
+                row_diff_column->load_fork_succ(config.infbase + kRowDiffForkSuccExt);
             }
         }
     }

From 6bdd9f63c3f97b458e4985ca3b62725e3c3922b1 Mon Sep 17 00:00:00 2001
From: Mikhail Karasikov <mike-karas@ya.ru>
Date: Wed, 21 Apr 2021 15:43:15 +0200
Subject: [PATCH 12/23] code negative values to pack into smaller integers,
 also other fixes

---
 .../int_matrix/row_diff/int_row_diff.hpp      | 43 ++++++++++++++++++-
 metagraph/src/annotation/row_diff_builder.cpp |  5 ++-
 2 files changed, 44 insertions(+), 4 deletions(-)

diff --git a/metagraph/src/annotation/int_matrix/row_diff/int_row_diff.hpp b/metagraph/src/annotation/int_matrix/row_diff/int_row_diff.hpp
index 5deb819e3f..5c52572d90 100644
--- a/metagraph/src/annotation/int_matrix/row_diff/int_row_diff.hpp
+++ b/metagraph/src/annotation/int_matrix/row_diff/int_row_diff.hpp
@@ -25,6 +25,24 @@ namespace matrix {
 const size_t RD_PATH_RESERVE_SIZE = 2;
 
 
+/**
+ * Convert deltas to positive integer for enable compression:
+ *      0  -> 0
+ *      -1 -> 1
+ *      1  -> 2
+ *      -2 -> 3
+ *      2  -> 4
+ *      ...
+ */
+inline uint64_t encode_diff(int64_t x) {
+    // TODO: skip 0?
+    return x >= 0 ? x * 2 : (-x) * 2 - 1;
+}
+
+inline int64_t decode_diff(uint64_t c) {
+    return !(c & 1) ? c / 2 : -((c + 1) / 2);
+}
+
 template <class BaseMatrix>
 class IntRowDiff : public binmat::IRowDiff, public IntMatrix {
   public:
@@ -60,6 +78,7 @@ class IntRowDiff : public binmat::IRowDiff, public IntMatrix {
     BaseMatrix& diffs() { return diffs_; }
 
   private:
+    static void decode_diffs(RowValues *diffs);
     static void add_diff(const RowValues &diff, RowValues *row);
 
     BaseMatrix diffs_;
@@ -107,6 +126,7 @@ IntMatrix::RowValues IntRowDiff<BaseMatrix>::get_row_values(Row row) const {
     assert(!fork_succ_.size() || fork_succ_.size() == graph_->get_boss().get_last().size());
 
     RowValues result = diffs_.get_row_values(row);
+    decode_diffs(&result);
     std::sort(result.begin(), result.end());
 
     uint64_t boss_edge = graph_->kmer_to_boss_index(
@@ -121,10 +141,15 @@ IntMatrix::RowValues IntRowDiff<BaseMatrix>::get_row_values(Row row) const {
                 graph_->boss_to_kmer_index(boss_edge));
 
         RowValues diff_row = diffs_.get_row_values(row);
+        decode_diffs(&diff_row);
         std::sort(diff_row.begin(), diff_row.end());
         add_diff(diff_row, &result);
     }
 
+    assert(std::all_of(result.begin(), result.end(),
+                       [](auto &p) { return p.second; }));
+    assert(std::all_of(result.begin(), result.end(),
+                       [](auto &p) { return (int64_t)p.second > 0; }));
     return result;
 }
 
@@ -200,6 +225,9 @@ IntRowDiff<BaseMatrix>::get_row_values(const std::vector<Row> &row_ids) const {
     node_to_rd = VectorMap<Row, size_t>();
 
     std::vector<RowValues> rd_rows = diffs_.get_row_values(rd_ids);
+    for (auto &row : rd_rows) {
+        decode_diffs(&row);
+    }
 
     rd_ids = std::vector<Row>();
 
@@ -215,6 +243,10 @@ IntRowDiff<BaseMatrix>::get_row_values(const std::vector<Row> &row_ids) const {
             // replace diff row with full reconstructed annotation
             rd_rows[*it] = result;
         }
+        assert(std::all_of(result.begin(), result.end(),
+                           [](auto &p) { return p.second; }));
+        assert(std::all_of(result.begin(), result.end(),
+                           [](auto &p) { return (int64_t)p.second > 0; }));
     }
 
     return rows;
@@ -235,6 +267,13 @@ void IntRowDiff<BaseMatrix>::serialize(std::ostream &out) const {
     diffs_.serialize(out);
 }
 
+template <class BaseMatrix>
+void IntRowDiff<BaseMatrix>::decode_diffs(RowValues *diffs) {
+    for (auto &[j, value] : *diffs) {
+        value = decode_diff(value);
+    }
+}
+
 template <class BaseMatrix>
 void IntRowDiff<BaseMatrix>::add_diff(const RowValues &diff, RowValues *row) {
     assert(std::is_sorted(row->begin(), row->end()));
@@ -256,8 +295,8 @@ void IntRowDiff<BaseMatrix>::add_diff(const RowValues &diff, RowValues *row) {
             result.push_back(*it2);
             ++it2;
         } else {
-            result.push_back(*it);
-            result.back().second += it2->second;
+            if (uint64_t sum = it->second + it2->second)
+                result.emplace_back(it->first, sum);
             ++it;
             ++it2;
         }
diff --git a/metagraph/src/annotation/row_diff_builder.cpp b/metagraph/src/annotation/row_diff_builder.cpp
index 1cba34f208..178d991e54 100644
--- a/metagraph/src/annotation/row_diff_builder.cpp
+++ b/metagraph/src/annotation/row_diff_builder.cpp
@@ -4,6 +4,7 @@
 #include <progress_bar.hpp>
 
 #include "annotation/binary_matrix/row_diff/row_diff.hpp"
+#include "annotation/int_matrix/row_diff/int_row_diff.hpp"
 #include "annotation/representation/annotation_matrix/static_annotators_def.hpp"
 #include "common/threads/threading.hpp"
 #include "common/elias_fano/elias_fano_merger.hpp"
@@ -946,7 +947,7 @@ void convert_batch_to_row_diff(const std::string &pred_succ_fprefix,
             if constexpr(with_counts) {
                 // diff counts may be negative, hence we need wider integers
                 counts[l_idx][j] = sdsl::int_vector<>(row_diff_bits[l_idx][j], 0,
-                                                      std::min(counts[l_idx][j].width() * 2, 64));
+                                                      std::min(counts[l_idx][j].width() + 1, 64));
             }
 
             auto call_ones = [&](const std::function<void(uint64_t)> &call) {
@@ -1001,7 +1002,7 @@ void convert_batch_to_row_diff(const std::string &pred_succ_fprefix,
                 elias_fano::merge_files<T>(filenames, [&](T v) {
                     call(utils::get_first(v));
                     if constexpr(with_counts) {
-                        counts[l_idx][j][rk] = v.second;
+                        counts[l_idx][j][rk] = matrix::encode_diff(v.second);
                     }
                     rk++;
                 });

From 12bc0aa88b85de68ae8e942897584332af3cb953 Mon Sep 17 00:00:00 2001
From: Mikhail Karasikov <mike-karas@ya.ru>
Date: Wed, 21 Apr 2021 16:12:26 +0200
Subject: [PATCH 13/23] added row-diff-int-brwt tests

---
 metagraph/integration_tests/base.py       | 108 ++++++++++++++++++---
 metagraph/integration_tests/test_query.py | 111 ++++------------------
 2 files changed, 109 insertions(+), 110 deletions(-)

diff --git a/metagraph/integration_tests/base.py b/metagraph/integration_tests/base.py
index ff41c0eea0..f28de6b371 100644
--- a/metagraph/integration_tests/base.py
+++ b/metagraph/integration_tests/base.py
@@ -17,7 +17,8 @@
                         'hashstr': '.hashstrdbg'}
 
 anno_file_extension = {'column': '.column.annodbg',
-                       'row': '.row.annodbg'}
+                       'row': '.row.annodbg',
+                       'row_diff': '.row_diff.annodbg'}
 
 GRAPH_TYPES = [graph_type for graph_type, _ in graph_file_extension.items()]
 
@@ -30,8 +31,8 @@ def setUpClass(cls):
         cls.tempdir = TemporaryDirectory()
 
     @staticmethod
-    def _get_stats(graph_filename):
-        stats_command = METAGRAPH + ' stats ' + graph_filename
+    def _get_stats(graph_path):
+        stats_command = METAGRAPH + ' stats ' + graph_path
         res = subprocess.run(stats_command.split(), stdout=PIPE, stderr=PIPE)
         return res
 
@@ -102,15 +103,92 @@ def _clean(graph, output, extra_params=''):
         assert res.returncode == 0
 
     @staticmethod
-    def _annotate_graph(input, graph_path, output, anno_repr):
-        annotate_command = '{exe} annotate --anno-header -i {graph} \
-                --anno-type {anno_repr} -o {outfile} -p {num_threads} {input}'.format(
-            exe=METAGRAPH,
-            num_threads=NUM_THREADS,
-            graph=graph_path,
-            anno_repr=anno_repr,
-            outfile=output,
-            input=input
-        )
-        res = subprocess.run([annotate_command], shell=True)
-        assert res.returncode == 0
+    def _annotate_graph(input, graph_path, output, anno_repr,
+                        separate=False, no_fork_opt=False, no_anchor_opt=False):
+        # if anno_repr.startswith('int_'):
+        #     # count k-mers
+        #     stats = TestingBase._get_stats(graph_path)
+        #     k = int(stats.stdout.decode().split('\n')[2].split(' ')[1])
+
+        #     tempdir = TemporaryDirectory()
+
+        #     command = f'{METAGRAPH} build -p {NUM_THREADS} --count-kmers \
+        #                 -k {k} -o {tempdir.name}/contigs {input}'
+
+        #     res = subprocess.run([command], shell=True, stdout=PIPE, stderr=PIPE)
+        #     assert res.returncode == 0
+
+        #     input = f'{tempdir.name}/contigs.fasta.gz'
+        #     command = f'{METAGRAPH} clean -p {NUM_THREADS} \
+        #                  --to-fasta -o {input} {tempdir.name}/contigs.dbg'
+
+        #     res = subprocess.run([command], shell=True, stdout=PIPE, stderr=PIPE)
+        #     assert res.returncode == 0
+
+        target_anno = anno_repr
+        if anno_repr in {'row_sparse'} or anno_repr.endswith('brwt') or anno_repr.startswith('row_diff'):
+            target_anno = anno_repr
+            anno_repr = 'column'
+        elif anno_repr in {'flat', 'rbfish'}:
+            target_anno = anno_repr
+            anno_repr = 'row'
+
+        command = f'{METAGRAPH} annotate -p {NUM_THREADS} --anno-header \
+                    -i {graph_path} --anno-type {anno_repr} \
+                    -o {output} {input}'
+
+        if with_counts := target_anno.endswith('int_brwt'):
+            command += ' --count-kmers'
+
+        res = subprocess.run([command], shell=True)
+        assert(res.returncode == 0)
+
+        if target_anno == anno_repr:
+            return
+
+        final_anno = target_anno
+        if final_anno.startswith('row_diff'):
+            target_anno = 'row_diff'
+
+        command = f'{METAGRAPH} transform_anno -p {NUM_THREADS} \
+                    --anno-type {target_anno} -o {output} \
+                    {output + anno_file_extension[anno_repr]}'
+
+        other_args = ' --count-kmers' if with_counts else ''
+
+        if target_anno == 'row_diff':
+            command += ' -i ' + graph_path
+
+        if not no_fork_opt:
+            if target_anno.startswith('row_diff'):
+                print('-- Building RowDiff without fork optimization...')
+            res = subprocess.run([command + other_args], shell=True)
+            assert(res.returncode == 0)
+
+        if target_anno == 'row_diff':
+            without_input_anno = command.split(' ')
+            without_input_anno.pop(-3)
+            without_input_anno = ' '.join(without_input_anno)
+            if not no_anchor_opt:
+                if separate:
+                    print('-- Building RowDiff succ/pred...')
+                    res = subprocess.run(['echo \"\" | ' + without_input_anno + ' --row-diff-stage 1'], shell=True)
+                    assert(res.returncode == 0)
+                res = subprocess.run([command + ' --row-diff-stage 1' + other_args], shell=True)
+                assert(res.returncode == 0)
+            if separate:
+                print('-- Assigning anchors...')
+                res = subprocess.run(['echo \"\" | ' + without_input_anno + ' --row-diff-stage 2'], shell=True)
+                assert(res.returncode == 0)
+            res = subprocess.run([command + ' --row-diff-stage 2' + other_args], shell=True)
+            assert(res.returncode == 0)
+
+            if final_anno != target_anno:
+                rd_type = 'column' if with_counts else 'row_diff'
+                command = f'{METAGRAPH} transform_anno --anno-type {final_anno} --greedy -o {output} ' \
+                                   f'-i {graph_path} -p {NUM_THREADS} {output}.{rd_type}.annodbg'
+                res = subprocess.run([command], shell=True)
+                assert (res.returncode == 0)
+                os.remove(output + anno_file_extension[rd_type])
+            else:
+                os.remove(output + anno_file_extension[anno_repr])
diff --git a/metagraph/integration_tests/test_query.py b/metagraph/integration_tests/test_query.py
index 0732ee4b1a..44ee0c19a0 100644
--- a/metagraph/integration_tests/test_query.py
+++ b/metagraph/integration_tests/test_query.py
@@ -7,20 +7,13 @@
 import glob
 import os
 from helpers import get_test_class_name
+from base import TestingBase, METAGRAPH, TEST_DATA_DIR, graph_file_extension
 
 
 """Test graph construction"""
 
-METAGRAPH = './metagraph'
 DNA_MODE = os.readlink(METAGRAPH).endswith("_DNA")
 PROTEIN_MODE = os.readlink(METAGRAPH).endswith("_Protein")
-TEST_DATA_DIR = os.path.dirname(os.path.realpath(__file__)) + '/../tests/data'
-
-graph_file_extension = {'succinct': '.dbg',
-                        'bitmap': '.bitmapdbg',
-                        'hash': '.orhashdbg',
-                        'hashfast': '.hashfastdbg',
-                        'hashstr': '.hashstrdbg'}
 
 anno_file_extension = {'column': '.column.annodbg',
                        'row': '.row.annodbg',
@@ -31,6 +24,7 @@
                        'rb_brwt': '.rb_brwt.annodbg',
                        'brwt': '.brwt.annodbg',
                        'int_brwt': '.int_brwt.annodbg',
+                       'row_diff_int_brwt': '.row_diff_int_brwt.annodbg',
                        'rbfish': '.rbfish.annodbg',
                        'flat': '.flat.annodbg'}
 
@@ -47,79 +41,6 @@ def product(graph_types, anno_types):
                 result.append((graph, anno))
     return result
 
-def build_annotation(graph_filename, input_fasta, anno_repr, output_filename,
-                     separate=False, no_fork_opt=False, no_anchor_opt=False):
-    target_anno = anno_repr
-    if anno_repr in {'rb_brwt', 'brwt', 'row_sparse'} or anno_repr.startswith('row_diff'):
-        target_anno = anno_repr
-        anno_repr = 'column'
-    elif anno_repr in {'flat', 'rbfish'}:
-        target_anno = anno_repr
-        anno_repr = 'row'
-
-    annotate_command = '{exe} annotate -p {num_threads} --anno-header -i {graph} \
-            --anno-type {anno_repr} -o {outfile} {input}'.format(
-        exe=METAGRAPH,
-        num_threads=NUM_THREADS,
-        graph=graph_filename,
-        anno_repr=anno_repr,
-        outfile=output_filename,
-        input=input_fasta
-    )
-    res = subprocess.run([annotate_command], shell=True)
-    assert(res.returncode == 0)
-
-    if target_anno == anno_repr:
-        return
-
-    final_anno = target_anno
-    if final_anno.startswith('row_diff'):
-        target_anno = 'row_diff'
-
-    annotate_command = '{exe} transform_anno -p {num_threads} \
-            --anno-type {target_anno} -o {outfile} {input}'.format(
-        exe=METAGRAPH,
-        num_threads=NUM_THREADS,
-        graph=graph_filename,
-        target_anno=target_anno,
-        outfile=output_filename,
-        input=output_filename + anno_file_extension[anno_repr]
-    )
-    if target_anno == 'row_diff':
-        annotate_command += ' -i ' + graph_filename
-
-    if not no_fork_opt:
-        print('-- Building RowDiff without fork optimization...')
-        res = subprocess.run([annotate_command], shell=True)
-        assert(res.returncode == 0)
-
-    if target_anno == 'row_diff':
-        without_input_anno = annotate_command.split(' ')
-        without_input_anno.pop(-3)
-        without_input_anno = ' '.join(without_input_anno)
-        if not no_anchor_opt:
-            if separate:
-                print('-- Building RowDiff succ/pred...')
-                res = subprocess.run(['echo \"\" | ' + without_input_anno + ' --row-diff-stage 1'], shell=True)
-                assert(res.returncode == 0)
-            res = subprocess.run([annotate_command + ' --row-diff-stage 1'], shell=True)
-            assert(res.returncode == 0)
-        if separate:
-            print('-- Assigning anchors...')
-            res = subprocess.run(['echo \"\" | ' + without_input_anno + ' --row-diff-stage 2'], shell=True)
-            assert(res.returncode == 0)
-        res = subprocess.run([annotate_command + ' --row-diff-stage 2'], shell=True)
-        assert(res.returncode == 0)
-
-        os.remove(output_filename + anno_file_extension[anno_repr])
-
-        if final_anno != target_anno:
-            annotate_command = f'{METAGRAPH} transform_anno --anno-type {final_anno} --greedy -o {output_filename} ' \
-                               f'-i {graph_filename} -p {NUM_THREADS} {output_filename}.row_diff.annodbg'
-            res = subprocess.run([annotate_command], shell=True)
-            assert (res.returncode == 0)
-            os.remove(output_filename + anno_file_extension['row_diff'])
-
 
 @parameterized_class(('graph_repr', 'anno_repr'),
     input_values=product(
@@ -130,7 +51,7 @@ def build_annotation(graph_filename, input_fasta, anno_repr, output_filename,
     ),
     class_name_func=get_test_class_name
 )
-class TestQuery(unittest.TestCase):
+class TestQuery(TestingBase):
     @classmethod
     def setUpClass(cls):
         cls.tempdir = TemporaryDirectory()
@@ -190,11 +111,11 @@ def check_suffix(anno_repr, suffix):
         cls.anno_repr, no_fork_opt = check_suffix(cls.anno_repr, '_no_fork_opt')
         cls.anno_repr, no_anchor_opt = check_suffix(cls.anno_repr, '_no_anchor_opt')
 
-        build_annotation(
-            cls.tempdir.name + '/graph' + graph_file_extension[cls.graph_repr],
+        cls._annotate_graph(
             TEST_DATA_DIR + '/transcripts_100.fa',
-            cls.anno_repr,
+            cls.tempdir.name + '/graph' + graph_file_extension[cls.graph_repr],
             cls.tempdir.name + '/annotation',
+            cls.anno_repr,
             separate,
             no_fork_opt,
             no_anchor_opt
@@ -591,7 +512,7 @@ def test_batch_query_with_tiny_batch(self):
     class_name_func=get_test_class_name
 )
 @unittest.skipIf(PROTEIN_MODE, "No canonical mode for Protein alphabets")
-class TestQueryCanonical(unittest.TestCase):
+class TestQueryCanonical(TestingBase):
     @classmethod
     def setUpClass(cls):
         cls.tempdir = TemporaryDirectory()
@@ -641,11 +562,11 @@ def setUpClass(cls):
             res = subprocess.run([convert_command], shell=True)
             assert(res.returncode == 0)
 
-        build_annotation(
-            cls.tempdir.name + '/graph' + graph_file_extension[cls.graph_repr],
+        cls._annotate_graph(
             TEST_DATA_DIR + '/transcripts_100.fa',
-            cls.anno_repr,
-            cls.tempdir.name + '/annotation'
+            cls.tempdir.name + '/graph' + graph_file_extension[cls.graph_repr],
+            cls.tempdir.name + '/annotation',
+            cls.anno_repr
         )
 
         # check annotation
@@ -775,7 +696,7 @@ def test_batch_query_with_tiny_batch(self):
     class_name_func=get_test_class_name
 )
 @unittest.skipIf(PROTEIN_MODE, "No canonical mode for Protein alphabets")
-class TestQueryPrimary(unittest.TestCase):
+class TestQueryPrimary(TestingBase):
     @classmethod
     def setUpClass(cls):
         cls.tempdir = TemporaryDirectory()
@@ -849,11 +770,11 @@ def setUpClass(cls):
             res = subprocess.run([convert_command], shell=True)
             assert(res.returncode == 0)
 
-        build_annotation(
-            cls.tempdir.name + '/graph' + graph_file_extension[cls.graph_repr],
+        cls._annotate_graph(
             TEST_DATA_DIR + '/transcripts_100.fa',
-            cls.anno_repr,
-            cls.tempdir.name + '/annotation'
+            cls.tempdir.name + '/graph' + graph_file_extension[cls.graph_repr],
+            cls.tempdir.name + '/annotation',
+            cls.anno_repr
         )
 
         # check annotation

From 8be155481dfbeb8a7df076ca413b1cb9754f28dd Mon Sep 17 00:00:00 2001
From: Mikhail Karasikov <mike-karas@ya.ru>
Date: Wed, 21 Apr 2021 16:14:39 +0200
Subject: [PATCH 14/23] cleanup

---
 metagraph/integration_tests/base.py | 20 --------------------
 1 file changed, 20 deletions(-)

diff --git a/metagraph/integration_tests/base.py b/metagraph/integration_tests/base.py
index f28de6b371..a9b90e3017 100644
--- a/metagraph/integration_tests/base.py
+++ b/metagraph/integration_tests/base.py
@@ -105,26 +105,6 @@ def _clean(graph, output, extra_params=''):
     @staticmethod
     def _annotate_graph(input, graph_path, output, anno_repr,
                         separate=False, no_fork_opt=False, no_anchor_opt=False):
-        # if anno_repr.startswith('int_'):
-        #     # count k-mers
-        #     stats = TestingBase._get_stats(graph_path)
-        #     k = int(stats.stdout.decode().split('\n')[2].split(' ')[1])
-
-        #     tempdir = TemporaryDirectory()
-
-        #     command = f'{METAGRAPH} build -p {NUM_THREADS} --count-kmers \
-        #                 -k {k} -o {tempdir.name}/contigs {input}'
-
-        #     res = subprocess.run([command], shell=True, stdout=PIPE, stderr=PIPE)
-        #     assert res.returncode == 0
-
-        #     input = f'{tempdir.name}/contigs.fasta.gz'
-        #     command = f'{METAGRAPH} clean -p {NUM_THREADS} \
-        #                  --to-fasta -o {input} {tempdir.name}/contigs.dbg'
-
-        #     res = subprocess.run([command], shell=True, stdout=PIPE, stderr=PIPE)
-        #     assert res.returncode == 0
-
         target_anno = anno_repr
         if anno_repr in {'row_sparse'} or anno_repr.endswith('brwt') or anno_repr.startswith('row_diff'):
             target_anno = anno_repr

From 1a919b06173cee75f4a7dafb356231a12d1e0c7a Mon Sep 17 00:00:00 2001
From: Mikhail Karasikov <mike-karas@ya.ru>
Date: Thu, 22 Apr 2021 11:08:51 +0200
Subject: [PATCH 15/23] minor

---
 metagraph/integration_tests/base.py | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/metagraph/integration_tests/base.py b/metagraph/integration_tests/base.py
index a9b90e3017..c655bf7055 100644
--- a/metagraph/integration_tests/base.py
+++ b/metagraph/integration_tests/base.py
@@ -117,7 +117,8 @@ def _annotate_graph(input, graph_path, output, anno_repr,
                     -i {graph_path} --anno-type {anno_repr} \
                     -o {output} {input}'
 
-        if with_counts := target_anno.endswith('int_brwt'):
+        with_counts = target_anno.endswith('int_brwt')
+        if with_counts:
             command += ' --count-kmers'
 
         res = subprocess.run([command], shell=True)

From 1eb2a68bc59b6e8b439a676584abb4a0584b1ea6 Mon Sep 17 00:00:00 2001
From: Mikhail Karasikov <mike-karas@ya.ru>
Date: Mon, 26 Apr 2021 10:53:06 +0200
Subject: [PATCH 16/23] skip empty columns, renamed relation counts -> values

---
 metagraph/src/annotation/row_diff_builder.cpp | 69 ++++++++++---------
 1 file changed, 36 insertions(+), 33 deletions(-)

diff --git a/metagraph/src/annotation/row_diff_builder.cpp b/metagraph/src/annotation/row_diff_builder.cpp
index 178d991e54..cf2e283ce0 100644
--- a/metagraph/src/annotation/row_diff_builder.cpp
+++ b/metagraph/src/annotation/row_diff_builder.cpp
@@ -645,8 +645,8 @@ void convert_batch_to_row_diff(const std::string &pred_succ_fprefix,
                                const std::string &row_reduction_fname,
                                uint64_t buf_size_bytes,
                                bool compute_row_reduction,
-                               bool with_counts) {
-    if (with_counts) {
+                               bool with_values) {
+    if (with_values) {
         convert_batch_to_row_diff<std::pair<uint64_t, uint64_t>>(
                 pred_succ_fprefix, source_files, col_out_dir, swap_dir,
                 row_reduction_fname, buf_size_bytes, compute_row_reduction);
@@ -657,7 +657,7 @@ void convert_batch_to_row_diff(const std::string &pred_succ_fprefix,
     }
 }
 
-// 'T' is either a row index, or a pair (row index, relation count)
+// 'T' is either a row index, or a pair (row index, value)
 template <typename T>
 void convert_batch_to_row_diff(const std::string &pred_succ_fprefix,
                                const std::vector<std::string> &source_files,
@@ -678,32 +678,35 @@ void convert_batch_to_row_diff(const std::string &pred_succ_fprefix,
         return;
     }
 
-    constexpr bool with_counts = utils::is_pair_v<T>;
-    std::vector<std::vector<sdsl::int_vector<>>> counts;
-    if (with_counts) {
-        counts.resize(source_files.size());
+    constexpr bool with_values = utils::is_pair_v<T>;
+    std::vector<std::vector<sdsl::int_vector<>>> values;
+    if (with_values) {
+        values.resize(source_files.size());
         #pragma omp parallel for num_threads(get_num_threads())
         for (size_t i = 0; i < source_files.size(); ++i) {
-            const auto &counts_fname = utils::remove_suffix(source_files[i],
+            if (!sources[i].num_labels())
+                continue;
+
+            const auto &values_fname = utils::remove_suffix(source_files[i],
                                                             ColumnCompressed<>::kExtension)
                                             + ColumnCompressed<>::kCountExtension;
-            std::ifstream instream(counts_fname, std::ios::binary);
-            if (!instream.good()) {
-                logger->error("Could not open file with relation counts {}", counts_fname);
+            std::ifstream instream(values_fname, std::ios::binary);
+            if (!instream) {
+                logger->error("Could not open file with column values {}", values_fname);
                 exit(1);
             }
 
-            counts[i].resize(sources[i].num_labels());
-            for (size_t j = 0; j < counts[i].size(); ++j) {
+            values[i].resize(sources[i].num_labels());
+            for (size_t j = 0; j < values[i].size(); ++j) {
                 try {
-                    counts[i][j].load(instream);
+                    values[i][j].load(instream);
                 } catch (...) {
-                    logger->error("Couldn't read relation counts from {}", counts_fname);
+                    logger->error("Couldn't read column values from {}", values_fname);
                     exit(1);
                 }
             }
         }
-        logger->trace("Done loading relation counts", counts.size());
+        logger->trace("Done loading column values");
     }
 
     anchor_bv_type anchor;
@@ -799,14 +802,14 @@ void convert_batch_to_row_diff(const std::string &pred_succ_fprefix,
     // buffer for writing previous block while populating next row_nbits_block
     std::vector<uint32_t> row_nbits_block_other;
 
-    // get bit at position |i| or its relation count
+    // get bit at position |i| or its value
     auto get_value = [&](const bit_vector &col,
                          size_t s, size_t j, uint64_t i) -> uint64_t {
-        if (!with_counts) {
+        if (!with_values) {
             return col[i];
         } else {
             if (uint64_t rk = col.conditional_rank1(i)) {
-                return counts[s][j][rk - 1];
+                return values[s][j][rk - 1];
             } else {
                 return 0;
             }
@@ -824,10 +827,10 @@ void convert_batch_to_row_diff(const std::string &pred_succ_fprefix,
                     const uint64_t *succ,
                     const uint64_t *pred_begin, const uint64_t *pred_end) {
 
-                // get bits for these positions (or relation counts, hence uint64_t)
+                // get bits for these positions (or values, hence uint64_t)
                 uint64_t curr_value;
-                if constexpr(with_counts) {
-                    curr_value = counts[source_idx][j][source_col.rank1(row_idx) - 1];
+                if constexpr(with_values) {
+                    curr_value = values[source_idx][j][source_col.rank1(row_idx) - 1];
                 } else {
                     curr_value = 1;
                 }
@@ -845,7 +848,7 @@ void convert_batch_to_row_diff(const std::string &pred_succ_fprefix,
                     if (is_anchor || succ_value != curr_value) {
                         // no reduction, we must keep the bit
                         auto &v = set_rows_fwd[source_idx][j];
-                        if constexpr(with_counts) {
+                        if constexpr(with_values) {
                             v.emplace_back(row_idx, curr_value - succ_value);
                         } else {
                             v.push_back(row_idx);
@@ -863,7 +866,7 @@ void convert_batch_to_row_diff(const std::string &pred_succ_fprefix,
                 for (const uint64_t *pred_p = pred_begin; pred_p < pred_end; ++pred_p) {
                     if (curr_value && !source_col[*pred_p] && (compute_row_reduction || !anchor[*pred_p])) {
                         auto &v = set_rows_bwd[source_idx][j];
-                        if constexpr(with_counts) {
+                        if constexpr(with_values) {
                             v.emplace_back(*pred_p, -curr_value);
                         } else {
                             v.push_back(*pred_p);
@@ -944,10 +947,10 @@ void convert_batch_to_row_diff(const std::string &pred_succ_fprefix,
         std::vector<std::unique_ptr<bit_vector>> columns(label_encoders[l_idx].size());
 
         for (size_t j = 0; j < label_encoders[l_idx].size(); ++j) {
-            if constexpr(with_counts) {
-                // diff counts may be negative, hence we need wider integers
-                counts[l_idx][j] = sdsl::int_vector<>(row_diff_bits[l_idx][j], 0,
-                                                      std::min(counts[l_idx][j].width() + 1, 64));
+            if constexpr(with_values) {
+                // diff values may be negative, hence we need wider integers
+                values[l_idx][j] = sdsl::int_vector<>(row_diff_bits[l_idx][j], 0,
+                                                      std::min(values[l_idx][j].width() + 1, 64));
             }
 
             auto call_ones = [&](const std::function<void(uint64_t)> &call) {
@@ -1001,8 +1004,8 @@ void convert_batch_to_row_diff(const std::string &pred_succ_fprefix,
                 uint64_t rk = 0;
                 elias_fano::merge_files<T>(filenames, [&](T v) {
                     call(utils::get_first(v));
-                    if constexpr(with_counts) {
-                        counts[l_idx][j][rk] = matrix::encode_diff(v.second);
+                    if constexpr(with_values) {
+                        values[l_idx][j][rk] = matrix::encode_diff(v.second);
                     }
                     rk++;
                 });
@@ -1016,7 +1019,7 @@ void convert_batch_to_row_diff(const std::string &pred_succ_fprefix,
 
         } else {
             auto fpath = col_out_dir/fs::path(source_files[l_idx]).filename();
-            if constexpr(with_counts) {
+            if constexpr(with_values) {
                 fpath.replace_extension().replace_extension(ColumnCompressed<>::kExtension);
                 ColumnCompressed<>(std::move(columns), label_encoders[l_idx]).serialize(fpath);
                 logger->trace("Serialized {}", fpath);
@@ -1024,10 +1027,10 @@ void convert_batch_to_row_diff(const std::string &pred_succ_fprefix,
                 fpath.replace_extension().replace_extension(ColumnCompressed<>::kCountExtension);
                 std::ofstream outstream(fpath, std::ios::binary);
                 for (size_t j = 0; j < label_encoders[l_idx].size(); ++j) {
-                    counts[l_idx][j].serialize(outstream);
+                    values[l_idx][j].serialize(outstream);
                 }
                 logger->trace("Serialized {}", fpath);
-                counts[l_idx].clear();
+                values[l_idx].clear();
 
             } else {
                 fpath.replace_extension().replace_extension(RowDiffColumnAnnotator::kExtension);

From 63aca3b88fa115407d83f918fe0659d419e5da2a Mon Sep 17 00:00:00 2001
From: Mikhail Karasikov <mike-karas@ya.ru>
Date: Mon, 26 Apr 2021 16:13:05 +0200
Subject: [PATCH 17/23] allow missing values for empty columns

---
 .../src/annotation/annotation_converters.cpp  | 21 ++++++++-----------
 1 file changed, 9 insertions(+), 12 deletions(-)

diff --git a/metagraph/src/annotation/annotation_converters.cpp b/metagraph/src/annotation/annotation_converters.cpp
index b082d25f26..34e46b7d59 100644
--- a/metagraph/src/annotation/annotation_converters.cpp
+++ b/metagraph/src/annotation/annotation_converters.cpp
@@ -1147,7 +1147,7 @@ void convert_to_row_diff(const std::vector<std::string> &files,
                          fs::path swap_dir,
                          RowDiffStage construction_stage,
                          fs::path count_vector_fname,
-                         bool with_counts) {
+                         bool with_values) {
     if (out_dir.empty())
         out_dir = "./";
 
@@ -1185,17 +1185,14 @@ void convert_to_row_diff(const std::vector<std::string> &files,
         for ( ; i < files.size(); ++i) {
             // also add some space for buffers for each column
             uint64_t file_size = fs::file_size(files[i]) + ROW_DIFF_BUFFER_BYTES;
-            if (with_counts) {
+            if (with_values) {
                 // also add k-mer counts
-                const auto &counts_fname
-                    = utils::remove_suffix(files[i], ColumnCompressed<>::kExtension)
-                                                + ColumnCompressed<>::kCountExtension;
-                file_size += fs::file_size(counts_fname);
-                if (!fs::exists(counts_fname)) {
-                    logger->warn("Could not find counts for annotation {}, skipped",
-                                 counts_fname);
-                    continue;
-                }
+                try {
+                    const auto &values_fname
+                        = utils::remove_suffix(files[i], ColumnCompressed<>::kExtension)
+                                                    + ColumnCompressed<>::kCountExtension;
+                    file_size += fs::file_size(values_fname);
+                } catch (...) {}
             }
             if (file_size > mem_bytes) {
                 logger->warn("Not enough memory to process {}, requires {} MB, skipped",
@@ -1228,7 +1225,7 @@ void convert_to_row_diff(const std::vector<std::string> &files,
             convert_batch_to_row_diff(graph_fname,
                     file_batch, out_dir, swap_dir, count_vector_fname, ROW_DIFF_BUFFER_BYTES,
                     construction_stage == RowDiffStage::COMPUTE_REDUCTION,
-                    with_counts);
+                    with_values);
         }
 
         logger->trace("Batch processed in {} sec", timer.elapsed());

From c546c9923f8c713843fba8fd5e5786b07fe499f6 Mon Sep 17 00:00:00 2001
From: Mikhail Karasikov <mike-karas@ya.ru>
Date: Mon, 26 Apr 2021 18:03:02 +0200
Subject: [PATCH 18/23] skip empty annotations

---
 metagraph/src/annotation/row_diff_builder.cpp | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/metagraph/src/annotation/row_diff_builder.cpp b/metagraph/src/annotation/row_diff_builder.cpp
index cf2e283ce0..60c0910123 100644
--- a/metagraph/src/annotation/row_diff_builder.cpp
+++ b/metagraph/src/annotation/row_diff_builder.cpp
@@ -1018,6 +1018,9 @@ void convert_batch_to_row_diff(const std::string &pred_succ_fprefix,
             diff_columns[l_idx] = std::move(columns);
 
         } else {
+            if (!columns.size())
+                continue;
+
             auto fpath = col_out_dir/fs::path(source_files[l_idx]).filename();
             if constexpr(with_values) {
                 fpath.replace_extension().replace_extension(ColumnCompressed<>::kExtension);

From fa8990ede809813e52c3efd08b4d7a740e26e19c Mon Sep 17 00:00:00 2001
From: Mikhail Karasikov <mike-karas@ya.ru>
Date: Wed, 28 Apr 2021 20:48:09 +0200
Subject: [PATCH 19/23] minor

---
 metagraph/src/annotation/annotation_converters.hpp | 4 ++--
 metagraph/src/annotation/row_diff_builder.hpp      | 2 +-
 2 files changed, 3 insertions(+), 3 deletions(-)

diff --git a/metagraph/src/annotation/annotation_converters.hpp b/metagraph/src/annotation/annotation_converters.hpp
index b670cce74f..e3e60570f5 100644
--- a/metagraph/src/annotation/annotation_converters.hpp
+++ b/metagraph/src/annotation/annotation_converters.hpp
@@ -104,7 +104,7 @@ void convert_to_row_annotator(const ColumnCompressed<Label> &annotator,
  * @param out_dir directory where the transformed columns will be dumped. Filenames are
  * kept, extension is changed from 'column.annodbg' to 'row_diff.annodbg'
  * @param swap_dir directory for temporary files
- * @param with_counts row-diff transform for k-mer counts
+ * @param with_values row-diff transform with k-mer counts/attributes
  */
 enum class RowDiffStage { COUNT_LABELS = 0, COMPUTE_REDUCTION, CONVERT };
 void convert_to_row_diff(const std::vector<std::string> &files,
@@ -115,7 +115,7 @@ void convert_to_row_diff(const std::vector<std::string> &files,
                          std::filesystem::path swap_dir,
                          RowDiffStage construction_stage,
                          std::filesystem::path count_vector_fname = "",
-                         bool with_counts = false);
+                         bool with_values = false);
 
 void convert_row_diff_to_col_compressed(const std::vector<std::string> &files,
                                         const std::string &outfbase);
diff --git a/metagraph/src/annotation/row_diff_builder.hpp b/metagraph/src/annotation/row_diff_builder.hpp
index fcbf046f38..c1b36517c5 100644
--- a/metagraph/src/annotation/row_diff_builder.hpp
+++ b/metagraph/src/annotation/row_diff_builder.hpp
@@ -35,7 +35,7 @@ void convert_batch_to_row_diff(const std::string &pred_succ_fprefix,
                                const std::string &row_reduction_fname,
                                uint64_t buf_size_bytes,
                                bool compute_row_reduction = true,
-                               bool with_counts = false);
+                               bool with_values = false);
 
 } // namespace annot
 } // namespace mtg

From 2020bcb4a012a297579167120066f3c9aa0bd795 Mon Sep 17 00:00:00 2001
From: Mikhail Karasikov <mike-karas@ya.ru>
Date: Wed, 28 Apr 2021 21:18:42 +0200
Subject: [PATCH 20/23] skip zero diffs, changed encoding

---
 .../annotation/int_matrix/row_diff/int_row_diff.hpp | 13 +++++++------
 metagraph/src/annotation/row_diff_builder.cpp       |  6 +++---
 2 files changed, 10 insertions(+), 9 deletions(-)

diff --git a/metagraph/src/annotation/int_matrix/row_diff/int_row_diff.hpp b/metagraph/src/annotation/int_matrix/row_diff/int_row_diff.hpp
index 5c52572d90..a58cbfe047 100644
--- a/metagraph/src/annotation/int_matrix/row_diff/int_row_diff.hpp
+++ b/metagraph/src/annotation/int_matrix/row_diff/int_row_diff.hpp
@@ -3,6 +3,7 @@
 
 #include <algorithm>
 #include <iostream>
+#include <cassert>
 #include <string>
 #include <vector>
 
@@ -27,20 +28,20 @@ const size_t RD_PATH_RESERVE_SIZE = 2;
 
 /**
  * Convert deltas to positive integer for enable compression:
- *      0  -> 0
+ *      0  -> X (not allowed, zero diffs must be skipped)
+ *      1  -> 0
  *      -1 -> 1
- *      1  -> 2
+ *      2  -> 2
  *      -2 -> 3
- *      2  -> 4
  *      ...
  */
 inline uint64_t encode_diff(int64_t x) {
-    // TODO: skip 0?
-    return x >= 0 ? x * 2 : (-x) * 2 - 1;
+    assert(x);
+    return x > 0 ? (x - 1) * 2 : (-x) * 2 - 1;
 }
 
 inline int64_t decode_diff(uint64_t c) {
-    return !(c & 1) ? c / 2 : -((c + 1) / 2);
+    return !(c & 1) ? c / 2 + 1 : -((c + 1) / 2);
 }
 
 template <class BaseMatrix>
diff --git a/metagraph/src/annotation/row_diff_builder.cpp b/metagraph/src/annotation/row_diff_builder.cpp
index 60c0910123..ce56d83188 100644
--- a/metagraph/src/annotation/row_diff_builder.cpp
+++ b/metagraph/src/annotation/row_diff_builder.cpp
@@ -843,9 +843,8 @@ void convert_batch_to_row_diff(const std::string &pred_succ_fprefix,
                 } else {
                     // add current bit if this node is an anchor
                     // or if the successor has zero diff
-                    bool is_anchor = anchor[row_idx];
-                    uint64_t succ_value = is_anchor ? 0 : get_value(source_col, source_idx, j, *succ);
-                    if (is_anchor || succ_value != curr_value) {
+                    uint64_t succ_value = anchor[row_idx] ? 0 : get_value(source_col, source_idx, j, *succ);
+                    if (succ_value != curr_value) {
                         // no reduction, we must keep the bit
                         auto &v = set_rows_fwd[source_idx][j];
                         if constexpr(with_values) {
@@ -1005,6 +1004,7 @@ void convert_batch_to_row_diff(const std::string &pred_succ_fprefix,
                 elias_fano::merge_files<T>(filenames, [&](T v) {
                     call(utils::get_first(v));
                     if constexpr(with_values) {
+                        assert(v.second && "zero diffs must have been skipped");
                         values[l_idx][j][rk] = matrix::encode_diff(v.second);
                     }
                     rk++;

From 0e87207090b2e7e554de62d6f5c9fd488fc3610c Mon Sep 17 00:00:00 2001
From: Mikhail Karasikov <mike-karas@ya.ru>
Date: Mon, 3 May 2021 11:46:10 +0200
Subject: [PATCH 21/23] added construction scripts

---
 projects/GTEx/counts.md      | 257 ++++++++++++++++++++++++++++++
 projects/kingsford/README.md | 299 +++++++++++++++++++++++++----------
 2 files changed, 475 insertions(+), 81 deletions(-)
 create mode 100644 projects/GTEx/counts.md

diff --git a/projects/GTEx/counts.md b/projects/GTEx/counts.md
new file mode 100644
index 0000000000..021eb8dc65
--- /dev/null
+++ b/projects/GTEx/counts.md
@@ -0,0 +1,257 @@
+```bash
+mkdir ~/metagenome/data/gtex_counts/build/nobackup;
+mkdir ~/metagenome/data/gtex_counts/build/nobackup/logs;
+mkdir ~/metagenome/data/gtex_counts/build/nobackup/graphs;
+mkdir ~/metagenome/data/gtex_counts/build/nobackup/unitigs;
+
+find /cluster/work/grlab_shared/GTEx/extract/dbGaP-9608/fastq -name "*_1.fastq.gz" \
+    > ~/metagenome/data/gtex_counts/build/nobackup/input_fastq.txt
+
+DIR=~/metagenome/data/gtex_counts/build/nobackup;
+METAGRAPH=~/projects/projects2014-metagenome/metagraph/build_test/metagraph;
+bsub -J "build_graphs[1-9858]" \
+     -o $DIR/logs/build_graphs.lsf \
+     -W 12:00 \
+     -n 8 -R "rusage[mem=10000] span[hosts=1]" \
+        "file=\\\$(sed -n \${LSB_JOBINDEX}p ${DIR}/input_fastq.txt); \
+        file=\\\${file%_1.fastq.gz}; \
+        name=\\\$(basename \\\${file}); \
+        /usr/bin/time -v $METAGRAPH build -v -k 31 \
+            --mode canonical \
+            --count-kmers --count-width 12 \
+            --mem-cap-gb 30 \
+            --disk-swap ~/metagenome/scratch/nobackup/stripe_1 \
+            -p 16 \
+            -o $DIR/graphs/\\\$name \
+            \\\${file}_1.fastq.gz \\\${file}_2.fastq.gz \
+        2>&1"
+
+DIR=~/metagenome/data/gtex_counts/build/nobackup;
+METAGRAPH=~/projects/projects2014-metagenome/metagraph/build_test/metagraph;
+bsub -J "extract_contigs[1-9858]" \
+     -w "exit(build_graphs[*])" \
+     -o $DIR/logs/extract_contigs.lsf \
+     -W 12:00 \
+     -n 8 -R "rusage[mem=10000] span[hosts=1]" \
+        "file=\\\$(sed -n \${LSB_JOBINDEX}p ${DIR}/input_fastq.txt); \
+        file=\\\${file%_1.fastq.gz}; \
+        name=\\\$(basename \\\${file}); \
+        /usr/bin/time -v $METAGRAPH clean \
+            --to-fasta --primary-kmers \
+            --prune-tips 62 \
+            --prune-unitigs 0 \
+            --fallback 2 \
+            -p 16 \
+            -o $DIR/unitigs/\\\$name \
+            $DIR/graphs/\\\$name.dbg \
+        2>&1 && rm \"$DIR/graphs/\\\$name.*\""
+
+DIR=~/metagenome/data/gtex_counts/build/nobackup;
+METAGRAPH=~/projects/projects2014-metagenome/metagraph/build_test/metagraph;
+bsub -J "build_graph" \
+     -o $DIR/logs/build.lsf \
+     -W 48:00 \
+     -n 36 -R "rusage[mem=19000] span[hosts=1]" \
+        "find $DIR/unitigs -name \"*.fasta.gz\" \
+        | /usr/bin/time -v $METAGRAPH build -v \
+            -k 31 \
+            --mode canonical \
+            --mem-cap-gb 100 \
+            --disk-swap ~/metagenome/scratch/nobackup/stripe_1 \
+            -p 72 \
+            -o $DIR/graph \
+        2>&1"
+
+DIR=~/metagenome/data/gtex_counts/build/nobackup;
+METAGRAPH=~/projects/projects2014-metagenome/metagraph/build_test/metagraph;
+bsub -J "extract_primary" \
+     -w "build_graph" \
+     -o $DIR/logs/extract_primary_contigs.lsf \
+     -W 48:00 \
+     -n 36 -R "rusage[mem=19000] span[hosts=1]" \
+        "/usr/bin/time -v $METAGRAPH transform -v \
+            --to-fasta \
+            --primary-kmers \
+            -p 72 \
+            -o $DIR/primary_contigs \
+            $DIR/graph.dbg \
+        2>&1"
+
+DIR=~/metagenome/data/gtex_counts/build/nobackup;
+METAGRAPH=~/projects/projects2014-metagenome/metagraph/build_test/metagraph;
+bsub -J "build_primary" \
+     -w "extract_primary" \
+     -o $DIR/logs/build_primary.lsf \
+     -W 48:00 \
+     -n 36 -R "rusage[mem=19000] span[hosts=1]" \
+        "/usr/bin/time -v $METAGRAPH build -v \
+            -k 31 \
+            --mode primary \
+            --mem-cap-gb 100 \
+            --disk-swap ~/metagenome/scratch/nobackup/stripe_1 \
+            -p 72 \
+            -o $DIR/graph_primary \
+            $DIR/primary_contigs.fasta.gz \
+        2>&1"
+
+DIR=~/metagenome/data/gtex_counts/build/nobackup;
+METAGRAPH=~/projects/projects2014-metagenome/metagraph/build_test/metagraph;
+find $DIR/unitigs -name "*.fasta.gz" | shuf > ${DIR}/unitigs.txt;
+
+bsub -J "build_clean_graphs[1-9759]" \
+     -o $DIR/logs/build_clean_graphs.lsf \
+     -W 4:00 \
+     -n 4 -R "rusage[mem=10000] span[hosts=1]" \
+        "file=\\\$(sed -n \${LSB_JOBINDEX}p ${DIR}/unitigs.txt); \
+        file=\\\${file%.fasta.gz}; \
+        name=\\\$(basename \\\${file}); \
+        /usr/bin/time -v $METAGRAPH build -v -k 31 \
+            --mode canonical \
+            --count-kmers --count-width 12 \
+            --mem-cap-gb 30 \
+            --disk-swap ~/metagenome/scratch/nobackup/stripe_1 \
+            -p 8 \
+            -o $DIR/graphs/\\\$name \
+            \\\${file}.fasta.gz \
+        2>&1"
+
+WINDOW_SIZE=31;
+DIR=~/metagenome/data/gtex_counts/build/nobackup/smoothing_${WINDOW_SIZE};
+mkdir ${DIR};
+mkdir ${DIR}/logs;
+mkdir ${DIR}/unitigs;
+
+DIR=~/metagenome/data/gtex_counts/build/nobackup/smoothing_${WINDOW_SIZE};
+METAGRAPH=~/projects/projects2014-metagenome/metagraph/build_test/metagraph;
+bsub -J "extract_contigs[1-9759]%500" \
+     -o $DIR/logs/extract_contigs.lsf \
+     -W 4:00 \
+     -n 4 -R "rusage[mem=10000] span[hosts=1]" \
+        "file=\\\$(sed -n \${LSB_JOBINDEX}p ${DIR}/../unitigs.txt); \
+        file=\\\${file%.fasta.gz}; \
+        name=\\\$(basename \\\${file}); \
+        /usr/bin/time -v $METAGRAPH clean \
+            --to-fasta --primary-kmers \
+            --smoothing-window ${WINDOW_SIZE} \
+            -p 8 \
+            -o $DIR/unitigs/\\\$name \
+            $DIR/../graphs/\\\${name}.dbg \
+        2>&1"
+
+DIR=~/metagenome/data/gtex_counts/build/nobackup/smoothing_${WINDOW_SIZE};
+cd $DIR;
+mkdir -p batches;
+cd batches;
+split -d -n r/40 <(find $DIR/unitigs -name "*.fasta.gz" | shuf);
+mkdir -p ${DIR}/columns;
+
+DIR=~/metagenome/data/gtex_counts/build/nobackup/smoothing_${WINDOW_SIZE};
+METAGRAPH=~/projects/projects2014-metagenome/metagraph/build_test/metagraph;
+for N in {0..39}; do
+    N=$(printf "%02d" $N);
+    list=x$N;
+    bsub -J "gtex_annotate_${WINDOW_SIZE}_${list}" \
+         -oo ${DIR}/logs/annotate_${list}.lsf \
+         -W 4:00 \
+         -n 18 -R "rusage[mem=1500] span[hosts=1]" \
+        "cat $DIR/batches/${list} \
+            | /usr/bin/time -v $METAGRAPH annotate \
+                -i $DIR/../graph_primary.dbg \
+                --anno-filename \
+                --separately \
+                --count-kmers --count-width 12 \
+                -o ${DIR}/columns \
+                -p 36 \
+                2>&1 | tee ${DIR}/logs/annotate_${list}.log"; \
+done
+
+
+DIR=~/metagenome/data/gtex_counts/build/nobackup/smoothing_${WINDOW_SIZE};
+mkdir $DIR/rd;
+mkdir $DIR/rd/rd_columns;
+ln -s $DIR/../graph_primary.dbg ${DIR}/rd/graph.dbg;
+
+DIR=~/metagenome/data/gtex_counts/build/nobackup/smoothing_${WINDOW_SIZE};
+METAGRAPH=~/projects/projects2014-metagenome/metagraph/build_test/metagraph;
+bsub -J "gtex_count_${WINDOW_SIZE}_rd_0" \
+     -w "gtex_annotate_${WINDOW_SIZE}_*" \
+     -oo ${DIR}/logs/count_rd_0.lsf \
+     -W 24:00 \
+     -n 36 -R "rusage[mem=19000] span[hosts=1]" \
+    "find ${DIR}/columns -name \"*.column.annodbg\" \
+        | /usr/bin/time -v $METAGRAPH transform_anno -v \
+            --anno-type row_diff --count-kmers \
+            --row-diff-stage 0 \
+            --mem-cap-gb 500 \
+            --disk-swap ~/metagenome/scratch/nobackup/stripe_1 \
+            -i ${DIR}/rd/graph.dbg \
+            -o ${DIR}/rd/rd_columns/out \
+            -p 72 \
+            2>&1 | tee ${DIR}/logs/count_rd_0.log";
+
+DIR=~/metagenome/data/gtex_counts/build/nobackup/smoothing_${WINDOW_SIZE};
+METAGRAPH=~/projects/projects2014-metagenome/metagraph/build_test/metagraph;
+bsub -J "gtex_count_${WINDOW_SIZE}_rd_1" \
+     -w "gtex_count_${WINDOW_SIZE}_rd_0" \
+     -oo ${DIR}/logs/count_rd_1.lsf \
+     -W 24:00 \
+     -n 36 -R "rusage[mem=19000] span[hosts=1]" \
+    "find ${DIR}/columns -name \"*.column.annodbg\" \
+        | /usr/bin/time -v $METAGRAPH transform_anno -v \
+            --anno-type row_diff --count-kmers \
+            --row-diff-stage 1 \
+            --mem-cap-gb 500 \
+            --disk-swap ~/metagenome/scratch/nobackup/stripe_1 \
+            -i ${DIR}/rd/graph.dbg \
+            -o ${DIR}/rd/rd_columns/out \
+            -p 72 \
+            2>&1 | tee ${DIR}/logs/count_rd_1.log";
+
+DIR=~/metagenome/data/gtex_counts/build/nobackup/smoothing_${WINDOW_SIZE};
+METAGRAPH=~/projects/projects2014-metagenome/metagraph/build_test/metagraph;
+bsub -J "gtex_count_${WINDOW_SIZE}_rd_2" \
+     -w "gtex_count_${WINDOW_SIZE}_rd_1" \
+     -oo ${DIR}/logs/count_rd_2.lsf \
+     -W 24:00 \
+     -n 36 -R "rusage[mem=19000] span[hosts=1]" \
+    "find ${DIR}/columns -name \"*.column.annodbg\" \
+        | /usr/bin/time -v $METAGRAPH transform_anno -v \
+            --anno-type row_diff --count-kmers \
+            --row-diff-stage 2 \
+            --mem-cap-gb 500 \
+            --disk-swap ~/metagenome/scratch/nobackup/stripe_1 \
+            -i ${DIR}/rd/graph.dbg \
+            -o ${DIR}/rd/rd_columns/out \
+            -p 72 \
+            2>&1 | tee ${DIR}/logs/count_rd_2.log";
+
+DIR=~/metagenome/data/gtex_counts/build/nobackup/smoothing_${WINDOW_SIZE};
+METAGRAPH=~/projects/projects2014-metagenome/metagraph/build_test/metagraph;
+bsub -J "gtex_count_${WINDOW_SIZE}_rd_brwt" \
+     -w "gtex_count_${WINDOW_SIZE}_rd_2" \
+     -oo ${DIR}/logs/count_rd_brwt.lsf \
+     -W 24:00 \
+     -n 36 -R "rusage[mem=19000] span[hosts=1]" \
+    "find ${DIR}/rd/rd_columns -name \"*.column.annodbg\" \
+        | /usr/bin/time -v $METAGRAPH transform_anno -v \
+            --anno-type row_diff_int_brwt \
+            --greedy --fast --subsample 1000000 \
+            -i ${DIR}/rd/graph.dbg \
+            -o ${DIR}/annotation \
+            -p 72 --parallel-nodes 10 \
+            2>&1 | tee ${DIR}/logs/count_rd_brwt.log";
+
+DIR=~/metagenome/data/gtex_counts/build/nobackup/smoothing_${WINDOW_SIZE};
+METAGRAPH=~/projects/projects2014-metagenome/metagraph/build_test/metagraph;
+bsub -J "gtex_count_${WINDOW_SIZE}_rd_brwt_relax" \
+     -w "gtex_count_${WINDOW_SIZE}_rd_brwt" \
+     -oo ${DIR}/logs/count_rd_brwt_relax.lsf \
+     -W 24:00 \
+     -n 16 -R "rusage[mem=20000] span[hosts=1]" \
+    "/usr/bin/time -v $METAGRAPH relax_brwt -v \
+            -p 24 \
+            --relax-arity 32 \
+            -o ${DIR}/annotation.relaxed \
+            ${DIR}/annotation.row_diff_int_brwt.annodbg \
+            2>&1 | tee ${DIR}/logs/count_rd_brwt_relax.log";
+```
diff --git a/projects/kingsford/README.md b/projects/kingsford/README.md
index 07da3da65d..531e37746b 100644
--- a/projects/kingsford/README.md
+++ b/projects/kingsford/README.md
@@ -4,61 +4,69 @@ find ~/metagenome/data/kingsford/kmc_20/ -name "*kmc_suf" > ~/metagenome/data/ki
 bsub -J "build_single[1-2652]%500" \
      -o ~/metagenome/data/kingsford/build_single.lsf \
      -W 4:00 \
-     -n 1 -R "rusage[mem=10000] span[hosts=1]" \
-        "file=\"\$(sed -n \${LSB_JOBINDEX}p ~/metagenome/data/kingsford/kingsford_list_kmc.txt)\"; \
+     -n 4 -R "rusage[mem=10000] span[hosts=1]" \
+        "file=\\\$(sed -n \${LSB_JOBINDEX}p ~/metagenome/data/kingsford/kingsford_list_kmc.txt); \
         /usr/bin/time -v ~/projects/projects2014-metagenome/metagraph/build_test/metagraph_DNA build -v \
             -k 20 \
             --mode canonical \
-            --count-kmers \
+            --count-kmers --count-width 32 \
             --mem-cap-gb 8 \
-            -p 1 \
-            -o ~/metagenome/data/kingsford/single_graphs/\$(basename \${file%.fasta.gz.kmc.kmc_suf}) \
-            \$file \
+            -p 8 \
+            -o ~/metagenome/data/kingsford/single_graphs/\\\$(basename \\\${file%.fasta.gz.kmc.kmc_suf}) \
+            \\\$file \
         2>&1"
 
-bsub -w "done(build_single[*])" \
-     -J "extract_contigs[1-2652]%500" \
-     -o ~/metagenome/data/kingsford/extract_contigs.lsf \
+mkdir ~/metagenome/data/kingsford/smoothing_1000;
+mkdir ~/metagenome/data/kingsford/smoothing_1000/logs;
+mkdir ~/metagenome/data/kingsford/smoothing_1000/unitigs;
+
+DIR=~/metagenome/data/kingsford/smoothing_1000;
+METAGRAPH=~/projects/projects2014-metagenome/metagraph/build_test2/metagraph;
+bsub -J "extract_contigs[1-2652]%500" \
+     -o $DIR/logs/extract_contigs.lsf \
      -W 4:00 \
-     -n 1 -R "rusage[mem=2000] span[hosts=1]" \
-        "file=\"\$(sed -n \${LSB_JOBINDEX}p ~/metagenome/data/kingsford/kingsford_list_kmc.txt)\"; \
-        file=\$(basename \${file%.fasta.gz.kmc.kmc_suf})
-        /usr/bin/time -v ~/projects/projects2014-metagenome/metagraph/build_test/metagraph_DNA clean -v \
+     -n 4 -R "rusage[mem=2000] span[hosts=1]" \
+        "file=\\\$(sed -n \${LSB_JOBINDEX}p ~/metagenome/data/kingsford/kingsford_list_kmc.txt); \
+        file=\\\$(basename \\\${file%.fasta.gz.kmc.kmc_suf})
+        /usr/bin/time -v $METAGRAPH clean -v \
             --to-fasta --primary-kmers \
-            -p 1 \
-            -o ~/metagenome/data/kingsford/single_graphs/\$file \
-            ~/metagenome/data/kingsford/single_graphs/\${file}.dbg \
+            --smoothing-window 1000 \
+            -p 8 \
+            -o $DIR/unitigs/\\\$file \
+            ~/metagenome/data/kingsford/single_graphs/\\\${file}.dbg \
         2>&1"
 
 bsub -J "kingsford_build" \
-     -oo ~/metagenome/data/kingsford/build_canonical.lsf \
+     -oo $DIR/logs/build_canonical.lsf \
      -W 4:00 \
      -n 36 -R "rusage[mem=3000] span[hosts=1]" \
-    "find ~/metagenome/data/kingsford/single_graphs/ -name \"*.fasta.gz\" \
-        | gtime -v ~/projects/projects2014-metagenome/metagraph/build_test/metagraph build -v \
+    "find $DIR/unitigs -name \"*.fasta.gz\" \
+        | /usr/bin/time -v $METAGRAPH build -v \
             -k 20 \
             --mode canonical \
             --mem-cap-gb 80 \
-            -p 36 \
-            -o ~/metagenome/data/kingsford/kingsford_canonical \
+            -p 72 \
+            -o $DIR/kingsford_canonical \
             2>&1"
 
 bsub -J "kingsford_to_primary" \
-     -oo ~/metagenome/data/kingsford/extract_primary_contigs.lsf \
+     -w "kingsford_build" \
+     -oo $DIR/logs/extract_primary_contigs.lsf \
      -W 4:00 \
      -n 36 -R "rusage[mem=2000] span[hosts=1]" \
-    "gtime -v ~/projects/projects2014-metagenome/metagraph/build_test/metagraph transform -v \
+    "/usr/bin/time -v $METAGRAPH transform -v \
             --to-fasta --primary-kmers \
-            -o ~/metagenome/data/kingsford/kingsford_primary \
-            ~/metagenome/data/kingsford/kingsford_canonical.dbg \
-            -p 36 \
+            -o $DIR/kingsford \
+            $DIR/kingsford_canonical.dbg \
+            -p 72 \
             2>&1"
 
 bsub -J "kingsford_build" \
+     -w "kingsford_to_primary" \
      -oo ~/metagenome/data/kingsford/build_primary.lsf \
      -W 4:00 \
      -n 36 -R "rusage[mem=1500] span[hosts=1]" \
-    "gtime -v ~/projects/projects2014-metagenome/metagraph/build_test/metagraph build -v \
+    "/usr/bin/time -v ~/projects/projects2014-metagenome/metagraph/build_test/metagraph build -v \
             -k 20 \
             --mode primary \
             --mem-cap-gb 40 \
@@ -68,76 +76,205 @@ bsub -J "kingsford_build" \
             2>&1"
 
 
-mkdir temp
-cd temp
-find ~/metagenome/data/kingsford/single_graphs/ -name "*.fasta.gz" > files_to_annotate.txt
-split -l 150 files_to_annotate.txt
+DIR=~/metagenome/data/kingsford/smoothing_1000;
+cd $DIR
+mkdir batches
+cd batches
+split -l 150 <(find $DIR/unitigs -name "*.fasta.gz" | shuf)
+mkdir ${DIR}/columns
 
+DIR=~/metagenome/data/kingsford/smoothing_1000;
+METAGRAPH=~/projects/projects2014-metagenome/metagraph/build_test/metagraph;
+cd $DIR/batches;
 for list in x*; do
     bsub -J "kingsford_annotate_${list}" \
-         -oo ~/metagenome/data/kingsford/annotation/annotate_primary_${list}.lsf \
+         -oo ${DIR}/logs/annotate_${list}.lsf \
          -W 4:00 \
-         -n 15 -R "rusage[mem=1000] span[hosts=1]" \
+         -n 18 -R "rusage[mem=1500] span[hosts=1]" \
         "cat ${list} \
-            | /usr/bin/time -v ~/projects/projects2014-metagenome/metagraph/build_test/metagraph annotate -v \
+            | /usr/bin/time -v $METAGRAPH annotate -v \
                 -i ~/metagenome/data/kingsford/kingsford_primary.dbg \
-                --mode canonical \
-                --count-kmers \
-                --parallel 15 \
                 --anno-filename \
                 --separately \
-                -o ~/metagenome/data/kingsford/annotation/columns_primary \
-                2>&1"; \
+                -o ${DIR}/columns \
+                -p 36 \
+                2>&1 | tee ${DIR}/logs/annotate_${list}.log"; \
 done
-cd ..
-
-find ~/metagenome/data/kingsford/annotation/columns_primary/ -name "*.column.annodbg" \
-    > ~/metagenome/data/kingsford/annotation/columns_primary.txt
 
-bsub -J "kingsford_build_rbbrwt_relax20_primary" \
-     -oo ~/metagenome/data/kingsford/annotation/build_rbbrwt_primary.lsf \
-     -W 24:00 \
-     -n 36 -R "rusage[mem=6300] span[hosts=1]" \
-    "cat ~/metagenome/data/kingsford/annotation/columns_primary.txt \
-        | gtime -v ~/projects/projects2014-metagenome/metagraph/build_test/metagraph transform_anno -v \
-            --anno-type rb_brwt --relax-arity 20 \
-            -o ~/metagenome/data/kingsford/annotation/annotation_primary \
-            -p 36 \
-            2>&1"
-
-bsub -J "kingsford_cluster_columns_primary" \
-     -oo ~/metagenome/data/kingsford/annotation/cluster_columns_primary.lsf \
-     -W 4:00 \
-     -n 36 -R "rusage[mem=2100] span[hosts=1]" \
-    "cat ~/metagenome/data/kingsford/annotation/columns_primary.txt \
-        | /usr/bin/time -v ~/projects/projects2014-metagenome/metagraph/build_test/metagraph transform_anno -v \
-            --linkage --greedy \
-            --subsample 100000000 \
-            -o ~/metagenome/data/kingsford/annotation/linkage_kingsford_primary.csv \
-            -p 36 \
-            2>&1";
 
-bsub -J "kingsford_build_brwt_primary" \
-     -oo ~/metagenome/data/kingsford/annotation/build_brwt_primary.lsf \
+DIR=~/metagenome/data/kingsford/annotation;
+METAGRAPH=~/projects/projects2014-metagenome/metagraph/build_test/metagraph;
+cd $DIR
+bsub -J "kingsford_brwt" \
+     -oo ${DIR}/logs/brwt.lsf \
      -W 4:00 \
      -n 36 -R "rusage[mem=1000] span[hosts=1]" \
-    "cat ~/metagenome/data/kingsford/annotation/columns_primary.txt \
-        | gtime -v ~/projects/projects2014-metagenome/metagraph/build_test/metagraph transform_anno -v \
+    "find ${DIR}/columns_primary -name \"*.column.annodbg\" \
+        | /usr/bin/time -v $METAGRAPH transform_anno -v \
             --anno-type brwt \
-            --linkage-file ~/metagenome/data/kingsford/annotation/linkage_kingsford_primary.csv \
-            -o ~/metagenome/data/kingsford/annotation/annotation_primary \
-            -p 36 --parallel-nodes 10 \
-            2>&1"
+            --greedy --fast --subsample 1000000 \
+            -o ${DIR}/annotation \
+            -p 72 --parallel-nodes 10 \
+            2>&1 | tee ${DIR}/logs/brwt.log"
 
-bsub -J "kingsford_relax_brwt_primary" \
-     -oo ~/metagenome/data/kingsford/annotation/relax_brwt_primary.lsf \
+bsub -J "kingsford_brwt_relax" \
+     -w "kingsford_brwt" \
+     -oo ${DIR}/logs/brwt_relax.lsf \
      -W 4:00 \
      -n 36 -R "rusage[mem=1000] span[hosts=1]" \
-    "gtime -v ~/projects/projects2014-metagenome/metagraph/build_test/metagraph relax_brwt -v \
-        -p 36 \
-        --relax-arity 16 \
-        -o ~/metagenome/data/kingsford/annotation/annotation_primary.relaxed \
-        ~/metagenome/data/kingsford/annotation/annotation_primary.brwt.annodbg \
-        2>&1"
+    "/usr/bin/time -v $METAGRAPH relax_brwt -v \
+        -p 72 \
+        --relax-arity 32 \
+        -o ${DIR}/annotation.relaxed \
+        ${DIR}/annotation.brwt.annodbg \
+        2>&1 | tee ${DIR}/logs/brwt_relax.log"
 ```
 
+
+
+
+
+## Index with k-mer counts
+
+```bash
+WINDOW_SIZE=1;
+mkdir ~/metagenome/data/kingsford/nobackup/smoothing_${WINDOW_SIZE};
+mkdir ~/metagenome/data/kingsford/nobackup/smoothing_${WINDOW_SIZE}/logs;
+mkdir ~/metagenome/data/kingsford/nobackup/smoothing_${WINDOW_SIZE}/unitigs;
+
+DIR=~/metagenome/data/kingsford/nobackup/smoothing_${WINDOW_SIZE};
+METAGRAPH=~/projects/projects2014-metagenome/metagraph/build_test/metagraph;
+bsub -J "extract_contigs[1-2652]%500" \
+     -o $DIR/logs/extract_contigs.lsf \
+     -W 4:00 \
+     -n 4 -R "rusage[mem=2000] span[hosts=1]" \
+        "file=\\\$(sed -n \${LSB_JOBINDEX}p ~/metagenome/data/kingsford/kingsford_list_kmc.txt); \
+        file=\\\$(basename \\\${file%.fasta.gz.kmc.kmc_suf})
+        /usr/bin/time -v $METAGRAPH clean \
+            --to-fasta --primary-kmers \
+            --smoothing-window ${WINDOW_SIZE} \
+            -p 8 \
+            -o $DIR/unitigs/\\\$file \
+            ~/metagenome/data/kingsford/single_graphs/\\\${file}.dbg \
+        2>&1"
+
+DIR=~/metagenome/data/kingsford/nobackup/smoothing_${WINDOW_SIZE};
+bsub -J "split" \
+     -w "extract_contigs" \
+     -o /dev/null -W 1:00 -n 1 -R "rusage[mem=1000]" \
+        "cd $DIR; \
+        mkdir -p batches; \
+        cd batches; \
+        split -d -n r/40 <(find $DIR/unitigs -name "*.fasta.gz" | shuf); \
+        mkdir -p ${DIR}/columns;";
+
+DIR=~/metagenome/data/kingsford/nobackup/smoothing_${WINDOW_SIZE};
+METAGRAPH=~/projects/projects2014-metagenome/metagraph/build_test/metagraph;
+for N in {0..39}; do
+    N=$(printf "%02d" $N);
+    list=x$N;
+    bsub -J "kingsford_annotate_${list}" \
+         -w "split" \
+         -oo ${DIR}/logs/annotate_${list}.lsf \
+         -W 4:00 \
+         -n 18 -R "rusage[mem=1500] span[hosts=1]" \
+        "cat $DIR/batches/${list} \
+            | /usr/bin/time -v $METAGRAPH annotate \
+                -i ~/metagenome/data/kingsford/kingsford_primary.dbg \
+                --anno-filename \
+                --separately \
+                --count-kmers --count-width 32 \
+                -o ${DIR}/columns \
+                -p 36 \
+                2>&1 | tee ${DIR}/logs/annotate_${list}.log"; \
+done
+
+DIR=~/metagenome/data/kingsford/nobackup/smoothing_${WINDOW_SIZE};
+mkdir $DIR/rd;
+mkdir $DIR/rd/rd_columns;
+ln -s ~/metagenome/data/kingsford/kingsford_primary.dbg ${DIR}/rd/graph.dbg;
+
+DIR=~/metagenome/data/kingsford/nobackup/smoothing_${WINDOW_SIZE};
+METAGRAPH=~/projects/projects2014-metagenome/metagraph/build_test/metagraph;
+bsub -J "kingsford_count_rd_0" \
+     -w "kingsford_annotate_*" \
+     -oo ${DIR}/logs/count_rd_0.lsf \
+     -W 24:00 \
+     -n 36 -R "rusage[mem=19000] span[hosts=1]" \
+    "find ${DIR}/columns -name \"*.column.annodbg\" \
+        | /usr/bin/time -v $METAGRAPH transform_anno -v \
+            --anno-type row_diff --count-kmers \
+            --row-diff-stage 0 \
+            --mem-cap-gb 500 \
+            --disk-swap ~/metagenome/scratch/nobackup/stripe_1 \
+            -i ${DIR}/rd/graph.dbg \
+            -o ${DIR}/rd/rd_columns/out \
+            -p 72 \
+            2>&1 | tee ${DIR}/logs/count_rd_0.log";
+
+DIR=~/metagenome/data/kingsford/nobackup/smoothing_${WINDOW_SIZE};
+METAGRAPH=~/projects/projects2014-metagenome/metagraph/build_test/metagraph;
+bsub -J "kingsford_count_rd_1" \
+     -w "kingsford_count_rd_0" \
+     -oo ${DIR}/logs/count_rd_1.lsf \
+     -W 24:00 \
+     -n 36 -R "rusage[mem=19000] span[hosts=1]" \
+    "find ${DIR}/columns -name \"*.column.annodbg\" \
+        | /usr/bin/time -v $METAGRAPH transform_anno -v \
+            --anno-type row_diff --count-kmers \
+            --row-diff-stage 1 \
+            --mem-cap-gb 500 \
+            --disk-swap ~/metagenome/scratch/nobackup/stripe_1 \
+            -i ${DIR}/rd/graph.dbg \
+            -o ${DIR}/rd/rd_columns/out \
+            -p 72 \
+            2>&1 | tee ${DIR}/logs/count_rd_1.log";
+
+DIR=~/metagenome/data/kingsford/nobackup/smoothing_${WINDOW_SIZE};
+METAGRAPH=~/projects/projects2014-metagenome/metagraph/build_test/metagraph;
+bsub -J "kingsford_count_rd_2" \
+     -w "kingsford_count_rd_1" \
+     -oo ${DIR}/logs/count_rd_2.lsf \
+     -W 24:00 \
+     -n 36 -R "rusage[mem=19000] span[hosts=1]" \
+    "find ${DIR}/columns -name \"*.column.annodbg\" \
+        | /usr/bin/time -v $METAGRAPH transform_anno -v \
+            --anno-type row_diff --count-kmers \
+            --row-diff-stage 2 \
+            --mem-cap-gb 500 \
+            --disk-swap ~/metagenome/scratch/nobackup/stripe_1 \
+            -i ${DIR}/rd/graph.dbg \
+            -o ${DIR}/rd/rd_columns/out \
+            -p 72 \
+            2>&1 | tee ${DIR}/logs/count_rd_2.log";
+
+DIR=~/metagenome/data/kingsford/nobackup/smoothing_${WINDOW_SIZE};
+METAGRAPH=~/projects/projects2014-metagenome/metagraph/build_test/metagraph;
+bsub -J "kingsford_count_rd_brwt" \
+     -w "kingsford_count_rd_2" \
+     -oo ${DIR}/logs/count_rd_brwt.lsf \
+     -W 24:00 \
+     -n 36 -R "rusage[mem=4000] span[hosts=1]" \
+    "find ${DIR}/rd/rd_columns -name \"*.column.annodbg\" \
+        | /usr/bin/time -v $METAGRAPH transform_anno -v \
+            --anno-type row_diff_int_brwt \
+            --greedy --fast --subsample 1000000 \
+            -i ${DIR}/rd/graph.dbg \
+            -o ${DIR}/annotation \
+            -p 72 --parallel-nodes 10 \
+            2>&1 | tee ${DIR}/logs/count_rd_brwt.log";
+
+DIR=~/metagenome/data/kingsford/nobackup/smoothing_${WINDOW_SIZE};
+METAGRAPH=~/projects/projects2014-metagenome/metagraph/build_test/metagraph;
+bsub -J "kingsford_count_rd_brwt_relax" \
+     -w "kingsford_count_rd_brwt" \
+     -oo ${DIR}/logs/count_rd_brwt_relax.lsf \
+     -W 24:00 \
+     -n 12 -R "rusage[mem=5000] span[hosts=1]" \
+    "/usr/bin/time -v $METAGRAPH relax_brwt -v \
+            -p 24 \
+            --relax-arity 32 \
+            -o ${DIR}/annotation.relaxed \
+            ${DIR}/annotation.row_diff_int_brwt.annodbg \
+            2>&1 | tee ${DIR}/logs/count_rd_brwt_relax.log";
+```

From ec4aa467064b31ea208d43ecde97135e4f504976 Mon Sep 17 00:00:00 2001
From: Mikhail Karasikov <mike-karas@ya.ru>
Date: Mon, 3 May 2021 15:30:59 +0200
Subject: [PATCH 22/23] added comment

---
 metagraph/src/annotation/annotation_converters.cpp | 6 +++++-
 1 file changed, 5 insertions(+), 1 deletion(-)

diff --git a/metagraph/src/annotation/annotation_converters.cpp b/metagraph/src/annotation/annotation_converters.cpp
index 34e46b7d59..07b1e58cdc 100644
--- a/metagraph/src/annotation/annotation_converters.cpp
+++ b/metagraph/src/annotation/annotation_converters.cpp
@@ -1192,7 +1192,11 @@ void convert_to_row_diff(const std::vector<std::string> &files,
                         = utils::remove_suffix(files[i], ColumnCompressed<>::kExtension)
                                                     + ColumnCompressed<>::kCountExtension;
                     file_size += fs::file_size(values_fname);
-                } catch (...) {}
+                } catch (...) {
+                    // Count vectors may be missing for empty annotations. If a count file
+                    // is missing for a non-empty annotation, the error will be thrown later
+                    // in convert_batch_to_row_diff, so we skip it here in any case.
+                }
             }
             if (file_size > mem_bytes) {
                 logger->warn("Not enough memory to process {}, requires {} MB, skipped",

From 64a80f8c9346b07a831395e371cb398927d7d610 Mon Sep 17 00:00:00 2001
From: Mikhail Karasikov <mike-karas@ya.ru>
Date: Mon, 3 May 2021 15:41:32 +0200
Subject: [PATCH 23/23] cleanup

---
 metagraph/src/annotation/int_matrix/row_diff/int_row_diff.hpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/metagraph/src/annotation/int_matrix/row_diff/int_row_diff.hpp b/metagraph/src/annotation/int_matrix/row_diff/int_row_diff.hpp
index a58cbfe047..3f5655d76d 100644
--- a/metagraph/src/annotation/int_matrix/row_diff/int_row_diff.hpp
+++ b/metagraph/src/annotation/int_matrix/row_diff/int_row_diff.hpp
@@ -37,7 +37,7 @@ const size_t RD_PATH_RESERVE_SIZE = 2;
  */
 inline uint64_t encode_diff(int64_t x) {
     assert(x);
-    return x > 0 ? (x - 1) * 2 : (-x) * 2 - 1;
+    return (std::abs(x) - 1) * 2 + (x < 0);
 }
 
 inline int64_t decode_diff(uint64_t c) {