use paryfor, update ips4o to avoid a sorting bug

Author: ekg

.gitmodules

@@ -34,3 +34,6 @@
 [submodule "deps/odgi"]
 	path = deps/odgi
 	url = https://github.com/vgteam/odgi.git
+[submodule "deps/paryfor"]
+	path = deps/paryfor
+	url = https://github.com/ekg/paryfor.git

CMakeLists.txt

@@ -165,6 +165,16 @@ ExternalProject_Add(cgranges
 ExternalProject_Get_property(cgranges SOURCE_DIR)
 set(cgranges_INCLUDE "${SOURCE_DIR}/cpp")
 
+# paryfor parallel_for
+ExternalProject_Add(paryfor
+  SOURCE_DIR "${CMAKE_SOURCE_DIR}/deps/paryfor"
+  UPDATE_COMMAND ""
+  INSTALL_COMMAND ""
+  BUILD_COMMAND ""
+  CONFIGURE_COMMAND "")
+ExternalProject_Get_property(paryfor SOURCE_DIR)
+set(paryfor_INCLUDE "${SOURCE_DIR}")
+
 add_subdirectory(deps/spoa EXCLUDE_FROM_ALL)
 set(spoa_INCLUDE "${CMAKE_SOURCE_DIR}/deps/spoa/include")
 
@@ -187,6 +197,7 @@ add_dependencies(smoothxg_objs tayweeargs)
 add_dependencies(smoothxg_objs gfakluge)
 add_dependencies(smoothxg_objs ips4o)
 add_dependencies(smoothxg_objs cgranges)
+add_dependencies(smoothxg_objs paryfor)
 
 set_target_properties(smoothxg_objs PROPERTIES POSITION_INDEPENDENT_CODE TRUE)
 
@@ -228,7 +239,8 @@ set(smoothxg_INCLUDES
   "${dynamic_INCLUDE}"
   "${hopscotch_map_INCLUDE}"
   "${sparsepp_INCLUDE}"
-  "${flathashmap_INCLUDE}")
+  "${flathashmap_INCLUDE}"
+  "${paryfor_INCLUDE}")
 
 set(smoothxg_LIBS
   "${sdsl-lite_LIB}/libsdsl.a"

deps/ips4o

@@ -115,63 +115,68 @@ odgi::graph_t smooth_and_lace(const xg::XG& graph,
     std::vector<path_position_range_t> path_mapping;
     std::vector<path_position_range_t> consensus_mapping;
     bool add_consensus = !consensus_base_name.empty();
-#pragma omp parallel for schedule(dynamic, 1)
-    for (uint64_t block_id = 0; block_id < blocks.size(); ++block_id) {
-        //for (auto& block : blocks) {
-        auto& block = blocks[block_id];
-        std::string consensus_name = consensus_base_name + std::to_string(block_id);
-        //std::cerr << "on block " << block_id+1 << " of " << blocks.size() << std::endl;
-        block_graphs[block_id] = smooth(graph, block, consensus_name);
-        auto& block_graph = block_graphs[block_id];
-        if (block_graph.get_node_count() > 0) {
-            //auto& block_graph = block_graphs.back();
-            // record the start and end paths
-            // nb: the path order is the same in the input block and output graph
-            uint64_t path_id = 0;
-            for (auto& path_range : block.path_ranges) {
-                auto path_handle = graph.get_path_handle_of_step(path_range.begin);
-                auto last_step = graph.get_previous_step(path_range.end);
-#pragma omp critical (path_mapping)
-                path_mapping.push_back({
-                        path_handle, // target path
-                        graph.get_position_of_step(path_range.begin), // start position
-                        (graph.get_position_of_step(last_step) // end position
-                         + graph.get_length(graph.get_handle_of_step(last_step))),
-                        path_range.begin,
-                        path_range.end,
-                        as_path_handle(++path_id),
-                        block_id
-                    });
-            }
-            // make the graph
+    std::mutex path_mapping_mutex, consensus_mapping_mutex;
+    paryfor::parallel_for<uint64_t>(
+        0, blocks.size(),
+        odgi::get_thread_count(),
+        [&](uint64_t block_id, int tid) {
+            auto& block = blocks[block_id];
+            std::string consensus_name = consensus_base_name + std::to_string(block_id);
+            //std::cerr << "on block " << block_id+1 << " of " << blocks.size() << std::endl;
+            auto& block_graph = block_graphs[block_id];
+            block_graph = smooth(graph, block, consensus_name);
+            if (block_graph.get_node_count() > 0) {
+                //auto& block_graph = block_graphs.back();
+                // record the start and end paths
+                // nb: the path order is the same in the input block and output graph
+                uint64_t path_id = 0;
+                for (auto& path_range : block.path_ranges) {
+                    auto path_handle = graph.get_path_handle_of_step(path_range.begin);
+                    auto last_step = graph.get_previous_step(path_range.end);
+                    {
+                        std::lock_guard<std::mutex> guard(path_mapping_mutex);
+                        path_mapping.push_back({
+                                path_handle, // target path
+                                graph.get_position_of_step(path_range.begin), // start position
+                                (graph.get_position_of_step(last_step) // end position
+                                 + graph.get_length(graph.get_handle_of_step(last_step))),
+                                path_range.begin,
+                                path_range.end,
+                                as_path_handle(++path_id),
+                                block_id
+                            });
+                    }
+                }
+                // make the graph
 
-            // record the consensus path
-            if (add_consensus) {
-                auto consensus_handle = block_graph.get_path_handle(consensus_name);
-                uint64_t path_end = 0;
-                step_handle_t empty_step;
-                as_integers(empty_step)[0] = 0;
-                as_integers(empty_step)[1] = 0;
-                block_graph.for_each_step_in_path(
-                    consensus_handle,
-                    [&](const step_handle_t& step) {
-                        path_end += block_graph.get_length(block_graph.get_handle_of_step(step));
-                    });
-#pragma omp critical (consensus_mapping)
-                consensus_mapping.push_back({
-                        as_path_handle(0), // consensus = 0 path handle
-                        0, // start position
-                        path_end, // end position
-                        empty_step,
-                        empty_step,
+                // record the consensus path
+                if (add_consensus) {
+                    auto consensus_handle = block_graph.get_path_handle(consensus_name);
+                    uint64_t path_end = 0;
+                    step_handle_t empty_step;
+                    as_integers(empty_step)[0] = 0;
+                    as_integers(empty_step)[1] = 0;
+                    block_graph.for_each_step_in_path(
                         consensus_handle,
-                        block_id
-                    });
+                        [&](const step_handle_t& step) {
+                            path_end += block_graph.get_length(block_graph.get_handle_of_step(step));
+                        });
+                    {
+                        std::lock_guard<std::mutex> guard(consensus_mapping_mutex);
+                        consensus_mapping.push_back({
+                                as_path_handle(0), // consensus = 0 path handle
+                                0, // start position
+                                path_end, // end position
+                                empty_step,
+                                empty_step,
+                                consensus_handle,
+                                block_id
+                            });
+                    }
+                }
             }
-            // increment our block id
-            //++block_id;
-        }
-    }
+        });
+
     // sort the path range mappings by path handle id, then start position
     // this will allow us to walk through them in order
     ips4o::parallel::sort(

deps/paryfor

@@ -6,6 +6,7 @@
 #include <algorithm>
 #include <cmath>
 #include <sstream>
+#include <mutex>
 #include "spoa/spoa.hpp"
 #include "xg.hpp"
 #include "blocks.hpp"
@@ -14,6 +15,7 @@
 #include "odgi/unchop.hpp"
 #include "odgi/topological_sort.hpp"
 #include "odgi/dna.hpp"
+#include "paryfor.hpp"
 
 namespace smoothxg {