hogwild path cover of the graph

src/algorithms/cover.cpp

@@ -383,5 +383,115 @@ namespace odgi {
             }
         }
 
+
+
+void hogwild_path_cover(handlegraph::MutablePathDeletableHandleGraph &graph,
+                        double target_depth,
+                        //bool write_node_coverages, std::string &node_coverages,
+                        const uint64_t& nthreads, const bool& ignore_paths, const bool& show_progress) {
+
+    // get depth -> graph depth
+    uint64_t node_count = graph.get_node_count();
+    uint64_t graph_bp = 0;
+    uint64_t step_count = 0;
+    // make a rank select dictionary over our sequence space
+    // to get a node randomly distributed in the total length of the graph
+    graph.for_each_handle(
+        [&](const handle_t& h) {
+            graph_bp += graph.get_length(h);
+            step_count += graph.get_step_count(h);
+        });
+    //sdsl::bit_vector graph_vec();
+    sdsl::bit_vector graph_bv(graph_bp);
+    graph_bp = 0;
+    graph.for_each_handle(
+        [&](const handle_t& h) {
+            graph_bv[graph_bp] = 1;
+            graph_bp += graph.get_length(h);
+        });
+    sdsl::bit_vector::rank_1_type graph_bv_rank;
+    sdsl::util::assign(graph_bv_rank, sdsl::bit_vector::rank_1_type(&graph_bv));
+    uint64_t target_step_count = node_count * target_depth - (ignore_paths ? 0 : step_count);
+    // run hogwild until our thread count is
+    std::unique_ptr<progress_meter::ProgressMeter> progress_meter;
+    if (show_progress) {
+        progress_meter = std::make_unique<progress_meter::ProgressMeter>(
+            target_step_count, "[odgi::hogwild_cover] covering the graph:");
+    }
+    // launch a thread to update the learning rate, count iterations, and decide when to stop
+    std::atomic<uint64_t> added_steps; added_steps.store(0);
+    std::atomic<uint64_t> added_paths; added_paths.store(0);
+    auto worker_lambda =
+        [&](uint64_t tid) {
+            // everyone tries to seed with their own random data
+            const std::uint64_t seed = 9399220 + tid;
+            XoshiroCpp::Xoshiro256Plus gen(seed); // a nice, fast PRNG
+
+            // we'll sample from all path steps
+            std::uniform_int_distribution<uint64_t> dis_graph_pos = std::uniform_int_distribution<uint64_t>(0, graph_bp-1);
+            std::uniform_int_distribution<uint64_t> flip(0, 1);
+            while (added_steps.load() < target_step_count) {
+                // create a path
+                std::stringstream ss;
+                ss << "cover_" << added_paths++;
+                path_handle_t path = graph.create_path_handle(ss.str());
+                // find a random handle
+                handle_t h = graph.get_handle(graph_bv_rank(dis_graph_pos(gen)), flip(gen));
+                uint64_t iter = 0;
+                handle_t lowest;
+                bool seen_low = false;
+                while (graph.get_step_count(h) > 0 && iter++ < 100) {
+                    h = graph.get_handle(graph_bv_rank(dis_graph_pos(gen)), flip(gen));
+                    if (!seen_low || graph.get_step_count(h) < graph.get_step_count(lowest)) {
+                        lowest = h;
+                        seen_low = true;
+                    }
+                }
+                if (seen_low) {
+                    h = lowest;
+                }
+                while (true) {
+                    graph.append_step(path, h);
+                    if (show_progress) progress_meter->increment(1);
+                    ++added_steps;
+                    if (added_steps.load() >= target_step_count) {
+                        break;
+                    }
+                    handle_t best_next;
+                    uint64_t lowest_cov = std::numeric_limits<uint64_t>::max();
+                    graph.follow_edges(
+                        h, false,
+                        [&](const handle_t& n) {
+                            uint64_t next_cov = graph.get_step_count(n);
+                            if (next_cov < lowest_cov) {
+                                best_next = n;
+                                lowest_cov = next_cov;
+                            }
+                        });
+                    if (lowest_cov < std::numeric_limits<uint64_t>::max()) {
+                        h = best_next;
+                    } else {
+                        break;
+                    }
+                }
+            }
+        };
+
+    std::vector<std::thread> workers;
+    workers.reserve(nthreads);
+    for (uint64_t t = 0; t < nthreads; ++t) {
+        workers.emplace_back(worker_lambda, t);
+    }
+
+    for (uint64_t t = 0; t < nthreads; ++t) {
+        workers[t].join();
+    }
+
+    if (show_progress) {
+        progress_meter->finish();
+    }
+}
+
+
     }
 }

src/algorithms/cover.hpp

@@ -3,14 +3,16 @@
 #include <unordered_map>
 #include <stack>
 #include <map>
-
+#include <sdsl/bit_vectors.hpp>
 #include <handlegraph/mutable_path_deletable_handle_graph.hpp>
 #include <handlegraph/util.hpp>
 #include "weakly_connected_components.hpp"
 #include <deps/ips4o/ips4o.hpp>
+#include "XoshiroCpp.hpp"
+#include "progress.hpp"
 
 namespace odgi {
-    namespace algorithms {
+namespace algorithms {
 
         using namespace handlegraph;
 
@@ -57,5 +59,11 @@ namespace odgi {
                         size_t min_node_coverage, size_t max_number_of_paths_generable,
                         bool write_node_coverages, std::string &node_coverages,
                         const uint64_t& nthreads, const bool& ignore_paths, const bool& show_progress);
-    }
-}
+
+void hogwild_path_cover(handlegraph::MutablePathDeletableHandleGraph &graph,
+                        double target_depth,
+                        //bool write_node_coverages, std::string &node_coverages,
+                        const uint64_t& nthreads, const bool& ignore_paths, const bool& show_progress);
+
+}
+}

src/subcommand/cover_main.cpp

@@ -21,6 +21,7 @@ namespace odgi {
         args::HelpFlag help(parser, "help", "display this help summary", {'h', "help"});
         args::ValueFlag<std::string> dg_in_file(parser, "FILE", "load the graph from this file", {'i', "idx"});
         args::ValueFlag<std::string> dg_out_file(parser, "FILE","store the graph with the generated paths in this file", {'o', "out"});
+        args::ValueFlag<double> hogwild_depth(parser, "DEPTH", "randomly cover the graph until it has reaches the given average DEPTH",{'H', "hogwild-depth"});
         args::ValueFlag<uint64_t> num_paths_per_component(parser, "N", "number of paths to generate per component",{'n', "num-paths-per-component"});
         args::ValueFlag<uint64_t> node_window_size(parser, "N","size of the node window to check each time a new path is extended (it has to be greater than or equal to 2)",{'k', "node-window-size"});
         args::ValueFlag<uint64_t> min_node_coverage(parser, "N","minimum node coverage to reach (it has to be greater than 0)",{'c', "min-node-coverage"});
@@ -91,32 +92,36 @@ namespace odgi {
             }
         }
 
-        uint64_t max_number_of_paths_generable = graph.get_node_count() * 5;
-        if (args::get(debug)){
-            if (_min_node_coverage) {
-                std::cerr << "[odgi::cover] there will be generated paths until the minimum node coverage is " << _min_node_coverage
-                          << ", or until the maximum number of allowed generated paths is reached ("
-                          << max_number_of_paths_generable << ")." << std::endl;
-            } else {
-                std::cerr << "[odgi::cover] there will be generated " << _num_paths_per_component << " paths per component."
-                          << std::endl;
-            }
-        }
-
         uint64_t num_threads = args::get(nthreads) ? args::get(nthreads) : 1;
 
-        std::string node_coverages;
-        algorithms::path_cover(graph, _num_paths_per_component, _node_window_size, _min_node_coverage,
-                               max_number_of_paths_generable,
-                               write_node_coverages, node_coverages,
-                               num_threads, args::get(ignore_paths), args::get(debug));
+        if (hogwild_depth) {
+            algorithms::hogwild_path_cover(graph, args::get(hogwild_depth), num_threads, args::get(ignore_paths), args::get(debug));
+        } else {
+            uint64_t max_number_of_paths_generable = graph.get_node_count() * 5;
+            if (args::get(debug)){
+                if (_min_node_coverage) {
+                    std::cerr << "[odgi::cover] there will be generated paths until the minimum node coverage is " << _min_node_coverage
+                              << ", or until the maximum number of allowed generated paths is reached ("
+                              << max_number_of_paths_generable << ")." << std::endl;
+                } else {
+                    std::cerr << "[odgi::cover] there will be generated " << _num_paths_per_component << " paths per component."
+                              << std::endl;
+                }
+            }
+
+            std::string node_coverages;
+            algorithms::path_cover(graph, _num_paths_per_component, _node_window_size, _min_node_coverage,
+                                   max_number_of_paths_generable,
+                                   write_node_coverages, node_coverages,
+                                   num_threads, args::get(ignore_paths), args::get(debug));
 
-        if (write_node_coverages) {
-            std::string covfile = args::get(write_node_coverages);
+            if (write_node_coverages) {
+                std::string covfile = args::get(write_node_coverages);
 
-            ofstream f(covfile.c_str());
-            f << node_coverages;
-            f.close();
+                ofstream f(covfile.c_str());
+                f << node_coverages;
+                f.close();
+            }
         }
 
         std::string outfile = args::get(dg_out_file);