Merge pull request #954 from EnergySystemsModellingLab/topo_sort_cycles

Solve investment order for commodity graphs with cycles

Author: tsmbland

src/graph.rs

@@ -2,16 +2,16 @@
 use crate::commodity::{CommodityID, CommodityMap, CommodityType};
 use crate::process::{ProcessID, ProcessMap};
 use crate::region::RegionID;
+use crate::simulation::investment::InvestmentSet;
 use crate::time_slice::{TimeSliceInfo, TimeSliceLevel, TimeSliceSelection};
 use crate::units::{Dimensionless, Flow};
-use anyhow::{Context, Result, anyhow, ensure};
+use anyhow::{Context, Result, ensure};
 use indexmap::IndexSet;
-use itertools::{Itertools, iproduct};
+use itertools::iproduct;
 use petgraph::Directed;
-use petgraph::algo::toposort;
+use petgraph::algo::{condensation, toposort};
 use petgraph::dot::Dot;
 use petgraph::graph::{EdgeReference, Graph};
-use petgraph::visit::EdgeFiltered;
 use std::collections::HashMap;
 use std::fmt::Display;
 use std::fs::File;
@@ -301,43 +301,64 @@ fn validate_commodities_graph(
 /// Performs topological sort on the commodity graph to get the ordering for investments
 ///
 /// The returned Vec only includes SVD and SED commodities.
-fn topo_sort_commodities(
+fn solve_investment_order(
     graph: &CommoditiesGraph,
     commodities: &CommodityMap,
-) -> Result<Vec<CommodityID>> {
-    // We only consider primary edges
-    let primary_graph =
-        EdgeFiltered::from_fn(graph, |edge| matches!(edge.weight(), GraphEdge::Primary(_)));
-
-    // Perform a topological sort on the graph
-    let order = toposort(&primary_graph, None).map_err(|cycle| {
-        let cycle_commodity = graph.node_weight(cycle.node_id()).unwrap().clone();
-        anyhow!("Cycle detected in commodity graph for commodity {cycle_commodity}")
-    })?;
-
-    // We return the order in reverse so that leaf-node commodities are solved first
-    // We also filter to only include SVD and SED commodities
-    let order = order
-        .iter()
-        .rev()
-        .filter_map(|node_idx| {
+) -> Vec<InvestmentSet> {
+    // Filter the graph to only include SVD/SED commodities and primary edges
+    let graph_filtered = graph.filter_map(
+        // Consider only SVD/SED commodities
+        |_, node_weight| {
             // Get the commodity for the node
-            let GraphNode::Commodity(commodity_id) = graph.node_weight(*node_idx).unwrap() else {
+            let GraphNode::Commodity(commodity_id) = node_weight else {
                 // Skip special nodes
                 return None;
             };
             let commodity = &commodities[commodity_id];
-
-            // Only include SVD and SED commodities
             matches!(
                 commodity.kind,
                 CommodityType::ServiceDemand | CommodityType::SupplyEqualsDemand
             )
-            .then(|| commodity_id.clone())
+            .then_some(node_weight)
+        },
+        // Consider only primary edges
+        |_, edge_weight| matches!(edge_weight, GraphEdge::Primary(_)).then_some(edge_weight),
+    );
+
+    // Condense strongly connected components
+    let condensed_graph = condensation(graph_filtered, true);
+
+    // Perform a topological sort on the condensed graph
+    // We can safely unwrap because `toposort` will only return an error in case of cycles, which
+    // should have been detected and compressed with `condensation`
+    let order = toposort(&condensed_graph, None).unwrap();
+
+    // Create investment sets (reverse topological order)
+    order
+        .iter()
+        .rev()
+        .filter_map(|node_idx| {
+            // Get set of commodity ID(s) for the node, referring back to `condensed_graph`
+            let commodities: Vec<CommodityID> = condensed_graph
+                .node_weight(*node_idx)
+                .unwrap()
+                .iter()
+                .filter_map(|node| match node {
+                    GraphNode::Commodity(id) => Some(id.clone()),
+                    _ => None,
+                })
+                .collect();
+
+            // Create investment set
+            // If a single commodity in the node this is `InvestmentSet::Single`, if multiple
+            // commodities this is `InvestmentSet::Cycle`
+            match commodities.as_slice() {
+                [] => None,
+                [only] => Some(InvestmentSet::Single(only.clone())),
+                _ => Some(InvestmentSet::Cycle(commodities)),
+            }
         })
-        .collect();
-
-    Ok(order)
+        .collect()
 }
 
 /// Builds base commodity graphs for each region and year
@@ -347,65 +368,42 @@ pub fn build_commodity_graphs_for_model(
     processes: &ProcessMap,
     region_ids: &IndexSet<RegionID>,
     years: &[u32],
-) -> Result<HashMap<(RegionID, u32), CommoditiesGraph>> {
-    let commodity_graphs: HashMap<(RegionID, u32), CommoditiesGraph> =
-        iproduct!(region_ids, years.iter())
-            .map(|(region_id, year)| {
-                let graph = create_commodities_graph_for_region_year(processes, region_id, *year);
-                ((region_id.clone(), *year), graph)
-            })
-            .collect();
-
-    Ok(commodity_graphs)
+) -> HashMap<(RegionID, u32), CommoditiesGraph> {
+    iproduct!(region_ids, years.iter())
+        .map(|(region_id, year)| {
+            let graph = create_commodities_graph_for_region_year(processes, region_id, *year);
+            ((region_id.clone(), *year), graph)
+        })
+        .collect()
 }
 
 /// Validates commodity graphs for the entire model.
 ///
-/// This function creates commodity flow graphs for each region/year combination in the model,
-/// validates the graph structure against commodity type rules, and determines the optimal
-/// investment order for commodities.
-///
 /// The validation process checks three time slice levels:
 /// - **Annual**: Validates annual-level commodities and processes
 /// - **Seasonal**: Validates seasonal-level commodities and processes for each season
 /// - **Day/Night**: Validates day/night-level commodities and processes for each time slice
 ///
 /// # Arguments
 ///
+/// * `commodity_graphs` - Commodity graphs for each region and year, outputted from `build_commodity_graphs_for_model`
 /// * `processes` - All processes in the model with their flows and activity limits
 /// * `commodities` - All commodities with their types and demand specifications
 /// * `region_ids` - Collection of regions to model
 /// * `years` - Years to analyse
 /// * `time_slice_info` - Time slice configuration (seasons, day/night periods)
 ///
-/// # Returns
-///
-/// A map from `(region, year)` to the ordered list of commodities for investment decisions. The
-/// ordering ensures that leaf-node commodities (those with no outgoing edges) are solved first.
-///
 /// # Errors
 ///
 /// Returns an error if:
-/// - Any commodity graph contains cycles
 /// - Commodity type rules are violated (e.g., SVD commodities being consumed)
 /// - Demand cannot be satisfied
 pub fn validate_commodity_graphs_for_model(
     commodity_graphs: &HashMap<(RegionID, u32), CommoditiesGraph>,
     processes: &ProcessMap,
     commodities: &CommodityMap,
     time_slice_info: &TimeSliceInfo,
-) -> Result<HashMap<(RegionID, u32), Vec<CommodityID>>> {
-    // Determine commodity ordering for each region and year
-    let commodity_order: HashMap<(RegionID, u32), Vec<CommodityID>> = commodity_graphs
-        .iter()
-        .map(|((region_id, year), graph)| -> Result<_> {
-            let order = topo_sort_commodities(graph, commodities).with_context(|| {
-                format!("Error validating commodity graph for {region_id} in {year}")
-            })?;
-            Ok(((region_id.clone(), *year), order))
-        })
-        .try_collect()?;
-
+) -> Result<()> {
     // Validate graphs at all time slice levels (taking into account process availability and demand)
     for ((region_id, year), base_graph) in commodity_graphs {
         for ts_level in TimeSliceLevel::iter() {
@@ -428,9 +426,33 @@ pub fn validate_commodity_graphs_for_model(
             }
         }
     }
+    Ok(())
+}
 
-    // If all the validation passes, return the commodity ordering
-    Ok(commodity_order)
+/// Determine commodity ordering for each region and year
+///
+/// # Arguments
+///
+/// * `commodity_graphs` - Commodity graphs for each region and year, outputted from `build_commodity_graphs_for_model`
+/// * `commodities` - All commodities with their types and demand specifications
+///
+/// # Returns
+///
+/// A map from `(region, year)` to the ordered list of commodities for investment decisions. The
+/// ordering ensures that leaf-node commodities (those with no outgoing edges) are solved first.
+pub fn solve_investment_order_for_model(
+    commodity_graphs: &HashMap<(RegionID, u32), CommoditiesGraph>,
+    commodities: &CommodityMap,
+) -> HashMap<(RegionID, u32), Vec<InvestmentSet>> {
+    commodity_graphs
+        .iter()
+        .map(|((region_id, year), graph)| {
+            (
+                (region_id.clone(), *year),
+                solve_investment_order(graph, commodities),
+            )
+        })
+        .collect()
 }
 
 /// Gets custom DOT attributes for edges in a commodity graph
@@ -473,7 +495,10 @@ mod tests {
     use std::rc::Rc;
 
     #[rstest]
-    fn test_topo_sort_linear_graph(sed_commodity: Commodity, svd_commodity: Commodity) {
+    fn test_solve_investment_order_linear_graph(
+        sed_commodity: Commodity,
+        svd_commodity: Commodity,
+    ) {
         // Create a simple linear graph: A -> B -> C
         let mut graph = Graph::new();
 
@@ -491,17 +516,18 @@ mod tests {
         commodities.insert("B".into(), Rc::new(sed_commodity));
         commodities.insert("C".into(), Rc::new(svd_commodity));
 
-        let result = topo_sort_commodities(&graph, &commodities).unwrap();
+        let result = solve_investment_order(&graph, &commodities);
 
         // Expected order: C, B, A (leaf nodes first)
+        // No cycles, so all investment sets should be `Single`
         assert_eq!(result.len(), 3);
-        assert_eq!(result[0], "C".into());
-        assert_eq!(result[1], "B".into());
-        assert_eq!(result[2], "A".into());
+        assert_eq!(result[0], InvestmentSet::Single("C".into()));
+        assert_eq!(result[1], InvestmentSet::Single("B".into()));
+        assert_eq!(result[2], InvestmentSet::Single("A".into()));
     }
 
     #[rstest]
-    fn test_topo_sort_cyclic_graph(sed_commodity: Commodity) {
+    fn test_solve_investment_order_cyclic_graph(sed_commodity: Commodity) {
         // Create a simple cyclic graph: A -> B -> A
         let mut graph = Graph::new();
 
@@ -517,11 +543,14 @@ mod tests {
         commodities.insert("A".into(), Rc::new(sed_commodity.clone()));
         commodities.insert("B".into(), Rc::new(sed_commodity));
 
-        // This should return an error due to the cycle
-        // The error message should flag commodity B
-        // Note: A is also involved in the cycle, but B is flagged as it is encountered first
-        let result = topo_sort_commodities(&graph, &commodities);
-        assert_error!(result, "Cycle detected in commodity graph for commodity B");
+        let result = solve_investment_order(&graph, &commodities);
+
+        // Should be a single `Cycle` investment set containing both commodities
+        assert_eq!(result.len(), 1);
+        assert_eq!(
+            result[0],
+            InvestmentSet::Cycle(vec!["A".into(), "B".into()])
+        );
     }
 
     #[rstest]
@@ -548,8 +577,7 @@ mod tests {
         graph.add_edge(node_c, node_d, GraphEdge::Demand);
 
         // Validate the graph at DayNight level
-        let result = validate_commodities_graph(&graph, &commodities, TimeSliceLevel::Annual);
-        assert!(result.is_ok());
+        assert!(validate_commodities_graph(&graph, &commodities, TimeSliceLevel::Annual).is_ok());
     }
 
     #[rstest]
@@ -574,8 +602,10 @@ mod tests {
         graph.add_edge(node_a, node_b, GraphEdge::Primary("process2".into()));
 
         // Validate the graph at DayNight level
-        let result = validate_commodities_graph(&graph, &commodities, TimeSliceLevel::DayNight);
-        assert_error!(result, "SVD commodity A cannot be an input to a process");
+        assert_error!(
+            validate_commodities_graph(&graph, &commodities, TimeSliceLevel::DayNight),
+            "SVD commodity A cannot be an input to a process"
+        );
     }
 
     #[rstest]
@@ -592,8 +622,10 @@ mod tests {
         graph.add_edge(node_a, node_b, GraphEdge::Demand);
 
         // Validate the graph at DayNight level
-        let result = validate_commodities_graph(&graph, &commodities, TimeSliceLevel::DayNight);
-        assert_error!(result, "SVD commodity A is demanded but has no producers");
+        assert_error!(
+            validate_commodities_graph(&graph, &commodities, TimeSliceLevel::DayNight),
+            "SVD commodity A is demanded but has no producers"
+        );
     }
 
     #[rstest]
@@ -611,9 +643,8 @@ mod tests {
         graph.add_edge(node_b, node_a, GraphEdge::Primary("process1".into()));
 
         // Validate the graph at DayNight level
-        let result = validate_commodities_graph(&graph, &commodities, TimeSliceLevel::DayNight);
         assert_error!(
-            result,
+            validate_commodities_graph(&graph, &commodities, TimeSliceLevel::DayNight),
             "SED commodity B may be consumed but has no producers"
         );
     }
@@ -639,9 +670,8 @@ mod tests {
         graph.add_edge(node_a, node_c, GraphEdge::Primary("process2".into()));
 
         // Validate the graph at DayNight level
-        let result = validate_commodities_graph(&graph, &commodities, TimeSliceLevel::DayNight);
         assert_error!(
-            result,
+            validate_commodities_graph(&graph, &commodities, TimeSliceLevel::DayNight),
             "OTH commodity A cannot have both producers and consumers"
         );
     }

src/input.rs

@@ -1,7 +1,8 @@
 //! Common routines for handling input data.
 use crate::asset::AssetPool;
 use crate::graph::{
-    CommoditiesGraph, build_commodity_graphs_for_model, validate_commodity_graphs_for_model,
+    CommoditiesGraph, build_commodity_graphs_for_model, solve_investment_order_for_model,
+    validate_commodity_graphs_for_model,
 };
 use crate::id::{HasID, IDLike};
 use crate::model::{Model, ModelParameters};
@@ -232,15 +233,17 @@ pub fn load_model<P: AsRef<Path>>(model_dir: P) -> Result<(Model, AssetPool)> {
     let assets = read_assets(model_dir.as_ref(), &agent_ids, &processes, &region_ids)?;
 
     // Build and validate commodity graphs for all regions and years
-    // This gives us the commodity order for each region/year which is passed to the model
-    let commodity_graphs = build_commodity_graphs_for_model(&processes, &region_ids, years)?;
-    let commodity_order = validate_commodity_graphs_for_model(
+    let commodity_graphs = build_commodity_graphs_for_model(&processes, &region_ids, years);
+    validate_commodity_graphs_for_model(
         &commodity_graphs,
         &processes,
         &commodities,
         &time_slice_info,
     )?;
 
+    // Solve investment order for each region/year
+    let investment_order = solve_investment_order_for_model(&commodity_graphs, &commodities);
+
     let model_path = model_dir
         .as_ref()
         .canonicalize()
@@ -253,7 +256,7 @@ pub fn load_model<P: AsRef<Path>>(model_dir: P) -> Result<(Model, AssetPool)> {
         processes,
         time_slice_info,
         regions,
-        commodity_order,
+        investment_order,
     };
     Ok((model, AssetPool::new(assets)))
 }
@@ -284,7 +287,7 @@ pub fn load_commodity_graphs<P: AsRef<Path>>(
         years,
     )?;
 
-    let commodity_graphs = build_commodity_graphs_for_model(&processes, &region_ids, years)?;
+    let commodity_graphs = build_commodity_graphs_for_model(&processes, &region_ids, years);
     Ok(commodity_graphs)
 }