summaryrefslogtreecommitdiff
path: root/networkx/algorithms
diff options
context:
space:
mode:
authorMatt Schwennesen <mjschwenne@gmail.com>2021-08-23 08:50:12 -0400
committerGitHub <noreply@github.com>2021-08-23 08:50:12 -0400
commite098dccd8be4e08497117fa39d66800907c2a932 (patch)
treef00ce92d55bcbe545436bf3adea44d52aa2cc6bd /networkx/algorithms
parent5867cbae81ca5e95e1b624a197c8cd6ddcb3fef8 (diff)
downloadnetworkx-e098dccd8be4e08497117fa39d66800907c2a932.tar.gz
GSoC Asadpour ATSP Implementation Pull Request (#4740)
* Add greedy algorithm for solving TSP Many problems of Combinational Optimization can be represented as graphs. These problems have enormous significance in many aspects of science, but there are not any algorithms to solve some of them in polynomial time. However many, heuristic and metaheuristic algorithms have been published over the past years in order to solve / approximate the solutions to these problems. The purpose of this commit is to add implementation of such algorithms for solve one of the most famous problems of Combinational Optimizations, Travelling Salesman Problem (TSP). A greedy algorithm has been implemented at the moment for this reason. "applications" package has been created which include modules that represent a problem. Each module contains several algorithms for solving the specific problem. At this commit, tsp.py module is added which contains greedy_tsp() function; a implementation of a greedy algorithm. * Fix example error * Trivial changes List of changes: Removal of unnesecary _is_weighted() function Improvements on documentation * Add applications package to setup.py file * Change output of greedy algorithm Algorithm's output is a list of nodes now * Add simulated annealing algorithm Add a metaheuristic local search algorithm for solving TSP * Minor changes * Fix example doc errors * Compatible with python 3 * Move tsp module to algorithms package * Code improvements * Handle small graphs and fix doc examples * Documentation changes and rename variables * Adds Threshold Accepting algorithm for TSP * Implemented maximal matching of minimal weight and created test suite. * Removed useless print * Implemented Christofides. * Coding was missing * Add more general traveling_salesman_problem using christofides Also reconfigure import structure and remove min_weight_matching from module since it is now in matching.py * Add new functions to the docs and minor typos * pep8 fixes * fix pep8 and change .gitignore * Add tests of the approximation namespace update docs in approximation/__init__.py * Fix is_matching to check if edges in G. Other tweaks: doc changes and put not_implemented_for on find_matching functions * Improve is_matching selfloop handling and expand tests * Move tsp to approximation directory. Apply black. * Move tsp tests to approximation tests folder * Attempt to bring tsp up to current code. * commit pep8 that my black didnt change, but pep8speaks did find. ?? * tweak a few things and run black * combine #4083 and #3585 into traveling_salesman.py * Match chistofides output to other tsp functions and adjust calling syntax of tests tweak docs tweak see also section * Put big-O complexity in in-line math env. Prevents sphinx from trying to do variable substitution between pipes. * Minor touchups to christofides docstring. * RST touchups to tsp module docstring. * Rm extra string from tsp module. * Docstring touchups for traveling_salesman_problem. * rst fixups for greedy_tsp docstring. * rst formatting for simulated annealing docstring. * More math in-lining for simulated annealing docstring. * rst and minor grammatical fixes to TA docstring. * Fix path-finding and test all methods for tsp function * the refactoring was incomplete. Now maybe is - Add tests of TSP with all methods. - Refactor tests to match simulated_annealing tests and threshold tests. - Unify treatment of weight so unweighted edges use default weight 1. weight now defaults to "weight" with a default value of 1. - Rename tolerance to max_iterations (tolerance is used for error bound) - Rename iterations to N_inner (each iteration takes this many inner loops) - Introduce idioms like `pairwise` and `cycle.copy()` (over cycle[:]) - Allow passthrough of method kwargs for traveling_salesman_problem Still need to: - add test of case where path is more than one edge less that cycle (incomplete_graph) - require cycle input (maybe make default list(G)??) - consider the complexity claims in the doc_strings * More api changes to TSP functions - `chritofides` now allows (and ignores) selfloops - `move` can be a function as well as "1-1" and "1-0" - `method` for traveling_salesman_problem must have 2 arguments instead of passing kwargs. User must "curry" to set parameters - changed doc_string typos in matching.py * Add test to check that cycle=False can remove many edges * Change init_cycle api to require input from user The idea is to make the user specify the initial cycle to start from rather than relying on the programmers default of a greedy algorithm. To easy usage, I check for a string "greedy" as a shortcut. * Update docs with more correct complexity info. * Check for complete graph now more efficient and selfloops ignored * merge is_matching changes * Stub for Asadpour. Needed to create GSoC PR * Update to integrate changes from main * Added function stubs and draft docstrings for the Asadpour algorithm * testing * I'm not entirly sure how the commit hook works... * Moved iterators into the correct files to maintain proper codebase visibility * Including Black reformat * Grabbing black reformats * Working on debugging ascent method plus black reformats * Ascent method terminating, but at non-optimal solution * minor edits * Fixed termination condition, still given non-optimal result * Minor bugfix, still non-optimal result * Fixed subtle bug in find_epsilon() * Cleaned code and tried something which didn't work * Modified the ArborescenceIterator to accept init partition * Black formats * Branch and bound returning optimal solution * Working Ascent method, code needs cleaning * black formatting changes * Performance tweaks and testing fractional answers * Fixed test bug, I hope * Asadpour output for ascent method * Fixed numpy imports crashing pypi tests * Removed branch and bound method. One unit test misbehaving * Added asymmetric fractional test for the ascent method * Removed printn statements and tweaked final test to be more asymmetric * Draft of spanning_tree_distribution * Black changes * Changed HK to only report on the support of the answer * Fixed contraction bug by changing to MultiGraph. Problem with prob > 1 * Black reformats * Fixed pypi test error * Further testing of dist fix * Can sample spanning trees * Developing test for sampling spanning tree * Changed sample_spanning_tree test to Chi squared test * Tweaked signifiance level * Found true minimum sample size * fixed typo * untested implementation of asadpour_tsp * Fixed issue reading flow_dict * Fixed runtime errors in asadpour_tsp * black reformats * Adding test cases * documentation update * Fixed rounding error with tests * One new test and check * Documentation update for the iterators * Attempting to fix class documentation * Reventing documentation changes * Update mst.py to accept suggestion Co-authored-by: Dan Schult <dschult@colgate.edu> * Update branchings.py accept doc string edit Co-authored-by: Dan Schult <dschult@colgate.edu> * Review suggestions from dshult * Cleaned code, merged functions if possible and opened partition functionality to all * Fixed pypi test error * Implemented review suggestions from rossbar * review edits added SpanningTreeIterator to algorithms/__init__.py * Update __init__.py ack, hasty / stupid change was meant to be a draft; github isn't letting me make a new branch to PR into this one * fixed misspelling of Kirchhoff * Implement suggestions from boothby Co-authored-by: Thodoris Sotiropoulos <theosotr@windowslive.com> Co-authored-by: Luca Cappelletti <cappelletti.luca94@gmail.com> Co-authored-by: Dan Schult <dschult@colgate.edu> Co-authored-by: Ross Barnowski <rossbar@berkeley.edu> Co-authored-by: Kelly Boothby <kelly.r.boothby@gmail.com> Co-authored-by: Kelly Boothby <boothby@dwavesys.com>
Diffstat (limited to 'networkx/algorithms')
-rw-r--r--networkx/algorithms/__init__.py1
-rw-r--r--networkx/algorithms/approximation/tests/test_traveling_salesman.py722
-rw-r--r--networkx/algorithms/approximation/traveling_salesman.py797
-rw-r--r--networkx/algorithms/minors/contraction.py8
-rw-r--r--networkx/algorithms/tree/branchings.py339
-rw-r--r--networkx/algorithms/tree/mst.py333
-rw-r--r--networkx/algorithms/tree/tests/test_branchings.py139
-rw-r--r--networkx/algorithms/tree/tests/test_mst.py102
8 files changed, 2327 insertions, 114 deletions
diff --git a/networkx/algorithms/__init__.py b/networkx/algorithms/__init__.py
index 93195475..2dc01011 100644
--- a/networkx/algorithms/__init__.py
+++ b/networkx/algorithms/__init__.py
@@ -120,6 +120,7 @@ from networkx.algorithms.tree.branchings import maximum_branching
from networkx.algorithms.tree.branchings import maximum_spanning_arborescence
from networkx.algorithms.tree.branchings import minimum_branching
from networkx.algorithms.tree.branchings import minimum_spanning_arborescence
+from networkx.algorithms.tree.branchings import ArborescenceIterator
from networkx.algorithms.tree.coding import *
from networkx.algorithms.tree.decomposition import *
from networkx.algorithms.tree.mst import *
diff --git a/networkx/algorithms/approximation/tests/test_traveling_salesman.py b/networkx/algorithms/approximation/tests/test_traveling_salesman.py
index f36981d4..db86d7b3 100644
--- a/networkx/algorithms/approximation/tests/test_traveling_salesman.py
+++ b/networkx/algorithms/approximation/tests/test_traveling_salesman.py
@@ -216,20 +216,8 @@ class TestSimulatedAnnealingTSP:
pytest.raises(nx.NetworkXError, self.tsp, self.UG, "weight")
def test_not_complete_graph(self):
- pytest.raises(
- nx.NetworkXError,
- self.tsp,
- self.incompleteUG,
- "greedy",
- source=0,
- )
- pytest.raises(
- nx.NetworkXError,
- self.tsp,
- self.incompleteDG,
- "greedy",
- source=0,
- )
+ pytest.raises(nx.NetworkXError, self.tsp, self.incompleteUG, "greedy", source=0)
+ pytest.raises(nx.NetworkXError, self.tsp, self.incompleteDG, "greedy", source=0)
def test_ignore_selfloops(self):
G = nx.complete_graph(5)
@@ -301,12 +289,7 @@ class TestThresholdAcceptingTSP(TestSimulatedAnnealingTSP):
# set threshold too low
initial_sol = ["D", "A", "B", "C", "D"]
cycle = self.tsp(
- self.DG,
- initial_sol,
- source="D",
- move="1-0",
- threshold=-3,
- seed=42,
+ self.DG, initial_sol, source="D", move="1-0", threshold=-3, seed=42
)
cost = sum(self.DG[n][nbr]["weight"] for n, nbr in pairwise(cycle))
assert cost > self.DG_cost
@@ -355,10 +338,7 @@ def test_TSP_weighted():
[6, 7, 8, 0, 1, 2, 3, 2, 1, 0, 8, 7, 6],
)
# path through all nodes
- expected_tourpaths = (
- [5, 6, 7, 8, 0, 1, 2, 3, 4],
- [4, 3, 2, 1, 0, 8, 7, 6, 5],
- )
+ expected_tourpaths = ([5, 6, 7, 8, 0, 1, 2, 3, 4], [4, 3, 2, 1, 0, 8, 7, 6, 5])
# Check default method
cycle = tsp(G, nodes=[3, 6], weight="weight")
@@ -402,3 +382,697 @@ def test_TSP_incomplete_graph_short_path():
path = nx_app.traveling_salesman_problem(G, cycle=False)
print(path)
assert len(path) == 13 and len(set(path)) == 12
+
+
+def test_held_karp_ascent():
+ """
+ Test the Held-Karp relaxation with the ascent method
+ """
+ import networkx.algorithms.approximation.traveling_salesman as tsp
+
+ np = pytest.importorskip("numpy")
+
+ # Adjacency matrix from page 1153 of the 1970 Held and Karp paper
+ # which have been edited to be directional, but also symmetric
+ G_array = np.array(
+ [
+ [0, 97, 60, 73, 17, 52],
+ [97, 0, 41, 52, 90, 30],
+ [60, 41, 0, 21, 35, 41],
+ [73, 52, 21, 0, 95, 46],
+ [17, 90, 35, 95, 0, 81],
+ [52, 30, 41, 46, 81, 0],
+ ]
+ )
+
+ solution_edges = [(1, 3), (2, 4), (3, 2), (4, 0), (5, 1), (0, 5)]
+
+ G = nx.from_numpy_array(G_array, create_using=nx.DiGraph)
+ opt_hk, z_star = tsp.held_karp_ascent(G)
+
+ # Check that the optimal weights are the same
+ assert round(opt_hk, 2) == 207.00
+ # Check that the z_stars are the same
+ solution = nx.DiGraph()
+ solution.add_edges_from(solution_edges)
+ assert nx.utils.edges_equal(z_star.edges, solution.edges)
+
+
+def test_ascent_fractional_solution():
+ """
+ Test the ascent method using a modified version of Figure 2 on page 1140
+ in 'The Traveling Salesman Problem and Minimum Spanning Trees' by Held and
+ Karp
+ """
+ import networkx.algorithms.approximation.traveling_salesman as tsp
+
+ np = pytest.importorskip("numpy")
+
+ # This version of Figure 2 has all of the edge weights multiplied by 100
+ # and is a complete directed graph with infinite edge weights for the
+ # edges not listed in the original graph
+ G_array = np.array(
+ [
+ [0, 100, 100, 100000, 100000, 1],
+ [100, 0, 100, 100000, 1, 100000],
+ [100, 100, 0, 1, 100000, 100000],
+ [100000, 100000, 1, 0, 100, 100],
+ [100000, 1, 100000, 100, 0, 100],
+ [1, 100000, 100000, 100, 100, 0],
+ ]
+ )
+
+ solution_z_star = {
+ (0, 1): 5 / 12,
+ (0, 2): 5 / 12,
+ (0, 5): 5 / 6,
+ (1, 0): 5 / 12,
+ (1, 2): 1 / 3,
+ (1, 4): 5 / 6,
+ (2, 0): 5 / 12,
+ (2, 1): 1 / 3,
+ (2, 3): 5 / 6,
+ (3, 2): 5 / 6,
+ (3, 4): 1 / 3,
+ (3, 5): 1 / 2,
+ (4, 1): 5 / 6,
+ (4, 3): 1 / 3,
+ (4, 5): 1 / 2,
+ (5, 0): 5 / 6,
+ (5, 3): 1 / 2,
+ (5, 4): 1 / 2,
+ }
+
+ G = nx.from_numpy_array(G_array, create_using=nx.DiGraph)
+ opt_hk, z_star = tsp.held_karp_ascent(G)
+
+ # Check that the optimal weights are the same
+ assert round(opt_hk, 2) == 303.00
+ # Check that the z_stars are the same
+ assert {key: round(z_star[key], 4) for key in z_star} == {
+ key: round(solution_z_star[key], 4) for key in solution_z_star
+ }
+
+
+def test_ascent_method_asymmetric():
+ """
+ Tests the ascent method using a truly asymmetric graph for which the
+ solution has been brute forced
+ """
+ import networkx.algorithms.approximation.traveling_salesman as tsp
+
+ np = pytest.importorskip("numpy")
+
+ G_array = np.array(
+ [
+ [0, 26, 63, 59, 69, 31, 41],
+ [62, 0, 91, 53, 75, 87, 47],
+ [47, 82, 0, 90, 15, 9, 18],
+ [68, 19, 5, 0, 58, 34, 93],
+ [11, 58, 53, 55, 0, 61, 79],
+ [88, 75, 13, 76, 98, 0, 40],
+ [41, 61, 55, 88, 46, 45, 0],
+ ]
+ )
+
+ solution_edges = [(0, 1), (1, 3), (3, 2), (2, 5), (5, 6), (4, 0), (6, 4)]
+
+ G = nx.from_numpy_array(G_array, create_using=nx.DiGraph)
+ opt_hk, z_star = tsp.held_karp_ascent(G)
+
+ # Check that the optimal weights are the same
+ assert round(opt_hk, 2) == 190.00
+ # Check that the z_stars match.
+ solution = nx.DiGraph()
+ solution.add_edges_from(solution_edges)
+ assert nx.utils.edges_equal(z_star.edges, solution.edges)
+
+
+def test_ascent_method_asymmetric_2():
+ """
+ Tests the ascent method using a truly asymmetric graph for which the
+ solution has been brute forced
+ """
+ import networkx.algorithms.approximation.traveling_salesman as tsp
+
+ np = pytest.importorskip("numpy")
+
+ G_array = np.array(
+ [
+ [0, 45, 39, 92, 29, 31],
+ [72, 0, 4, 12, 21, 60],
+ [81, 6, 0, 98, 70, 53],
+ [49, 71, 59, 0, 98, 94],
+ [74, 95, 24, 43, 0, 47],
+ [56, 43, 3, 65, 22, 0],
+ ]
+ )
+
+ solution_edges = [(0, 5), (5, 4), (1, 3), (3, 0), (2, 1), (4, 2)]
+
+ G = nx.from_numpy_array(G_array, create_using=nx.DiGraph)
+ opt_hk, z_star = tsp.held_karp_ascent(G)
+
+ # Check that the optimal weights are the same
+ assert round(opt_hk, 2) == 144.00
+ # Check that the z_stars match.
+ solution = nx.DiGraph()
+ solution.add_edges_from(solution_edges)
+ assert nx.utils.edges_equal(z_star.edges, solution.edges)
+
+
+def test_held_karp_ascent_asymmetric_3():
+ """
+ Tests the ascent method using a truly asymmetric graph with a fractional
+ solution for which the solution has been brute forced.
+
+ In this graph their are two different optimal, integral solutions (which
+ are also the overall atsp solutions) to the Held Karp relaxation. However,
+ this particular graph has two different tours of optimal value and the
+ possible solutions in the held_karp_ascent function are not stored in an
+ ordered data structure.
+ """
+ import networkx.algorithms.approximation.traveling_salesman as tsp
+
+ np = pytest.importorskip("numpy")
+
+ G_array = np.array(
+ [
+ [0, 1, 5, 2, 7, 4],
+ [7, 0, 7, 7, 1, 4],
+ [4, 7, 0, 9, 2, 1],
+ [7, 2, 7, 0, 4, 4],
+ [5, 5, 4, 4, 0, 3],
+ [3, 9, 1, 3, 4, 0],
+ ]
+ )
+
+ solution1_edges = [(0, 3), (1, 4), (2, 5), (3, 1), (4, 2), (5, 0)]
+
+ solution2_edges = [(0, 3), (3, 1), (1, 4), (4, 5), (2, 0), (5, 2)]
+
+ G = nx.from_numpy_array(G_array, create_using=nx.DiGraph)
+ opt_hk, z_star = tsp.held_karp_ascent(G)
+
+ assert round(opt_hk, 2) == 13.00
+ # Check that the z_stars are the same
+ solution1 = nx.DiGraph()
+ solution1.add_edges_from(solution1_edges)
+ solution2 = nx.DiGraph()
+ solution2.add_edges_from(solution2_edges)
+ assert nx.utils.edges_equal(z_star.edges, solution1.edges) or nx.utils.edges_equal(
+ z_star.edges, solution2.edges
+ )
+
+
+def test_held_karp_ascent_fractional_asymmetric():
+ """
+ Tests the ascent method using a truly asymmetric graph with a fractional
+ solution for which the solution has been brute forced
+ """
+ import networkx.algorithms.approximation.traveling_salesman as tsp
+
+ np = pytest.importorskip("numpy")
+
+ G_array = np.array(
+ [
+ [0, 100, 150, 100000, 100000, 1],
+ [150, 0, 100, 100000, 1, 100000],
+ [100, 150, 0, 1, 100000, 100000],
+ [100000, 100000, 1, 0, 150, 100],
+ [100000, 2, 100000, 100, 0, 150],
+ [2, 100000, 100000, 150, 100, 0],
+ ]
+ )
+
+ solution_z_star = {
+ (0, 1): 5 / 12,
+ (0, 2): 5 / 12,
+ (0, 5): 5 / 6,
+ (1, 0): 5 / 12,
+ (1, 2): 5 / 12,
+ (1, 4): 5 / 6,
+ (2, 0): 5 / 12,
+ (2, 1): 5 / 12,
+ (2, 3): 5 / 6,
+ (3, 2): 5 / 6,
+ (3, 4): 5 / 12,
+ (3, 5): 5 / 12,
+ (4, 1): 5 / 6,
+ (4, 3): 5 / 12,
+ (4, 5): 5 / 12,
+ (5, 0): 5 / 6,
+ (5, 3): 5 / 12,
+ (5, 4): 5 / 12,
+ }
+
+ G = nx.from_numpy_array(G_array, create_using=nx.DiGraph)
+ opt_hk, z_star = tsp.held_karp_ascent(G)
+
+ # Check that the optimal weights are the same
+ assert round(opt_hk, 2) == 304.00
+ # Check that the z_stars are the same
+ assert {key: round(z_star[key], 4) for key in z_star} == {
+ key: round(solution_z_star[key], 4) for key in solution_z_star
+ }
+
+
+def test_spanning_tree_distribution():
+ """
+ Test that we can create an exponential distribution of spanning trees such
+ that the probability of each tree is proportional to the product of edge
+ weights.
+
+ Results of this test have been confirmed with hypothesis testing from the
+ created distribution.
+
+ This test uses the symmetric, fractional Held Karp solution.
+ """
+ import networkx.algorithms.approximation.traveling_salesman as tsp
+
+ pytest.importorskip("numpy")
+
+ z_star = {
+ (0, 1): 5 / 12,
+ (0, 2): 5 / 12,
+ (0, 5): 5 / 6,
+ (1, 0): 5 / 12,
+ (1, 2): 1 / 3,
+ (1, 4): 5 / 6,
+ (2, 0): 5 / 12,
+ (2, 1): 1 / 3,
+ (2, 3): 5 / 6,
+ (3, 2): 5 / 6,
+ (3, 4): 1 / 3,
+ (3, 5): 1 / 2,
+ (4, 1): 5 / 6,
+ (4, 3): 1 / 3,
+ (4, 5): 1 / 2,
+ (5, 0): 5 / 6,
+ (5, 3): 1 / 2,
+ (5, 4): 1 / 2,
+ }
+
+ solution_gamma = {
+ (0, 1): -0.6383,
+ (0, 2): -0.6827,
+ (0, 5): 0,
+ (1, 2): -1.0781,
+ (1, 4): 0,
+ (2, 3): 0,
+ (5, 3): -0.2820,
+ (5, 4): -0.3327,
+ (4, 3): -0.9927,
+ }
+
+ # The undirected support of z_star
+ G = nx.MultiGraph()
+ for u, v in z_star:
+ if (u, v) in G.edges or (v, u) in G.edges:
+ continue
+ G.add_edge(u, v)
+
+ gamma = tsp.spanning_tree_distribution(G, z_star)
+
+ assert {key: round(gamma[key], 4) for key in gamma} == solution_gamma
+
+
+def test_sample_spanning_tree():
+ """
+ Using a fixed seed, sample one tree for repeatability.
+ """
+ import networkx.algorithms.approximation.traveling_salesman as tsp
+ from math import exp
+
+ pytest.importorskip("numpy")
+
+ gamma = {
+ (0, 1): -0.6383,
+ (0, 2): -0.6827,
+ (0, 5): 0,
+ (1, 2): -1.0781,
+ (1, 4): 0,
+ (2, 3): 0,
+ (5, 3): -0.2820,
+ (5, 4): -0.3327,
+ (4, 3): -0.9927,
+ }
+
+ # The undirected support of gamma
+ G = nx.Graph()
+ for u, v in gamma:
+ G.add_edge(u, v, lambda_key=exp(gamma[(u, v)]))
+
+ solution_edges = [(2, 3), (3, 4), (0, 5), (5, 4), (4, 1)]
+ solution = nx.Graph()
+ solution.add_edges_from(solution_edges)
+
+ sampled_tree = tsp.sample_spanning_tree(G, "lambda_key", 42)
+
+ assert nx.utils.edges_equal(solution.edges, sampled_tree.edges)
+
+
+def test_sample_spanning_tree_large_sample():
+ """
+ Sample a single spanning tree from the distribution created in the last test
+ """
+ import networkx.algorithms.approximation.traveling_salesman as tsp
+ from math import exp
+ from random import Random
+
+ pytest.importorskip("numpy")
+ stats = pytest.importorskip("scipy.stats")
+
+ gamma = {
+ (0, 1): -0.6383,
+ (0, 2): -0.6827,
+ (0, 5): 0,
+ (1, 2): -1.0781,
+ (1, 4): 0,
+ (2, 3): 0,
+ (5, 3): -0.2820,
+ (5, 4): -0.3327,
+ (4, 3): -0.9927,
+ }
+
+ # The undirected support of gamma
+ G = nx.Graph()
+ for u, v in gamma:
+ G.add_edge(u, v, lambda_key=exp(gamma[(u, v)]))
+
+ # Find the multiplicative weight for each tree.
+ total_weight = 0
+ tree_expected = {}
+ for t in nx.SpanningTreeIterator(G):
+ # Find the multiplicative weight of the spanning tree
+ weight = 1
+ for u, v, d in t.edges(data="lambda_key"):
+ weight *= d
+ tree_expected[t] = weight
+ total_weight += weight
+
+ # Assert that every tree has an entry in the expected distribution
+ assert len(tree_expected) == 75
+
+ # Set the sample size and then calculate the expected number of times we
+ # expect to see each tree. This test uses a near minimum sample size where
+ # the most unlikely tree has an expected frequency of 5.15.
+ # (Minimum required is 5)
+ #
+ # Here we also initialize the tree_actual dict so that we know the keys
+ # match between the two. We will later take advantage of the fact that since
+ # python 3.7 dict order is guaranteed so the expected and actual data will
+ # have the same order.
+ sample_size = 1200
+ tree_actual = {}
+ for t in tree_expected:
+ tree_expected[t] = (tree_expected[t] / total_weight) * sample_size
+ tree_actual[t] = 0
+
+ # Sample the spanning trees
+ #
+ # Assert that they are actually trees and record which of the 75 trees we
+ # have sampled.
+ #
+ # For repeatability, we want to take advantage of the decorators in NetworkX
+ # to randomly sample the same sample each time. However, if we pass in a
+ # constant seed to sample_spanning_tree we will get the same tree each time.
+ # Instead, we can create our own random number generator with a fixed seed
+ # and pass those into sample_spanning_tree.
+ rng = Random(37)
+ for _ in range(sample_size):
+ sampled_tree = tsp.sample_spanning_tree(G, "lambda_key", rng)
+ assert nx.is_tree(sampled_tree)
+
+ for t in tree_expected:
+ if nx.utils.edges_equal(t.edges, sampled_tree.edges):
+ tree_actual[t] += 1
+
+ # Conduct a Chi squared test to see if the actual distribution matches the
+ # expected one at an alpha = 0.05 significance level.
+ #
+ # H_0: The distribution of trees in tree_actual matches the normalized product
+ # of the edge weights in the tree.
+ #
+ # H_a: The distribution of trees in tree_actual follows some other
+ # distribution of spanning trees.
+ chisq, p = stats.chisquare(list(tree_actual.values()), list(tree_expected.values()))
+
+ # Assert that p is greater than the significance level so that we do not
+ # reject the null hypothesis
+ assert not p < 0.05
+
+
+def test_asadpour_tsp():
+ """
+ Test the complete asadpour tsp algorithm with the fractional, symmetric
+ Held Karp solution. This test also uses an incomplete graph as input.
+ """
+ # This version of Figure 2 has all of the edge weights multiplied by 100
+ # and the 0 weight edges have a weight of 1.
+ pytest.importorskip("numpy")
+
+ edge_list = [
+ (0, 1, 100),
+ (0, 2, 100),
+ (0, 5, 1),
+ (1, 2, 100),
+ (1, 4, 1),
+ (2, 3, 1),
+ (3, 4, 100),
+ (3, 5, 100),
+ (4, 5, 100),
+ (1, 0, 100),
+ (2, 0, 100),
+ (5, 0, 1),
+ (2, 1, 100),
+ (4, 1, 1),
+ (3, 2, 1),
+ (4, 3, 100),
+ (5, 3, 100),
+ (5, 4, 100),
+ ]
+
+ G = nx.DiGraph()
+ G.add_weighted_edges_from(edge_list)
+
+ def fixed_asadpour(G, weight):
+ return nx_app.asadpour_atsp(G, weight, 19)
+
+ tour = nx_app.traveling_salesman_problem(G, weight="weight", method=fixed_asadpour)
+
+ # Check that the returned list is a valid tour. Because this is an
+ # incomplete graph, the conditions are not as strict. We need the tour to
+ #
+ # Start and end at the same node
+ # Pass through every vertex at least once
+ # Have a total cost at most ln(6) / ln(ln(6)) = 3.0723 times the optimal
+ #
+ # For the second condition it is possible to have the tour pass through the
+ # same vertex more then. Imagine that the tour on the complete version takes
+ # an edge not in the original graph. In the output this is substituted with
+ # the shortest path between those vertices, allowing vertices to appear more
+ # than once.
+ #
+ # However, we are using a fixed random number generator so we know what the
+ # expected tour is.
+ expected_tours = [[1, 4, 5, 0, 2, 3, 2, 1], [3, 2, 0, 1, 4, 5, 3]]
+
+ assert tour in expected_tours
+
+
+def test_asadpour_real_world():
+ """
+ This test uses airline prices between the six largest cities in the US.
+
+ * New York City -> JFK
+ * Los Angeles -> LAX
+ * Chicago -> ORD
+ * Houston -> IAH
+ * Phoenix -> PHX
+ * Philadelphia -> PHL
+
+ Flight prices from August 2021 using Delta or American airlines to get
+ nonstop flight. The brute force solution found the optimal tour to cost $872
+
+ This test also uses the `source` keyword argument to ensure that the tour
+ always starts at city 0.
+ """
+ np = pytest.importorskip("numpy")
+
+ G_array = np.array(
+ [
+ # JFK LAX ORD IAH PHX PHL
+ [0, 243, 199, 208, 169, 183], # JFK
+ [277, 0, 217, 123, 127, 252], # LAX
+ [297, 197, 0, 197, 123, 177], # ORD
+ [303, 169, 197, 0, 117, 117], # IAH
+ [257, 127, 160, 117, 0, 319], # PHX
+ [183, 332, 217, 117, 319, 0], # PHL
+ ]
+ )
+
+ node_map = {0: "JFK", 1: "LAX", 2: "ORD", 3: "IAH", 4: "PHX", 5: "PHL"}
+
+ expected_tours = [
+ ["JFK", "LAX", "PHX", "ORD", "IAH", "PHL", "JFK"],
+ ["JFK", "ORD", "PHX", "LAX", "IAH", "PHL", "JFK"],
+ ]
+
+ G = nx.from_numpy_array(G_array, create_using=nx.DiGraph)
+ nx.relabel_nodes(G, node_map, copy=False)
+
+ def fixed_asadpour(G, weight):
+ return nx_app.asadpour_atsp(G, weight, 37, source="JFK")
+
+ tour = nx_app.traveling_salesman_problem(G, weight="weight", method=fixed_asadpour)
+
+ assert tour in expected_tours
+
+
+def test_asadpour_real_world_path():
+ """
+ This test uses airline prices between the six largest cities in the US. This
+ time using a path, not a cycle.
+
+ * New York City -> JFK
+ * Los Angeles -> LAX
+ * Chicago -> ORD
+ * Houston -> IAH
+ * Phoenix -> PHX
+ * Philadelphia -> PHL
+
+ Flight prices from August 2021 using Delta or American airlines to get
+ nonstop flight. The brute force solution found the optimal tour to cost $872
+ """
+ np = pytest.importorskip("numpy")
+
+ G_array = np.array(
+ [
+ # JFK LAX ORD IAH PHX PHL
+ [0, 243, 199, 208, 169, 183], # JFK
+ [277, 0, 217, 123, 127, 252], # LAX
+ [297, 197, 0, 197, 123, 177], # ORD
+ [303, 169, 197, 0, 117, 117], # IAH
+ [257, 127, 160, 117, 0, 319], # PHX
+ [183, 332, 217, 117, 319, 0], # PHL
+ ]
+ )
+
+ node_map = {0: "JFK", 1: "LAX", 2: "ORD", 3: "IAH", 4: "PHX", 5: "PHL"}
+
+ expected_paths = [
+ ["ORD", "PHX", "LAX", "IAH", "PHL", "JFK"],
+ ["JFK", "PHL", "IAH", "ORD", "PHX", "LAX"],
+ ]
+
+ G = nx.from_numpy_array(G_array, create_using=nx.DiGraph)
+ nx.relabel_nodes(G, node_map, copy=False)
+
+ def fixed_asadpour(G, weight):
+ return nx_app.asadpour_atsp(G, weight, 56)
+
+ path = nx_app.traveling_salesman_problem(
+ G, weight="weight", cycle=False, method=fixed_asadpour
+ )
+
+ assert path in expected_paths
+
+
+def test_asadpour_disconnected_graph():
+ """
+ Test that the proper exception is raised when asadpour_atsp is given an
+ disconnected graph.
+ """
+
+ G = nx.complete_graph(4, create_using=nx.DiGraph)
+ # have to set edge weights so that if the exception is not raised, the
+ # function will complete and we will fail the test
+ nx.set_edge_attributes(G, 1, "weight")
+ G.add_node(5)
+
+ pytest.raises(nx.NetworkXError, nx_app.asadpour_atsp, G)
+
+
+def test_asadpour_incomplete_graph():
+ """
+ Test that the proper exception is raised when asadpour_atsp is given an
+ incomplete graph
+ """
+
+ G = nx.complete_graph(4, create_using=nx.DiGraph)
+ # have to set edge weights so that if the exception is not raised, the
+ # function will complete and we will fail the test
+ nx.set_edge_attributes(G, 1, "weight")
+ G.remove_edge(0, 1)
+
+ pytest.raises(nx.NetworkXError, nx_app.asadpour_atsp, G)
+
+
+def test_asadpour_empty_graph():
+ """
+ Test the asadpour_atsp function with an empty graph
+ """
+ G = nx.DiGraph()
+
+ pytest.raises(nx.NetworkXError, nx_app.asadpour_atsp, G)
+
+
+def test_asadpour_integral_held_karp():
+ """
+ This test uses an integral held karp solution and the held karp function
+ will return a graph rather than a dict, bypassing most of the asadpour
+ algorithm.
+
+ At first glance, this test probably doesn't look like it ensures that we
+ skip the rest of the asadpour algorithm, but it does. We are not fixing a
+ see for the random number generator, so if we sample any spanning trees
+ the approximation would be different basically every time this test is
+ executed but it is not since held karp is deterministic and we do not
+ reach the portion of the code with the dependence on random numbers.
+ """
+ np = pytest.importorskip("numpy")
+
+ G_array = np.array(
+ [
+ [0, 26, 63, 59, 69, 31, 41],
+ [62, 0, 91, 53, 75, 87, 47],
+ [47, 82, 0, 90, 15, 9, 18],
+ [68, 19, 5, 0, 58, 34, 93],
+ [11, 58, 53, 55, 0, 61, 79],
+ [88, 75, 13, 76, 98, 0, 40],
+ [41, 61, 55, 88, 46, 45, 0],
+ ]
+ )
+
+ G = nx.from_numpy_array(G_array, create_using=nx.DiGraph)
+
+ for _ in range(2):
+ tour = nx_app.traveling_salesman_problem(G, method=nx_app.asadpour_atsp)
+
+ assert [1, 3, 2, 5, 2, 6, 4, 0, 1] == tour
+
+
+def test_directed_tsp_impossible():
+ """
+ Test the asadpour algorithm with a graph without a hamiltonian circuit
+ """
+ pytest.importorskip("numpy")
+
+ # In this graph, once we leave node 0 we cannot return
+ edges = [
+ (0, 1, 10),
+ (0, 2, 11),
+ (0, 3, 12),
+ (1, 2, 4),
+ (1, 3, 6),
+ (2, 1, 3),
+ (2, 3, 2),
+ (3, 1, 5),
+ (3, 2, 1),
+ ]
+
+ G = nx.DiGraph()
+ G.add_weighted_edges_from(edges)
+
+ pytest.raises(nx.NetworkXError, nx_app.traveling_salesman_problem, G)
diff --git a/networkx/algorithms/approximation/traveling_salesman.py b/networkx/algorithms/approximation/traveling_salesman.py
index dc6964c8..598d8be5 100644
--- a/networkx/algorithms/approximation/traveling_salesman.py
+++ b/networkx/algorithms/approximation/traveling_salesman.py
@@ -12,6 +12,7 @@ Categories of algorithms which are implemented:
- Greedy
- Simulated Annealing (SA)
- Threshold Accepting (TA)
+- Asadpour Asymmetric Traveling Salesman Algorithm
The Travelling Salesman Problem tries to find, given the weight
(distance) between all points where a salesman has to visit, the
@@ -33,12 +34,15 @@ important in operations research and theoretical computer science.
http://en.wikipedia.org/wiki/Travelling_salesman_problem
"""
import math
+
import networkx as nx
+
from networkx.utils import py_random_state, not_implemented_for, pairwise
__all__ = [
"traveling_salesman_problem",
"christofides",
+ "asadpour_atsp",
"greedy_tsp",
"simulated_annealing_tsp",
"threshold_accepting_tsp",
@@ -203,14 +207,15 @@ def traveling_salesman_problem(G, weight="weight", nodes=None, cycle=True, metho
graph using the all-pairs shortest_paths between nodes in `nodes`.
Edge weights in the new graph are the lengths of the paths
between each pair of nodes in the original graph.
- Second, an algorithm (default: `christofides`) is used to approximate
- the minimal Hamiltonian cycle on this new graph. The available
- algorithms are:
+ Second, an algorithm (default: `christofides` for undirected and
+ `asadpour_atsp` for directed) is used to approximate the minimal Hamiltonian
+ cycle on this new graph. The available algorithms are:
- christofides
- greedy_tsp
- simulated_annealing_tsp
- threshold_accepting_tsp
+ - asadpour_atsp
Once the Hamiltonian Cycle is found, this function post-processes to
accommodate the structure of the original graph. If `cycle` is ``False``,
@@ -221,7 +226,7 @@ def traveling_salesman_problem(G, weight="weight", nodes=None, cycle=True, metho
Parameters
----------
G : NetworkX graph
- Undirected possibly weighted graph
+ A possibly weighted graph
nodes : collection of nodes (default=G.nodes)
collection (list, set, etc.) of nodes to visit
@@ -255,9 +260,16 @@ def traveling_salesman_problem(G, weight="weight", nodes=None, cycle=True, metho
Returns
-------
list
- List of nodes in `G` along a path with a 3/2-approximation of the minimal
+ List of nodes in `G` along a path with an approximation of the minimal
path through `nodes`.
+
+ Raises
+ ------
+ NetworkXError
+ If `G` is a directed graph it has to be strongly connected or the
+ complete version cannot be generated.
+
Examples
--------
>>> tsp = nx.approximation.traveling_salesman_problem
@@ -280,11 +292,7 @@ def traveling_salesman_problem(G, weight="weight", nodes=None, cycle=True, metho
"""
if method is None:
if G.is_directed():
-
- def threshold_tsp(G, weight):
- return threshold_accepting_tsp(G, "greedy", weight)
-
- method = threshold_tsp
+ method = asadpour_atsp
else:
method = christofides
if nodes is None:
@@ -296,7 +304,13 @@ def traveling_salesman_problem(G, weight="weight", nodes=None, cycle=True, metho
dist[n] = d
path[n] = p
- GG = nx.Graph()
+ if G.is_directed():
+ # If the graph is not strongly connected, raise an exception
+ if not nx.is_strongly_connected(G):
+ raise nx.NetworkXError("G is not strongly connected")
+ GG = nx.DiGraph()
+ else:
+ GG = nx.Graph()
for u in nodes:
for v in nodes:
if u == v:
@@ -306,8 +320,6 @@ def traveling_salesman_problem(G, weight="weight", nodes=None, cycle=True, metho
if not cycle:
# find and remove the biggest edge
- biggest_edge = None
- length_biggest = float("-inf")
(u, v) = max(pairwise(best_GG), key=lambda x: dist[x[0]][x[1]])
pos = best_GG.index(u) + 1
while best_GG[pos] != v:
@@ -321,6 +333,764 @@ def traveling_salesman_problem(G, weight="weight", nodes=None, cycle=True, metho
return best_path
+@not_implemented_for("undirected")
+def asadpour_atsp(G, weight="weight", seed=None, source=None):
+ """
+ Returns an approximate solution to the traveling salesman problem.
+
+ This approximate solution is one of the best known approximations for the
+ asymmetric traveling salesman problem developed by Asadpour et al,
+ [1]_. The algorithm first solves the Held-Karp relaxation to find a lower
+ bound for the weight of the cycle. Next, it constructs an exponential
+ distribution of undirected spanning trees where the probability of an
+ edge being in the tree corresponds to the weight of that edge using a
+ maximum entropy rounding scheme. Next we sample that distribution
+ $2 \\lceil \\ln n \\rceil$ times and save the minimum sampled tree once the
+ direction of the arcs is added back to the edges. Finally, we augment
+ then short circuit that graph to find the approximate tour for the
+ salesman.
+
+ Parameters
+ ----------
+ G : nx.DiGraph
+ The graph should be a complete weighted directed graph. The
+ distance between all paris of nodes should be included and the triangle
+ inequality should hold. That is, the direct edge between any two nodes
+ should be the path of least cost.
+
+ weight : string, optional (default="weight")
+ Edge data key corresponding to the edge weight.
+ If any edge does not have this attribute the weight is set to 1.
+
+ seed : integer, random_state, or None (default)
+ Indicator of random number generation state.
+ See :ref:`Randomness<randomness>`.
+
+ source : node label (default=`None`)
+ If given, return the cycle starting and ending at the given node.
+
+ Returns
+ -------
+ cycle : list of nodes
+ Returns the cycle (list of nodes) that a salesman can follow to minimize
+ the total weight of the trip.
+
+ Raises
+ ------
+ NetworkXError
+ If `G` is not complete or has less than two nodes, the algorithm raises
+ an exception.
+
+ NetworkXError
+ If 'source` is not `None` and is not a node in `G`, the algorithm raises
+ an exception.
+
+ NetworkXNotImplemented
+ If `G` is an undirected graph.
+
+ References
+ ----------
+ .. [1] A. Asadpour, M. X. Goemans, A. Madry, S. O. Gharan, and A. Saberi,
+ An o(log n/log log n)-approximation algorithm for the asymmetric
+ traveling salesman problem, Operations research, 65 (2017),
+ pp. 1043–1061
+
+ Examples
+ --------
+ >>> import networkx as nx
+ >>> import networkx.algorithms.approximation as approx
+ >>> G = nx.complete_graph(3, create_using=nx.DiGraph)
+ >>> nx.set_edge_attributes(G, {(0, 1): 2, (1, 2): 2, (2, 0): 2, (0, 2): 1, (2, 1): 1, (1, 0): 1}, "weight")
+ >>> tour = approx.asadpour_atsp(G,source=0)
+ >>> tour
+ [0, 2, 1, 0]
+ """
+ from math import log as ln
+ from math import exp
+ from math import ceil
+
+ # Check that G is a complete graph
+ N = len(G) - 1
+ if N < 2:
+ raise nx.NetworkXError("G must have at least two nodes")
+ # This check ignores selfloops which is what we want here.
+ if any(len(nbrdict) - (n in nbrdict) != N for n, nbrdict in G.adj.items()):
+ raise nx.NetworkXError("G is not a complete DiGraph")
+ # Check that the source vertex, if given, is in the graph
+ if source is not None and source not in G.nodes:
+ raise nx.NetworkXError("Given source node not in G.")
+
+ opt_hk, z_star = held_karp_ascent(G, weight)
+
+ # Test to see if the ascent method found an integer solution or a fractional
+ # solution. If it is integral then z_star is a nx.Graph, otherwise it is
+ # a dict
+ if not isinstance(z_star, dict):
+ # Here we are using the shortcutting method to go from the list of edges
+ # returned from eularian_circuit to a list of nodes
+ return _shortcutting(nx.eulerian_circuit(z_star, source=source))
+
+ # Create the undirected support of z_star
+ z_support = nx.MultiGraph()
+ for u, v in z_star:
+ if (u, v) not in z_support.edges:
+ edge_weight = min(G[u][v][weight], G[v][u][weight])
+ z_support.add_edge(u, v, **{weight: edge_weight})
+
+ # Create the exponential distribution of spanning trees
+ gamma = spanning_tree_distribution(z_support, z_star)
+
+ # Write the lambda values to the edges of z_support
+ z_support = nx.Graph(z_support)
+ lambda_dict = {(u, v): exp(gamma[(u, v)]) for u, v in z_support.edges()}
+ nx.set_edge_attributes(z_support, lambda_dict, "lambda_key")
+ del gamma, lambda_dict
+
+ # Sample 2 * ceil( ln(n) ) spanning trees and record the minimum one
+ minimum_sampled_tree = None
+ minimum_sampled_tree_weight = math.inf
+ for _ in range(2 * ceil(ln(G.number_of_nodes()))):
+ sampled_tree = sample_spanning_tree(z_support, "lambda_key", seed)
+ sampled_tree_weight = sampled_tree.size(weight)
+ if sampled_tree_weight < minimum_sampled_tree_weight:
+ minimum_sampled_tree = sampled_tree.copy()
+ minimum_sampled_tree_weight = sampled_tree_weight
+
+ # Orient the edges in that tree to keep the cost of the tree the same.
+ t_star = nx.MultiDiGraph()
+ for u, v, d in minimum_sampled_tree.edges(data=weight):
+ if d == G[u][v][weight]:
+ t_star.add_edge(u, v, **{weight: d})
+ else:
+ t_star.add_edge(v, u, **{weight: d})
+
+ # Find the node demands needed to neutralize the flow of t_star in G
+ node_demands = {n: t_star.out_degree(n) - t_star.in_degree(n) for n in t_star}
+ nx.set_node_attributes(G, node_demands, "demand")
+
+ # Find the min_cost_flow
+ flow_dict = nx.min_cost_flow(G, "demand")
+
+ # Build the flow into t_star
+ for source, values in flow_dict.items():
+ for target in values:
+ if (source, target) not in t_star.edges and values[target] > 0:
+ # IF values[target] > 0 we have to add that many edges
+ for _ in range(values[target]):
+ t_star.add_edge(source, target)
+
+ # Return the shortcut eulerian circuit
+ circuit = nx.eulerian_circuit(t_star, source=source)
+ return _shortcutting(circuit)
+
+
+def held_karp_ascent(G, weight="weight"):
+ """
+ Minimizes the Held-Karp relaxation of the TSP for `G`
+
+ Solves the Held-Karp relaxation of the input complete digraph and scales
+ the output solution for use in the Asadpour [1]_ ASTP algorithm.
+
+ The Held-Karp relaxation defines the lower bound for solutions to the
+ ATSP, although it does return a fractional solution. This is used in the
+ Asadpour algorithm as an initial solution which is later rounded to a
+ integral tree within the spanning tree polytopes. This function solves
+ the relaxation with the branch and bound method in [2]_.
+
+ Parameters
+ ----------
+ G : nx.DiGraph
+ The graph should be a complete weighted directed graph.
+ The distance between all paris of nodes should be included.
+
+ weight : string, optional (default="weight")
+ Edge data key corresponding to the edge weight.
+ If any edge does not have this attribute the weight is set to 1.
+
+ Returns
+ -------
+ OPT : float
+ The cost for the optimal solution to the Held-Karp relaxation
+ z : dict or nx.Graph
+ A symmetrized and scaled version of the optimal solution to the
+ Held-Karp relaxation for use in the Asadpour algorithm.
+
+ If an integral solution is found, then that is an optimal solution for
+ the ATSP problem and that is returned instead.
+
+ References
+ ----------
+ .. [1] A. Asadpour, M. X. Goemans, A. Madry, S. O. Gharan, and A. Saberi,
+ An o(log n/log log n)-approximation algorithm for the asymmetric
+ traveling salesman problem, Operations research, 65 (2017),
+ pp. 1043–1061
+
+ .. [2] M. Held, R. M. Karp, The traveling-salesman problem and minimum
+ spanning trees, Operations Research, 1970-11-01, Vol. 18 (6),
+ pp.1138-1162
+ """
+ import numpy as np
+ import scipy.optimize as optimize
+
+ def k_pi():
+ """
+ Find the set of minimum 1-Arborescences for G at point pi.
+
+ Returns
+ -------
+ Set
+ The set of minimum 1-Arborescences
+ """
+ # Create a copy of G without vertex 1.
+ G_1 = G.copy()
+ minimum_1_arborescences = set()
+ minimum_1_arborescence_weight = math.inf
+
+ # node is node '1' in the Held and Karp paper
+ n = next(G.__iter__())
+ G_1.remove_node(n)
+
+ # Iterate over the spanning arborescences of the graph until we know
+ # that we have found the minimum 1-arborescences. My proposed strategy
+ # is to find the most extensive root to connect to from 'node 1' and
+ # the least expensive one. We then iterate over arborescences until
+ # the cost of the basic arborescence is the cost of the minimum one
+ # plus the difference between the most and least expensive roots,
+ # that way the cost of connecting 'node 1' will by definition not by
+ # minimum
+ min_root = {"node": None, weight: math.inf}
+ max_root = {"node": None, weight: -math.inf}
+ for u, v, d in G.edges(n, data=True):
+ if d[weight] < min_root[weight]:
+ min_root = {"node": v, weight: d[weight]}
+ if d[weight] > max_root[weight]:
+ max_root = {"node": v, weight: d[weight]}
+
+ min_in_edge = min(G.in_edges(n, data=True), key=lambda x: x[2][weight])
+ min_root[weight] = min_root[weight] + min_in_edge[2][weight]
+ max_root[weight] = max_root[weight] + min_in_edge[2][weight]
+
+ min_arb_weight = math.inf
+ for arb in nx.ArborescenceIterator(G_1):
+ arb_weight = arb.size(weight)
+ if min_arb_weight == math.inf:
+ min_arb_weight = arb_weight
+ elif arb_weight > min_arb_weight + max_root[weight] - min_root[weight]:
+ break
+ # We have to pick the root node of the arborescence for the out
+ # edge of the first vertex as that is the only node without an
+ # edge directed into it.
+ for N, deg in arb.in_degree:
+ if deg == 0:
+ # root found
+ arb.add_edge(n, N, **{weight: G[n][N][weight]})
+ arb_weight += G[n][N][weight]
+ break
+
+ # We can pick the minimum weight in-edge for the vertex with
+ # a cycle. If there are multiple edges with the same, minimum
+ # weight, We need to add all of them.
+ #
+ # Delete the edge (N, v) so that we cannot pick it.
+ edge_data = G[N][n]
+ G.remove_edge(N, n)
+ min_weight = min(G.in_edges(n, data=weight), key=lambda x: x[2])[2]
+ min_edges = [
+ (u, v, d) for u, v, d in G.in_edges(n, data=weight) if d == min_weight
+ ]
+ for u, v, d in min_edges:
+ new_arb = arb.copy()
+ new_arb.add_edge(u, v, **{weight: d})
+ new_arb_weight = arb_weight + d
+ # Check to see the weight of the arborescence, if it is a
+ # new minimum, clear all of the old potential minimum
+ # 1-arborescences and add this is the only one. If its
+ # weight is above the known minimum, do not add it.
+ if new_arb_weight < minimum_1_arborescence_weight:
+ minimum_1_arborescences.clear()
+ minimum_1_arborescence_weight = new_arb_weight
+ # We have a 1-arborescence, add it to the set
+ if new_arb_weight == minimum_1_arborescence_weight:
+ minimum_1_arborescences.add(new_arb)
+ G.add_edge(N, n, **edge_data)
+
+ return minimum_1_arborescences
+
+ def direction_of_ascent():
+ """
+ Find the direction of ascent at point pi.
+
+ See [1]_ for more information.
+
+ Returns
+ -------
+ dict
+ A mapping from the nodes of the graph which represents the direction
+ of ascent.
+
+ References
+ ----------
+ .. [1] M. Held, R. M. Karp, The traveling-salesman problem and minimum
+ spanning trees, Operations Research, 1970-11-01, Vol. 18 (6),
+ pp.1138-1162
+ """
+ # 1. Set d equal to the zero n-vector.
+ d = {}
+ for n in G:
+ d[n] = 0
+ del n
+ # 2. Find a 1-Aborescence T^k such that k is in K(pi, d).
+ minimum_1_arborescences = k_pi()
+ while True:
+ # Reduce K(pi) to K(pi, d)
+ # Find the arborescence in K(pi) which increases the lest in
+ # direction d
+ min_k_d_weight = math.inf
+ min_k_d = None
+ for arborescence in minimum_1_arborescences:
+ weighted_cost = 0
+ for n, deg in arborescence.degree:
+ weighted_cost += d[n] * (deg - 2)
+ if weighted_cost < min_k_d_weight:
+ min_k_d_weight = weighted_cost
+ min_k_d = arborescence
+
+ # 3. If sum of d_i * v_{i, k} is greater than zero, terminate
+ if min_k_d_weight > 0:
+ return d, min_k_d
+ # 4. d_i = d_i + v_{i, k}
+ for n, deg in min_k_d.degree:
+ d[n] += deg - 2
+ # Check that we do not need to terminate because the direction
+ # of ascent does not exist. This is done with linear
+ # programming.
+ c = np.full(len(minimum_1_arborescences), -1, dtype=int)
+ a_eq = np.empty((len(G) + 1, len(minimum_1_arborescences)), dtype=int)
+ b_eq = np.zeros(len(G) + 1, dtype=int)
+ b_eq[len(G)] = 1
+ arb_count = 0
+ for arborescence in minimum_1_arborescences:
+ n_count = len(G) - 1
+ for n, deg in arborescence.degree:
+ a_eq[n_count][arb_count] = deg - 2
+ n_count -= 1
+ a_eq[len(G)][arb_count] = 1
+ arb_count += 1
+ program_result = optimize.linprog(c, A_eq=a_eq, b_eq=b_eq)
+ bool_result = program_result.x >= 0
+ if program_result.status == 0 and np.sum(bool_result) == len(
+ minimum_1_arborescences
+ ):
+ # There is no direction of ascent
+ return None, minimum_1_arborescences
+
+ # 5. GO TO 2
+
+ def find_epsilon(k, d):
+ """
+ Given the direction of ascent at pi, find the maximum distance we can go
+ in that direction.
+
+ Parameters
+ ----------
+ k_xy : set
+ The set of 1-arborescences which have the minimum rate of increase
+ in the direction of ascent
+
+ d : dict
+ The direction of ascent
+
+ Returns
+ -------
+ float
+ The distance we can travel in direction `d`
+ """
+ min_epsilon = math.inf
+ for e_u, e_v, e_w in G.edges(data=weight):
+ if (e_u, e_v) in k.edges:
+ continue
+ # Now, I have found a condition which MUST be true for the edges to
+ # be a valid substitute. The edge in the graph which is the
+ # substitute is the one with the same terminated end. This can be
+ # checked rather simply.
+ #
+ # Find the edge within k which is the substitute. Because k is a
+ # 1-arborescence, we know that they is only one such edges
+ # leading into every vertex.
+ if len(k.in_edges(e_v, data=weight)) > 1:
+ raise Exception
+ sub_u, sub_v, sub_w = next(k.in_edges(e_v, data=weight).__iter__())
+ k.add_edge(e_u, e_v, **{weight: e_w})
+ k.remove_edge(sub_u, sub_v)
+ if (
+ max(d for n, d in k.in_degree()) <= 1
+ and len(G) == k.number_of_edges()
+ and nx.is_weakly_connected(k)
+ ):
+ # Ascent method calculation
+ if d[sub_u] == d[e_u] or sub_w == e_w:
+ # Revert to the original graph
+ k.remove_edge(e_u, e_v)
+ k.add_edge(sub_u, sub_v, **{weight: sub_w})
+ continue
+ epsilon = (sub_w - e_w) / (d[e_u] - d[sub_u])
+ if 0 < epsilon < min_epsilon:
+ min_epsilon = epsilon
+ # Revert to the original graph
+ k.remove_edge(e_u, e_v)
+ k.add_edge(sub_u, sub_v, **{weight: sub_w})
+
+ return min_epsilon
+
+ # I have to know that the elements in pi correspond to the correct elements
+ # in the direction of ascent, even if the node labels are not integers.
+ # Thus, I will use dictionaries to made that mapping.
+ pi_dict = {}
+ for n in G:
+ pi_dict[n] = 0
+ del n
+ original_edge_weights = {}
+ for u, v, d in G.edges(data=True):
+ original_edge_weights[(u, v)] = d[weight]
+ dir_ascent, k_d = direction_of_ascent()
+ while dir_ascent is not None:
+ max_distance = find_epsilon(k_d, dir_ascent)
+ for n, v in dir_ascent.items():
+ pi_dict[n] += max_distance * v
+ for u, v, d in G.edges(data=True):
+ d[weight] = original_edge_weights[(u, v)] + pi_dict[u]
+ dir_ascent, k_d = direction_of_ascent()
+ # k_d is no longer an individual 1-arborescence but rather a set of
+ # minimal 1-arborescences at the maximum point of the polytope and should
+ # be reflected as such
+ k_max = k_d
+
+ # Search for a cycle within k_max. If a cycle exists, return it as the
+ # solution
+ for k in k_max:
+ if len([n for n in k if k.degree(n) == 2]) == G.order():
+ # Tour found
+ return k.size(weight), k
+
+ # Write the original edge weights back to G and every member of k_max at
+ # the maximum point. Also average the number of times that edge appears in
+ # the set of minimal 1-arborescences.
+ x_star = {}
+ size_k_max = len(k_max)
+ for u, v, d in G.edges(data=True):
+ edge_count = 0
+ d[weight] = original_edge_weights[(u, v)]
+ for k in k_max:
+ if (u, v) in k.edges():
+ edge_count += 1
+ k[u][v][weight] = original_edge_weights[(u, v)]
+ x_star[(u, v)] = edge_count / size_k_max
+ # Now symmetrize the edges in x_star and scale them according to (5) in
+ # reference [1]
+ z_star = {}
+ scale_factor = (G.order() - 1) / G.order()
+ for u, v in x_star.keys():
+ frequency = x_star[(u, v)] + x_star[(v, u)]
+ if frequency > 0:
+ z_star[(u, v)] = scale_factor * frequency
+ del x_star
+ # Return the optimal weight and the z dict
+ return next(k_max.__iter__()).size(weight), z_star
+
+
+def total_spanning_tree_weight(G, weight=None):
+ """
+ Apply Kirchhoff's Tree Matrix Theorem a graph in order to find the total
+ weight of all spanning trees.
+
+ The theorem states that the determinant of any cofactor of the Laplacian
+ matrix of a graph is the number of spanning trees in the graph. For a
+ weighted Laplacian matrix, it is the sum across all spanning trees of the
+ multiplicative weight of each tree. That is, the weight of each tree is the
+ product of its edge weights.
+
+ Parameters
+ ----------
+ G : NetworkX Graph
+ The graph to use Kirchhoff's theorem on.
+
+ weight : string or None
+ The key for the edge attribute holding the edge weight. If `None`, then
+ each edge is assumed to have a weight of 1.
+
+ Returns
+ -------
+ float
+ The sum of the total multiplicative weight for all spanning trees in the
+ graph.
+ """
+ import numpy as np
+
+ G_laplacian = nx.laplacian_matrix(G, weight=weight).toarray()
+ # Determinant ignoring first row and column
+ return abs(np.linalg.det(G_laplacian[1:, 1:]))
+
+
+def spanning_tree_distribution(G, z):
+ """
+ Find the asadpour exponential distribution of spanning trees.
+
+ Solves the Maximum Entropy Convex Program in the Asadpour algorithm [1]_
+ using the approach in section 7 to build an exponential distribution of
+ undirected spanning trees.
+
+ This algorithm ensures that the probability of any edge in a spanning
+ tree is proportional to the sum of the probabilities of the tress
+ containing that edge over the sum of the probabilities of all spanning
+ trees of the graph.
+
+ Parameters
+ ----------
+ G : nx.MultiGraph
+ The undirected support graph for the Held Karp relaxation
+
+ z : dict
+ The output of `held_karp_ascent()`, a scaled version of the Held-Karp
+ solution.
+
+ Returns
+ -------
+ gamma : dict
+ The probability distribution which approximately preserves the marginal
+ probabilities of `z`.
+ """
+ from math import exp
+ from math import log as ln
+
+ def q(e):
+ """
+ The value of q(e) is described in the Asadpour paper is "the
+ probability that edge e will be included in a spanning tree T that is
+ chosen with probability proportional to exp(gamma(T))" which
+ basically means that it is the total probability of the edge appearing
+ across the whole distribution.
+
+ Parameters
+ ----------
+ e : tuple
+ The `(u, v)` tuple describing the edge we are interested in
+
+ Returns
+ -------
+ float
+ The probability that a spanning tree chosen according to the
+ current values of gamma will include edge `e`.
+ """
+ # Create the laplacian matrices
+ for u, v, d in G.edges(data=True):
+ d[lambda_key] = exp(gamma[(u, v)])
+ G_Kirchhoff = total_spanning_tree_weight(G, lambda_key)
+ G_e = nx.contracted_edge(G, e, self_loops=False)
+ G_e_Kirchhoff = total_spanning_tree_weight(G_e, lambda_key)
+
+ # Multiply by the weight of the contracted edge since it is not included
+ # in the total weight of the contracted graph.
+ return exp(gamma[(e[0], e[1])]) * G_e_Kirchhoff / G_Kirchhoff
+
+ # initialize gamma to the zero dict
+ gamma = {}
+ for u, v, _ in G.edges:
+ gamma[(u, v)] = 0
+
+ # set epsilon
+ EPSILON = 0.2
+
+ # pick an edge attribute name that is unlikely to be in the graph
+ lambda_key = "spanning_tree_distribution's secret attribute name for lambda"
+
+ while True:
+ # We need to know that know that no values of q_e are greater than
+ # (1 + epsilon) * z_e, however changing one gamma value can increase the
+ # value of a different q_e, so we have to complete the for loop without
+ # changing anything for the condition to be meet
+ in_range_count = 0
+ # Search for an edge with q_e > (1 + epsilon) * z_e
+ for u, v in gamma:
+ e = (u, v)
+ q_e = q(e)
+ z_e = z[e]
+ if q_e > (1 + EPSILON) * z_e:
+ delta = ln(
+ (q_e * (1 - (1 + EPSILON / 2) * z_e))
+ / ((1 - q_e) * (1 + EPSILON / 2) * z_e)
+ )
+ gamma[e] -= delta
+ # Check that delta had the desired effect
+ new_q_e = q(e)
+ desired_q_e = (1 + EPSILON / 2) * z_e
+ if round(new_q_e, 8) != round(desired_q_e, 8):
+ raise nx.NetworkXError(
+ f"Unable to modify probability for edge ({u}, {v})"
+ )
+ else:
+ in_range_count += 1
+ # Check if the for loop terminated without changing any gamma
+ if in_range_count == len(gamma):
+ break
+
+ # Remove the new edge attributes
+ for _, _, d in G.edges(data=True):
+ if lambda_key in d:
+ del d[lambda_key]
+
+ return gamma
+
+
+@py_random_state(2)
+def sample_spanning_tree(G, lambda_key, seed=None):
+ """
+ Sample a spanning tree using the edges weights of the graph.
+
+ The edge weights are multiplicative, so the probability of each tree is
+ proportional to the product of edge weights.
+
+ The algorithm itself uses algorithm A8 in [1]_ .
+
+ We 'shuffle' the edges in the graph, and then probabilistically
+ determine weather to add the edge conditioned on all of the previous
+ edges which where added to the tree. Probabilities are calculated using
+ Kirchhoff's Matrix Tree Theorem and a weighted Laplacian matrix.
+
+ At each iteration, we contract the edges we have decided to include in the
+ sampled tree and delete those which we have decided not to include.
+
+ Parameters
+ ----------
+ G : nx.Graph
+ An undirected version of the original graph.
+
+ lambda_key : string
+ The edge key for the edge attribute holding edge weight.
+
+ seed : integer, random_state, or None (default)
+ Indicator of random number generation state.
+ See :ref:`Randomness<randomness>`.
+
+ Returns
+ -------
+ nx.Graph
+ A spanning tree using the distribution defined by `gamma`.
+
+ References
+ ----------
+ .. [1] V. Kulkarni, Generating random combinatorial objects, Journal of
+ algorithms, 11 (1990), pp. 185–207
+ """
+
+ def find_node(merged_nodes, n):
+ """
+ We can think of clusters of contracted nodes as having one
+ representative in the graph. Each node which is not in merged_nodes
+ is still its own representative. Since a representative can be later
+ contracted, we need to recursively search though the dict to find
+ the final representative, but once we know it we can use path
+ compression to speed up the access of the representative for next time.
+
+ Parameters
+ ----------
+ merged_nodes : dict
+ The dict storing the mapping from node to representative
+ n
+ The node whose representative we seek
+
+ Returns
+ -------
+ The representative of the `n`
+ """
+ if n not in merged_nodes:
+ return n
+ else:
+ rep = find_node(merged_nodes, merged_nodes[n])
+ merged_nodes[n] = rep
+ return rep
+
+ def prepare_graph():
+ """
+ For the graph `G`, remove all edges not in the set `V` and then
+ contract all edges in the set `U`.
+
+ Returns
+ -------
+ A copy of `G` which has had all edges not in `V` removed and all edges
+ in `U` contracted.
+ """
+
+ # The result is a MultiGraph version of G so that parallel edges are
+ # allowed during edge contraction
+ result = nx.MultiGraph(incoming_graph_data=G)
+
+ # Remove all edges not in V
+ edges_to_remove = set(result.edges()).difference(V)
+ result.remove_edges_from(edges_to_remove)
+
+ # Contract all edges in U
+ #
+ # Imagine that you have two edges to contract and they share an
+ # endpoint like this:
+ # [0] ----- [1] ----- [2]
+ # If we contract (0, 1) first, the contraction function will always
+ # delete the second node it is passed so the resulting graph would be
+ # [0] ----- [2]
+ # and edge (1, 2) no longer exists but (0, 2) would need to be contracted
+ # in its place now. That is why I use the below dict as a merge-find
+ # data structure with path compression to track how the nodes are merged.
+ merged_nodes = {}
+
+ for u, v in U:
+ u_rep = find_node(merged_nodes, u)
+ v_rep = find_node(merged_nodes, v)
+ # We cannot contract a node with itself
+ if u_rep == v_rep:
+ continue
+ nx.contracted_nodes(result, u_rep, v_rep, self_loops=False, copy=False)
+ merged_nodes[v_rep] = u_rep
+
+ return merged_nodes, result
+
+ U = set()
+ V = set(G.edges())
+ shuffled_edges = list(G.edges())
+ seed.shuffle(shuffled_edges)
+
+ for u, v in shuffled_edges:
+ node_map, prepared_G = prepare_graph()
+ G_total_tree_weight = total_spanning_tree_weight(prepared_G, lambda_key)
+ # Add the edge to U so that we can compute the total tree weight
+ # assuming we include that edge
+ U.add((u, v))
+ # Now, if (u, v) cannot exist in G because it is fully contracted out
+ # of existence, then it by definition cannot influence G_e's Kirchhoff
+ # value. But, we also cannot pick it.
+ _, prepared_G_e = prepare_graph()
+ rep_edge = (find_node(node_map, u), find_node(node_map, v))
+ # Check to see if the 'representative edge' for the current edge is
+ # in prepared_G. If so, then we can pick it.
+ if rep_edge in prepared_G.edges:
+ G_e_total_tree_weight = total_spanning_tree_weight(prepared_G_e, lambda_key)
+ else:
+ G_e_total_tree_weight = 0.0
+ z = seed.uniform(0.0, 1.0)
+ # This will be useful if I move this random spanning tree method to
+ # the boarder NetworkX library
+ e_weight = G[u][v][lambda_key] if lambda_key is not None else 1
+ if z > e_weight * G_e_total_tree_weight / G_total_tree_weight:
+ # Remove the edge from U since we did not decide to include it in
+ # the sampled spanning tree. Also remove the edge from V because if
+ # we did not decide to include it we must reject it.
+ U.remove((u, v))
+ V.remove((u, v))
+ # If we decide to keep an edge, it may complete the spanning tree.
+ elif len(U) == G.number_of_nodes() - 1:
+ spanning_tree = nx.Graph()
+ spanning_tree.add_edges_from(U)
+ return spanning_tree
+
+
def greedy_tsp(G, weight="weight", source=None):
"""Return a low cost cycle starting at `source` and its cost.
@@ -682,7 +1452,6 @@ def threshold_accepting_tsp(
move : "1-1" or "1-0" or function, optional (default="1-1")
Indicator of what move to use when finding new trial solutions.
Strings indicate two special built-in moves:
-
- "1-1": 1-1 exchange which transposes the position
of two elements of the current solution.
The function called is :func:`swap_two_nodes`.
diff --git a/networkx/algorithms/minors/contraction.py b/networkx/algorithms/minors/contraction.py
index 199c2272..04574bc3 100644
--- a/networkx/algorithms/minors/contraction.py
+++ b/networkx/algorithms/minors/contraction.py
@@ -534,7 +534,7 @@ def contracted_nodes(G, u, v, self_loops=True, copy=True):
identified_nodes = contracted_nodes
-def contracted_edge(G, edge, self_loops=True):
+def contracted_edge(G, edge, self_loops=True, copy=True):
"""Returns the graph that results from contracting the specified edge.
Edge contraction identifies the two endpoints of the edge as a single node
@@ -555,6 +555,10 @@ def contracted_edge(G, edge, self_loops=True):
endpoints of `edge` in `G` become self-loops on the new node in the
returned graph.
+ copy : Boolean (default True)
+ If this is True, a the contraction will be performed on a copy of `G`,
+ otherwise the contraction will happen in place.
+
Returns
-------
Networkx graph
@@ -596,4 +600,4 @@ def contracted_edge(G, edge, self_loops=True):
u, v = edge[:2]
if not G.has_edge(u, v):
raise ValueError(f"Edge {edge} does not exist in graph G; cannot contract it")
- return contracted_nodes(G, u, v, self_loops=self_loops)
+ return contracted_nodes(G, u, v, self_loops=self_loops, copy=copy)
diff --git a/networkx/algorithms/tree/branchings.py b/networkx/algorithms/tree/branchings.py
index c68564e0..410be310 100644
--- a/networkx/algorithms/tree/branchings.py
+++ b/networkx/algorithms/tree/branchings.py
@@ -26,17 +26,17 @@ This implementation is based on:
# pages={109-122},
# language={English}
# }
-
-
import string
+from dataclasses import dataclass, field
+from enum import Enum
from operator import itemgetter
+from queue import PriorityQueue
import networkx as nx
from networkx.utils import py_random_state
from .recognition import is_arborescence, is_branching
-
__all__ = [
"branching_weight",
"greedy_branching",
@@ -44,6 +44,7 @@ __all__ = [
"minimum_branching",
"maximum_spanning_arborescence",
"minimum_spanning_arborescence",
+ "ArborescenceIterator",
"Edmonds",
]
@@ -239,6 +240,7 @@ def get_path(G, u, v):
"""
nodes = nx.shortest_path(G, u, v)
+
# We are guaranteed that there is only one edge connected every node
# in the shortest path.
@@ -255,7 +257,23 @@ def get_path(G, u, v):
class Edmonds:
"""
- Edmonds algorithm for finding optimal branchings and spanning arborescences.
+ Edmonds algorithm [1]_ for finding optimal branchings and spanning
+ arborescences.
+
+ This algorithm can find both minimum and maximum spanning arborescences and
+ branchings.
+
+ Notes
+ -----
+ While this algorithm can find a minimum branching, since it isn't required
+ to be spanning, the minimum branching is always from the set of negative
+ weight edges which is most likely the empty set for most graphs.
+
+ References
+ ----------
+ .. [1] J. Edmonds, Optimum Branchings, Journal of Research of the National
+ Bureau of Standards, 1967, Vol. 71B, p.233-240,
+ https://archive.org/details/jresv71Bn4p233
"""
@@ -272,7 +290,7 @@ class Edmonds:
# sure that our node names do not conflict with the real node names.
self.template = random_string(seed=seed) + "_{0}"
- def _init(self, attr, default, kind, style, preserve_attrs, seed):
+ def _init(self, attr, default, kind, style, preserve_attrs, seed, partition):
if kind not in KINDS:
raise nx.NetworkXException("Unknown value for `kind`.")
@@ -305,6 +323,9 @@ class Edmonds:
for key, (u, v, data) in enumerate(self.G_original.edges(data=True)):
d = {attr: trans(data.get(attr, default))}
+ if data.get(partition) is not None:
+ d[partition] = data.get(partition)
+
if preserve_attrs:
for (d_k, d_v) in data.items():
if d_k != attr:
@@ -330,11 +351,12 @@ class Edmonds:
# A list of lists of edge indexes. Each list is a circuit for graph G^i.
# Note the edge list will not, in general, be a circuit in graph G^0.
self.circuits = []
- # Stores the index of the minimum edge in the circuit found in G^i and B^i.
- # The ordering of the edges seems to preserve the weight ordering from G^0.
- # So even if the circuit does not form a circuit in G^0, it is still true
- # that the minimum edge of the circuit in G^i is still the minimum edge
- # in circuit G^0 (depsite their weights being different).
+ # Stores the index of the minimum edge in the circuit found in G^i
+ # and B^i. The ordering of the edges seems to preserve the weight
+ # ordering from G^0. So even if the circuit does not form a circuit
+ # in G^0, it is still true that the minimum edge of the circuit in
+ # G^i is still the minimum edge in circuit G^0 (despite their weights
+ # being different).
self.minedge_circuit = []
def find_optimum(
@@ -344,6 +366,7 @@ class Edmonds:
kind="max",
style="branching",
preserve_attrs=False,
+ partition=None,
seed=None,
):
"""
@@ -367,6 +390,9 @@ class Edmonds:
preserve_attrs : bool
If True, preserve the other edge attributes of the original
graph (that are not the one passed to `attr`)
+ partition : str
+ The edge attribute holding edge partition data. Used in the
+ spanning arborescence iterator.
seed : integer, random_state, or None (default)
Indicator of random number generation state.
See :ref:`Randomness<randomness>`.
@@ -377,7 +403,7 @@ class Edmonds:
The branching.
"""
- self._init(attr, default, kind, style, preserve_attrs, seed)
+ self._init(attr, default, kind, style, preserve_attrs, seed, partition)
uf = self.uf
# This enormous while loop could use some refactoring...
@@ -392,14 +418,27 @@ class Edmonds:
"""
Find the edge directed toward v with maximal weight.
+ If an edge partition exists in this graph, return the included edge
+ if it exists and no not return any excluded edges. There can only
+ be one included edge for each vertex otherwise the edge partition is
+ empty.
"""
edge = None
weight = -INF
for u, _, key, data in G.in_edges(v, data=True, keys=True):
+ # Skip excluded edges
+ if data.get(partition) == nx.EdgePartition.EXCLUDED:
+ continue
new_weight = data[attr]
+ # Return the included edge
+ if data.get(partition) == nx.EdgePartition.INCLUDED:
+ weight = new_weight
+ edge = (u, v, key, new_weight, data)
+ return edge, weight
+ # Find the best open edge
if new_weight > weight:
weight = new_weight
- edge = (u, v, key, new_weight)
+ edge = (u, v, key, new_weight, data)
return edge, weight
@@ -452,7 +491,7 @@ class Edmonds:
# from the paper does hold. We need to store the circuit
# for future reference.
Q_nodes, Q_edges = get_path(B, v, u)
- Q_edges.append(edge[2])
+ Q_edges.append(edge[2]) # Edge key
else:
# Then B with the edge is still a branching and condition
# (a) from the paper does not hold.
@@ -468,6 +507,8 @@ class Edmonds:
# print(f"Edge is acceptable: {acceptable}")
if acceptable:
dd = {attr: weight}
+ if edge[4].get(partition) is not None:
+ dd[partition] = edge[4].get(partition)
B.add_edge(u, v, edge[2], **dd)
G[u][v][edge[2]][self.candidate_attr] = True
uf.union(u, v)
@@ -486,8 +527,12 @@ class Edmonds:
Q_incoming_weight = {}
for edge_key in Q_edges:
u, v, data = B.edge_index[edge_key]
+ # We cannot remove an included edges, even if it is
+ # the minimum edge in the circuit
w = data[attr]
Q_incoming_weight[v] = w
+ if data.get(partition) == nx.EdgePartition.INCLUDED:
+ continue
if w < minweight:
minweight = w
minedge = edge_key
@@ -638,26 +683,47 @@ class Edmonds:
return H
-def maximum_branching(G, attr="weight", default=1, preserve_attrs=False):
+def maximum_branching(
+ G, attr="weight", default=1, preserve_attrs=False, partition=None
+):
ed = Edmonds(G)
B = ed.find_optimum(
- attr, default, kind="max", style="branching", preserve_attrs=preserve_attrs
+ attr,
+ default,
+ kind="max",
+ style="branching",
+ preserve_attrs=preserve_attrs,
+ partition=partition,
)
return B
-def minimum_branching(G, attr="weight", default=1, preserve_attrs=False):
+def minimum_branching(
+ G, attr="weight", default=1, preserve_attrs=False, partition=None
+):
ed = Edmonds(G)
B = ed.find_optimum(
- attr, default, kind="min", style="branching", preserve_attrs=preserve_attrs
+ attr,
+ default,
+ kind="min",
+ style="branching",
+ preserve_attrs=preserve_attrs,
+ partition=partition,
)
return B
-def maximum_spanning_arborescence(G, attr="weight", default=1, preserve_attrs=False):
+def maximum_spanning_arborescence(
+ G, attr="weight", default=1, preserve_attrs=False, partition=None
+):
ed = Edmonds(G)
B = ed.find_optimum(
- attr, default, kind="max", style="arborescence", preserve_attrs=preserve_attrs
+ attr,
+ default,
+ kind="max",
+ style="arborescence",
+ preserve_attrs=preserve_attrs,
+ partition=partition,
)
if not is_arborescence(B):
msg = "No maximum spanning arborescence in G."
@@ -665,10 +731,17 @@ def maximum_spanning_arborescence(G, attr="weight", default=1, preserve_attrs=Fa
return B
-def minimum_spanning_arborescence(G, attr="weight", default=1, preserve_attrs=False):
+def minimum_spanning_arborescence(
+ G, attr="weight", default=1, preserve_attrs=False, partition=None
+):
ed = Edmonds(G)
B = ed.find_optimum(
- attr, default, kind="min", style="arborescence", preserve_attrs=preserve_attrs
+ attr,
+ default,
+ kind="min",
+ style="arborescence",
+ preserve_attrs=preserve_attrs,
+ partition=partition,
)
if not is_arborescence(B):
msg = "No minimum spanning arborescence in G."
@@ -691,6 +764,10 @@ default : float
preserve_attrs : bool
If True, preserve the other attributes of the original graph (that are not
passed to `attr`)
+partition : str
+ The key for the edge attribute containing the partition
+ data on the graph. Edges can be included, excluded or open using the
+ `EdgePartition` enum.
Returns
-------
@@ -724,3 +801,223 @@ maximum_spanning_arborescence.__doc__ = docstring_arborescence.format(
minimum_spanning_arborescence.__doc__ = docstring_arborescence.format(
kind="minimum", style="spanning arborescence"
)
+
+
+class ArborescenceIterator:
+ """
+ Iterate over all spanning arborescences of a graph in either increasing or
+ decreasing cost.
+
+ Notes
+ -----
+ This iterator uses the partition scheme from [1]_ (included edges,
+ excluded edges and open edges). It generates minimum spanning
+ arborescences using a modified Edmonds' Algorithm which respects the
+ partition of edges. For arborescences with the same weight, ties are
+ broken arbitrarily.
+
+ References
+ ----------
+ .. [1] G.K. Janssens, K. Sörensen, An algorithm to generate all spanning
+ trees in order of increasing cost, Pesquisa Operacional, 2005-08,
+ Vol. 25 (2), p. 219-229,
+ https://www.scielo.br/j/pope/a/XHswBwRwJyrfL88dmMwYNWp/?lang=en
+ """
+
+ @dataclass(order=True)
+ class Partition:
+ """
+ This dataclass represents a partition and stores a dict with the edge
+ data and the weight of the minimum spanning arborescence of the
+ partition dict.
+ """
+
+ mst_weight: float
+ partition_dict: dict = field(compare=False)
+
+ def __copy__(self):
+ return ArborescenceIterator.Partition(
+ self.mst_weight, self.partition_dict.copy()
+ )
+
+ def __init__(self, G, weight="weight", minimum=True, init_partition=None):
+ """
+ Initialize the iterator
+
+ Parameters
+ ----------
+ G : nx.DiGraph
+ The directed graph which we need to iterate trees over
+
+ weight : String, default = "weight"
+ The edge attribute used to store the weight of the edge
+
+ minimum : bool, default = True
+ Return the trees in increasing order while true and decreasing order
+ while false.
+
+ init_partition : tuple, default = None
+ In the case that certain edges have to be included or excluded from
+ the arborescences, `init_partition` should be in the form
+ `(included_edges, excluded_edges)` where each edges is a
+ `(u, v)`-tuple inside an iterable such as a list or set.
+
+ """
+ self.G = G.copy()
+ self.weight = weight
+ self.minimum = minimum
+ self.method = (
+ minimum_spanning_arborescence if minimum else maximum_spanning_arborescence
+ )
+ # Randomly create a key for an edge attribute to hold the partition data
+ self.partition_key = (
+ "ArborescenceIterators super secret partition attribute name"
+ )
+ if init_partition is not None:
+ partition_dict = {}
+ for e in init_partition[0]:
+ partition_dict[e] = nx.EdgePartition.INCLUDED
+ for e in init_partition[1]:
+ partition_dict[e] = nx.EdgePartition.EXCLUDED
+ self.init_partition = ArborescenceIterator.Partition(0, partition_dict)
+ else:
+ self.init_partition = None
+
+ def __iter__(self):
+ """
+ Returns
+ -------
+ ArborescenceIterator
+ The iterator object for this graph
+ """
+ self.partition_queue = PriorityQueue()
+ self._clear_partition(self.G)
+
+ # Write the initial partition if it exists.
+ if self.init_partition is not None:
+ self._write_partition(self.init_partition)
+
+ mst_weight = self.method(
+ self.G,
+ self.weight,
+ partition=self.partition_key,
+ preserve_attrs=True,
+ ).size(weight=self.weight)
+
+ self.partition_queue.put(
+ self.Partition(
+ mst_weight if self.minimum else -mst_weight,
+ dict()
+ if self.init_partition is None
+ else self.init_partition.partition_dict,
+ )
+ )
+
+ return self
+
+ def __next__(self):
+ """
+ Returns
+ -------
+ (multi)Graph
+ The spanning tree of next greatest weight, which ties broken
+ arbitrarily.
+ """
+ if self.partition_queue.empty():
+ del self.G, self.partition_queue
+ raise StopIteration
+
+ partition = self.partition_queue.get()
+ self._write_partition(partition)
+ next_arborescence = self.method(
+ self.G,
+ self.weight,
+ partition=self.partition_key,
+ preserve_attrs=True,
+ )
+ self._partition(partition, next_arborescence)
+
+ self._clear_partition(next_arborescence)
+ return next_arborescence
+
+ def _partition(self, partition, partition_arborescence):
+ """
+ Create new partitions based of the minimum spanning tree of the
+ current minimum partition.
+
+ Parameters
+ ----------
+ partition : Partition
+ The Partition instance used to generate the current minimum spanning
+ tree.
+ partition_arborescence : nx.Graph
+ The minimum spanning arborescence of the input partition.
+ """
+ # create two new partitions with the data from the input partition dict
+ p1 = self.Partition(0, partition.partition_dict.copy())
+ p2 = self.Partition(0, partition.partition_dict.copy())
+ for e in partition_arborescence.edges:
+ # determine if the edge was open or included
+ if e not in partition.partition_dict:
+ # This is an open edge
+ p1.partition_dict[e] = nx.EdgePartition.EXCLUDED
+ p2.partition_dict[e] = nx.EdgePartition.INCLUDED
+
+ self._write_partition(p1)
+ try:
+ p1_mst = self.method(
+ self.G,
+ self.weight,
+ partition=self.partition_key,
+ preserve_attrs=True,
+ )
+
+ p1_mst_weight = p1_mst.size(weight=self.weight)
+ p1.mst_weight = p1_mst_weight if self.minimum else -p1_mst_weight
+ self.partition_queue.put(p1.__copy__())
+ except nx.NetworkXException:
+ pass
+
+ p1.partition_dict = p2.partition_dict.copy()
+
+ def _write_partition(self, partition):
+ """
+ Writes the desired partition into the graph to calculate the minimum
+ spanning tree. Also, if one incoming edge is included, mark all others
+ as excluded so that if that vertex is merged during Edmonds' algorithm
+ we cannot still pick another of that vertex's included edges.
+
+ Parameters
+ ----------
+ partition : Partition
+ A Partition dataclass describing a partition on the edges of the
+ graph.
+ """
+ for u, v, d in self.G.edges(data=True):
+ if (u, v) in partition.partition_dict:
+ d[self.partition_key] = partition.partition_dict[(u, v)]
+ else:
+ d[self.partition_key] = nx.EdgePartition.OPEN
+
+ for n in self.G:
+ included_count = 0
+ excluded_count = 0
+ for u, v, d in self.G.in_edges(nbunch=n, data=True):
+ if d.get(self.partition_key) == nx.EdgePartition.INCLUDED:
+ included_count += 1
+ elif d.get(self.partition_key) == nx.EdgePartition.EXCLUDED:
+ excluded_count += 1
+ # Check that if there is an included edges, all other incoming ones
+ # are excluded. If not fix it!
+ if included_count == 1 and excluded_count != self.G.in_degree(n) - 1:
+ for u, v, d in self.G.in_edges(nbunch=n, data=True):
+ if d.get(self.partition_key) != nx.EdgePartition.INCLUDED:
+ d[self.partition_key] = nx.EdgePartition.EXCLUDED
+
+ def _clear_partition(self, G):
+ """
+ Removes partition data from the graph
+ """
+ for u, v, d in G.edges(data=True):
+ if self.partition_key in d:
+ del d[self.partition_key]
diff --git a/networkx/algorithms/tree/mst.py b/networkx/algorithms/tree/mst.py
index bad8c132..5bf193fc 100644
--- a/networkx/algorithms/tree/mst.py
+++ b/networkx/algorithms/tree/mst.py
@@ -2,10 +2,13 @@
Algorithms for calculating min/max spanning trees/forests.
"""
+from dataclasses import dataclass, field
+from enum import Enum
from heapq import heappop, heappush
from operator import itemgetter
from itertools import count
from math import isnan
+from queue import PriorityQueue
import networkx as nx
from networkx.utils import UnionFind, not_implemented_for
@@ -15,9 +18,27 @@ __all__ = [
"maximum_spanning_edges",
"minimum_spanning_tree",
"maximum_spanning_tree",
+ "partition_spanning_tree",
+ "EdgePartition",
+ "SpanningTreeIterator",
]
+class EdgePartition(Enum):
+ """
+ An enum to store the state of an edge partition. The enum is written to the
+ edges of a graph before being pasted to `kruskal_mst_edges`. Options are:
+
+ - EdgePartition.OPEN
+ - EdgePartition.INCLUDED
+ - EdgePartition.EXCLUDED
+ """
+
+ OPEN = 0
+ INCLUDED = 1
+ EXCLUDED = 2
+
+
@not_implemented_for("multigraph")
def boruvka_mst_edges(
G, minimum=True, weight="weight", keys=False, data=True, ignore_nan=False
@@ -116,9 +137,10 @@ def boruvka_mst_edges(
def kruskal_mst_edges(
- G, minimum, weight="weight", keys=True, data=True, ignore_nan=False
+ G, minimum, weight="weight", keys=True, data=True, ignore_nan=False, partition=None
):
- """Iterate over edges of a Kruskal's algorithm min/max spanning tree.
+ """
+ Iterate over edge of a Kruskal's algorithm min/max spanning tree.
Parameters
----------
@@ -144,40 +166,64 @@ def kruskal_mst_edges(
If a NaN is found as an edge weight normally an exception is raised.
If `ignore_nan is True` then that edge is ignored instead.
+ partition : string (default: None)
+ The name of the edge attribute holding the partition data, if it exists.
+ Partition data is written to the edges using the `EdgePartition` enum.
+ If a partition exists, all included edges and none of the excluded edges
+ will appear in the final tree. Open edges may or may not be used.
+
+ Yields
+ ------
+ edge tuple
+ The edges as discovered by Kruskal's method. Each edge can
+ take the following forms: `(u, v)`, `(u, v, d)` or `(u, v, k, d)`
+ depending on the `key` and `data` parameters
"""
subtrees = UnionFind()
if G.is_multigraph():
edges = G.edges(keys=True, data=True)
+ else:
+ edges = G.edges(data=True)
- def filter_nan_edges(edges=edges, weight=weight):
- sign = 1 if minimum else -1
- for u, v, k, d in edges:
- wt = d.get(weight, 1) * sign
- if isnan(wt):
- if ignore_nan:
- continue
- msg = f"NaN found as an edge weight. Edge {(u, v, k, d)}"
- raise ValueError(msg)
- yield wt, u, v, k, d
+ """
+ Sort the edges of the graph with respect to the partition data.
+ Edges are returned in the following order:
+
+ * Included edges
+ * Open edges from smallest to largest weight
+ * Excluded edges
+ """
+ included_edges = []
+ open_edges = []
+ for e in edges:
+ d = e[-1]
+ wt = d.get(weight, 1)
+ if isnan(wt):
+ if ignore_nan:
+ continue
+ raise ValueError(f"NaN found as an edge weight. Edge {e}")
+
+ edge = (wt,) + e
+ if d.get(partition) == EdgePartition.INCLUDED:
+ included_edges.append(edge)
+ elif d.get(partition) == EdgePartition.EXCLUDED:
+ continue
+ else:
+ open_edges.append(edge)
+ if minimum:
+ sorted_open_edges = sorted(open_edges, key=itemgetter(0))
else:
- edges = G.edges(data=True)
+ sorted_open_edges = sorted(open_edges, key=itemgetter(0), reverse=True)
- def filter_nan_edges(edges=edges, weight=weight):
- sign = 1 if minimum else -1
- for u, v, d in edges:
- wt = d.get(weight, 1) * sign
- if isnan(wt):
- if ignore_nan:
- continue
- msg = f"NaN found as an edge weight. Edge {(u, v, d)}"
- raise ValueError(msg)
- yield wt, u, v, d
+ # Condense the lists into one
+ included_edges.extend(sorted_open_edges)
+ sorted_edges = included_edges
+ del open_edges, sorted_open_edges, included_edges
- edges = sorted(filter_nan_edges(), key=itemgetter(0))
# Multigraphs need to handle edge keys in addition to edge data.
if G.is_multigraph():
- for wt, u, v, k, d in edges:
+ for wt, u, v, k, d in sorted_edges:
if subtrees[u] != subtrees[v]:
if keys:
if data:
@@ -191,12 +237,12 @@ def kruskal_mst_edges(
yield u, v
subtrees.union(u, v)
else:
- for wt, u, v, d in edges:
+ for wt, u, v, d in sorted_edges:
if subtrees[u] != subtrees[v]:
if data:
- yield (u, v, d)
+ yield u, v, d
else:
- yield (u, v)
+ yield u, v
subtrees.union(u, v)
@@ -550,6 +596,69 @@ def minimum_spanning_tree(G, weight="weight", algorithm="kruskal", ignore_nan=Fa
return T
+def partition_spanning_tree(
+ G, minimum=True, weight="weight", partition="partition", ignore_nan=False
+):
+ """
+ Find a spanning tree while respecting a partition of edges.
+
+ Edges can be flagged as either `INLCUDED` which are required to be in the
+ returned tree, `EXCLUDED`, which cannot be in the returned tree and `OPEN`.
+
+ This is used in the SpanningTreeIterator to create new partitions following
+ the algorithm of Sörensen and Janssens [1]_.
+
+ Parameters
+ ----------
+ G : undirected graph
+ An undirected graph.
+
+ minimum : bool (default: True)
+ Determines whether the returned tree is the minimum spanning tree of
+ the partition of the maximum one.
+
+ weight : str
+ Data key to use for edge weights.
+
+ partition : str
+ The key for the edge attribute containing the partition
+ data on the graph. Edges can be included, excluded or open using the
+ `EdgePartition` enum.
+
+ ignore_nan : bool (default: False)
+ If a NaN is found as an edge weight normally an exception is raised.
+ If `ignore_nan is True` then that edge is ignored instead.
+
+
+ Returns
+ -------
+ G : NetworkX Graph
+ A minimum spanning tree using all of the included edges in the graph and
+ none of the excluded edges.
+
+ References
+ ----------
+ .. [1] G.K. Janssens, K. Sörensen, An algorithm to generate all spanning
+ trees in order of increasing cost, Pesquisa Operacional, 2005-08,
+ Vol. 25 (2), p. 219-229,
+ https://www.scielo.br/j/pope/a/XHswBwRwJyrfL88dmMwYNWp/?lang=en
+ """
+ edges = kruskal_mst_edges(
+ G,
+ minimum,
+ weight,
+ keys=True,
+ data=True,
+ ignore_nan=ignore_nan,
+ partition=partition,
+ )
+ T = G.__class__() # Same graph class as G
+ T.graph.update(G.graph)
+ T.add_nodes_from(G.nodes.items())
+ T.add_edges_from(edges)
+ return T
+
+
def maximum_spanning_tree(G, weight="weight", algorithm="kruskal", ignore_nan=False):
"""Returns a maximum spanning tree or forest on an undirected graph `G`.
@@ -610,3 +719,171 @@ def maximum_spanning_tree(G, weight="weight", algorithm="kruskal", ignore_nan=Fa
T.add_nodes_from(G.nodes.items())
T.add_edges_from(edges)
return T
+
+
+class SpanningTreeIterator:
+ """
+ Iterate over all spanning trees of a graph in either increasing or
+ decreasing cost.
+
+ Notes
+ -----
+ This iterator uses the partition scheme from [1]_ (included edges,
+ excluded edges and open edges) as well as a modified Kruskal's Algorithm
+ to generate minimum spanning trees which respect the partition of edges.
+ For spanning trees with the same weight, ties are broken arbitrarily.
+
+ References
+ ----------
+ .. [1] G.K. Janssens, K. Sörensen, An algorithm to generate all spanning
+ trees in order of increasing cost, Pesquisa Operacional, 2005-08,
+ Vol. 25 (2), p. 219-229,
+ https://www.scielo.br/j/pope/a/XHswBwRwJyrfL88dmMwYNWp/?lang=en
+ """
+
+ @dataclass(order=True)
+ class Partition:
+ """
+ This dataclass represents a partition and stores a dict with the edge
+ data and the weight of the minimum spanning tree of the partition dict.
+ """
+
+ mst_weight: float
+ partition_dict: dict = field(compare=False)
+
+ def __copy__(self):
+ return SpanningTreeIterator.Partition(
+ self.mst_weight, self.partition_dict.copy()
+ )
+
+ def __init__(self, G, weight="weight", minimum=True, ignore_nan=False):
+ """
+ Initialize the iterator
+
+ Parameters
+ ----------
+ G : nx.Graph
+ The directed graph which we need to iterate trees over
+
+ weight : String, default = "weight"
+ The edge attribute used to store the weight of the edge
+
+ minimum : bool, default = True
+ Return the trees in increasing order while true and decreasing order
+ while false.
+
+ ignore_nan : bool, default = False
+ If a NaN is found as an edge weight normally an exception is raised.
+ If `ignore_nan is True` then that edge is ignored instead.
+ """
+ self.G = G.copy()
+ self.weight = weight
+ self.minimum = minimum
+ self.ignore_nan = ignore_nan
+ # Randomly create a key for an edge attribute to hold the partition data
+ self.partition_key = (
+ "SpanningTreeIterators super secret partition attribute name"
+ )
+
+ def __iter__(self):
+ """
+ Returns
+ -------
+ SpanningTreeIterator
+ The iterator object for this graph
+ """
+ self.partition_queue = PriorityQueue()
+ self._clear_partition(self.G)
+ mst_weight = partition_spanning_tree(
+ self.G, self.minimum, self.weight, self.partition_key, self.ignore_nan
+ ).size(weight=self.weight)
+
+ self.partition_queue.put(
+ self.Partition(mst_weight if self.minimum else -mst_weight, dict())
+ )
+
+ return self
+
+ def __next__(self):
+ """
+ Returns
+ -------
+ (multi)Graph
+ The spanning tree of next greatest weight, which ties broken
+ arbitrarily.
+ """
+ if self.partition_queue.empty():
+ del self.G, self.partition_queue
+ raise StopIteration
+
+ partition = self.partition_queue.get()
+ self._write_partition(partition)
+ next_tree = partition_spanning_tree(
+ self.G, self.minimum, self.weight, self.partition_key, self.ignore_nan
+ )
+ self._partition(partition, next_tree)
+
+ self._clear_partition(next_tree)
+ return next_tree
+
+ def _partition(self, partition, partition_tree):
+ """
+ Create new partitions based of the minimum spanning tree of the
+ current minimum partition.
+
+ Parameters
+ ----------
+ partition : Partition
+ The Partition instance used to generate the current minimum spanning
+ tree.
+ partition_tree : nx.Graph
+ The minimum spanning tree of the input partition.
+ """
+ # create two new partitions with the data from the input partition dict
+ p1 = self.Partition(0, partition.partition_dict.copy())
+ p2 = self.Partition(0, partition.partition_dict.copy())
+ for e in partition_tree.edges:
+ # determine if the edge was open or included
+ if e not in partition.partition_dict:
+ # This is an open edge
+ p1.partition_dict[e] = EdgePartition.EXCLUDED
+ p2.partition_dict[e] = EdgePartition.INCLUDED
+
+ self._write_partition(p1)
+ p1_mst = partition_spanning_tree(
+ self.G,
+ self.minimum,
+ self.weight,
+ self.partition_key,
+ self.ignore_nan,
+ )
+ p1_mst_weight = p1_mst.size(weight=self.weight)
+ if nx.is_connected(p1_mst):
+ p1.mst_weight = p1_mst_weight if self.minimum else -p1_mst_weight
+ self.partition_queue.put(p1.__copy__())
+ p1.partition_dict = p2.partition_dict.copy()
+
+ def _write_partition(self, partition):
+ """
+ Writes the desired partition into the graph to calculate the minimum
+ spanning tree.
+
+ Parameters
+ ----------
+ partition : Partition
+ A Partition dataclass describing a partition on the edges of the
+ graph.
+ """
+ for u, v, d in self.G.edges(data=True):
+ if (u, v) in partition.partition_dict:
+ d[self.partition_key] = partition.partition_dict[(u, v)]
+ else:
+ d[self.partition_key] = EdgePartition.OPEN
+
+ def _clear_partition(self, G):
+ """
+ Removes partition data from the graph
+ """
+ for u, v, d in G.edges(data=True):
+ if self.partition_key in d:
+ del d[self.partition_key]
diff --git a/networkx/algorithms/tree/tests/test_branchings.py b/networkx/algorithms/tree/tests/test_branchings.py
index 7f5e9ccd..4da563e5 100644
--- a/networkx/algorithms/tree/tests/test_branchings.py
+++ b/networkx/algorithms/tree/tests/test_branchings.py
@@ -1,10 +1,11 @@
+import math
+
import pytest
np = pytest.importorskip("numpy")
import networkx as nx
-
from networkx.algorithms.tree import branchings
from networkx.algorithms.tree import recognition
@@ -17,16 +18,18 @@ from networkx.algorithms.tree import recognition
# fmt: off
G_array = np.array([
# 0 1 2 3 4 5 6 7 8
- [0, 0, 12, 0, 12, 0, 0, 0, 0], # 0
- [4, 0, 0, 0, 0, 13, 0, 0, 0], # 1
- [0, 17, 0, 21, 0, 12, 0, 0, 0], # 2
- [5, 0, 0, 0, 17, 0, 18, 0, 0], # 3
- [0, 0, 0, 0, 0, 0, 0, 12, 0], # 4
- [0, 0, 0, 0, 0, 0, 14, 0, 12], # 5
- [0, 0, 21, 0, 0, 0, 0, 0, 15], # 6
- [0, 0, 0, 19, 0, 0, 15, 0, 0], # 7
- [0, 0, 0, 0, 0, 0, 0, 18, 0], # 8
+ [0, 0, 12, 0, 12, 0, 0, 0, 0], # 0
+ [4, 0, 0, 0, 0, 13, 0, 0, 0], # 1
+ [0, 17, 0, 21, 0, 12, 0, 0, 0], # 2
+ [5, 0, 0, 0, 17, 0, 18, 0, 0], # 3
+ [0, 0, 0, 0, 0, 0, 0, 12, 0], # 4
+ [0, 0, 0, 0, 0, 0, 14, 0, 12], # 5
+ [0, 0, 21, 0, 0, 0, 0, 0, 15], # 6
+ [0, 0, 0, 19, 0, 0, 15, 0, 0], # 7
+ [0, 0, 0, 0, 0, 0, 0, 18, 0], # 8
], dtype=int)
+
+
# fmt: on
@@ -143,11 +146,6 @@ def assert_equal_branchings(G1, G2, attr="weight", default=1):
e1 = sorted_edges(G1, attr, default)
e2 = sorted_edges(G2, attr, default)
- # If we have an exception, let's see the edges.
- print(e1)
- print(e2)
- print
-
for a, b in zip(e1, e2):
assert a[:2] == b[:2]
np.testing.assert_almost_equal(a[2], b[2])
@@ -413,7 +411,6 @@ def test_edge_attribute_preservation_normal_graph():
def test_edge_attribute_preservation_multigraph():
-
# Test that edge attributes are preserved when finding an optimum graph
# using the Edmonds class for multigraphs.
G = nx.MultiGraph()
@@ -449,3 +446,113 @@ def test_edge_attribute_discard():
edge_dict = B[0][1]
with pytest.raises(KeyError):
_ = edge_dict["otherattr"]
+
+
+def test_partition_spanning_arborescence():
+ """
+ Test that we can generate minimum spanning arborescences which respect the
+ given partition.
+ """
+ G = nx.from_numpy_array(G_array, create_using=nx.DiGraph)
+ G[3][0]["partition"] = nx.EdgePartition.EXCLUDED
+ G[2][3]["partition"] = nx.EdgePartition.INCLUDED
+ G[7][3]["partition"] = nx.EdgePartition.EXCLUDED
+ G[0][2]["partition"] = nx.EdgePartition.EXCLUDED
+ G[6][2]["partition"] = nx.EdgePartition.INCLUDED
+
+ actual_edges = [
+ (0, 4, 12),
+ (1, 0, 4),
+ (1, 5, 13),
+ (2, 3, 21),
+ (4, 7, 12),
+ (5, 6, 14),
+ (5, 8, 12),
+ (6, 2, 21),
+ ]
+
+ B = branchings.minimum_spanning_arborescence(G, partition="partition")
+ assert_equal_branchings(build_branching(actual_edges), B)
+
+
+def test_arborescence_iterator_min():
+ """
+ Tests the arborescence iterator.
+
+ A brute force method found 680 arboresecences in this graph.
+ This test will not verify all of them individually, but will check two
+ things
+
+ * The iterator returns 680 arboresecences
+ * The weight of the arborescences is non-strictly increasing
+
+ for more information please visit
+ https://mjschwenne.github.io/2021/06/10/implementing-the-iterators.html
+ """
+ G = nx.from_numpy_array(G_array, create_using=nx.DiGraph)
+
+ arborescence_count = 0
+ arborescence_weight = -math.inf
+ for B in branchings.ArborescenceIterator(G):
+ arborescence_count += 1
+ new_arborescence_weight = B.size(weight="weight")
+ assert new_arborescence_weight >= arborescence_weight
+ arborescence_weight = new_arborescence_weight
+
+ assert arborescence_count == 680
+
+
+def test_arborescence_iterator_max():
+ """
+ Tests the arborescence iterator.
+
+ A brute force method found 680 arboresecences in this graph.
+ This test will not verify all of them individually, but will check two
+ things
+
+ * The iterator returns 680 arboresecences
+ * The weight of the arborescences is non-strictly decreasing
+
+ for more information please visit
+ https://mjschwenne.github.io/2021/06/10/implementing-the-iterators.html
+ """
+ G = nx.from_numpy_array(G_array, create_using=nx.DiGraph)
+
+ arborescence_count = 0
+ arborescence_weight = math.inf
+ for B in branchings.ArborescenceIterator(G, minimum=False):
+ arborescence_count += 1
+ new_arborescence_weight = B.size(weight="weight")
+ assert new_arborescence_weight <= arborescence_weight
+ arborescence_weight = new_arborescence_weight
+
+ assert arborescence_count == 680
+
+
+def test_arborescence_iterator_initial_partition():
+ """
+ Tests the arborescence iterator with three included edges and three excluded
+ in the initial partition.
+
+ A brute force method similar to the one used in the above tests found that
+ there are 16 arborescences which contain the included edges and not the
+ excluded edges.
+ """
+ G = nx.from_numpy_array(G_array, create_using=nx.DiGraph)
+ included_edges = [(1, 0), (5, 6), (8, 7)]
+ excluded_edges = [(0, 2), (3, 6), (1, 5)]
+
+ arborescence_count = 0
+ arborescence_weight = -math.inf
+ for B in branchings.ArborescenceIterator(
+ G, init_partition=(included_edges, excluded_edges)
+ ):
+ arborescence_count += 1
+ new_arborescence_weight = B.size(weight="weight")
+ assert new_arborescence_weight >= arborescence_weight
+ arborescence_weight = new_arborescence_weight
+ for e in included_edges:
+ assert e in B.edges
+ for e in excluded_edges:
+ assert e not in B.edges
+ assert arborescence_count == 16
diff --git a/networkx/algorithms/tree/tests/test_mst.py b/networkx/algorithms/tree/tests/test_mst.py
index 5253a29f..6eebe0cb 100644
--- a/networkx/algorithms/tree/tests/test_mst.py
+++ b/networkx/algorithms/tree/tests/test_mst.py
@@ -13,20 +13,17 @@ def test_unknown_algorithm():
class MinimumSpanningTreeTestBase:
"""Base class for test classes for minimum spanning tree algorithms.
-
This class contains some common tests that will be inherited by
subclasses. Each subclass must have a class attribute
:data:`algorithm` that is a string representing the algorithm to
run, as described under the ``algorithm`` keyword argument for the
:func:`networkx.minimum_spanning_edges` function. Subclasses can
then implement any algorithm-specific tests.
-
"""
def setup_method(self, method):
"""Creates an example graph and stores the expected minimum and
maximum spanning tree edges.
-
"""
# This stores the class attribute `algorithm` in an instance attribute.
self.algo = self.algorithm
@@ -208,7 +205,6 @@ class MinimumSpanningTreeTestBase:
class TestBoruvka(MinimumSpanningTreeTestBase):
"""Unit tests for computing a minimum (or maximum) spanning tree
using Borůvka's algorithm.
-
"""
algorithm = "boruvka"
@@ -216,7 +212,6 @@ class TestBoruvka(MinimumSpanningTreeTestBase):
def test_unicode_name(self):
"""Tests that using a Unicode string can correctly indicate
Borůvka's algorithm.
-
"""
edges = nx.minimum_spanning_edges(self.G, algorithm="borůvka")
# Edges from the spanning edges functions don't come in sorted
@@ -231,7 +226,6 @@ class MultigraphMSTTestBase(MinimumSpanningTreeTestBase):
def test_multigraph_keys_min(self):
"""Tests that the minimum spanning edges of a multigraph
preserves edge keys.
-
"""
G = nx.MultiGraph()
G.add_edge(0, 1, key="a", weight=2)
@@ -243,7 +237,6 @@ class MultigraphMSTTestBase(MinimumSpanningTreeTestBase):
def test_multigraph_keys_max(self):
"""Tests that the maximum spanning edges of a multigraph
preserves edge keys.
-
"""
G = nx.MultiGraph()
G.add_edge(0, 1, key="a", weight=2)
@@ -256,7 +249,6 @@ class MultigraphMSTTestBase(MinimumSpanningTreeTestBase):
class TestKruskal(MultigraphMSTTestBase):
"""Unit tests for computing a minimum (or maximum) spanning tree
using Kruskal's algorithm.
-
"""
algorithm = "kruskal"
@@ -265,7 +257,6 @@ class TestKruskal(MultigraphMSTTestBase):
class TestPrim(MultigraphMSTTestBase):
"""Unit tests for computing a minimum (or maximum) spanning tree
using Prim's algorithm.
-
"""
algorithm = "prim"
@@ -283,3 +274,96 @@ class TestPrim(MultigraphMSTTestBase):
G.add_edge(0, 1, key="b", weight=1)
T = nx.maximum_spanning_tree(G)
assert edges_equal([(0, 1, 2)], list(T.edges(data="weight")))
+
+
+class TestSpanningTreeIterator:
+ """
+ Tests the spanning tree iterator on the example graph in the 2005 Sörensen
+ and Janssens paper An Algorithm to Generate all Spanning Trees of a Graph in
+ Order of Increasing Cost
+ """
+
+ def setup(self):
+ # Original Graph
+ edges = [(0, 1, 5), (1, 2, 4), (1, 4, 6), (2, 3, 5), (2, 4, 7), (3, 4, 3)]
+ self.G = nx.Graph()
+ self.G.add_weighted_edges_from(edges)
+ # List of lists of spanning trees in increasing order
+ self.spanning_trees = [
+ # 1, MST, cost = 17
+ [
+ (0, 1, {"weight": 5}),
+ (1, 2, {"weight": 4}),
+ (2, 3, {"weight": 5}),
+ (3, 4, {"weight": 3}),
+ ],
+ # 2, cost = 18
+ [
+ (0, 1, {"weight": 5}),
+ (1, 2, {"weight": 4}),
+ (1, 4, {"weight": 6}),
+ (3, 4, {"weight": 3}),
+ ],
+ # 3, cost = 19
+ [
+ (0, 1, {"weight": 5}),
+ (1, 4, {"weight": 6}),
+ (2, 3, {"weight": 5}),
+ (3, 4, {"weight": 3}),
+ ],
+ # 4, cost = 19
+ [
+ (0, 1, {"weight": 5}),
+ (1, 2, {"weight": 4}),
+ (2, 4, {"weight": 7}),
+ (3, 4, {"weight": 3}),
+ ],
+ # 5, cost = 20
+ [
+ (0, 1, {"weight": 5}),
+ (1, 2, {"weight": 4}),
+ (1, 4, {"weight": 6}),
+ (2, 3, {"weight": 5}),
+ ],
+ # 6, cost = 21
+ [
+ (0, 1, {"weight": 5}),
+ (1, 4, {"weight": 6}),
+ (2, 4, {"weight": 7}),
+ (3, 4, {"weight": 3}),
+ ],
+ # 7, cost = 21
+ [
+ (0, 1, {"weight": 5}),
+ (1, 2, {"weight": 4}),
+ (2, 3, {"weight": 5}),
+ (2, 4, {"weight": 7}),
+ ],
+ # 8, cost = 23
+ [
+ (0, 1, {"weight": 5}),
+ (1, 4, {"weight": 6}),
+ (2, 3, {"weight": 5}),
+ (2, 4, {"weight": 7}),
+ ],
+ ]
+
+ def test_minimum_spanning_tree_iterator(self):
+ """
+ Tests that the spanning trees are correctly returned in increasing order
+ """
+ tree_index = 0
+ for tree in nx.SpanningTreeIterator(self.G):
+ actual = sorted(tree.edges(data=True))
+ assert edges_equal(actual, self.spanning_trees[tree_index])
+ tree_index += 1
+
+ def test_maximum_spanning_tree_iterator(self):
+ """
+ Tests that the spanning trees are correctly returned in decreasing order
+ """
+ tree_index = 7
+ for tree in nx.SpanningTreeIterator(self.G, minimum=False):
+ actual = sorted(tree.edges(data=True))
+ assert edges_equal(actual, self.spanning_trees[tree_index])
+ tree_index -= 1
View Lorry Controller Status