diff --git a/docs/source/api.rst b/docs/source/api.rst index 1e3bb187c..2715a955c 100644 --- a/docs/source/api.rst +++ b/docs/source/api.rst @@ -4,9 +4,15 @@ API Reference .. automodule:: gurobi_optimods.bipartite_matching :members: maximum_bipartite_matching +.. automodule:: gurobi_optimods.max_flow + :members: max_flow + .. automodule:: gurobi_optimods.min_cost_flow :members: min_cost_flow_pandas, min_cost_flow_networkx, min_cost_flow_scipy +.. automodule:: gurobi_optimods.min_cut + :members: min_cut, MinCutResult + .. automodule:: gurobi_optimods.mwis :members: maximum_weighted_independent_set diff --git a/docs/source/gallery.rst b/docs/source/gallery.rst index bb0b76725..c12ef84ad 100644 --- a/docs/source/gallery.rst +++ b/docs/source/gallery.rst @@ -18,6 +18,12 @@ The OptiMods Gallery :text-align: center :img-top: mods/icons/lad-regression.png + .. grid-item-card:: Maximum-Flow/Mininum-Cut + :link: mods/max-flow-min-cut + :link-type: doc + :text-align: center + :img-top: mods/figures/max-flow-min-cut.png + .. grid-item-card:: Minimum-Cost Flow :link: mods/min-cost-flow :link-type: doc @@ -68,6 +74,7 @@ The OptiMods Gallery mods/bipartite-matching mods/lad-regression mods/min-cost-flow + mods/max-flow-min-cut mods/mwis mods/opf/opf mods/portfolio diff --git a/docs/source/mods/figures/max-flow-min-cut.png b/docs/source/mods/figures/max-flow-min-cut.png new file mode 100644 index 000000000..cd5a20e45 Binary files /dev/null and b/docs/source/mods/figures/max-flow-min-cut.png differ diff --git a/docs/source/mods/figures/max-flow.png b/docs/source/mods/figures/max-flow.png new file mode 100644 index 000000000..91439a130 Binary files /dev/null and b/docs/source/mods/figures/max-flow.png differ diff --git a/docs/source/mods/figures/min-cut.png b/docs/source/mods/figures/min-cut.png new file mode 100644 index 000000000..e2a95c717 Binary files /dev/null and b/docs/source/mods/figures/min-cut.png differ diff --git a/docs/source/mods/max-flow-min-cut.rst b/docs/source/mods/max-flow-min-cut.rst new file mode 100644 index 000000000..5dfe95c53 --- /dev/null +++ b/docs/source/mods/max-flow-min-cut.rst @@ -0,0 +1,310 @@ +Maximum Flow / Minimum Cut +========================== + +The maximum flow problem finds the total flow that can be routed through a +capacitated network from a given source vertex to a given sink vertex. The +minimum cut problem finds a set of edges in the same graph with minimal combined +capacity that, when removed, would disconnect the source from the sink. The two +problems are related through duality: the maximum flow is equal to the capacity +of the minimum cut. + +The first algorithm proposed to solve this problem was the Ford-Fulkerson +algorithm :footcite:p:`ford1956maximal`. :footcite:t:`goldberg1988new` later +proposed the famous push-relabel algorithm, and more recently, +:footcite:t:`orlin2013max` and other authors have proposed polynomial-time +algorithms. The maximum flow problem is a special case of the minimum cost flow +problem, so it can also be solved efficiently using the network simplex +algorithm :footcite:p:`caliskan2012faster` making a linear programming (LP) +solver a suitable approach. + +Problem Specification +--------------------- + +Let :math:`G` be a graph with set of vertices :math:`V` and edges :math:`E`. +Each edge :math:`(i,j)\in E` has capacity :math:`B_{ij}\in\mathbb{R}`. Given a +source vertex :math:`s\in V` and a sink vertex :math:`t\in V` the two problems +can be stated as follows: + +- Maximum Flow: Find the flow with maximum value from source to sink such that + the flow is capacity feasible. +- Minimum Cut: Find the set of edges which disconnects the source and sink such + that the total capacity of the cut is minimised. + +.. dropdown:: Background: Optimization Model + + Let us define a set of continuous variables :math:`x_{ij}` to represent + the amount of non-negative (:math:`\geq 0`) flow going through an edge + :math:`(i,j) \in E`. + + The mathematical formulation for the maximum flow can be stated as: + + .. math:: + + \begin{alignat}{2} + \max \quad & \sum_{j \in \delta^+(s)} x_{sj} \\ + \mbox{s.t.} \quad & \sum_{j \in \delta^-(i)} x_{ji} = \sum_{j \in \delta^+(i)} x_{ij} & \quad\forall i \in V\setminus\{s,t\} \\ + & 0 \leq x_{ij} \le B_{ij} & \forall (i, j) \in E \\ + \end{alignat} + + where :math:`\delta^+(\cdot)` (:math:`\delta^-(\cdot)`) denotes the + outgoing (incoming) neighbours of a given vertex. + + The objective maximises the total outgoing flow from the source vertex. The + first set of constraint ensure flow balance for all vertices. That is, for a + given node that is neither the source or the sink, the outgoing flow must be + equal to the incoming flow. The last set of constraints ensure + non-negativity of the variables and that the capacity of each edge is not + exceeded. + + The minimum cut problem formulation is the dual of the maximum flow + formulation. The minimum cut is found by solving the maximum flow problem + and extracting the constraint dual values. Specifically, the edges in the + cutset are those whose capacity constraints have non-zero dual values. The + graph partitions formed by removing the edges in the cutset can then be + round by following the predecessor and successor vertices of these edges. + +| + +Code and Inputs +--------------- + +This Mod accepts input graphs of any of the following types: + +* pandas: using a ``pd.DataFrame``; +* NetworkX: using a ``nx.DiGraph`` or ``nx.Graph``; +* SciPy.sparse: using a ``sp.sparray`` matrix. + +An example of these inputs with their respective requirements is shown below. + +.. tabs:: + + .. group-tab:: pandas + + .. doctest:: pandas + :options: +NORMALIZE_WHITESPACE + + >>> from gurobi_optimods import datasets + >>> edge_data, _ = datasets.simple_graph_pandas() + >>> edge_data[["capacity"]] + capacity + source target + 0 1 2 + 2 2 + 1 3 1 + 2 3 1 + 4 2 + 3 5 2 + 4 5 2 + + The ``edge_data`` DataFrame is indexed by ``source`` and ``target`` nodes + and contains columns labelled ``capacity`` with the edge attributes. + + .. group-tab:: NetworkX + + .. doctest:: networkx + :options: +NORMALIZE_WHITESPACE + + >>> from gurobi_optimods import datasets + >>> G = datasets.simple_graph_networkx() + >>> for u, v, capacity in G.edges.data("capacity"): + ... print(f"{u} -> {v}: {capacity = }") + 0 -> 1: capacity = 2 + 0 -> 2: capacity = 2 + 1 -> 3: capacity = 1 + 2 -> 3: capacity = 1 + 2 -> 4: capacity = 2 + 3 -> 5: capacity = 2 + 4 -> 5: capacity = 2 + + Edges have attributes ``capacity``. + + .. group-tab:: scipy.sparse + + .. doctest:: scipy + :options: +NORMALIZE_WHITESPACE + + >>> from gurobi_optimods import datasets + >>> G, capacities, _, _ = datasets.simple_graph_scipy() + >>> G.data = capacities.data # Copy capacity data + >>> G + <5x6 sparse array of type '' + with 7 stored elements in COOrdinate format> + >>> print(G) + (0, 1) 2 + (0, 2) 2 + (1, 3) 1 + (2, 3) 1 + (2, 4) 2 + (3, 5) 2 + (4, 5) 2 + + We only need the adjacency matrix for the graph (as a sparse array) where + each each entry contains the capacity of the edge. + +| + +Solution +-------- + +Maximum Flow +^^^^^^^^^^^^ + +Let us use the data to solve the maximum flow problem. + +.. tabs:: + + .. group-tab:: pandas + + .. doctest:: pandas + :options: +NORMALIZE_WHITESPACE + + >>> from gurobi_optimods.max_flow import max_flow + >>> obj, flow = max_flow(edge_data, 0, 5, verbose=False) # Find max-flow between nodes 0 and 5 + >>> obj + 3.0 + >>> flow + source target + 0 1 1.0 + 2 2.0 + 1 3 1.0 + 2 3 1.0 + 4 1.0 + 3 5 2.0 + 4 5 1.0 + Name: flow, dtype: float64 + + The ``max_flow`` function returns the cost of the solution as well + as ``pd.Series`` with the flow per edge. Similarly as the input + DataFrame the resulting series is indexed by ``source`` and ``target``. + In this case, the resulting maximum flow has value 3. + + .. group-tab:: NetworkX + + .. doctest:: networkx + :options: +NORMALIZE_WHITESPACE + + >>> from gurobi_optimods.max_flow import max_flow + >>> obj, sol = max_flow(G, 0, 5, verbose=False) + >>> obj + 3.0 + >>> type(sol) + + >>> list(sol.edges(data=True)) + [(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})] + + The ``max_flow`` function returns the cost of the solution + as well as a dictionary indexed by edge with the non-zero flow. + + .. group-tab:: scipy.sparse + + .. doctest:: scipy + :options: +NORMALIZE_WHITESPACE + + >>> from gurobi_optimods.max_flow import max_flow + >>> obj, sol = max_flow(G, 0, 5, verbose=False) + >>> obj + 3.0 + >>> sol + <5x6 sparse array of type '' + with 6 stored elements in COOrdinate format> + >>> print(sol) + (0, 1) 1.0 + (0, 2) 2.0 + (1, 3) 1.0 + (2, 4) 2.0 + (3, 5) 1.0 + (4, 5) 2.0 + + The ``max_flow`` function returns the value of the maximum flow as well a + sparse array with the amount of non-zero flow in each edge in the + solution. + +The solution for this example is shown in the figure below. The edge labels +denote the edge capacity and resulting flow: :math:`x^*_{ij}/B_{ij}`. All +edges in the maximum flow solution carry some flow, totalling at 3.0 at the +sink. + +.. image:: figures/max-flow.png + :width: 600 + :alt: Maximum Flow Solution. + :align: center + +Minimum Cut +^^^^^^^^^^^ + +Let us use the data to solve the minimum cut problem. + +.. tabs:: + + .. group-tab:: pandas + + .. doctest:: pandas + :options: +NORMALIZE_WHITESPACE + + >>> from gurobi_optimods.min_cut import min_cut + >>> res = min_cut(edge_data, 0, 5, verbose=False) + >>> res + MinCutResult(cut_value=3.0, partition=({0, 1}, {2, 3, 4, 5}), cutset={(0, 2), (1, 3)}) + >>> res.cut_value + 3.0 + >>> res.partition + ({0, 1}, {2, 3, 4, 5}) + >>> res.cutset + {(0, 2), (1, 3)} + + The ``min_cut`` function returns a ``MinCutResult`` which contains the + cutset value, the partition of the nodes and the edges in the cutset. + + + .. group-tab:: NetworkX + + .. doctest:: networkx + :options: +NORMALIZE_WHITESPACE + + >>> from gurobi_optimods.min_cut import min_cut + >>> res = min_cut(G, 0, 5, verbose=False) + >>> res + MinCutResult(cut_value=3.0, partition=({0, 1}, {2, 3, 4, 5}), cutset={(0, 2), (1, 3)}) + >>> res.cut_value + 3.0 + >>> res.partition + ({0, 1}, {2, 3, 4, 5}) + >>> res.cutset + {(0, 2), (1, 3)} + + The ``min_cut`` function returns a ``MinCutResult`` which contains the + cutset value, the partition of the nodes and the edges in the cutset. + + .. group-tab:: scipy.sparse + + .. doctest:: scipy + :options: +NORMALIZE_WHITESPACE + + >>> from gurobi_optimods.min_cut import min_cut + >>> res = min_cut(G, 0, 5, verbose=False) + >>> res + MinCutResult(cut_value=3.0, partition=({0, 1}, {2, 3, 4, 5}), cutset={(0, 2), (1, 3)}) + >>> res.cut_value + 3.0 + >>> res.partition + ({0, 1}, {2, 3, 4, 5}) + >>> res.cutset + {(0, 2), (1, 3)} + + The ``min_cut`` function returns a ``MinCutResult`` which contains the + cutset value, the partition of the nodes and the edges in the cutset. + +The solution for the minimum cut problem is shown in the figure below. The edges +in the cutset are shown in blue (with their capacity values shown), and the +nodes in the partitions are shown in blue (nodes 0 and 1) and in green (nodes 2, +3, 4 and 5). The capacity of the minimum cut is :math:`B_{0,2}+B_{1,3}=2+1=3` +which is also the value of the maximum flow. We can also see that if we remove +the blue edges we would be left with a disconnected graph with the two +partitions. + +.. image:: figures/min-cut.png + :width: 600 + :alt: Minimum Cut solution. + :align: center + +.. footbibliography:: diff --git a/docs/source/refs/graphs.bib b/docs/source/refs/graphs.bib index 240343913..418cf2aa6 100644 --- a/docs/source/refs/graphs.bib +++ b/docs/source/refs/graphs.bib @@ -43,3 +43,43 @@ @article{kovacs2015minimum year={2015}, publisher={Taylor \& Francis} } + +@article{ford1956maximal, + title={Maximal flow through a network}, + author={Ford, Lester Randolph and Fulkerson, Delbert R}, + journal={Canadian journal of Mathematics}, + volume={8}, + pages={399--404}, + year={1956}, + publisher={Cambridge University Press} +} + +@article{goldberg1988new, + title={A new approach to the maximum-flow problem}, + author={Goldberg, Andrew V and Tarjan, Robert E}, + journal={Journal of the ACM (JACM)}, + volume={35}, + number={4}, + pages={921--940}, + year={1988}, + publisher={ACM New York, NY, USA} +} + +@inproceedings{orlin2013max, + title={Max flows in O(nm) time, or better}, + author={Orlin, James B}, + booktitle={Proceedings of the forty-fifth annual ACM symposium on Theory of computing}, + pages={765--774}, + year={2013} +} + +@article{caliskan2012faster, + title={A faster polynomial algorithm for the constrained maximum flow problem}, + author={Çalışkan, Cenk}, + journal={Computers \& operations research}, + volume={39}, + number={11}, + pages={2634--2641}, + year={2012}, + publisher={Elsevier} +} diff --git a/src/gurobi_optimods/max_flow.py b/src/gurobi_optimods/max_flow.py new file mode 100644 index 000000000..48a3fa1ac --- /dev/null +++ b/src/gurobi_optimods/max_flow.py @@ -0,0 +1,127 @@ +""" +Maximum Flow +------------ +""" + +import logging + +import numpy as np +import pandas as pd +import scipy.sparse as sp + +try: + import networkx as nx +except ImportError: + nx = None + +from gurobi_optimods.min_cost_flow import ( + min_cost_flow_networkx, + min_cost_flow_pandas, + min_cost_flow_scipy, +) + +logger = logging.getLogger(__name__) + + +def max_flow(graph, source: int, sink: int, **kwargs): + """Solve the maximum flow problem for a given graph. + + Parameters + ---------- + graph : spmatrix or Graph or DataFrame + A graph, specified either as a scipy.sparse adjacency matrix, networkx + graph, or pandas dataframe + source : int or str + The source (or origin) node for the maximum flow. + sink : int or str + The sink (or destination) node for the maximum flow. + + Returns + ------- + flow: float + Maximum feasible flow through the network. + subgraph: spmatrix or Graph or DataFrame + A subgraph of the original graph specifying the flow. + """ + if isinstance(graph, sp.spmatrix): + return _max_flow_scipy(graph, source, sink, **kwargs) + elif isinstance(graph, pd.DataFrame): + return _max_flow_pandas(graph, source, sink, **kwargs) + elif nx is not None and isinstance(graph, nx.Graph): + return _max_flow_networkx(graph, source, sink, **kwargs) + else: + raise ValueError(f"Unknown graph type: {type(graph)}") + + +def _remove_dummy_edge(graph, source, sink): + if isinstance(graph, sp.spmatrix): + graph = graph.tolil() + graph = graph[:-1, :] + elif isinstance(graph, pd.Series) or isinstance(graph, pd.DataFrame): + graph.drop((sink, source), inplace=True) + elif nx is not None and isinstance(graph, nx.Graph): + graph.remove_edge(sink, source) + return graph + + +def _max_flow_pandas(arc_data, source, sink, **kwargs): + f, t = arc_data.index.names + arc_data["cost"] = [0] * len(arc_data) + # Find maximum flow through (minumum of sum of all outgoing/incoming + # capacities at the source/sink) + max_flow = min( + arc_data.loc[([source], slice(None)),]["capacity"].sum(), + arc_data.loc[(slice(None), [sink]),]["capacity"].sum(), + ) + # Add dummy edge + arc_data = pd.concat( + [ + arc_data, + pd.DataFrame( + [{f: sink, t: source, "capacity": max_flow, "cost": -1}] + ).set_index([f, t]), + ] + ) + demand_data = pd.DataFrame([{"node": source, "demand": 0}]).set_index("node") + # Solve + obj, flow = min_cost_flow_pandas(arc_data, demand_data, **kwargs) + arc_data = _remove_dummy_edge(arc_data, source, sink) + flow = _remove_dummy_edge(flow, source, sink) + return -obj, flow + + +def _max_flow_scipy(G, source, sink, **kwargs): + # Create new matrix with dummy edge (sink, source) + max_flow = G.tolil()[[0], :].sum() + data = np.append(G.data, max_flow) + from_arc = np.append(G.row, sink) + to_arc = np.append(G.col, source) + G = sp.coo_array((data, (from_arc, to_arc)), dtype=float) + + capacities = sp.coo_array(G) + + costs = np.zeros(G.row.shape, dtype=float) + costs[-1] = -1 + costs = sp.coo_array((costs, (G.row, G.col)), dtype=float) + demands = np.zeros(G.shape[1], dtype=float) + # Solve + obj, flow = min_cost_flow_scipy(G, capacities, costs, demands, **kwargs) + G = _remove_dummy_edge(G, source, sink) + flow = _remove_dummy_edge(flow, source, sink) + return -obj, sp.coo_array(flow) + + +def _max_flow_networkx(G, source, sink, **kwargs): + nx.set_edge_attributes(G, 0, "cost") + nx.set_node_attributes(G, 0, "demand") + max_flow = 0 + for j in G.successors(source): + # G.edges[source, j]["cost"] = -1 + max_flow += G.edges[source, j]["capacity"] + G.add_edge(sink, source, cost=-1, capacity=max_flow) + # Set cost attribute to -1 for all outgoing edges from source + obj, flow_graph = min_cost_flow_networkx(G, **kwargs) + G = _remove_dummy_edge(G, source, sink) + if len(flow_graph.edges) > 0: + flow_graph = _remove_dummy_edge(flow_graph, source, sink) + return -obj, flow_graph diff --git a/src/gurobi_optimods/min_cut.py b/src/gurobi_optimods/min_cut.py new file mode 100644 index 000000000..2d8d5e777 --- /dev/null +++ b/src/gurobi_optimods/min_cut.py @@ -0,0 +1,292 @@ +""" +Minimum Cut +----------- +""" + +import logging +from dataclasses import dataclass + +import gurobipy as gp +import gurobipy_pandas as gppd +import numpy as np +import pandas as pd +import scipy.sparse as sp +from gurobipy import GRB + +try: + import networkx as nx +except ImportError: + nx = None + +from gurobi_optimods.max_flow import _remove_dummy_edge +from gurobi_optimods.utils import optimod + +logger = logging.getLogger(__name__) + + +@dataclass +class MinCutResult: + """ + Solution to a Minimum-cut problem. + + Attributes + ---------- + cut: float + Cut value of the minimum cut. + partition: tuple with sets + Partition of size 2 with cut sets. + cutset: set of tuple + Cutset with edges. + """ + + cut_value: float + partition: tuple + cutset: set + + +@optimod() +def min_cut(graph, source, sink, *, create_env): + """Solve the minimum cut problem for a given graph. + + Parameters + ---------- + graph : spmatrix or Graph or DataFrame + A graph, specified either as a scipy.sparse adjacency matrix, networkx + graph, or pandas dataframe + source : int or str + The source (or origin) node for the cutset. + sink : int or str + The sink (or destination) node for the cutset. + + Returns + ------- + min_cut_result: MinCutResult + A dataclass containing the cut value, and set of nodes and edges in the + minimum cut. + """ + if isinstance(graph, sp.spmatrix): + return _min_cut_scipy(graph, source, sink, create_env) + elif isinstance(graph, pd.DataFrame): + return _min_cut_pandas(graph, source, sink, create_env) + elif nx is not None and isinstance(graph, nx.Graph): + return _min_cut_networkx(graph, source, sink, create_env) + else: + raise ValueError(f"Unknown graph type: {type(graph)}") + + +def _min_cut_pandas(arc_data, source, sink, create_env): + f, t = arc_data.index.names + arc_data["cost"] = [0] * len(arc_data) + # Create dummy edge to find maximum flow through (minimum of sum of all + # outgoing/incoming capacities at the source/sink) + max_flow = min( + arc_data.loc[([source], slice(None)),]["capacity"].sum(), + arc_data.loc[(slice(None), [sink]),]["capacity"].sum(), + ) + arc_data = pd.concat( + [ + arc_data, + pd.DataFrame( + [{f: sink, t: source, "capacity": max_flow, "cost": 1}] + ).set_index([f, t]), + ] + ) + + with create_env() as env, gp.Model(env=env) as model: + # Solve max-flow problem + model.ModelSense = GRB.MAXIMIZE + arc_df = arc_data.gppd.add_vars(model, obj="cost", name="flow") + balance_df = pd.DataFrame( + { + "inflow": arc_df["flow"].groupby(t).sum(), + "outflow": arc_df["flow"].groupby(f).sum(), + } + ).gppd.add_constrs(model, "inflow == outflow", name="balance") + capacity_constrs = gppd.add_constrs( + model, + arc_df["flow"], + GRB.LESS_EQUAL, + arc_df["capacity"], + name="capacity", + ) + logger.info( + f"Solving min-cut problem with {len(balance_df)} nodes and " + f"{len(arc_data)-1} edges" + ) + model.optimize() + _remove_dummy_edge(arc_data, source, sink) + + if model.Status in [GRB.INFEASIBLE, GRB.INF_OR_UNBD]: + raise ValueError("Unsatisfiable flows") + + # Construct partition and cutset + cap_pi = capacity_constrs.gppd.Pi + + # Find edges in the cutset (excluding the dummy (sink, source) edge. + cutset = set([a for a in arc_data.index if cap_pi[a] > 1e-3 if a[0] != sink]) + + if len(cutset) == 0: # No arc in the cutset + return MinCutResult(0.0, (set(), set()), set()) + + p1 = set() + queue = [source] + while len(queue) > 0: + node = queue.pop() + p1.add(node) + # Add successors of `node` that are not in the cutset + queue.extend( + [ + a[1] + for a in arc_data.loc[([node], slice(None)),].index + if a not in cutset and a[1] not in p1 and a[1] not in queue + ] + ) + p2 = set([n for n in balance_df.index if n not in p1]) + return MinCutResult(model.ObjVal, (p1, p2), cutset) + + +def _min_cut_scipy(G, source, sink, create_env): + # Create new matrix with dummy edge (sink, source) + original_shape = G.shape + max_flow = G.tolil()[[0], :].sum() + + data = np.append(G.data, max_flow) + + from_arc = np.append(G.row, sink) + to_arc = np.append(G.col, source) + G = sp.coo_array((data, (from_arc, to_arc)), dtype=float) + + capacities = data + + costs = np.zeros(G.row.shape, dtype=float) + costs[-1] = 1 + demands = np.zeros(G.shape[1], dtype=float) + + G = G.tocoo() + + # Create incidence matrix from edge lists. + indices = np.column_stack((from_arc, to_arc)).reshape(-1, order="C") + indptr = np.arange(0, 2 * from_arc.shape[0] + 2, 2) + ones = np.ones(from_arc.shape) + data = np.column_stack((ones * -1.0, ones)).reshape(-1, order="C") + + A = sp.csc_array((data, indices, indptr)) + + logger.info( + f"Solving min-cut problem with {A.shape[0]} nodes and " f"{A.shape[1]-1} edges" + ) + + with create_env() as env, gp.Model(env=env) as model: + # Solve max-flow problem + model.ModelSense = GRB.MAXIMIZE + x = model.addMVar(A.shape[1], lb=0, obj=costs, name="x") + cap = model.addConstr(x <= capacities, name="capacity") + model.addMConstr(A, x, GRB.EQUAL, demands, name="flow") + model.optimize() + _remove_dummy_edge(G, source, sink) + + if model.Status in [GRB.INFEASIBLE, GRB.INF_OR_UNBD]: + raise ValueError("Unsatisfiable flows") + + # Construct partition and cutset + cap_pi = cap.Pi + + p1 = set() + cutset = set( + [ + (i, j) + for n, (i, j) in enumerate(zip(G.row, G.col)) + if cap_pi[n] > 1e-3 and i != sink + ] + ) + if len(cutset) == 0: # No arc in the cutset + return MinCutResult(0.0, (set(), set()), set()) + + queue = [source] + while len(queue) > 0: + node = queue.pop() + p1.add(node) + row = G.getrow(node) + # Add successors of `node` that are not in the cutset + queue.extend( + [ + j + for j in row.tocoo().col + if (node, j) not in cutset and j not in p1 and j not in queue + ] + ) + p2 = set([n for n in range(G.shape[1]) if n not in p1]) + return MinCutResult(model.ObjVal, (p1, p2), cutset) + + +def _min_cut_networkx(G, source, sink, create_env): + logger.info( + f"Solving min-cut problem with {len(G.nodes)} nodes and " + f"{len(G.edges)} edges" + ) + nx.set_edge_attributes(G, 0, "cost") + max_flow = 0 + for j in G.successors(source): + max_flow += G.edges[source, j]["capacity"] + G.add_edge(sink, source, cost=1, capacity=max_flow) + with create_env() as env, gp.Model(env=env) as model: + model.ModelSense = GRB.MAXIMIZE + edges, capacities, costs = gp.multidict( + {(i, j): [d["capacity"], d["cost"]] for i, j, d in G.edges(data=True)} + ) + nodes = list(G.nodes(data=True)) + x = { + (i, j): model.addVar(obj=costs[i, j], name=f"flow[{i},{j}]") + for i, j in edges + } + + 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)) + ), + name=f"flow_balance[{n}]", + ) + for n, data in nodes + } + model.update() + capacity_constrs = { + (i, j): model.addConstr( + x[i, j] <= capacities[i, j], name=f"capacity[{i},{j}]" + ) + for i, j in edges + } + + model.optimize() + _remove_dummy_edge(G, source, sink) + + if model.Status in [GRB.INFEASIBLE, GRB.INF_OR_UNBD]: + raise ValueError("Unsatisfiable flows") + + # Construct partition and cutset + cutset = set( + [ + (i, j) + for (i, j) in edges + if capacity_constrs[i, j].Pi > 1e-3 and i != sink + ] + ) + if len(cutset) == 0: + return MinCutResult(0.0, (set(), set()), set()) + + p1 = set() + queue = [source] + while len(queue) > 0: + node = queue.pop() + p1.add(node) + # Add successors of `node` that are not in the cutset + queue.extend( + [ + j + for j in G.successors(node) + if (node, j) not in cutset and j not in p1 and j not in queue + ] + ) + p2 = set([n for n in G.nodes if n not in p1]) + return MinCutResult(model.ObjVal, (p1, p2), cutset) diff --git a/tests/test_max_flow.py b/tests/test_max_flow.py new file mode 100644 index 000000000..536343c2f --- /dev/null +++ b/tests/test_max_flow.py @@ -0,0 +1,174 @@ +import unittest + +import numpy as np + +try: + import networkx as nx +except ImportError: + nx = None + +import gurobi_optimods.datasets as datasets +from gurobi_optimods.max_flow import max_flow + +from .test_graph_utils import ( + check_solution_networkx, + check_solution_pandas, + check_solution_scipy, +) +from .test_min_cost_flow import ( + load_graph2_networkx, + load_graph2_pandas, + load_graph2_scipy, +) + + +class TestMaxFlow(unittest.TestCase): + def setUp(self): + self.expected_max_flow = 3.0 + + def test_pandas(self): + edge_data, _ = datasets.simple_graph_pandas() + 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, + } + self.assertTrue(check_solution_pandas(sol, [candidate, candidate2])) + + def test_empty_pandas(self): + edge_data, _ = datasets.simple_graph_pandas() + edge_data["capacity"] = [0] * len(edge_data) + obj, sol = max_flow(edge_data, 0, 5) + self.assertEqual(obj, 0.0) + + def test_scipy(self): + G, capacity, _, _ = datasets.simple_graph_scipy() + G.data = capacity.data + obj, sol = max_flow(G, 0, 5) + self.assertEqual(obj, self.expected_max_flow) + expected = np.array( + [ + [0, 1, 2, 0, 0, 0], + [0, 0, 0, 1, 0, 0], + [0, 0, 0, 0, 2, 0], + [0, 0, 0, 0, 0, 1], + [0, 0, 0, 0, 0, 2], + ] + ) + self.assertTrue(check_solution_scipy(sol, [expected])) + + def test_empty_scipy(self): + G, capacity, _, _ = datasets.simple_graph_scipy() + # Use vsmall values to set the capacity to nearly zero + G.data = np.repeat(1e-20, len(G.data)) + obj, sol = max_flow(G, 0, 5) + self.assertEqual(obj, 0.0) + + @unittest.skipIf(nx is None, "networkx is not installed") + def test_networkx(self): + 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}, + } + self.assertTrue(check_solution_networkx(sol, [candidate, candidate2])) + + @unittest.skipIf(nx is None, "networkx is not installed") + def test_empty_networkx(self): + G = datasets.simple_graph_networkx() + nx.set_edge_attributes(G, 0, "capacity") + obj, sol = max_flow(G, 0, 5) + self.assertEqual(obj, 0.0) + + +class TestMaxFlow2(unittest.TestCase): + def setUp(self): + self.expected_max_flow = 23.0 + + def test_pandas(self): + edge_data, _ = load_graph2_pandas() + 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, + } + self.assertTrue(check_solution_pandas(sol, [candidate])) + + def test_scipy(self): + G, capacity, _, _ = load_graph2_scipy() + G.data = capacity.data + obj, sol = max_flow(G, 0, 4) + self.assertEqual(obj, self.expected_max_flow) + expected = np.array( + [ + [0, 15, 8, 0, 0], + [0, 0, 1, 4, 10], + [0, 0, 0, 9, 0], + [0, 0, 0, 0, 13], + ] + ) + self.assertTrue(check_solution_scipy(sol, [expected])) + + @unittest.skipIf(nx is None, "networkx is not installed") + def test_networkx(self): + 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}, + } + self.assertTrue(check_solution_networkx(sol, [candidate, candidate2])) diff --git a/tests/test_min_cut.py b/tests/test_min_cut.py new file mode 100644 index 000000000..35692c523 --- /dev/null +++ b/tests/test_min_cut.py @@ -0,0 +1,173 @@ +import unittest + +import numpy as np +import pandas as pd + +try: + import networkx as nx +except ImportError: + nx = None + +import gurobi_optimods.datasets as datasets +from gurobi_optimods.datasets import ( + _convert_pandas_to_digraph, + _convert_pandas_to_scipy, +) +from gurobi_optimods.min_cut import min_cut + +from .test_min_cost_flow import ( + load_graph2_networkx, + load_graph2_pandas, + load_graph2_scipy, +) + + +class TestMinCut(unittest.TestCase): + def setUp(self): + self.expected_partition = (set({0, 1}), set({2, 3, 4, 5})) + self.expected_cut_value = 3.0 + self.expected_cut_set = set([(0, 2), (1, 3)]) + + def test_pandas(self): + edge_data, node_data = datasets.simple_graph_pandas() + res = min_cut(edge_data, 0, 5) + cut_value, partition, cutset = res.cut_value, res.partition, res.cutset + self.assertEqual(cut_value, self.expected_cut_value) + self.assertEqual(partition[0], self.expected_partition[0]) + self.assertEqual(partition[1], self.expected_partition[1]) + self.assertEqual(cutset, self.expected_cut_set) + + def test_empty_pandas(self): + edge_data, _ = datasets.simple_graph_pandas() + edge_data["capacity"] = [0] * len(edge_data) + res = min_cut(edge_data, 0, 5) + cut_value, partition, cutset = res.cut_value, res.partition, res.cutset + self.assertEqual(cut_value, 0.0) + self.assertEqual(partition, (set(), set())) + self.assertEqual(cutset, set()) + + @unittest.skipIf(nx is None, "networkx is not installed") + def test_networkx(self): + G = datasets.simple_graph_networkx() + res = min_cut(G, 0, 5) + cut_value, partition, cutset = res.cut_value, res.partition, res.cutset + self.assertEqual(cut_value, self.expected_cut_value) + self.assertEqual(partition[0], self.expected_partition[0]) + self.assertEqual(partition[1], self.expected_partition[1]) + self.assertEqual(cutset, self.expected_cut_set) + + @unittest.skipIf(nx is None, "networkx is not installed") + def test_empty_networkx(self): + G = datasets.simple_graph_networkx() + nx.set_edge_attributes(G, 0, "capacity") + res = min_cut(G, 0, 5) + cut_value, partition, cutset = res.cut_value, res.partition, res.cutset + self.assertEqual(cut_value, 0.0) + self.assertEqual(partition, (set(), set())) + self.assertEqual(cutset, set()) + + def test_scipy(self): + G, capacity, _, _ = datasets.simple_graph_scipy() + G.data = capacity.data + res = min_cut(G, 0, 5) + cut_value, partition, cutset = res.cut_value, res.partition, res.cutset + self.assertEqual(cut_value, self.expected_cut_value) + self.assertEqual(partition[0], self.expected_partition[0]) + self.assertEqual(partition[1], self.expected_partition[1]) + self.assertEqual(cutset, self.expected_cut_set) + + def test_empty_scipy(self): + G, capacity, _, _ = datasets.simple_graph_scipy() + G.data = np.repeat(1e-20, len(G.data)) + res = min_cut(G, 0, 5) + cut_value, partition, cutset = res.cut_value, res.partition, res.cutset + self.assertEqual(cut_value, 0.0) + self.assertEqual(partition, (set(), set())) + self.assertEqual(cutset, set()) + + +class TestMinCut2(unittest.TestCase): + def setUp(self): + self.expected_partition = (set({0}), set({1, 2, 3, 4})) + self.expected_cut_value = 23.0 + self.expected_cut_set = set([(0, 1), (0, 2)]) + + def test_pandas(self): + edge_data, _ = load_graph2_pandas() + res = min_cut(edge_data, 0, 4) + cut_value, partition, cutset = res.cut_value, res.partition, res.cutset + self.assertEqual(cut_value, self.expected_cut_value) + self.assertEqual(partition[0], self.expected_partition[0]) + self.assertEqual(partition[1], self.expected_partition[1]) + self.assertEqual(cutset, self.expected_cut_set) + + @unittest.skipIf(nx is None, "networkx is not installed") + def test_networkx(self): + G = load_graph2_networkx() + res = min_cut(G, 0, 4) + cut_value, partition, cutset = res.cut_value, res.partition, res.cutset + self.assertEqual(cut_value, self.expected_cut_value) + self.assertEqual(partition[0], self.expected_partition[0]) + self.assertEqual(partition[1], self.expected_partition[1]) + self.assertEqual(cutset, self.expected_cut_set) + + def test_scipy(self): + G, capacity, _, _ = load_graph2_scipy() + G.data = capacity.data + res = min_cut(G, 0, 4) + cut_value, partition, cutset = res.cut_value, res.partition, res.cutset + self.assertEqual(cut_value, self.expected_cut_value) + self.assertEqual(partition[0], self.expected_partition[0]) + self.assertEqual(partition[1], self.expected_partition[1]) + self.assertEqual(cutset, self.expected_cut_set) + + +class TestMinCut3(unittest.TestCase): + def setUp(self): + self.expected_partition = (set({3, 1, 0}), set({4, 6, 5, 2})) + self.expected_cut_value = 3.0 + self.expected_cut_set = set([(0, 2), (3, 6)]) + self.arc_data = pd.DataFrame( + [ + {"source": 0, "target": 1, "capacity": 5.0}, + {"source": 0, "target": 2, "capacity": 1.0}, + {"source": 1, "target": 3, "capacity": 3.0}, + {"source": 2, "target": 3, "capacity": 5.0}, + {"source": 2, "target": 4, "capacity": 4.0}, + {"source": 4, "target": 5, "capacity": 2.0}, + {"source": 3, "target": 6, "capacity": 2.0}, + {"source": 5, "target": 6, "capacity": 3.0}, + ] + ).set_index(["source", "target"]) + + def test_pandas(self): + res = min_cut(self.arc_data, 0, 6) + cut_value, partition, cutset = res.cut_value, res.partition, res.cutset + self.assertEqual(cut_value, self.expected_cut_value) + self.assertEqual(partition[0], self.expected_partition[0]) + self.assertEqual(partition[1], self.expected_partition[1]) + self.assertEqual(cutset, self.expected_cut_set) + + @unittest.skipIf(nx is None, "networkx is not installed") + def test_networkx(self): + res = min_cut( + _convert_pandas_to_digraph(self.arc_data, None, demand=False), 0, 6 + ) + cut_value, partition, cutset = res.cut_value, res.partition, res.cutset + self.assertEqual(cut_value, self.expected_cut_value) + self.assertEqual(partition[0], self.expected_partition[0]) + self.assertEqual(partition[1], self.expected_partition[1]) + self.assertEqual(cutset, self.expected_cut_set) + + def test_scipy(self): + G, capacity, _, _ = _convert_pandas_to_scipy( + self.arc_data, None, cost=False, demand=False + ) + + G.data = capacity.data + res = min_cut(G, 0, 6) + cut_value, partition, cutset = res.cut_value, res.partition, res.cutset + self.assertEqual(cut_value, self.expected_cut_value) + self.assertEqual(partition[0], self.expected_partition[0]) + self.assertEqual(partition[1], self.expected_partition[1]) + self.assertEqual(cutset, self.expected_cut_set)