Attempt to correctly partition in the presence of a terminal nexus with no defined outflow destination (NGWPC-6553)

include/core/network.hpp

@@ -132,7 +132,7 @@ namespace network {
   /**
    * @brief A type for holding a pair of network::NetworkIndexT iterators,  pair.first=begin,  pair.second=end
    * 
-   * these indicies should be const to the caller
+   * these indices should be const to the caller
    * 
    */
   using IndexPair = std::pair< NetworkIndexT::const_iterator, NetworkIndexT::const_iterator>;
@@ -164,7 +164,7 @@ namespace network {
          * @param features a geojson::GeoJSON collection of features to add as nodes to the graph
          * @param link_key the property to read from features to determine edge linking, i.e. 'toid'
          */
-        Network( geojson::GeoJSON features, std::string* link_key);
+        Network( geojson::GeoJSON features, std::string const* link_key);
 
         /**
          * @brief Destroy the Network object
@@ -330,7 +330,7 @@ namespace network {
          * @brief Initializes the head/tailwater iterators after the underlying graph is constructed.
          * 
          */
-        void init_indicies();
+        void init_indices();
 
         /**
          * @brief Vector of topologically sorted features

include/geojson/FeatureCollection.hpp

@@ -166,7 +166,7 @@ namespace geojson {
 
             void set(const std::string& key, JSONProperty& property);
 
-            void add_feature(Feature feature, std::string *id = nullptr);
+            void add_feature(Feature feature, std::string const*id = nullptr);
 
             /**
              * Add a reference to a feature by id
@@ -189,9 +189,9 @@ namespace geojson {
             */
             void update_ids(const std::string& alt_id = {});
 
-            int link_features_from_property(std::string* from_property = nullptr, std::string* to_property = nullptr);
+            int link_features_from_property(std::string const* from_property = nullptr, std::string const* to_property = nullptr);
 
-            int link_features_from_attribute(std::string* from_attribute = nullptr, std::string* to_attribute = nullptr);
+            int link_features_from_attribute(std::string const* from_attribute = nullptr, std::string const* to_attribute = nullptr);
 
             /* Untested, excluded for now in favor of filter constructor above.
             template<typename C>
@@ -224,4 +224,4 @@ namespace geojson {
     };
 }
 
-#endif // GEOJSON_FEATURE_COLLECTION_H
+#endif // GEOJSON_FEATURE_COLLECTION_H

src/core/HY_Features.cpp

@@ -71,4 +71,4 @@ HY_Features::HY_Features(network::Network network, std::shared_ptr<Formulation_M
         }
       }
 
-}
+}

src/core/network.cpp

@@ -30,6 +30,8 @@ Network::Network( geojson::GeoJSON fabric ){
   std::string feature_id, downstream_id;
   Graph::vertex_descriptor v1, v2;
 
+  const std::string terminal = "wb-TERMINAL_SENTINEL-";
+
   for(auto& feature: *fabric)
   {
     feature_id = feature->get_id();
@@ -59,7 +61,17 @@ Network::Network( geojson::GeoJSON fabric ){
     }
 
     //Add the downstream features/edges
-    for( auto& downstream: feature->destination_features() )
+    auto const& destination_features = feature->destination_features();
+    if (destination_features.empty()) {
+        //std::cerr << "Adding sentinel for " << feature_id << std::endl;
+        auto name = terminal + feature_id;
+        auto v_sentinel = add_vertex(name, this->graph);
+        this->descriptor_map.emplace( name, v_sentinel );
+        this->layer_map.emplace(name, DEFAULT_LAYER_ID);
+        add_edge(v1, v_sentinel, this->graph);
+    }
+
+    for( auto& downstream: destination_features )
     {
       downstream_id = downstream->get_id();
       if( this->descriptor_map.find(downstream_id) != this->descriptor_map.end() )
@@ -77,10 +89,10 @@ Network::Network( geojson::GeoJSON fabric ){
       //std::cout<<"Added edge: "<<feature_id<<" -> "<<downstream_id<<std::endl;
     }
   }
-  init_indicies();
+  init_indices();
 }
 
-void Network::init_indicies(){
+void Network::init_indices(){
 
   Graph::vertex_iterator begin, end;
   boost::tie(begin, end) = boost::vertices(this->graph);
@@ -100,43 +112,45 @@ void Network::init_indicies(){
                    boost::vertex_index_map(get(boost::vertex_index, this->graph)));
 }
 
-Network::Network( geojson::GeoJSON features, std::string* link_key ){
-
+Network::Network( geojson::GeoJSON features, std::string const* link_key ){
   std::string feature_id, downstream_id;
   Graph::vertex_descriptor v1, v2;
+  const std::string terminal = "wb-TERMINAL_SENTINEL-";
 
   //TODO ensure all features are the same logical HY_Features type?
-  for(auto& feature: *features)
-  {
+  for (auto& feature: *features) {
     feature_id = feature->get_id();
-    if( this->descriptor_map.find( feature_id ) == this->descriptor_map.end() )
-    {
+    if (this->descriptor_map.find( feature_id ) == this->descriptor_map.end()) {
       //Haven't visited this feature yet, add it to graph
       //add vertex to graph
       v1 = add_vertex( feature_id, this->graph );
       this->descriptor_map.emplace( feature_id, v1 );
-    }
-    else{
+    } else {
       v1 = this->descriptor_map[ feature_id ];
     }
 
-      if (link_key != nullptr and feature->has_property(*link_key)) {
-
-          downstream_id = feature->get_property(*link_key).as_string();
-          if( this->descriptor_map.find(downstream_id) != this->descriptor_map.end() )
-          {
-            v2 = this->descriptor_map[ downstream_id ];
-          }
-          else {
-            v2 = add_vertex( downstream_id, this->graph );
-            this->descriptor_map.emplace( downstream_id, v2 );
-          }
-            add_edge(v1, v2, this->graph);
+    if (link_key == nullptr)
+      continue;
+
+    if (feature->has_property(*link_key)) {
+      downstream_id = feature->get_property(*link_key).as_string();
+      std::cout << feature_id << " -> " << downstream_id << std::endl;
+      if (this->descriptor_map.find(downstream_id) != this->descriptor_map.end()) {
+        v2 = this->descriptor_map[ downstream_id ];
+      } else {
+        v2 = add_vertex( downstream_id, this->graph );
+        this->descriptor_map.emplace( downstream_id, v2 );
       }
+    } else {
+      std::cout << feature_id << " getting sentinel added\n";
+      downstream_id = terminal + feature_id;
+      v2 = add_vertex( downstream_id, this->graph );
+      this->descriptor_map.emplace( downstream_id, v2 );
+    }
+    add_edge(v1, v2, this->graph);
   }
 
-  init_indicies();
-
+  init_indices();
 }
 
 NetworkIndexT::const_reverse_iterator Network::begin(){

src/geojson/FeatureCollection.cpp

@@ -129,7 +129,7 @@ void FeatureCollection::set(const std::string& key, JSONProperty& property) {
     foreign_members.emplace(key, property);
 }
 
-int FeatureCollection::link_features_from_property(std::string* from_property, std::string* to_property) {
+int FeatureCollection::link_features_from_property(std::string const* from_property, std::string const* to_property) {
     int links_found = 0;
 
     for (Feature feature : features) {
@@ -155,7 +155,7 @@ int FeatureCollection::link_features_from_property(std::string* from_property, s
     return links_found;
 }
 
-int FeatureCollection::link_features_from_attribute(std::string* from_attribute, std::string* to_attribute) {
+int FeatureCollection::link_features_from_attribute(std::string const* from_attribute, std::string const* to_attribute) {
     int links_found = 0;
 
     for (Feature feature : features) {
@@ -198,7 +198,7 @@ int FeatureCollection::link_features_from_attribute(std::string* from_attribute,
     return links_found;
 }
 
-void FeatureCollection::add_feature(Feature feature, std::string *id) {
+void FeatureCollection::add_feature(Feature feature, std::string const* id) {
     features.push_back(feature);
 
     if (id != nullptr) {

src/partitionGenerator.cpp

@@ -110,12 +110,13 @@ void write_remote_connections(const PartitionVSet& catchment_part, const Partiti
  * @param num_catchments
  * @param catchment_part 
  */
-void generate_partitions(network::Network& network, const int& num_partitions, const int& num_catchments, PartitionVSet& catchment_part,
-     PartitionVSet& nexus_part)
+void generate_partitions(network::Network& network, const int& num_partitions, PartitionVSet& catchment_part, PartitionVSet& nexus_part)
 {
+    auto catchments = network.filter("cat", network::SortOrder::TransposedDepthFirstPreorder);
+
     int partition = 0;
     int counter = 0;
-    int total = num_catchments;
+    int total = size(catchments);
     int partition_size = total/num_partitions;
     int partition_size_norm = partition_size;
     int remainder;
@@ -134,9 +135,7 @@ void generate_partitions(network::Network& network, const int& num_partitions, c
     nexus_set.reserve(partition_size);
     std::string part_id, partition_str;
 
-    std::string up_nexus;
-    std::string down_nexus;
-    for(const auto& catchment : network.filter("cat", network::SortOrder::TransposedDepthFirstPreorder)){
+    for(const auto& catchment : catchments){
             if (partition < remainder)
                 partition_size = partition_size_plus1;
             else
@@ -150,9 +149,7 @@ void generate_partitions(network::Network& network, const int& num_partitions, c
             }
             if(nexus_set.size() == 0){
                 partgen_ss <<"Error: Catchment "<<catchment<<" has no destination nexus.\n";
-                LOG(partgen_ss.str(), LogLevel::FATAL); partgen_ss.str("");
-
-                exit(1);
+                LOG(partgen_ss.str(), LogLevel::WARNING); partgen_ss.str("");
             }
             for( auto upstream : network.get_origination_ids(catchment) ){
                 nexus_set.emplace(upstream);
@@ -161,18 +158,14 @@ void generate_partitions(network::Network& network, const int& num_partitions, c
 
             //keep track of all the features in this partition
             catchment_set.emplace(catchment);
+            LOG(catchment + " placed in partition " + std::to_string(partition), LogLevel::DEBUG);
+
             counter++;
             if(counter == partition_size)
             {
                 //partgen_ss<<"nexus "<<nexus<<" is remote DOWN on partition "<<partition<<std::endl;
                 //FIXME partitioning shouldn't have to assume dendritic network
                 std::vector<std::string> destinations = network.get_destination_ids(catchment);
-                if(destinations.size() == 0){
-                    partgen_ss <<"Error: Catchment "<<catchment<<" has no destination nexus.\n";
-                    LOG(partgen_ss.str(), LogLevel::FATAL); partgen_ss.str("");
-                    exit(1);
-                }
-                down_nexus = destinations[0];
 
                 part_id = std::to_string(partition);  // Is id used?
                 partition_str = std::to_string(partition);
@@ -183,15 +176,10 @@ void generate_partitions(network::Network& network, const int& num_partitions, c
                 catchment_set.clear();
                 nexus_set.clear();
 
-                partition_str = std::to_string(partition);
-
                 partition++;
                 counter = 0;
                 //partgen_ss<<"\nnexus "<<nexus<<" is remote UP on partition "<<partition<<std::endl;
 
-                //this nexus overlaps partitions
-                //Handled above by ensure all up/down stream nexuses are recorded
-                up_nexus = down_nexus;
                 //partgen_ss<<"\nin partition "<<partition<<":"<<std::endl;
             }
     }
@@ -474,12 +462,9 @@ int main(int argc, char* argv[])
 
     std::string link_key = "toid";
 
-    Network catchment_network(catchment_collection, &link_key);
     //Assumes dendritic, can add check in network if needed.
     PartitionVSet catchment_part, nexus_part;
 
-    //catchment_network.print_network();
-
     //build the remote connections from network
     // read the nexus hydrofabric, reuse the catchments
     geojson::GeoJSON global_nexus_collection;
@@ -521,13 +506,23 @@ int main(int argc, char* argv[])
     Network global_network(global_nexus_collection);
 
     //Generate the partitioning
-    generate_partitions(global_network, num_partitions, num_catchments, catchment_part, nexus_part);
+    generate_partitions(global_network, num_partitions, catchment_part, nexus_part);
 
     //global_network.print_network();
 
     //The container holding all remote_connections
     std::vector<RemoteConnectionVec> remote_connections_vec;
 
+    for (int i = 0; i < num_partitions; ++i) {
+        partgen_ss << "Partition " << i << " catchments: " << catchment_part[i].size() << "\n";
+        for (auto& c : catchment_part[i])
+            partgen_ss << c << std::endl;
+        partgen_ss << "nexuses " << nexus_part[i].size() << "\n";
+        for (auto& n : nexus_part[i])
+            partgen_ss << n << std::endl;
+        LOG(partgen_ss.str(), LogLevel::DEBUG); partgen_ss.str("");
+    }
+
     int total_remotes = 0;
     // loop over all partitions by partition id
     for (int ipart=0; ipart < catchment_part.size(); ++ipart)