Repository: Gurobi/gurobi-optimods

Issue:
Two edges between two nodes not properly calculated
**Describe the bug**
When solving a min cost flow problem, and if two or more edges in the same direction exist between two nodes, only the last defined edge use used. This causes the wrong solution to be found in many cases.

Partially this is because DiGraph does not support this feature, but MultiDiGraph does, which is not supported by gurobi optimods. Pandas also has a similar problem because the from and to columns as set as keys.

**To Reproduce**
```python
from gurobi_optimods import datasets
from gurobi_optimods.min_cost_flow import min_cost_flow_networkx
from networkx import DiGraph

G: DiGraph = datasets.simple_graph_networkx()
G.add_edge(2, 3, capacity=10, cost=-100)
G.add_edge(2, 3, capacity=1, cost=100)

obj, sol = min_cost_flow_networkx(G)
print(list(sol.edges(data=True)))
```

The first edge that is added is not used because it is overwritten by the second edge added.

**Expected behavior**
All three edges would be considered for flow rather than just the last flow defined for that path.

**Additional context**
gurobi-optimods v2.1.0
python v3.10.13

Discussion:
@jsfreischuetz: Draft of one possible solution with #156 @simonbowly: Thanks for raising this @jsfreischuetz. Definitely would make sense to handle these cases, I'll take a look at your PR.

Solution patch:
diff --git a/src/gurobi_optimods/datasets.py b/src/gurobi_optimods/datasets.py
index 862107f..e37d903 100644
--- a/src/gurobi_optimods/datasets.py
+++ b/src/gurobi_optimods/datasets.py
@@ -137,14 +137,21 @@ def load_portfolio():
 
 
 def _convert_pandas_to_digraph(
-    edge_data, node_data, capacity=True, cost=True, demand=True
+    edge_data,
+    node_data,
+    capacity=True,
+    cost=True,
+    demand=True,
+    use_multigraph=False,
 ):
     """
-    Convert from a pandas DataFrame to a networkx.DiGraph with the appropriate
+    Convert from a pandas DataFrame to a networkx.MultiDiGraph with the appropriate
     attributes. For edges: `capacity`, and `cost`. For nodes: `demand`.
     """
+    graph_type = nx.MultiDiGraph if use_multigraph else nx.DiGraph
+
     G = nx.from_pandas_edgelist(
-        edge_data.reset_index(), create_using=nx.DiGraph(), edge_attr=True
+        edge_data.reset_index(), create_using=graph_type(), edge_attr=True
     )
     if demand:
         for i, d in node_data.iterrows():
diff --git a/src/gurobi_optimods/min_cost_flow.py b/src/gurobi_optimods/min_cost_flow.py
index 39267b2..00a4f59 100644
--- a/src/gurobi_optimods/min_cost_flow.py
+++ b/src/gurobi_optimods/min_cost_flow.py
@@ -61,14 +61,17 @@ def min_cost_flow_pandas(
     with create_env() as env, gp.Model(env=env) as model:
         model.ModelSense = GRB.MINIMIZE
 
+        source_label, target_label = arc_data.index.names
+
+        arc_data = arc_data.reset_index()
+
         arc_df = arc_data.gppd.add_vars(model, ub="capacity", obj="cost", name="flow")
 
-        source_label, target_label = arc_data.index.names
         balance_df = (
             pd.DataFrame(
                 {
-                    "inflow": arc_df["flow"].groupby(target_label).sum(),
-                    "outflow": arc_df["flow"].groupby(source_label).sum(),
+                    "inflow": arc_df.groupby(target_label)["flow"].sum(),
+                    "outflow": arc_df.groupby(source_label)["flow"].sum(),
                     "demand": demand_data["demand"],
                 }
             )
@@ -84,6 +87,7 @@ def min_cost_flow_pandas(
         if model.Status in [GRB.INFEASIBLE, GRB.INF_OR_UNBD]:
             raise ValueError("Unsatisfiable flows")
 
+        arc_df = arc_df.set_index([source_label, target_label])
         return model.ObjVal, arc_df["flow"].gppd.X
 
 
@@ -171,27 +175,34 @@ def min_cost_flow_networkx(G, *, create_env):
     logger.info(
         f"Solving min-cost flow with {len(G.nodes)} nodes and {len(G.edges)} edges"
     )
-
     with create_env() as env, gp.Model(env=env) as model:
+        use_multigraph = isinstance(G, nx.MultiGraph)
+
+        G = nx.MultiDiGraph(G)
+
         edges, capacities, costs = gp.multidict(
-            {(i, j): [d["capacity"], d["cost"]] for i, j, d in G.edges(data=True)}
+            {
+                (i, j, idx): [d["capacity"], d["cost"]]
+                for i, j, idx, d in G.edges(data=True, keys=True)
+            }
         )
+
         nodes = list(G.nodes(data=True))
         x = {
-            (i, j): model.addVar(
-                name=f"flow[{i},{j}]",
-                ub=capacities[i, j],
-                obj=costs[i, j],
+            (i, j, idx): model.addVar(
+                name=f"flow[{i},{j},{idx}]",
+                ub=capacities[i, j, idx],
+                obj=costs[i, j, idx],
             )
-            for i, j in edges
+            for i, j, idx in edges
         }
 
         flow_constrs = {}
         for n, data in nodes:
             flow_constrs[n] = model.addConstr(
                 (
-                    gp.quicksum(x[j, n] for j in G.predecessors(n))
-                    - gp.quicksum(x[n, j] for j in G.successors(n))
+                    gp.quicksum(x[ie] for ie in G.in_edges(n, keys=True))
+                    - gp.quicksum(x[oe] for oe in G.out_edges(n, keys=True))
                     == data["demand"]
                 ),
                 name=f"flow_balance[{n}]",
@@ -203,7 +214,7 @@ def min_cost_flow_networkx(G, *, create_env):
             raise ValueError("Unsatisfiable flows")
 
         # Create a new Graph with selected edges in the matching
-        resulting_flow = nx.DiGraph()
+        resulting_flow = nx.MultiDiGraph() if use_multigraph else nx.DiGraph()
         resulting_flow.add_nodes_from(nodes)
         resulting_flow.add_edges_from(
             [(edge[0], edge[1], {"flow": v.X}) for edge, v in x.items() if v.X > 0.1]


Test patch:
diff --git a/tests/test_graph_utils.py b/tests/test_graph_utils.py
index 3b49e04..b3c9e25 100644
--- a/tests/test_graph_utils.py
+++ b/tests/test_graph_utils.py
@@ -3,25 +3,24 @@ import numpy as np
 
 def check_solution_pandas(solution, candidates):
     # Checks whether the solution (`pd.Series`) matches any of the list of
-    # candidates (containing `dict`)
-    if any(solution.to_dict() == c for c in candidates):
-        return True
-    return False
+    # candidates. Each candidate is a list of tuples ((i, j), v) tuples,
+    # compare with the solution in sorted order.
+    solution_list = sorted(solution.items())
+    return any(solution_list == sorted(candidate) for candidate in candidates)
 
 
 def check_solution_scipy(solution, candidates):
     # Checks whether the solution (`sp.sparray`) matches any of the list of
     # candidates (containing `np.ndarray`)
     arr = solution.toarray()
-    if any(np.array_equal(arr, c) for c in candidates):
-        return True
-    return False
+    return any(np.array_equal(arr, c) for c in candidates)
 
 
 def check_solution_networkx(solution, candidates):
     # Checks whether the solution (`nx.DiGraph`) matches any of the list of
     # candidates (containing tuples dict `{(i, j): data}`)
-    sol_dict = {(i, j): d for i, j, d in solution.edges(data=True)}
-    if any(sol_dict == c for c in candidates):
-        return True
-    return False
+    solution_list = sorted(
+        [((i, j), data["flow"]) for i, j, data in solution.edges(data=True)],
+        key=str,
+    )
+    return any(solution_list == sorted(candidate, key=str) for candidate in candidates)
diff --git a/tests/test_max_flow.py b/tests/test_max_flow.py
index 536343c..15afbc0 100644
--- a/tests/test_max_flow.py
+++ b/tests/test_max_flow.py
@@ -31,23 +31,23 @@ class TestMaxFlow(unittest.TestCase):
         obj, sol = max_flow(edge_data, 0, 5)
         sol = sol[sol > 0]
         self.assertEqual(obj, self.expected_max_flow)
-        candidate = {
-            (0, 1): 1.0,
-            (0, 2): 2.0,
-            (1, 3): 1.0,
-            (2, 3): 1.0,
-            (2, 4): 1.0,
-            (3, 5): 2.0,
-            (4, 5): 1.0,
-        }
-        candidate2 = {
-            (0, 1): 1.0,
-            (0, 2): 2.0,
-            (1, 3): 1.0,
-            (2, 4): 2.0,
-            (3, 5): 1.0,
-            (4, 5): 2.0,
-        }
+        candidate = [
+            ((0, 1), 1.0),
+            ((0, 2), 2.0),
+            ((1, 3), 1.0),
+            ((2, 3), 1.0),
+            ((2, 4), 1.0),
+            ((3, 5), 2.0),
+            ((4, 5), 1.0),
+        ]
+        candidate2 = [
+            ((0, 1), 1.0),
+            ((0, 2), 2.0),
+            ((1, 3), 1.0),
+            ((2, 4), 2.0),
+            ((3, 5), 1.0),
+            ((4, 5), 2.0),
+        ]
         self.assertTrue(check_solution_pandas(sol, [candidate, candidate2]))
 
     def test_empty_pandas(self):
@@ -84,23 +84,23 @@ class TestMaxFlow(unittest.TestCase):
         G = datasets.simple_graph_networkx()
         obj, sol = max_flow(G, 0, 5)
         self.assertEqual(obj, self.expected_max_flow)
-        candidate = {
-            (0, 1): {"flow": 1},
-            (0, 2): {"flow": 2},
-            (1, 3): {"flow": 1},
-            (2, 4): {"flow": 2},
-            (3, 5): {"flow": 1},
-            (4, 5): {"flow": 2},
-        }
-        candidate2 = {
-            (0, 1): {"flow": 1.0},
-            (0, 2): {"flow": 2.0},
-            (1, 3): {"flow": 1.0},
-            (2, 3): {"flow": 1.0},
-            (2, 4): {"flow": 1.0},
-            (3, 5): {"flow": 2.0},
-            (4, 5): {"flow": 1.0},
-        }
+        candidate = [
+            ((0, 1), 1),
+            ((0, 2), 2),
+            ((1, 3), 1),
+            ((2, 4), 2),
+            ((3, 5), 1),
+            ((4, 5), 2),
+        ]
+        candidate2 = [
+            ((0, 1), 1.0),
+            ((0, 2), 2.0),
+            ((1, 3), 1.0),
+            ((2, 3), 1.0),
+            ((2, 4), 1.0),
+            ((3, 5), 2.0),
+            ((4, 5), 1.0),
+        ]
         self.assertTrue(check_solution_networkx(sol, [candidate, candidate2]))
 
     @unittest.skipIf(nx is None, "networkx is not installed")
@@ -120,16 +120,16 @@ class TestMaxFlow2(unittest.TestCase):
         obj, sol = max_flow(edge_data, 0, 4)
         sol = sol[sol > 0]
         self.assertEqual(obj, self.expected_max_flow)
-        candidate = {
-            (0, 1): 15.0,
-            (0, 2): 8.0,
-            (1, 3): 4.0,
-            (1, 2): 1.0,
-            (1, 4): 10.0,
-            (2, 3): 4.0,
-            (2, 4): 5.0,
-            (3, 4): 8.0,
-        }
+        candidate = [
+            ((0, 1), 15.0),
+            ((0, 2), 8.0),
+            ((1, 3), 4.0),
+            ((1, 2), 1.0),
+            ((1, 4), 10.0),
+            ((2, 3), 4.0),
+            ((2, 4), 5.0),
+            ((3, 4), 8.0),
+        ]
         self.assertTrue(check_solution_pandas(sol, [candidate]))
 
     def test_scipy(self):
@@ -152,23 +152,23 @@ class TestMaxFlow2(unittest.TestCase):
         G = load_graph2_networkx()
         obj, sol = max_flow(G, 0, 4)
         self.assertEqual(obj, self.expected_max_flow)
-        candidate = {
-            (0, 1): {"flow": 15.0},
-            (0, 2): {"flow": 8.0},
-            (1, 3): {"flow": 4.0},
-            (1, 2): {"flow": 1.0},
-            (1, 4): {"flow": 10.0},
-            (2, 3): {"flow": 4.0},
-            (2, 4): {"flow": 5.0},
-            (3, 4): {"flow": 8.0},
-        }
-        candidate2 = {
-            (0, 1): {"flow": 15},
-            (0, 2): {"flow": 8},
-            (1, 2): {"flow": 1},
-            (1, 3): {"flow": 4},
-            (1, 4): {"flow": 10},
-            (2, 3): {"flow": 9},
-            (3, 4): {"flow": 13},
-        }
+        candidate = [
+            ((0, 1), 15.0),
+            ((0, 2), 8.0),
+            ((1, 3), 4.0),
+            ((1, 2), 1.0),
+            ((1, 4), 10.0),
+            ((2, 3), 4.0),
+            ((2, 4), 5.0),
+            ((3, 4), 8.0),
+        ]
+        candidate2 = [
+            ((0, 1), 15),
+            ((0, 2), 8),
+            ((1, 2), 1),
+            ((1, 3), 4),
+            ((1, 4), 10),
+            ((2, 3), 9),
+            ((3, 4), 13),
+        ]
         self.assertTrue(check_solution_networkx(sol, [candidate, candidate2]))
diff --git a/tests/test_min_cost_flow.py b/tests/test_min_cost_flow.py
index e4df9f5..eb45e58 100644
--- a/tests/test_min_cost_flow.py
+++ b/tests/test_min_cost_flow.py
@@ -41,6 +41,20 @@ node_data2 = """
 4,15
 """
 
+edge_data3 = """
+source,target,capacity,cost
+0,1,15,4
+0,2,8,4
+1,3,4,2
+1,2,20,2
+1,4,10,6
+2,3,15,1
+2,4,5,3
+3,4,20,2
+2,3,1,-100
+2,3,1,10
+"""
+
 
 def load_graph2_pandas():
     return (
@@ -59,13 +73,33 @@ def load_graph2_scipy():
     return datasets._convert_pandas_to_scipy(edge_data, node_data)
 
 
+def load_graph3_pandas():
+    return (
+        pd.read_csv(io.StringIO(edge_data3)).set_index(["source", "target"]),
+        pd.read_csv(io.StringIO(node_data2), index_col=0),
+    )
+
+
+def load_graph3_networkx(use_multigraph):
+    edge_data, node_data = load_graph3_pandas()
+    return datasets._convert_pandas_to_digraph(
+        edge_data, node_data, use_multigraph=use_multigraph
+    )
+
+
 class TestMinCostFlow(unittest.TestCase):
     def test_pandas(self):
         edge_data, node_data = datasets.simple_graph_pandas()
         cost, sol = mcf.min_cost_flow_pandas(edge_data, node_data)
         sol = sol[sol > 0]
         self.assertEqual(cost, 31)
-        candidate = {(0, 1): 1.0, (0, 2): 1.0, (1, 3): 1.0, (2, 4): 2.0, (4, 5): 2.0}
+        candidate = [
+            ((0, 1), 1.0),
+            ((0, 2), 1.0),
+            ((1, 3), 1.0),
+            ((2, 4), 2.0),
+            ((4, 5), 2.0),
+        ]
         self.assertIsInstance(sol, pd.Series)
         self.assertTrue(check_solution_pandas(sol, [candidate]))
 
@@ -97,13 +131,13 @@ class TestMinCostFlow(unittest.TestCase):
         G = datasets.simple_graph_networkx()
         cost, sol = mcf.min_cost_flow_networkx(G)
         self.assertEqual(cost, 31)
-        expected = {
-            (0, 1): {"flow": 1.0},
-            (0, 2): {"flow": 1.0},
-            (1, 3): {"flow": 1.0},
-            (2, 4): {"flow": 2.0},
-            (4, 5): {"flow": 2.0},
-        }
+        expected = [
+            ((0, 1), 1.0),
+            ((0, 2), 1.0),
+            ((1, 3), 1.0),
+            ((2, 4), 2.0),
+            ((4, 5), 2.0),
+        ]
         self.assertIsInstance(sol, nx.Graph)
         self.assertTrue(check_solution_networkx(sol, [expected]))
 
@@ -113,13 +147,13 @@ class TestMinCostFlow(unittest.TestCase):
         G = nx.relabel_nodes(G, {0: "s", 5: "t"})
         cost, sol = mcf.min_cost_flow_networkx(G)
         self.assertEqual(cost, 31)
-        expected = {
-            ("s", 1): {"flow": 1.0},
-            ("s", 2): {"flow": 1.0},
-            (1, 3): {"flow": 1.0},
-            (2, 4): {"flow": 2.0},
-            (4, "t"): {"flow": 2.0},
-        }
+        expected = [
+            (("s", 1), 1.0),
+            (("s", 2), 1.0),
+            ((1, 3), 1.0),
+            ((2, 4), 2.0),
+            ((4, "t"), 2.0),
+        ]
         self.assertIsInstance(sol, nx.Graph)
         self.assertTrue(check_solution_networkx(sol, [expected]))
 
@@ -130,24 +164,24 @@ class TestMinCostFlow2(unittest.TestCase):
         cost, sol = mcf.min_cost_flow_pandas(edge_data, node_data)
         sol = sol[sol > 0]
         self.assertEqual(cost, 150)
-        candidate = {
-            (0, 1): 12.0,
-            (0, 2): 8.0,
-            (1, 3): 4.0,
-            (1, 2): 8.0,
-            (2, 3): 15.0,
-            (2, 4): 1.0,
-            (3, 4): 14.0,
-        }
-        candidate2 = {
-            (0, 1): 12.0,
-            (0, 2): 8.0,
-            (1, 3): 4.0,
-            (1, 2): 8.0,
-            (2, 3): 11.0,
-            (2, 4): 5.0,
-            (3, 4): 10.0,
-        }
+        candidate = [
+            ((0, 1), 12.0),
+            ((0, 2), 8.0),
+            ((1, 3), 4.0),
+            ((1, 2), 8.0),
+            ((2, 3), 15.0),
+            ((2, 4), 1.0),
+            ((3, 4), 14.0),
+        ]
+        candidate2 = [
+            ((0, 1), 12.0),
+            ((0, 2), 8.0),
+            ((1, 3), 4.0),
+            ((1, 2), 8.0),
+            ((2, 3), 11.0),
+            ((2, 4), 5.0),
+            ((3, 4), 10.0),
+        ]
         self.assertTrue(check_solution_pandas(sol, [candidate, candidate2]))
 
     def test_scipy(self):
@@ -169,24 +203,24 @@ class TestMinCostFlow2(unittest.TestCase):
         G = load_graph2_networkx()
         cost, sol = mcf.min_cost_flow_networkx(G)
         self.assertEqual(cost, 150)
-        candidate = {
-            (0, 1): {"flow": 12.0},
-            (0, 2): {"flow": 8.0},
-            (1, 2): {"flow": 8.0},
-            (1, 3): {"flow": 4.0},
-            (2, 3): {"flow": 11.0},
-            (2, 4): {"flow": 5.0},
-            (3, 4): {"flow": 10.0},
-        }
-        candidate2 = {
-            (0, 1): {"flow": 12.0},
-            (0, 2): {"flow": 8.0},
-            (1, 3): {"flow": 4.0},
-            (1, 2): {"flow": 8.0},
-            (2, 3): {"flow": 15.0},
-            (2, 4): {"flow": 1.0},
-            (3, 4): {"flow": 14.0},
-        }
+        candidate = [
+            ((0, 1), 12.0),
+            ((0, 2), 8.0),
+            ((1, 2), 8.0),
+            ((1, 3), 4.0),
+            ((2, 3), 11.0),
+            ((2, 4), 5.0),
+            ((3, 4), 10.0),
+        ]
+        candidate2 = [
+            ((0, 1), 12.0),
+            ((0, 2), 8.0),
+            ((1, 3), 4.0),
+            ((1, 2), 8.0),
+            ((2, 3), 15.0),
+            ((2, 4), 1.0),
+            ((3, 4), 14.0),
+        ]
         self.assertTrue(check_solution_networkx(sol, [candidate, candidate2]))
 
     @unittest.skipIf(nx is None, "networkx is not installed")
@@ -195,22 +229,60 @@ class TestMinCostFlow2(unittest.TestCase):
         G = nx.relabel_nodes(G, {0: "s", 4: "t"})
         cost, sol = mcf.min_cost_flow_networkx(G)
         self.assertEqual(cost, 150)
-        candidate = {
-            ("s", 1): {"flow": 12.0},
-            ("s", 2): {"flow": 8.0},
-            (1, 2): {"flow": 8.0},
-            (1, 3): {"flow": 4.0},
-            (2, 3): {"flow": 11.0},
-            (2, "t"): {"flow": 5.0},
-            (3, "t"): {"flow": 10.0},
-        }
+        candidate = [
+            (("s", 1), 12.0),
+            (("s", 2), 8.0),
+            ((1, 2), 8.0),
+            ((1, 3), 4.0),
+            ((2, 3), 11.0),
+            ((2, "t"), 5.0),
+            ((3, "t"), 10.0),
+        ]
         candidate2 = {
-            ("s", 1): {"flow": 12.0},
-            ("s", 2): {"flow": 8.0},
-            (1, 3): {"flow": 4.0},
-            (1, 2): {"flow": 8.0},
-            (2, 3): {"flow": 15.0},
-            (2, "t"): {"flow": 1.0},
-            (3, "t"): {"flow": 14.0},
+            (("s", 1), 12.0),
+            (("s", 2), 8.0),
+            ((1, 3), 4.0),
+            ((1, 2), 8.0),
+            ((2, 3), 15.0),
+            ((2, "t"), 1.0),
+            ((3, "t"), 14.0),
         }
         self.assertTrue(check_solution_networkx(sol, [candidate, candidate2]))
+
+
+class TestMinCostFlow3(unittest.TestCase):
+    def test_pandas(self):
+        edge_data, node_data = load_graph3_pandas()
+        cost, sol = mcf.min_cost_flow_pandas(edge_data, node_data)
+        sol = sol[sol > 0]
+        self.assertEqual(cost, 49.0)
+
+        candidate = [
+            ((0, 1), 12.0),
+            ((0, 2), 8.0),
+            ((1, 3), 4.0),
+            ((1, 2), 8.0),
+            ((2, 3), 10.0),
+            ((2, 4), 5.0),
+            ((3, 4), 10.0),
+            ((2, 3), 1.0),
+        ]
+        self.assertTrue(check_solution_pandas(sol, [candidate]))
+
+    @unittest.skipIf(nx is None, "networkx is not installed")
+    def test_networkx(self):
+        G = load_graph3_networkx(use_multigraph=True)
+        cost, sol = mcf.min_cost_flow_networkx(G)
+        self.assertEqual(cost, 49.0)
+        candidate = [
+            ((0, 1), 12.0),
+            ((0, 2), 8.0),
+            ((1, 3), 4.0),
+            ((1, 2), 8.0),
+            ((2, 3), 10.0),
+            ((2, 3), 1.0),
+            ((2, 4), 5.0),
+            ((3, 4), 10.0),
+        ]
+
+        self.assertTrue(check_solution_networkx(sol, [candidate]))