Hydrofabric network cycle detection

include/core/network.hpp

@@ -137,6 +137,86 @@ namespace network {
    */
   using IndexPair = std::pair< NetworkIndexT::const_iterator, NetworkIndexT::const_iterator>;
 
+  /**
+   * @brief Graph visitor which detects and stores cycles in the network
+   * 
+   */
+  struct detect_cycles : public boost::dfs_visitor<>
+  {
+    using colormap = std::map<Graph::vertex_descriptor, boost::default_color_type>;
+    colormap vertex_coloring;
+
+    using edgeColorMap = std::map<Graph::edge_descriptor, boost::default_color_type>;
+    edgeColorMap  edge_coloring;
+    using vertex_t = Graph::vertex_descriptor;
+
+    //Delete the default constructor so we don't get a nullptr for cycles
+    detect_cycles() = delete;
+    /**
+     * @brief Construct a new detect_cycle object holding a reference to a vector to store cycle pairs
+     * 
+     * @param out pointer to a vector to hold vertex pairs that complete cycles
+     */
+    detect_cycles(std::vector<std::pair<Graph::vertex_descriptor, Graph::vertex_descriptor>>* out):cycles(out){}
+
+    /**
+     * @brief Mark target verticies of visited edge
+     * 
+     * @tparam Edge 
+     * @tparam Graph 
+     * @param e 
+     * @param g 
+     */
+    /*
+    template <class Edge, class Graph>
+    void tree_edge(Edge e, const Graph& g) {
+      //Start by pushing the first vertex (source)
+      if (vertexVisited.empty()) {
+          vertexVisited.push(boost::source(e, g));
+      }
+      //Push every assoicated target vertex of the edge
+      vertexVisited.push(boost::target(e, g));
+    }*/
+
+    /**
+      * @brief When a back edge is encountered, we have a cycle
+      * 
+      * @tparam Edge 
+      * @tparam Graph 
+      * @param e 
+      * @param g 
+      */
+    template <class Edge, class Graph>
+    void back_edge(Edge e, const Graph& g) {
+      // std::cout << source(e, g)
+      //   << " -- "
+      //   << target(e, g) << "\n";
+      // std::cout<< get(boost::vertex_name, g)[ source(e, g) ] <<
+      //   " -- " << get(boost::vertex_name, g)[ target(e,g)  ]<<"\n";
+      vertex_t v2;
+      //At this point, vertexVisited is every vertice which has been visted
+      //up to the point of finding this back edge (so each is part of the total cycle)
+      /*while ( vertexVisited.top() != boost::target(e, g) )
+      {
+          //std::cout << " Cycle middle=" << vertexVisited.top() << std::endl;
+          v2 = vertexVisited.top();
+          vertexVisited.pop();
+      }*/
+      //Reset the verticies stack to the last
+      //vertexVisited.push(v2);
+      if(cycles != nullptr){
+        cycles->push_back(std::make_pair(source(e,g), target(e,g)));
+      }
+    }
+    //Currently not used, but could be used to report the reset of the 
+    //verticies involved in the cycle if needed
+    //(would need become a pointer similar to cycles if consumed/reported externally)
+    //std::stack<Graph::vertex_descriptor> vertexVisited;
+    //Since visitors are passed by value, we need an external reference
+    //to store the cycle in
+    std::vector<std::pair<Graph::vertex_descriptor, Graph::vertex_descriptor>>* cycles;
+  };
+
     /**
      * @brief A lightweight, graph based index of hydrologic features.
      * 

src/core/network.cpp

@@ -3,6 +3,8 @@
 #include <stdexcept>
 #include <boost/graph/reverse_graph.hpp>
 #include <boost/graph/graph_utility.hpp>
+#include <boost/graph/undirected_dfs.hpp>
+#include<boost/graph/properties.hpp>
 
 using namespace network;
 
@@ -96,6 +98,20 @@ void Network::init_indicies(){
     }
   }
 
+  //Do a cycle check and report cycles in the topology if they exist
+  std::vector<std::pair<Graph::vertex_descriptor, Graph::vertex_descriptor>> cycles;
+  detect_cycles vis(&cycles);
+  //Run an undirected, depth first search applying the loop detector visitor
+  //cycles added as vertex pairs for the closing edge of the cycle
+  boost::undirected_dfs(this->graph, vis, make_assoc_property_map(vis.vertex_coloring), make_assoc_property_map(vis.edge_coloring));
+  if(cycles.size() > 0){
+    std::string msg = "The following features form a cycle:\n";
+    for(auto p : cycles){
+      msg += "\t"+get_id(p.first)+" --> "+get_id(p.second)+"\n";
+    }
+    throw std::domain_error("Cycle(s) detected in hydrofabric\n"+msg);
+  }
+
   boost::topological_sort(this->graph, std::back_inserter(this->topo_order),
                    boost::vertex_index_map(get(boost::vertex_index, this->graph)));
 }

test/core/NetworkTests.cpp

@@ -145,6 +145,52 @@ class Network_Test2 : public Network_Test, public ::testing::Test{
   }
 };
 
+class Network_Test_Cycle : public Network_Test, public ::testing::Test{
+public:
+  Network_Test_Cycle(){}
+  void SetUp(){
+    /**
+     * The following create a topology that looks like
+     *    ------------------------
+     *  /                          \
+     * |                            \
+     * cat-0            cat-3 -      |
+     *      \                   \    | 
+     *        nex 0 -> cat-2 ->  nex-1
+     *      /                   /
+     * cat-1            cat-4 -
+     * 
+     * This contains a cycle from nex-1 back to cat-0
+     * 
+     */
+    this->add_catchment("cat-0", "nex-0");
+    this->add_catchment("cat-1", "nex-0");
+    this->add_nexus("nex-0", "cat-2");
+    this->add_catchment("cat-2", "nex-1");
+    this->add_catchment("cat-3", "nex-1");
+    this->add_catchment("cat-4", "nex-1");
+    this->add_nexus("nex-1", "cat-0");
+    //Constructing this network SHOULD fail due to the cycle
+  }
+};
+
+//! Test that a network with cycles throws an error
+TEST_F(Network_Test_Cycle, TestCycleError){
+
+  EXPECT_THROW( {
+    try{
+      n =  Network(this->get_fabric());
+    }catch(std::domain_error &e){
+      std::cout<<e.what();
+      EXPECT_STREQ(e.what(), "Cycle(s) detected in hydrofabric\n"
+                             "The following features form a cycle:\n"
+                             "\tnex-1 --> cat-0\n");
+      throw;
+    }
+  }, std::domain_error);
+
+}
+
 //! Test that a network can be created.
 TEST_P(Network_Test1, TestNetworkConstructionNumberOfNodes)
 {