diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml index 1168bd9..f1df83a 100644 --- a/.github/workflows/test.yml +++ b/.github/workflows/test.yml @@ -26,7 +26,7 @@ jobs: - name: Install dependencies run: | python -m pip install --upgrade pip - pip install flake8 pytest + pip install flake8 pytest pytest-benchmark if [ -f requirements.txt ]; then pip install -r requirements.txt; fi - name: Lint with flake8 run: | diff --git a/examples/g2o_experiment.py b/examples/g2o_experiment.py index 81c65b8..9b9c240 100644 --- a/examples/g2o_experiment.py +++ b/examples/g2o_experiment.py @@ -3,14 +3,12 @@ import numpy as np import networkx as nx from timeit import default_timer as timer -from pose_graph_utils import read_g2o_file, plot_poses, rpm_to_mac, RelativePoseMeasurement, poses_ate_tran, poses_rpe_rot +from pose_graph_utils import split_edges, read_g2o_file, plot_poses, rpm_to_mac, RelativePoseMeasurement, poses_ate_tran, poses_rpe_rot # MAC requirements -from mac import MAC -from mac.baseline import NaiveGreedy -from mac.greedy_eig import GreedyEig -from mac.greedy_esp import GreedyESP -from mac.utils import split_edges, Edge, round_madow +from mac.solvers import MAC, NaiveGreedy +from mac.utils.graphs import Edge +from mac.utils.rounding import round_madow import matplotlib.pyplot as plt plt.rcParams['text.usetex'] = True @@ -267,14 +265,14 @@ def to_sesync_format(measurements): # pass # Make a MAC Solver - mac = MAC(odom_edges, lc_edges, num_poses, use_cache=True, fiedler_method="tracemin_cholesky") + mac = MAC(odom_edges, lc_edges, num_poses, fiedler_method="tracemin_cholesky") # Make a Naive Solver naive = NaiveGreedy(lc_edges) # Make a GreedyEig Solver if run_greedy: - # greedy_eig = GreedyEig(rpm_to_mac(odom_measurements), rpm_to_mac(lc_measurements), num_poses) + from mac.solvers.greedy_esp import GreedyESP greedy_esp = GreedyESP(odom_edges, lc_edges, num_poses, lazy=True) ############################# @@ -318,7 +316,7 @@ def to_sesync_format(measurements): # Solve the relaxed maximum algebraic connectivity augmentation problem. start = timer() - result, unrounded, upper, rtime = mac.fw_subset(w_init, num_lc, max_iters=20, rounding="nearest", return_rounding_time=True) + result, unrounded, upper, rtime = mac.solve(num_lc, w_init, max_iters=20, rounding="nearest", return_rounding_time=True, use_cache=True) end = timer() solve_time = end - start times.append(solve_time) @@ -337,14 +335,8 @@ def to_sesync_format(measurements): # point solution every time. madow_times.append(solve_time + (end - start) - rtime) - # Solve the relaxed maximum algebraic connectivity augmentation problem. + # Solve the greedy k-edge selection problem if run_greedy: - # start = timer() - # greedy_eig_result, _ = greedy_eig.subset(num_lc) - # end = timer() - # greedy_eig_times.append(end - start) - # greedy_eig_results.append(greedy_eig_result) - num_lcs = [int(pct_lc * len(lc_measurements)) for pct_lc in percent_lc] greedy_esp_results, _, greedy_esp_times = greedy_esp.subsets_lazy(num_lcs, verbose=True) pass diff --git a/examples/g2o_mac_cache.py b/examples/g2o_mac_cache.py deleted file mode 100644 index ea3653c..0000000 --- a/examples/g2o_mac_cache.py +++ /dev/null @@ -1,148 +0,0 @@ -import sys -import random -import numpy as np -import networkx as nx -from timeit import default_timer as timer -from pose_graph_utils import read_g2o_file, plot_poses, rpm_to_mac, RelativePoseMeasurement - -# MAC requirements -from mac import MAC -from mac.baseline import NaiveGreedy -from mac.utils import split_edges, Edge, round_madow - -import matplotlib.pyplot as plt -plt.rcParams['text.usetex'] = True -plt.rcParams.update({'font.size': 16}) - -if __name__ == '__main__': - if len(sys.argv) < 2: - print(f"Usage: {sys.argv[0]} [.g2o file] [optional: --run-greedy]") - sys.exit() - pass - - run_greedy = False - if len(sys.argv) > 2: - if sys.argv[2] == "--run-greedy": - run_greedy = True - pass - else: - print(f"Unknown argument: {sys.argv[2]}") - print(f"Usage: {sys.argv[0]} [.g2o file] [optional: --run-greedy]") - sys.exit() - pass - pass - - dataset_name = sys.argv[1].split('/')[-1].split('.')[0] - print(f"Loading dataset: {dataset_name}") - - # Load a g2o file - print("Reading g2o file") - start = timer() - measurements, num_poses = read_g2o_file(sys.argv[1]) - end = timer() - print("Success! elapsed time: ", (end - start)) - - # Split measurements into odom and loop closures - odom_measurements, lc_measurements = split_edges(measurements) - - # Convert measurements to MAC edge format - odom_edges = rpm_to_mac(odom_measurements) - lc_edges = rpm_to_mac(lc_measurements) - - # Print dataset stats - print(f"Loaded {len(measurements)} total measurements with: ") - print(f"\t {len(odom_measurements)} base (odometry) measurements and") - print(f"\t {len(lc_measurements)} candidate (loop closure) measurements") - - # Make a naive solver for use as an initialization - naive = NaiveGreedy(lc_edges) - - # Make a MAC Solver - mac = MAC(odom_edges, lc_edges, num_poses, use_cache=False) - mac_cache = MAC(odom_edges, lc_edges, num_poses, fiedler_method="tracemin_cholesky", use_cache=True) - - ############################# - # Running the tests! - ############################# - - # Test between 100% and 0% loop closures - # NOTE: If running greedy, these must be in increasing order! - percent_lc = [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0] - - # Containers for MAC results - results = [] - unrounded_results = [] - upper_bounds = [] - times = [] - - cache_results = [] - cache_unrounded_results = [] - cache_upper_bounds = [] - cache_times = [] - - for pct_lc in percent_lc: - num_lc = int(pct_lc * len(lc_measurements)) - print("Num LC to accept: ", num_lc) - - # Compute a solution using the naive method. This serves both as a - # baseline and as a sparse initializer for our method. - naive_result = naive.subset(num_lc) - - w_init = naive_result - - # Solve the relaxed maximum algebraic connectivity augmentation problem. - start = timer() - result, unrounded, upper = mac.fw_subset(w_init, num_lc, max_iters=20, rounding="nearest") - end = timer() - solve_time = end - start - times.append(solve_time) - results.append(result) - upper_bounds.append(upper) - unrounded_results.append(unrounded) - - start = timer() - result_cache, unrounded_cache, upper_cache = mac_cache.fw_subset(w_init, num_lc, max_iters=20, rounding="nearest", fallback=False) - end = timer() - solve_time = end - start - cache_times.append(solve_time) - cache_results.append(result_cache) - cache_upper_bounds.append(upper_cache) - cache_unrounded_results.append(unrounded_cache) - - # Verify that the results match - # assert np.allclose(result, result_cache), "Results should be the same" - # assert np.allclose(unrounded, unrounded_cache), "Unrounded results should be the same" - # assert np.allclose(upper, upper_cache), "Upper bounds should be the same" - pass - - # plot connectivity vs. percent_lc - our_objective_vals = [mac.evaluate_objective(result) for result in results] - unrounded_objective_vals = [mac.evaluate_objective(unrounded) for unrounded in unrounded_results] - cache_objective_vals = [mac.evaluate_objective(cache_result) for cache_result in cache_results] - cache_unrounded_objective_vals = [mac.evaluate_objective(cache_unrounded) for cache_unrounded in cache_unrounded_results] - - plt.plot(100.0*np.array(percent_lc), our_objective_vals, label='MAC (Ours) Nearest') - plt.plot(100.0*np.array(percent_lc), upper_bounds, label='Dual Upper Bound', linestyle='--', color='C0') - plt.fill_between(100.0*np.array(percent_lc), our_objective_vals, upper_bounds, alpha=0.2, label='Suboptimality Gap') - plt.plot(100.0*np.array(percent_lc), unrounded_objective_vals, label='Unrounded', c='C2') - - plt.plot(100.0*np.array(percent_lc), cache_objective_vals, label='MAC (Cache)', marker='x', color='C4') - plt.plot(100.0*np.array(percent_lc), cache_upper_bounds, label='Upper Bound (Cache)', linestyle='--', color='C4') - plt.fill_between(100.0*np.array(percent_lc), cache_objective_vals, cache_upper_bounds, alpha=0.1, label='Suboptimality Gap (Cache)', color='C5') - plt.plot(100.0*np.array(percent_lc), cache_unrounded_objective_vals, label='Unrounded (Cache)', c='C5') - - plt.ylabel(r'Algebraic Connectivity $\lambda_2$') - plt.xlabel(r'\% Edges Added') - plt.legend() - plt.savefig(f"alg_conn_{dataset_name}.png", dpi=600, bbox_inches='tight') - plt.show() - - plt.figure() - plt.plot(100.0*np.array(percent_lc), times, label='MAC (No Cache)') - plt.plot(100.0*np.array(percent_lc), cache_times, label='MAC (Cache)', marker='x', color='C4') - plt.xlim([0.0, 90.0]) - plt.ylabel(r'Time (s)') - plt.xlabel(r'\% Edges Added') - plt.savefig(f"comp_time_cache_compare_{dataset_name}.png", dpi=600, bbox_inches='tight') - plt.legend() - plt.show() diff --git a/examples/greedy_eig_test.py b/examples/greedy_eig_test.py deleted file mode 100644 index f9430a5..0000000 --- a/examples/greedy_eig_test.py +++ /dev/null @@ -1,61 +0,0 @@ -import numpy as np -import networkx as nx -import matplotlib.pyplot as plt -from mac.utils import select_edges, split_edges, Edge, nx_to_mac, mac_to_nx -from mac.baseline import NaiveGreedy -from mac.greedy_eig import GreedyEig -from mac.mac import MAC - -n = 20 -p = 0.99 -G = nx.erdos_renyi_graph(n, p) - -# Add a chain -for i in range(n-1): - if G.has_edge(i+1, i): - G.remove_edge(i+1, i) - if not G.has_edge(i, i+1): - G.add_edge(i, i+1) - -print(G) -nx.draw(G) -plt.show() - -# Ensure G is connected before proceeding -assert(nx.is_connected(G)) - -measurements = nx_to_mac(G) - -# Split chain and non-chain parts -fixed_meas, candidate_meas = split_edges(measurements) - -pct_candidates = 0.1 -num_candidates = int(pct_candidates * len(candidate_meas)) -mac = MAC(fixed_meas, candidate_meas, n) -naive_greedy_eig = GreedyEig(fixed_meas, candidate_meas, n) - -w_init = np.zeros(len(candidate_meas)) -w_init[:num_candidates] = 1.0 - -result, rounded, upper = mac.fw_subset(w_init, num_candidates, max_iters=200) -init_selected = select_edges(candidate_meas, w_init) -selected = select_edges(candidate_meas, result) - -init_selected_G = mac_to_nx(fixed_meas + init_selected) -selected_G = mac_to_nx(fixed_meas + selected) - -w_greedy_eig, selected_eig = naive_greedy_eig.subset(num_candidates) -selected_G_eig = mac_to_nx(fixed_meas + selected_eig) - -print(f"lambda2 Random: {mac.evaluate_objective(w_init)}") -print(f"lambda2 Greedy: {mac.evaluate_objective(w_greedy_eig)}") -print(f"lambda2 Ours: {mac.evaluate_objective(result)}") - -plt.subplot(131) -nx.draw(init_selected_G) -plt.subplot(132) -nx.draw(selected_G_eig) -plt.subplot(133) -nx.draw(selected_G) -plt.show() - diff --git a/examples/greedy_esp_example.py b/examples/greedy_esp_example.py deleted file mode 100644 index ba23e07..0000000 --- a/examples/greedy_esp_example.py +++ /dev/null @@ -1,50 +0,0 @@ -import networkx as nx -from mac.mac import MAC -from mac.greedy_esp import GreedyESP -from mac.utils import split_edges, nx_to_mac, mac_to_nx, get_edge_selection_as_binary_mask -import matplotlib.pyplot as plt -import time - -if __name__ == "__main__": - # set up graph in expected format - n = 200 - p = 0.30 - G = nx.erdos_renyi_graph(n, p, seed=42) - - # Add a chain - for i in range(n-1): - if G.has_edge(i+1, i): - G.remove_edge(i+1, i) - if not G.has_edge(i, i+1): - G.add_edge(i, i+1) - - edges = nx_to_mac(G) - fixed_edges, candidate_edges = split_edges(edges) - print(f"Graph has {len(edges)} edges") - print(f"Fixed edges: {len(fixed_edges)}, candidate edges: {len(candidate_edges)}") - mac = MAC(fixed_edges, candidate_edges, n) - - # build the problem - time_start = time.perf_counter() - edge_selection_problem = GreedyESP(fixed_edges, candidate_edges, n) - - # solve the problem - num_edges_to_select = 100 - selected_edges_mask, selected_edges = edge_selection_problem.subset(num_edges_to_select) - - print(f"Time taken: {time.perf_counter() - time_start}") - - # check that the solution is valid - assert len(selected_edges) == num_edges_to_select - assert len(set(selected_edges).intersection(set(fixed_edges))) == 0 - assert len(set(selected_edges).intersection(set(candidate_edges))) == num_edges_to_select - - # print the value of the objective function - print(f"lambda2 GreedyESP: {mac.evaluate_objective(selected_edges_mask)}") - - # plot the solution - all_edges_in_sparsified_graph = fixed_edges + selected_edges - selected_graph = mac_to_nx(all_edges_in_sparsified_graph) - plt.figure() - nx.draw(selected_graph, with_labels=True) - plt.show() diff --git a/examples/petersen_graph_sparsification.py b/examples/petersen_graph_sparsification.py index 57c80c8..cc56d94 100644 --- a/examples/petersen_graph_sparsification.py +++ b/examples/petersen_graph_sparsification.py @@ -1,62 +1,64 @@ import numpy as np import networkx as nx import matplotlib.pyplot as plt -from mac.utils import select_edges, split_edges, nx_to_mac, mac_to_nx -from mac.baseline import NaiveGreedy -from mac.greedy_eig import GreedyEig -from mac.greedy_esp import GreedyESP -from mac.greedy_esp_minimal import MinimalGreedyESP -from mac.mac import MAC +from mac.utils.graphs import select_edges +from mac.utils.conversions import nx_to_mac, mac_to_nx +from mac.solvers.greedy_eig import GreedyEig +from mac.solvers.greedy_esp import GreedyESP +from mac.solvers.mac import MAC plt.rcParams['text.usetex'] = True G = nx.petersen_graph() n = len(G.nodes()) -# Add a chain -for i in range(n-1): - if G.has_edge(i+1, i): - G.remove_edge(i+1, i) - if not G.has_edge(i, i+1): - G.add_edge(i, i+1) - print(G) -pos = nx.shell_layout(G, nlist=[range(5,10), range(5)]) +pos = nx.shell_layout(G, nlist=[range(5,10), range(5)], rotate=np.pi/32) nx.draw(G, pos=pos) plt.show() # Ensure G is connected before proceeding assert(nx.is_connected(G)) -edges = nx_to_mac(G) +# Split the graph into a tree part and a "loop" part +spanning_tree = nx.minimum_spanning_tree(G) +loop_graph = nx.difference(G, spanning_tree) + +nx.draw(spanning_tree, pos=pos) +plt.show() -# Split chain and non-chain parts -fixed_meas, candidate_meas = split_edges(edges) +fixed_edges = nx_to_mac(spanning_tree) +candidate_edges = nx_to_mac(loop_graph) -pct_candidates = 0.4 -num_candidates = int(pct_candidates * len(candidate_meas)) -mac = MAC(fixed_meas, candidate_meas, n) -greedy_eig = GreedyEig(fixed_meas, candidate_meas, n) -greedy_esp = GreedyESP(fixed_meas, candidate_meas, n) +pct_candidates = 0.6 +num_candidates = int(pct_candidates * len(candidate_edges)) +print(num_candidates) +mac = MAC(fixed_edges, candidate_edges, n) +greedy_eig = GreedyEig(fixed_edges, candidate_edges, n) +greedy_esp = GreedyESP(fixed_edges, candidate_edges, n) -w_init = np.zeros(len(candidate_meas)) -w_init[:num_candidates] = 1.0 +w_init = np.zeros(len(candidate_edges)) +# Sample num_candidates indices at random without replacement +rng = np.random.default_rng() +random_selection = rng.choice(len(candidate_edges), size=num_candidates, replace=False) +w_init[random_selection] = 1.0 +print(w_init) -result, unrounded, upper = mac.fw_subset(w_init, num_candidates, max_iters=100) +result, unrounded, upper = mac.solve(num_candidates, w_init, max_iters=100, rounding="nearest") greedy_eig_result, _ = greedy_eig.subset(num_candidates) greedy_esp_result, _ = greedy_esp.subset(num_candidates) -init_selected = select_edges(candidate_meas, w_init) -init_selected_G = mac_to_nx(fixed_meas + init_selected) +init_selected = select_edges(candidate_edges, w_init) +init_selected_G = mac_to_nx(fixed_edges + init_selected) -greedy_eig_selected = select_edges(candidate_meas, greedy_eig_result) -greedy_eig_selected_G = mac_to_nx(fixed_meas + greedy_eig_selected) +greedy_eig_selected = select_edges(candidate_edges, greedy_eig_result) +greedy_eig_selected_G = mac_to_nx(fixed_edges + greedy_eig_selected) -greedy_esp_selected = select_edges(candidate_meas, greedy_esp_result) -greedy_esp_selected_G = mac_to_nx(fixed_meas + greedy_esp_selected) +greedy_esp_selected = select_edges(candidate_edges, greedy_esp_result) +greedy_esp_selected_G = mac_to_nx(fixed_edges + greedy_esp_selected) -selected = select_edges(candidate_meas, result) -selected_G = mac_to_nx(fixed_meas + selected) +selected = select_edges(candidate_edges, result) +selected_G = mac_to_nx(fixed_edges + selected) print(f"lambda2 Random: {mac.evaluate_objective(w_init)}") print(f"lambda2 Ours: {mac.evaluate_objective(result)}") diff --git a/examples/pose_graph_utils.py b/examples/pose_graph_utils.py index 7aa8cb2..948277d 100644 --- a/examples/pose_graph_utils.py +++ b/examples/pose_graph_utils.py @@ -4,7 +4,7 @@ from scipy.sparse import csr_matrix, coo_matrix from typing import List, Tuple, Union import networkx as nx -from mac.utils import Edge +from mac.utils.graphs import Edge from evo.core.trajectory import PoseTrajectory3D from evo.core import sync @@ -15,6 +15,36 @@ RelativePoseMeasurement = namedtuple('RelativePoseMeasurement', ['i', 'j', 't', 'R', 'kappa', 'tau']) +def split_edges(edges: List[Edge]) -> Tuple[List[Edge], List[Edge]]: + """Splits list of edges into a "fixed" chain part and a set of candidate + loops. + + Args: + edges: List of edges. + + Returns: + A tuple containing two lists: + fixed: edges where |i - j| = 1, and + candidates: edges where |i - j| != 1 + + This is particularly useful for pose-graph SLAM applications where the + fixed edges correspond to an odometry chain and the candidate edges + correspond to loop closures. + + """ + chain_edges = [] + loop_edges = [] + for edge in edges: + id1 = edge.i + id2 = edge.j + if abs(id2 - id1) > 1: + loop_edges.append(edge) + else: + chain_edges.append(edge) + + return chain_edges, loop_edges + + def rot2D_from_theta(theta: float) -> np.ndarray: """ Simply builds a 2D rotation matrix: diff --git a/examples/random_graph_sparsification.py b/examples/random_graph_sparsification.py index fcac901..598e5fa 100644 --- a/examples/random_graph_sparsification.py +++ b/examples/random_graph_sparsification.py @@ -1,12 +1,12 @@ import numpy as np import networkx as nx import matplotlib.pyplot as plt -from mac.utils import select_edges, split_edges, Edge, nx_to_mac, mac_to_nx -from mac.baseline import NaiveGreedy -from mac.mac import MAC +from mac.utils.graphs import select_edges, Edge +from mac.utils.conversions import nx_to_mac, mac_to_nx +from mac.solvers import MAC n = 20 -p = 0.99 +p = 0.6 G = nx.erdos_renyi_graph(n, p) # Add a chain @@ -16,7 +16,10 @@ if not G.has_edge(i, i+1): G.add_edge(i, i+1) -print(G) +# Split the graph into a tree part and a "loop" part +spanning_tree = nx.minimum_spanning_tree(G) +loop_graph = nx.difference(G, spanning_tree) + nx.draw(G) plt.title("Original Graph") plt.show() @@ -24,31 +27,30 @@ # Ensure G is connected before proceeding assert(nx.is_connected(G)) -edges = nx_to_mac(G) - -# Split chain and non-chain parts -fixed_edges, candidate_edges = split_edges(edges) +fixed_edges = nx_to_mac(spanning_tree) +candidate_edges = nx_to_mac(loop_graph) -pct_candidates = 0.1 +pct_candidates = 0.2 num_candidates = int(pct_candidates * len(candidate_edges)) mac = MAC(fixed_edges, candidate_edges, n) w_init = np.zeros(len(candidate_edges)) w_init[:num_candidates] = 1.0 +np.random.shuffle(w_init) -result, rounded, upper = mac.fw_subset(w_init, num_candidates, max_iters=50) +result, rounded, upper = mac.solve(num_candidates, w_init, max_iters=100, rounding="madow", use_cache=True) init_selected = select_edges(candidate_edges, w_init) selected = select_edges(candidate_edges, result) init_selected_G = mac_to_nx(fixed_edges + init_selected) selected_G = mac_to_nx(fixed_edges + selected) -print(f"lambda2 Random: {mac.evaluate_objective(w_init)}") +print(f"lambda2 Initial: {mac.evaluate_objective(w_init)}") print(f"lambda2 Ours: {mac.evaluate_objective(result)}") plt.subplot(121) nx.draw(init_selected_G) -plt.title(f"Random Selection\n$\lambda_2$ = {mac.evaluate_objective(w_init):.2f}") +plt.title(f"Initial Selection\n$\lambda_2$ = {mac.evaluate_objective(w_init):.2f}") plt.subplot(122) nx.draw(selected_G) plt.title(f"Ours\n$\lambda_2$ = {mac.evaluate_objective(result):.2f}") diff --git a/mac/__init__.py b/mac/__init__.py index 1775a6d..e69de29 100644 --- a/mac/__init__.py +++ b/mac/__init__.py @@ -1 +0,0 @@ -from .mac import MAC diff --git a/mac/fiedler.py b/mac/fiedler.py deleted file mode 100644 index bd07d46..0000000 --- a/mac/fiedler.py +++ /dev/null @@ -1,95 +0,0 @@ -import numpy as np - -import networkx as nx -import networkx.linalg as la - -# TRACEMIN-Fiedler imports -import scipy as sp - -def find_fiedler_pair(L, method='tracemin_lu', tol=1e-8): - """ - Compute the second smallest eigenvalue of L and corresponding - eigenvector using `method` and tolerance `tol`. - - w: An element of [0,1]^m; this is the edge selection to use - method: Any method supported by NetworkX for computing algebraic - connectivity. See: - https://networkx.org/documentation/stable/reference/generated/networkx.linalg.algebraicconnectivity.algebraic_connectivity.html - - tol: Numerical tolerance for eigenvalue computation - - returns a tuple (lambda_2(L), v_2(L)) containing the Fiedler - value and corresponding vector. - - """ - assert method != 'lobpcg', "lobpcg is not currently supported" # LOBPCG not supported at the moment - - seed = np.random.RandomState(7) - q = min(4, L.shape[0] - 1) - X = np.asarray(seed.normal(size=(q, L.shape[0]))).T - if method == 'tracemin_cholesky': - from mac.cholesky_utils import tracemin_fiedler_cholesky - sigma, X = tracemin_fiedler_cholesky(L=L, X=X, normalized=False, tol=tol) - else: - sigma, X = la.algebraicconnectivity._tracemin_fiedler(L=L, X=X, normalized=False, tol=tol, method=method) - - return (sigma[0], X[:, 0]) - -def find_fiedler_pair_and_eigvecs(L, X=None, method='tracemin_lu', tol=1e-8, seed=None): - """ - Compute the second smallest eigenvalue of L and corresponding - eigenvector using `method` and tolerance `tol`. - - w: An element of [0,1]^m; this is the edge selection to use - method: Any method supported by NetworkX for computing algebraic - connectivity. See: - https://networkx.org/documentation/stable/reference/generated/networkx.linalg.algebraicconnectivity.algebraic_connectivity.html - - X: Initial guess for eigenvectors. If None, use random initial guess. - - tol: Numerical tolerance for eigenvalue computation - - returns a tuple (lambda_2(L), v_2(L)) containing the Fiedler - value and corresponding vector. - - """ - if seed is None: - seed = np.random.RandomState(7) - q = min(4, L.shape[0] - 1) - if X is None: - random_init = np.asarray(seed.normal(size=(q, L.shape[0]))).T - X = random_init - - # Check if X is a valid initial guess - assert X.shape[0] == L.shape[0] - assert X.shape[1] == q - - if method == 'tracemin_cholesky': - from mac.cholesky_utils import tracemin_fiedler_cholesky - sigma, X = tracemin_fiedler_cholesky(L=L, X=X, normalized=False, tol=tol) - else: - sigma, X = la.algebraicconnectivity._tracemin_fiedler(L=L, X=X, normalized=False, tol=tol, method=method) - - return (sigma[0], X[:, 0], X) - -def grad_from_fiedler(fiedler_vec, edge_list, weights): - """ - Compute a (super)gradient of the algebraic connectivity with respect to w - from the Fiedler vector. - - fiedler_vec: Eigenvector of the Laplacian corresponding to the second eigenvalue evaluated at a point 'w'. - edge_list: np.array of shape (m,2) containing the edge list of the graph. - weights: np.array of weights for each edge in the graph. - - returns grad F(w) from equation (8) of our paper: https://arxiv.org/pdf/2203.13897.pdf. - """ - grad = np.zeros(len(weights)) - - for k in range(len(weights)): - edge = edge_list[k] # get edge (i,j) - v_i = fiedler_vec[edge[0]] - v_j = fiedler_vec[edge[1]] - weight_k = weights[k] - kdelta = weight_k * (v_i - v_j) - grad[k] = kdelta * (v_i - v_j) - return grad diff --git a/mac/frankwolfe.py b/mac/frankwolfe.py deleted file mode 100644 index cb31d54..0000000 --- a/mac/frankwolfe.py +++ /dev/null @@ -1,183 +0,0 @@ -""" -Utilities for maximizing concave functions over simple feasible sets with Frank-Wolfe -""" - -from mac.utils import round_nearest -import numpy as np - -def frank_wolfe(initial, - problem, - solve_lp, - stepsize=None, - maxiter=50, - relative_duality_gap_tol=1e-5, - grad_norm_tol=1e-10, - verbose=False): - """Frank-Wolfe algorithm for maximizing a concave function. - - Parameters - ---------- - initial : array-like - Initial guess for the minimizer. - problem : callable - Function returning a tuple (f, gradf) where f is the objective function. - solve_lp : callable - Solve the Frank-Wolfe LP subproblem over the feasible set. - stepsize : callable, optional - Step size function. - maxiter : int, optional - Maximum number of iterations to perform. - relative_duality_gap_tol : float, optional - Tolerance for the *relative* duality gap (i.e., the duality gap divided - by the objective value). If the relative duality gap is less than this - tolerance, the algorithm terminates. - grad_norm_tol : float, optional - Tolerance for the norm of the gradient. If the norm of the gradient is - less than this tolerance, the algorithm terminates. - verbose : bool, optional - Whether to print the iteration number and the objective value. - - Returns - ------- - x : array-like - The minimizer. - - """ - if stepsize is None: - stepsize = lambda x, g, s, k: naive_stepsize(k) - - x = initial - u = float("inf") - for i in range(maxiter): - # Compute objective value and a (super)-gradient. - f, gradf = problem(x) - - # Solve the direction-finding subproblem by maximizing the linear - # approximation of f at x over the feasible set. - s = solve_lp(gradf) - - # Compute dual upper bound from linear approximation. - u = min(u, f + gradf @ (s - x)) - - # If the gradient norm is sufficiently small, we are done. - if np.linalg.norm(gradf) < grad_norm_tol: - if verbose: - print("Gradient norm is approximately 0. Found optimal solution") - return x, u - - # If the *relative* duality gap is sufficiently small, we are done. - if (u - f) / f < relative_duality_gap_tol: - if verbose: - print("Duality gap tolerance reached, found optimal solution") - return x, u - - x = x + stepsize(x, gradf, s, i) * (s - x) - pass - if verbose: - print("Reached maximum number of iterations, returning best solution") - return x, u - -def frank_wolfe_with_recycling(initial, - problem, - solve_lp, - stepsize=None, - maxiter=50, - relative_duality_gap_tol=1e-5, - grad_norm_tol=1e-10, - verbose=False): - """Frank-Wolfe algorithm for maximizing a concave function. Recycles problem data. - - Parameters - ---------- - initial : array-like - Initial guess for the minimizer. - problem : callable - Function returning a tuple (f, gradf, problem_data) where f is the - objective function. Takes as argument a point x and an object - problem_data - solve_lp : callable - Solve the Frank-Wolfe LP subproblem over the feasible set. - stepsize : callable, optional - Step size function. - maxiter : int, optional - Maximum number of iterations to perform. - relative_duality_gap_tol : float, optional - Tolerance for the *relative* duality gap (i.e., the duality gap divided - by the objective value). If the relative duality gap is less than this - tolerance, the algorithm terminates. - grad_norm_tol : float, optional - Tolerance for the norm of the gradient. If the norm of the gradient is - less than this tolerance, the algorithm terminates. - verbose : bool, optional - Whether to print the iteration number and the objective value. - - Returns - ------- - x : array-like - The minimizer. - - """ - if stepsize is None: - stepsize = lambda x, g, s, k: naive_stepsize(k) - - x = initial - u = float("inf") - problem_data = None - for i in range(maxiter): - # Compute objective value and a (super)-gradient. - f, gradf, problem_data = problem(x, problem_data) - - # Solve the direction-finding subproblem by maximizing the linear - # approximation of f at x over the feasible set. - s = solve_lp(gradf) - - # Compute dual upper bound from linear approximation. - u = min(u, f + gradf @ (s - x)) - - # If the gradient norm is sufficiently small, we are done. - if np.linalg.norm(gradf) < grad_norm_tol: - if verbose: - print("Gradient norm is approximately 0. Found optimal solution") - return x, u - - # If the *relative* duality gap is sufficiently small, we are done. - if (u - f) / f < relative_duality_gap_tol: - if verbose: - print("Duality gap tolerance reached, found optimal solution") - return x, u - - x = x + stepsize(x, gradf, s, i) * (s - x) - pass - if verbose: - print("Reached maximum number of iterations, returning best solution") - return x, u - -def solve_subset_box_lp(g, k): - """ - Maximize the objective function - g^T x - subject to - 0 <= x <= 1; ||x||_0 <= k - - The closed-form solution to this problem is given by a vector with 1's in - the positions of the largest k elements of g. - """ - return round_nearest(g, k) - -def solve_box_lp(g): - """ - Maximize the objective function - g^T x - subject to - 0 <= x <= 1 - - The closed-form solution to this problem is given by a vector with 1's in - the positions of all nonnegative elements of g. - - """ - solution = np.zeros_like(w) - solution[g >= 0.0] = 1.0 - return solution - -def naive_stepsize(k): - return 2.0 / (k + 2.0) diff --git a/mac/greedy_eig_minimal.py b/mac/greedy_eig_minimal.py deleted file mode 100644 index e115a34..0000000 --- a/mac/greedy_eig_minimal.py +++ /dev/null @@ -1,145 +0,0 @@ -from mac.utils import * -import numpy as np - -import networkx as nx -import networkx.linalg as la - -class MinimalGreedyEig: - def __init__(self, odom_measurements, lc_measurements, num_poses): - self.L_odom = weight_graph_lap_from_edge_list(odom_measurements, num_poses) - self.num_poses = num_poses - self.laplacian_e_list = [] - self.weights = [] - self.edge_list = [] - - for meas in lc_measurements: - laplacian_e = weight_graph_lap_from_edge_list([meas], num_poses) - self.laplacian_e_list.append(laplacian_e) - self.weights.append(meas.weight) - self.edge_list.append((meas.i,meas.j)) - - self.laplacian_e_list = np.array(self.laplacian_e_list) - self.weights = np.array(self.weights) - self.edge_list = np.array(self.edge_list) - - def find_fiedler_pair(self, L, method='tracemin_lu', tol=1e-8): - """ - Compute the second smallest eigenvalue of L and corresponding - eigenvector using `method` and tolerance `tol`. - - w: An element of [0,1]^m; this is the edge selection to use - method: Any method supported by NetworkX for computing algebraic - connectivity. See: - https://networkx.org/documentation/stable/reference/generated/networkx.linalg.algebraicconnectivity.algebraic_connectivity.html - - tol: Numerical tolerance for eigenvalue computation - - returns a tuple (lambda_2(L), v_2(L)) containing the Fiedler - value and corresponding vector. - - """ - find_fiedler_func = la.algebraicconnectivity._get_fiedler_func(method) - x = None - output = find_fiedler_func(L, x=x, normalized=False, tol=tol, seed=np.random.RandomState(7)) - return output - - def combined_laplacian(self, w, tol=1e-10): - """ - Construct the combined Laplacian (fixed edges plus candidate edges weighted by w). - - w: An element of [0,1]^m; this is the edge selection to use - tol: Tolerance for edges that are numerically zero. This improves speed - in situations where edges are not *exactly* zero, but close enough that - they have almost no influence on the graph. - - returns the matrix L(w) - """ - idx = np.where(w > tol) - prod = w[idx]*self.weights[idx] - C1 = weight_graph_lap_from_edges(self.edge_list[idx], prod, self.num_poses) - C = self.L_odom + C1 - return C - - def grad_from_fiedler(self, fiedler_vec): - """ - Compute a (super)gradient of the algebraic connectivity with respect to w - from the Fiedler vector. - - fiedler_vec: Eigenvector of the Laplacian corresponding to the second eigenvalue. - - returns grad F(w) from equation (8) of our paper: https://arxiv.org/pdf/2203.13897.pdf. - """ - grad = np.zeros(len(self.weights)) - - for k in range(len(self.weights)): - edge = self.edge_list[k] # get edge (i,j) - v_i = fiedler_vec[edge[0]] - v_j = fiedler_vec[edge[1]] - weight_k = self.weights[k] - kdelta = weight_k * (v_i - v_j) - grad[k] = kdelta * (v_i - v_j) - return grad - - def subset(self, k, save_intermediate=False): - # Initialize w to be all zeros - solution = np.zeros(len(self.weights)) - - # Compute the Fiedler vector and value for the fixed subgraph - l2, fiedler_vec = self.find_fiedler_pair(self.L_odom) - grad = self.grad_from_fiedler(fiedler_vec) - solution_l2 = l2 - solution_grad = grad - selected_edges = [] - for i in range(k): - # Placeholders to keep track of best measurement - best_idx = -1 - best_l2 = 0 - best_grad = np.zeros(len(self.weights)) - # Loop over all unselected measurements to find new best - # measurement to add - for j in range(len(self.weights)): - # If measurement j is already selected, skip it - if solution[j] > 0: - continue - - # Construct a test vector with edge j added - w = np.copy(solution) - w[j] = 1 - - # Compute a linear approximation to L2(w) evaluated at the - # current solution. The value of this linear approximation at - # w[j]=1 gives an upper bound on L2(w) at this point. This - # means that if u is less than the best l2, there's no way this - # edge will be the best, and we can safely discard it. - # u = best_l2 + solution_grad @ (w - solution) - u = solution_l2 + solution_grad[j] - if u < best_l2: - continue - - # If we made it here, we need to compute the true L2(w) and - # gradient. - L = self.combined_laplacian(w) - l2, fiedler_vec = self.find_fiedler_pair(L) - grad = self.grad_from_fiedler(fiedler_vec) - - # add a numerical tolerance for the comparison so we - # deterministically tie-break weighted effective resistances by - # selecting the first edge with the max weighted reff. This was - # actually a reason why some of the tests were failing - tol = 1e-8 - if l2 > best_l2 + tol: - best_idx = j - best_l2 = l2 - best_grad = grad - # If best_idx is still -1, something went terribly wrong, or there - # are no measurements - assert(best_idx != -1) - solution[best_idx] = 1 - solution_l2 = best_l2 - solution_grad = best_grad - best_edge_ij = self.edge_list[best_idx] - best_edge = Edge(best_edge_ij[0], best_edge_ij[1], self.weights[best_idx]) - selected_edges.append(best_edge) - - return solution, selected_edges - diff --git a/mac/greedy_esp_minimal.py b/mac/greedy_esp_minimal.py deleted file mode 100644 index 948989d..0000000 --- a/mac/greedy_esp_minimal.py +++ /dev/null @@ -1,108 +0,0 @@ -import numpy as np -from typing import List, Tuple -from mac.utils import weight_graph_lap_from_edge_list, weight_graph_lap_from_edges, Edge - -def incidence_vector(eij, num_nodes: int): - incidence_vec = np.zeros(num_nodes) - i = eij[0] - j = eij[1] - - incidence_vec[i] = 1 - incidence_vec[j] = -1 - return incidence_vec - -class MinimalGreedyESP: - def __init__(self, fixed_edges: List[Edge], candidate_edges: List[Edge], num_poses, use_reduced_laplacian=False): - self.L_odom = weight_graph_lap_from_edge_list(fixed_edges, num_poses) - self.num_poses = num_poses - self.laplacian_e_list = [] - self.edge_list = np.array([(e[0], e[1]) for e in candidate_edges]) - self.weights = np.array([e[2] for e in candidate_edges]) - self.use_reduced_laplacian = use_reduced_laplacian - - for edge in candidate_edges: - laplacian_e = weight_graph_lap_from_edge_list([edge], num_poses) - self.laplacian_e_list.append(laplacian_e) - - self.laplacian_e_list = np.array(self.laplacian_e_list) - - def combined_laplacian(self, w, tol=1e-10): - """ - Construct the combined Laplacian (fixed edges plus candidate edges weighted by w). - - w: An element of [0,1]^m; this is the edge selection to use - tol: Tolerance for edges that are numerically zero. This improves speed - in situations where edges are not *exactly* zero, but close enough that - they have almost no influence on the graph. - - returns the matrix L(w) - """ - idx = np.where(w > tol) - prod = w[idx]*self.weights[idx] - C1 = weight_graph_lap_from_edges(self.edge_list[idx], prod, self.num_poses) - C = self.L_odom + C1 - return C - - def subset(self, k: int) -> Tuple[np.ndarray, List[Edge]]: - """ - Apply GreedyEsp to select k edges from the candidate edges. - - Args: - k (int): the number of edges to select - - Returns: - np.ndarray: the selected edges as a boolean vector - List[Edge]: the selected edges - - """ - solution = np.zeros(len(self.weights)) - selected_edges: List[Edge] = [] - - for _ in range(k): - # Placeholders to keep track of best measurement - best_idx = -1 - best_weighted_reff = 0 - best_edge = None - # Loop over all unselected measurements to find new best - # measurement to add - for j in range(len(self.weights)): - # If measurement j is already selected, skip it - if solution[j] > 0: - continue - # Test solution - w = np.copy(solution) - curr_edge_ij = self.edge_list[j] - auv = incidence_vector(curr_edge_ij, self.num_poses).reshape((self.num_poses, 1)) - - reff = self.compute_reff(w, auv) - weighted_reff = self.weights[j] * reff - - # add a numerical tolerance for the comparison so we - # deterministically tie-break weighted effective resistances by - # selecting the first edge with the max weighted reff. This was - # actually a reason why some of the tests were failing - tol = 1e-8 - if weighted_reff > best_weighted_reff + tol: - best_idx = j - best_weighted_reff = weighted_reff - best_edge = Edge(curr_edge_ij[0], curr_edge_ij[1], self.weights[j]) - # If best_idx is still -1, something went terribly wrong, or there - # are no measurements - assert(best_idx != -1) - assert best_edge is not None - solution[best_idx] = 1 - selected_edges.append(best_edge) - return solution, selected_edges - - def compute_reff(self, w: np.ndarray, auv: np.ndarray): - """ - Compute the effective resistance of the selected edges - """ - L = self.combined_laplacian(w) - if self.use_reduced_laplacian: - auv = auv[1:] - Linv = np.linalg.inv(L.todense()[1:,1:]) - else: - Linv = np.linalg.pinv(L.todense()) - reff = (auv.T @ Linv @ auv)[0,0] - return reff diff --git a/mac/optimization/constraints.py b/mac/optimization/constraints.py new file mode 100644 index 0000000..41232cf --- /dev/null +++ b/mac/optimization/constraints.py @@ -0,0 +1,37 @@ +""" +Implements several linear program "oracles." + +Specifically, we provide utilities for computing the optimal solutions to +linear programs of the form: max_x , s.t. x in C, for a variety of +compact convex sets C. + +""" +from mac.utils.rounding import round_nearest +import numpy as np + +def solve_subset_box_lp(g, k): + """ + Maximize the objective function + g^T x + subject to + 0 <= x <= 1; ||x||_0 <= k + + The closed-form solution to this problem is given by a vector with 1's in + the positions of the largest k elements of g. + """ + return round_nearest(g, k) + +def solve_box_lp(g): + """ + Maximize the objective function + g^T x + subject to + 0 <= x <= 1 + + The closed-form solution to this problem is given by a vector with 1's in + the positions of all nonnegative elements of g. + + """ + solution = np.zeros_like(g) + solution[g > 0.0] = 1.0 + return solution diff --git a/mac/optimization/frankwolfe.py b/mac/optimization/frankwolfe.py new file mode 100644 index 0000000..db34673 --- /dev/null +++ b/mac/optimization/frankwolfe.py @@ -0,0 +1,80 @@ +""" +Utilities for maximizing concave functions over simple feasible sets with Frank-Wolfe +""" + +import numpy as np + +def naive_stepsize(k): + return 2.0 / (k + 2.0) + +def frank_wolfe(initial, + problem, + solve_lp, + stepsize=None, + maxiter=50, + relative_duality_gap_tol=1e-5, + grad_norm_tol=1e-10, + verbose=False): + """Frank-Wolfe algorithm for maximizing a concave function. + + Parameters + ---------- + initial : array-like + Initial guess for the minimizer. + problem : callable + Function returning a tuple (f, gradf) where f is the objective function. + solve_lp : callable + Solve the Frank-Wolfe LP subproblem over the feasible set. + stepsize : callable, optional + Step size function. + maxiter : int, optional + Maximum number of iterations to perform. + relative_duality_gap_tol : float, optional + Tolerance for the *relative* duality gap (i.e., the duality gap divided + by the objective value). If the relative duality gap is less than this + tolerance, the algorithm terminates. + grad_norm_tol : float, optional + Tolerance for the norm of the gradient. If the norm of the gradient is + less than this tolerance, the algorithm terminates. + verbose : bool, optional + Whether to print the iteration number and the objective value. + + Returns + ------- + x : array-like + The minimizer. + + """ + if stepsize is None: + stepsize = lambda x, g, s, k: naive_stepsize(k) + + x = initial + u = float("inf") + for i in range(maxiter): + # Compute objective value and a (super)-gradient. + f, gradf = problem(x) + + # Solve the direction-finding subproblem by maximizing the linear + # approximation of f at x over the feasible set. + s = solve_lp(gradf) + + # Compute dual upper bound from linear approximation. + u = min(u, f + gradf @ (s - x)) + + # If the gradient norm is sufficiently small, we are done. + if np.linalg.norm(gradf) < grad_norm_tol: + if verbose: + print("Gradient norm is approximately 0. Found optimal solution") + return x, u + + # If the *relative* duality gap is sufficiently small, we are done. + if (u - f) < relative_duality_gap_tol * abs(f): + if verbose: + print("Duality gap tolerance reached, found optimal solution") + return x, u + + x = x + stepsize(x, gradf, s, i) * (s - x) + if verbose: + print("Reached maximum number of iterations, returning best solution") + return x, u + diff --git a/mac/solvers/__init__.py b/mac/solvers/__init__.py new file mode 100644 index 0000000..b54ea33 --- /dev/null +++ b/mac/solvers/__init__.py @@ -0,0 +1,2 @@ +from .mac import MAC +from .baseline import NaiveGreedy diff --git a/mac/baseline.py b/mac/solvers/baseline.py similarity index 100% rename from mac/baseline.py rename to mac/solvers/baseline.py diff --git a/mac/greedy_eig.py b/mac/solvers/greedy_eig.py similarity index 98% rename from mac/greedy_eig.py rename to mac/solvers/greedy_eig.py index aa3cdbe..4108871 100644 --- a/mac/greedy_eig.py +++ b/mac/solvers/greedy_eig.py @@ -1,10 +1,10 @@ -from mac.utils import * +from mac.utils.graphs import * import numpy as np import networkx as nx import networkx.linalg as la -from mac.cholesky_utils import CholeskyFiedlerSolver +from mac.utils.cholesky import CholeskyFiedlerSolver class GreedyEig: def __init__(self, odom_measurements, lc_measurements, num_poses): diff --git a/mac/greedy_esp.py b/mac/solvers/greedy_esp.py similarity index 98% rename from mac/greedy_esp.py rename to mac/solvers/greedy_esp.py index eb6ed16..7565e33 100644 --- a/mac/greedy_esp.py +++ b/mac/solvers/greedy_esp.py @@ -34,11 +34,10 @@ from scipy.sparse import csc_matrix from typing import List, Set, Tuple, Union, Optional import math -from mac.utils import Edge, set_incidence_vector_for_edge_inplace -from mac.cholesky_utils import ( +from mac.utils.graphs import * +from mac.utils.cholesky import ( update_cholesky_factorization_inplace, get_cholesky_forward_solve, - weight_reduced_graph_lap_from_edge_list, ) def compute_weighted_effective_resistances( @@ -272,7 +271,7 @@ def subset_lazy(self, k: int, verbose: bool = False) -> Tuple[np.ndarray, List[E """A convenience function for subsets_lazy that only returns the result for a single budget. See subsets_lazy for more details. """ - results, selected_edges, times = self.subsets_lazy([k], return_times=return_time, verbose=verbose) + results, selected_edges, times = self.subsets_lazy([k], verbose=verbose) res = results[0] time = times[0] return res, selected_edges, time diff --git a/mac/mac.py b/mac/solvers/mac.py similarity index 50% rename from mac/mac.py rename to mac/solvers/mac.py index d8c819b..4f786f9 100644 --- a/mac/mac.py +++ b/mac/solvers/mac.py @@ -1,20 +1,26 @@ -from mac.utils import * -import mac.fiedler as fiedler -import mac.frankwolfe as fw -from timeit import default_timer as timer - import numpy as np import networkx as nx import networkx.linalg as la -import matplotlib.pyplot as plt - -from scipy.sparse import csc_matrix, csr_matrix +from dataclasses import dataclass +from typing import Optional from timeit import default_timer as timer +from mac.utils.graphs import * +from mac.utils.rounding import * + +import mac.utils.fiedler as fiedler +import mac.optimization.frankwolfe as fw +import mac.optimization.constraints as constraints + class MAC: + @dataclass + class Cache: + """Problem data cache""" + Q: Optional[np.ndarray] = None + def __init__(self, fixed_edges, candidate_edges, num_nodes, - fiedler_method='tracemin_lu', use_cache=False, fiedler_tol=1e-8, + fiedler_method='tracemin_lu', fiedler_tol=1e-8, min_selection_weight_tol=1e-10): """Parameters ---------- @@ -29,31 +35,34 @@ def __init__(self, fixed_edges, candidate_edges, num_nodes, 'tracemin_lu', 'tracemin_cholesky'. Default is 'tracemin_lu'. Using the 'tracemin_cholesky' method is faster but requires SuiteSparse to be installed. - use_cache : bool, optional - Whether to cache the Fiedler vector and the gradient. fiedler_tol : float, optional Tolerance for computing the Fiedler vector and corresponding eigenvalue. min_edge_selection_tol : float, optional Tolerance for the minimum edge selection weight. Default is 1e-10. """ - if (num_nodes == 0): - assert(len(fixed_edges) == len(candidate_edges) == 0) + # Check that we at least *could* have a spanning tree in the set + # {fixed_edges U candidate_edges} This does not guarantee that a + # spanning tree exists, but it's a good basic test. + num_edges = len(fixed_edges) + len(candidate_edges) + assert (num_nodes - 1) <= num_edges + + # We also check that there aren't "too many" edges. The number of edges + # in the complete graph K(n) is equal to n * (n - 1) / 2, so we cannot + # possibly have more edges than this. + assert num_edges <= 0.5 * num_nodes * (num_nodes - 1) + + # Pre-compute the Laplacian for the subgraph comprised of the "fixed edges". self.L_fixed = weight_graph_lap_from_edge_list(fixed_edges, num_nodes) self.num_nodes = num_nodes - self.laplacian_e_list = [] + self.weights = [] self.edge_list = [] - for edge in candidate_edges: - laplacian_e = weight_graph_lap_from_edge_list([edge], num_nodes) - self.laplacian_e_list.append(laplacian_e) self.weights.append(edge.weight) self.edge_list.append((edge.i, edge.j)) - self.laplacian_e_list = np.array(self.laplacian_e_list) self.weights = np.array(self.weights) self.edge_list = np.array(self.edge_list) - self.use_cache = use_cache # Configuration for Fiedler vector computation self.fiedler_method = fiedler_method @@ -62,95 +71,75 @@ def __init__(self, fixed_edges, candidate_edges, num_nodes, # Truncate edges with selection weights below this threshold self.min_selection_weight_tol = min_selection_weight_tol - def combined_laplacian(self, w): - """ - Construct the combined Laplacian (fixed edges plus candidate edges weighted by w). - - w: An element of [0,1]^m; this is the edge selection to use - tol: Tolerance for edges that are numerically zero. This improves speed - in situations where edges are not *exactly* zero, but close enough that - they have almost no influence on the graph. + def laplacian(self, x): + """Construct the combined Laplacian (fixed edges plus candidate edges weighted by x). + x: An element of [0,1]^m; this is the edge selection to use - returns the matrix L(w) - """ - idx = np.where(w > self.min_selection_weight_tol) - prod = w[idx]*self.weights[idx] - C1 = weight_graph_lap_from_edges(self.edge_list[idx], prod, self.num_nodes) - C = self.L_fixed + C1 - return C + The tolerance parameter `min_selection_weight_tol` is used to prune out + edges that are numerically zero. This improves speed in situations + where edges are not *exactly* zero, but close enough that they have + almost no influence on the graph. - def find_fiedler_pair(self, w): + returns the matrix L(x) """ - Compute the second smallest eigenvalue of L(w) and corresponding - eigenvector using `method` and tolerance `tol`. This is just a helper - that constructs L(w) and calls `fiedler.find_fiedler_pair` on the - resulting matrix, passing along the arguments. - - w: An element of [0,1]^m; this is the edge selection to use - method: Any method supported by NetworkX for computing algebraic - connectivity. See: - https://networkx.org/documentation/stable/reference/generated/networkx.linalg.algebraicconnectivity.algebraic_connectivity.html - - returns a tuple (lambda_2(L(w)), v_2(L(w))) containing the Fiedler - value and corresponding vector. + idx = np.where(x > self.min_selection_weight_tol) + prod = x[idx]*self.weights[idx] + L_candidate = weight_graph_lap_from_edges(self.edge_list[idx], prod, self.num_nodes) + L_x = self.L_fixed + L_candidate + return L_x + def evaluate_objective(self, x): """ - L = self.combined_laplacian(w) - if L.shape[0] == 0: - # If the graph is empty, then the Fiedler value is 0 and the - # Fiedler vector is empty. - return 0.0, np.array([]) - return fiedler.find_fiedler_pair(L, self.fiedler_method, tol=self.fiedler_tol) - - def evaluate_objective(self, w): - """ - Compute lambda_2(L(w)) where L(w) is the Laplacian with edge i weighted - by w_i and lambda_2 is the second smallest eigenvalue (this is the + Compute lambda_2(L(x)) where L(x) is the Laplacian with edge i weighted + by x_i*weight_i and lambda_2 is the second smallest eigenvalue (this is the algebraic connectivity). - w: Weights for each candidate edge (does not include fixed edges) + x: Weights for each candidate edge (does not include fixed edges) - returns F(w) = lambda_2(L(w)). + returns F(x) = lambda_2(L(x)). """ - return self.find_fiedler_pair(w)[0] + return fiedler.find_fiedler_pair(L=self.laplacian(x), + method=self.fiedler_method, tol=self.fiedler_tol)[0] - def grad_from_fiedler(self, fiedler_vec): - """ - Compute a (super)gradient of the algebraic connectivity with respect to w - from the Fiedler vector. + def problem(self, x, cache=None): + """Compute the algebraic connectivity of L(x) and a (super)gradient of the + algebraic connectivity with respect to x. - fiedler_vec: Eigenvector of the Laplacian corresponding to the second eigenvalue. + x: Weights for each candidate edge (does not include fixed edges) + cache: Mutable `Cache` object. If a `Cache` object is provided in the `cache` field, it will be used + and updated, but not explicitly returned. Rather, it will be updated directly. - returns grad F(w) from equation (8) of our paper: https://arxiv.org/pdf/2203.13897.pdf. + returns x, grad F(x). """ - return fiedler.grad_from_fiedler(fiedler_vec, self.edge_list, self.weights) - - def problem(self, x): - f, fiedler_vec = self.find_fiedler_pair(x) - gradf = self.grad_from_fiedler(fiedler_vec) - return (f, gradf) - - def problem_with_recycling(self, x, Q=None): - """ - Q is a set of recycled eigenvectors. - NOTE: experimental - """ - f, fiedler_vec, Q = fiedler.find_fiedler_pair_and_eigvecs(self.combined_laplacian(x), X=Q, method=self.fiedler_method, tol=self.fiedler_tol) - gradf = self.grad_from_fiedler(fiedler_vec) - return f, gradf, Q - - def fw_subset(self, w_init, k, rounding="nearest", fallback=False, - max_iters=5, relative_duality_gap_tol=1e-4, - grad_norm_tol=1e-8, random_rounding_max_iters=1, verbose=False, return_rounding_time=False): + Q = None if cache is None else cache.Q + f, fiedler_vec, Qnew = fiedler.find_fiedler_pair(L=self.laplacian(x), X=Q) + + gradf = np.zeros(len(self.weights)) + for k in range(len(self.weights)): + edge = self.edge_list[k] # get edge (i,j) + v_i = fiedler_vec[edge[0]] + v_j = fiedler_vec[edge[1]] + weight_k = self.weights[k] + kdelta = weight_k * (v_i - v_j) + gradf[k] = kdelta * (v_i - v_j) + + if cache is not None: + cache.Q = Q + return f, gradf + + def solve(self, k, x_init=None, rounding="nearest", fallback=False, + max_iters=5, relative_duality_gap_tol=1e-4, + grad_norm_tol=1e-8, random_rounding_max_iters=1, + verbose=False, return_rounding_time=False, use_cache=False): """Use the Frank-Wolfe method to solve the subset selection problem,. Parameters ---------- - w_init : Array-like - Initial weights for the candidate edges, must satisfy 0 <= w_i <= 1, |w| <= k. This - is the starting point for the Frank-Wolfe algorithm. TODO(kevin): make optional k : int Number of edges to select. + x_init : optional, array-like + Initial weights for the candidate edges, must satisfy 0 <= w_i <= 1, |w| <= k. This + is the starting point for the Frank-Wolfe algorithm. TODO(kevin): make optional rounding : str, optional Rounding method to use. Options are "nearest" (default) and "madow" (a random rounding procedure). @@ -190,27 +179,25 @@ def fw_subset(self, w_init, k, rounding="nearest", fallback=False, return result, result, self.evaluate_objective(np.ones(len(self.weights))) - assert(len(w_init) == len(self.weights)) + # TODO handle case where x is none + assert(len(x_init) == len(self.weights)) # Solution for the direction-finding subproblem - solve_lp = lambda g: fw.solve_subset_box_lp(g, k) + solve_lp = lambda g: constraints.solve_subset_box_lp(g, k) + + # Set up problem to use cache (or not) + cache = None + if use_cache: + cache = MAC.Cache() + problem = lambda x: self.problem(x, cache=cache) # Run Frank-Wolfe to solve the relaxation of subset constrained # algebraic connectivity maximization - if self.use_cache: - w, u = fw.frank_wolfe_with_recycling(initial=w_init, - problem=self.problem_with_recycling, - solve_lp=solve_lp, - maxiter=max_iters, - relative_duality_gap_tol=relative_duality_gap_tol, - grad_norm_tol=grad_norm_tol, - verbose=verbose) - else: - w, u = fw.frank_wolfe(initial=w_init, problem=self.problem, - solve_lp=solve_lp, maxiter=max_iters, - relative_duality_gap_tol=relative_duality_gap_tol, - grad_norm_tol=grad_norm_tol, - verbose=verbose) + w, u = fw.frank_wolfe(initial=x_init, problem=problem, + solve_lp=solve_lp, maxiter=max_iters, + relative_duality_gap_tol=relative_duality_gap_tol, + grad_norm_tol=grad_norm_tol, + verbose=verbose) start = timer() if rounding == "madow": @@ -222,7 +209,7 @@ def fw_subset(self, w_init, k, rounding="nearest", fallback=False, rounding_time = end - start if fallback: - init_f = self.evaluate_objective(w_init) + init_f = self.evaluate_objective(x_init) rounded_f = self.evaluate_objective(rounded) # If the rounded solution is worse than the initial solution, then diff --git a/mac/utils.py b/mac/utils.py deleted file mode 100644 index 11b156c..0000000 --- a/mac/utils.py +++ /dev/null @@ -1,335 +0,0 @@ -import numpy as np -from collections import namedtuple -import matplotlib.pyplot as plt -from scipy.sparse import csr_matrix, coo_matrix, csc_matrix -import networkx as nx -import math -from typing import List, Union, Tuple - -# Define Edge container -Edge = namedtuple("Edge", ["i", "j", "weight"]) - -def round_nearest(w, k, weights=None, break_ties_decimal_tol=None): - """ - Round a solution w to the relaxed problem, i.e. w \in [0,1]^m, |w| = k to a - solution to the original problem with w_i \in {0,1}. Ties between edges - are broken based on the original weight (all else being equal, we - prefer edges with larger weight). - - w: A solution in the feasible set for the relaxed problem - weights: The original weights of the edges to use as a tiebreaker - k: The number of edges to select - break_ties_decimal_tol: tolerance for determining floating point equality of weights w. If two selection weights are equal to this many decimal places, we break the tie based on the original weight. - - returns w': A solution in the feasible set for the original problem - """ - if weights is None or break_ties_decimal_tol is None: - # If there are no tiebreakers, just set the top k elements of w to 1, - # and the rest to 0 - idx = np.argpartition(w, -k)[-k:] - rounded = np.zeros(len(w)) - if k > 0: - rounded[idx] = 1.0 - return rounded - - # If there are tiebreakers, we truncate the selection weights to the - # specified number of decimal places, and then break ties based on the - # original weights - truncated_w = w.round(decimals=break_ties_decimal_tol) - zipped_vals = np.array( - [(truncated_w[i], weights[i]) for i in range(len(w))], - dtype=[("w", "float"), ("weight", "float")], - ) - idx = np.argpartition(zipped_vals, -k, order=["w", "weight"])[-k:] - rounded = np.zeros(len(w)) - if k > 0: - rounded[idx] = 1.0 - return rounded - -def round_random(w, k): - """ - Round a solution w to the relaxed problem, i.e. w \in [0,1]^m, |w| = k to - one with hard edge constraints and satisfying the constraint that the - expected number of selected edges is equal to k. - - w: A solution in the feasible set for the relaxed problem - k: The number of edges to select _in expectation_ - - returns w': A solution containing hard edge selections with an expected - number of selected edges equal to k. - """ - x = np.zeros(len(w)) - for i in range(len(w)): - r = np.random.rand() - if w[i] > r: - x[i] = 1.0 - return x - -def round_madow(w, k, seed=None, value_fn=None, max_iters=1): - if value_fn is None or max_iters == 1: - return round_madow_base(w, k, seed) - - best_x = None - best_val = -np.inf - for i in range(max_iters): - x = round_madow_base(w, k, seed) - val = value_fn(x) - if val > best_val: - best_val = val - best_x = x - return best_x - - -def round_madow_base(w, k, seed=None): - """ - Use Madow rounding - """ - if seed is None: - u = np.random.rand() - else: - u = seed.rand() - x = np.zeros(len(w)) - # pi = np.zeros(len(w) + 1) - pi = np.zeros(len(w)) - sumw = np.cumsum(w) - pi[1:] = sumw[:-1] - for i in range(k): - total = u + i - x[np.where((pi <= total) & (total < sumw))] = 1.0 - # for j in range(len(w)): - # if (x[j] != 1) and (pi[j] <= total) and (total < pi[j+1]): - # x[j] = 1.0 - # break - - assert np.sum(x) == k, f"Error: {np.sum(x)} != {k}" - return x - -def nx_to_mac(G: nx.Graph) -> List[Edge]: - """Returns the list of edges in the graph G - - Args: - G (nx.Graph): the graph - - Returns: - List[Edge]: the list of edges - """ - edges = [] - for nxedge in G.edges(): - edge = Edge(nxedge[0], nxedge[1], 1.0) - edges.append(edge) - return edges - - -def mac_to_nx(edges: List[Edge]) -> nx.Graph: - """returns the graph corresponding to the list of edges - - Args: - edges (List[Edge]): the list of edges - - Returns: - nx.Graph: the networkx graph - """ - G = nx.Graph() - for edge in edges: - G.add_edge(edge.i, edge.j, weight=edge.weight) - return G - - -def weight_graph_lap_from_edge_list(edges: List[Edge], num_nodes: int) -> csr_matrix: - """Returns the (sparse) weighted graph Laplacian matrix from the list of edges - - Args: - edges (List[Edge]): the list of edges - num_nodes (int): the number of variables - - Returns: - csr_matrix: the weighted graph Laplacian matrix - """ - # Preallocate triplets - rows = [] - cols = [] - data = [] - for edge in edges: - # Diagonal elem (u,u) - rows.append(edge.i) - cols.append(edge.i) - data.append(edge.weight) - - # Diagonal elem (v,v) - rows.append(edge.j) - cols.append(edge.j) - data.append(edge.weight) - - # Off diagonal (u,v) - rows.append(edge.i) - cols.append(edge.j) - data.append(-edge.weight) - - # Off diagonal (v,u) - rows.append(edge.j) - cols.append(edge.i) - data.append(-edge.weight) - - return csr_matrix(coo_matrix((data, (rows, cols)), shape=[num_nodes, num_nodes])) - - -def weight_reduced_graph_lap_from_edge_list( - edges: List[Edge], num_nodes: int -) -> csr_matrix: - graph_lap = weight_graph_lap_from_edge_list(edges, num_nodes) - return graph_lap[1:, 1:] - - -def weight_graph_lap_from_edges( - edges: List[Tuple[int, int]], weights: List[int], num_poses: int -) -> csr_matrix: - """Returns the sparse rotational weighted graph Laplacian matrix from a list - of edges and edge weights - - Args: - edges (List[Tuple[int, int]]): the list of edge indices (i,j) - weights (List[float]): the list of edge weights - num_poses (int): the number of poses - - Returns: - csr_matrix: the rotational weighted graph Laplacian matrix - """ - assert len(edges) == len(weights) - # Preallocate triplets - rows = [] - cols = [] - data = [] - for i in range(len(edges)): - # Diagonal elem (u,u) - rows.append(edges[i, 0]) - cols.append(edges[i, 0]) - data.append(weights[i]) - - # Diagonal elem (v,v) - rows.append(edges[i, 1]) - cols.append(edges[i, 1]) - data.append(weights[i]) - - # Off diagonal (u,v) - rows.append(edges[i, 0]) - cols.append(edges[i, 1]) - data.append(-weights[i]) - - # Off diagonal (v,u) - rows.append(edges[i, 1]) - cols.append(edges[i, 0]) - data.append(-weights[i]) - - return csr_matrix(coo_matrix((data, (rows, cols)), shape=[num_poses, num_poses])) - - -def split_edges(edges: List[Edge]) -> Tuple[List[Edge], List[Edge]]: - """Splits list of edges into a "fixed" chain part and a set of candidate - loops. - - Args: - edges: List of edges. - - Returns: - A tuple containing two lists: - fixed: edges where |i - j| = 1, and - candidates: edges where |i - j| != 1 - - This is particularly useful for pose-graph SLAM applications where the - fixed edges correspond to an odometry chain and the candidate edges - correspond to loop closures. - - """ - chain_edges = [] - loop_edges = [] - for edge in edges: - id1 = edge.i - id2 = edge.j - if abs(id2 - id1) > 1: - loop_edges.append(edge) - else: - chain_edges.append(edge) - - return chain_edges, loop_edges - - -def select_edges(edges, w): - """ - Select the subset of edges from `edges` with weight equal to one in - `w`. - """ - assert len(edges) == len(w), f"Selection mask length {len(w)} does not match number of edges {len(edges)}" - edges_out = [] - for i, edge in enumerate(edges): - if w[i] == 1.0: - edges_out.append(edge) - return edges_out - - -def get_incidence_vector(eij: Union[Edge, Tuple[int, int]], num_nodes: int): - """Returns the incidence vector for the edge eij - - Args: - eij (Union[Edge, Tuple[int, int]]): the edge - num_nodes (int): the number of nodes - - Returns: - np.ndarray: the incidence vector - """ - incidence_vec = np.zeros(num_nodes) - i = eij[0] - j = eij[1] - - incidence_vec[i] = 1 - incidence_vec[j] = -1 - return incidence_vec - - -def set_incidence_vector_for_edge_inplace( - auv_vec: np.ndarray, edge: Union[Edge, Tuple[int, int]], num_nodes: int -) -> None: - """Modifies the passed in auv vector to be the correct values for the - given edge. - - NOTE: Assumes that the edge indices are in the range [0, num_nodes) and - that we are getting the incidence vector for a reduced Laplacian - (i.e. indices are shifted by -1). - - Args: - auv_vec (np.ndarray): the auv vector to modify - edge (Union[Edge, Tuple[int, int]]): the edge - num_nodes (int): the number of nodes - - """ - assert len(auv_vec) == num_nodes - 1 - auv_vec.fill(0) - i = edge[0] - 1 - j = edge[1] - 1 - if i >= 0: - auv_vec[i] = 1 - if j >= 0: - auv_vec[j] = -1 - - -def get_edge_selection_as_binary_mask( - edges: List[Edge], selected_edges: List[Edge] -) -> np.ndarray: - """ - Returns a binary mask of the selected edges. - - Args: - edges: List of edges. - selected_edges: List of selected edges. - - Returns: - A binary mask of the selected edges. - """ - assert len(edges) >= len( - selected_edges - ), "The number of selected edges cannot be greater than the total number of edges." - mask = np.zeros(len(edges)) - for i, edge in enumerate(edges): - if edge in selected_edges: - mask[i] = 1.0 - return mask - diff --git a/mac/cholesky_utils.py b/mac/utils/cholesky.py similarity index 99% rename from mac/cholesky_utils.py rename to mac/utils/cholesky.py index 223646a..34cecb6 100644 --- a/mac/cholesky_utils.py +++ b/mac/utils/cholesky.py @@ -1,5 +1,6 @@ +import math import numpy as np -from mac.utils import * +from mac.utils.graphs import * from scipy.sparse import csr_matrix, coo_matrix, csc_matrix from sksparse.cholmod import cholesky, Factor, analyze, cholesky_AAt, CholmodNotPositiveDefiniteError @@ -107,7 +108,7 @@ class _CholeskySolver: """Cholesky factorization. To solve Ax = b: - solver = _LUSolver(A) + solver = _CholeskySolver(A) x = solver.solve(b) optional argument `tol` on solve method is ignored but included diff --git a/mac/utils/conversions.py b/mac/utils/conversions.py new file mode 100644 index 0000000..e134e68 --- /dev/null +++ b/mac/utils/conversions.py @@ -0,0 +1,49 @@ +""" +Utilities for converting between graph types. +""" +import networkx as nx +from typing import List + +from .graphs import Edge + +def nx_to_mac(G: nx.Graph) -> List[Edge]: + """Returns the list of edges in the graph G + + Args: + G (nx.Graph): the graph + + Returns: + List[Edge]: the list of edges + """ + edges = [] + for nxedge in G.edges(): + i = nxedge[0] + j = nxedge[1] + data = G.get_edge_data(i, j) + weight = 1.0 + if "weight" in data: + weight = data["weight"] + if i < j: + edge = Edge(i, j, weight) + else: + edge = Edge(j, i, weight) + edges.append(edge) + return edges + + +def mac_to_nx(edges: List[Edge]) -> nx.Graph: + """returns the graph corresponding to the list of edges + + Args: + edges (List[Edge]): the list of edges + + Returns: + nx.Graph: the networkx graph + """ + G = nx.Graph() + for edge in edges: + if edge.i < edge.j: + G.add_edge(edge.i, edge.j, weight=edge.weight) + else: + G.add_edge(edge.j, edge.i, weight=edge.weight) + return G diff --git a/mac/utils/fiedler.py b/mac/utils/fiedler.py new file mode 100644 index 0000000..98539ce --- /dev/null +++ b/mac/utils/fiedler.py @@ -0,0 +1,45 @@ +import numpy as np + +import networkx as nx +import networkx.linalg as la + +# TRACEMIN-Fiedler imports +import scipy as sp + +def find_fiedler_pair(L, X=None, method='tracemin_lu', tol=1e-8, seed=None): + """ + Compute the second smallest eigenvalue of L and corresponding + eigenvector using `method` and tolerance `tol`. + + w: An element of [0,1]^m; this is the edge selection to use + method: Any method supported by NetworkX for computing algebraic + connectivity. See: + https://networkx.org/documentation/stable/reference/generated/networkx.linalg.algebraicconnectivity.algebraic_connectivity.html + + X: Initial guess for eigenvectors. If None, use random initial guess. + + tol: Numerical tolerance for eigenvalue computation + + returns a tuple (lambda_2(L), v_2(L)) containing the Fiedler + value and corresponding vector. + + """ + if seed is None: + seed = np.random.RandomState(7) + q = min(4, L.shape[0] - 1) + if X is None: + random_init = np.asarray(seed.normal(size=(q, L.shape[0]))).T + X = random_init + + # Check if X is a valid initial guess + assert X.shape[0] == L.shape[0] + assert X.shape[1] == q + + if method == 'tracemin_cholesky': + from mac.utils.cholesky import tracemin_fiedler_cholesky + sigma, X = tracemin_fiedler_cholesky(L=L, X=X, normalized=False, tol=tol) + else: + sigma, X = la.algebraicconnectivity._tracemin_fiedler(L=L, X=X, normalized=False, tol=tol, method=method) + + return (sigma[0], X[:, 0], X) + diff --git a/mac/utils/graphs.py b/mac/utils/graphs.py new file mode 100644 index 0000000..0fb9cb0 --- /dev/null +++ b/mac/utils/graphs.py @@ -0,0 +1,179 @@ +""" +Types and utilities for working with graphs. +""" + +import numpy as np +from typing import List, Union, Tuple +from collections import namedtuple +from scipy.sparse import csr_matrix, coo_matrix, csc_matrix + +# Define Edge container +Edge = namedtuple("Edge", ["i", "j", "weight"]) + +def weight_graph_lap_from_edge_list(edges: List[Edge], num_nodes: int) -> csr_matrix: + """Returns the (sparse) weighted graph Laplacian matrix from the list of edges + + Args: + edges (List[Edge]): the list of edges + num_nodes (int): the number of variables + + Returns: + csr_matrix: the weighted graph Laplacian matrix + """ + # Preallocate triplets + rows = [] + cols = [] + data = [] + for edge in edges: + # Diagonal elem (u,u) + rows.append(edge.i) + cols.append(edge.i) + data.append(edge.weight) + + # Diagonal elem (v,v) + rows.append(edge.j) + cols.append(edge.j) + data.append(edge.weight) + + # Off diagonal (u,v) + rows.append(edge.i) + cols.append(edge.j) + data.append(-edge.weight) + + # Off diagonal (v,u) + rows.append(edge.j) + cols.append(edge.i) + data.append(-edge.weight) + + return csr_matrix(coo_matrix((data, (rows, cols)), shape=[num_nodes, num_nodes])) + + +def weight_reduced_graph_lap_from_edge_list( + edges: List[Edge], num_nodes: int +) -> csr_matrix: + graph_lap = weight_graph_lap_from_edge_list(edges, num_nodes) + return graph_lap[1:, 1:] + + +def weight_graph_lap_from_edges( + edges: List[Tuple[int, int]], weights: List[int], num_nodes: int +) -> csr_matrix: + """Returns the weighted graph Laplacian matrix from a list + of edges and edge weights + + Args: + edges (List[Tuple[int, int]]): the list of edge indices (i,j) + weights (List[float]): the list of edge weights + num_nodes (int): the number of nodes + + Returns: + csr_matrix: the weighted graph Laplacian matrix + """ + assert len(edges) == len(weights) + # Preallocate triplets + rows = [] + cols = [] + data = [] + for i in range(len(edges)): + # Diagonal elem (u,u) + rows.append(edges[i, 0]) + cols.append(edges[i, 0]) + data.append(weights[i]) + + # Diagonal elem (v,v) + rows.append(edges[i, 1]) + cols.append(edges[i, 1]) + data.append(weights[i]) + + # Off diagonal (u,v) + rows.append(edges[i, 0]) + cols.append(edges[i, 1]) + data.append(-weights[i]) + + # Off diagonal (v,u) + rows.append(edges[i, 1]) + cols.append(edges[i, 0]) + data.append(-weights[i]) + + return csr_matrix(coo_matrix((data, (rows, cols)), shape=[num_nodes, num_nodes])) + + +def select_edges(edges, w): + """ + Select the subset of edges from `edges` with weight equal to one in + `w`. + """ + assert len(edges) == len(w), f"Selection mask length {len(w)} does not match number of edges {len(edges)}" + edges_out = [] + for i, edge in enumerate(edges): + if w[i] == 1.0: + edges_out.append(edge) + return edges_out + + +def get_incidence_vector(eij: Union[Edge, Tuple[int, int]], num_nodes: int): + """Returns the incidence vector for the edge eij + + Args: + eij (Union[Edge, Tuple[int, int]]): the edge + num_nodes (int): the number of nodes + + Returns: + np.ndarray: the incidence vector + """ + incidence_vec = np.zeros(num_nodes) + i = eij[0] + j = eij[1] + + incidence_vec[i] = 1 + incidence_vec[j] = -1 + return incidence_vec + + +def set_incidence_vector_for_edge_inplace( + auv_vec: np.ndarray, edge: Union[Edge, Tuple[int, int]], num_nodes: int +) -> None: + """Modifies the passed in auv vector to be the correct values for the + given edge. + + NOTE: Assumes that the edge indices are in the range [0, num_nodes) and + that we are getting the incidence vector for a reduced Laplacian + (i.e. indices are shifted by -1). + + Args: + auv_vec (np.ndarray): the auv vector to modify + edge (Union[Edge, Tuple[int, int]]): the edge + num_nodes (int): the number of nodes + + """ + assert len(auv_vec) == num_nodes - 1 + auv_vec.fill(0) + i = edge[0] - 1 + j = edge[1] - 1 + if i >= 0: + auv_vec[i] = 1 + if j >= 0: + auv_vec[j] = -1 + + +def get_edge_selection_as_binary_mask( + edges: List[Edge], selected_edges: List[Edge] +) -> np.ndarray: + """ + Returns a binary mask of the selected edges. + + Args: + edges: List of edges. + selected_edges: List of selected edges. + + Returns: + A binary mask of the selected edges. + """ + assert len(edges) >= len( + selected_edges + ), "The number of selected edges cannot be greater than the total number of edges." + mask = np.zeros(len(edges)) + for i, edge in enumerate(edges): + if edge in selected_edges: + mask[i] = 1.0 + return mask diff --git a/mac/utils/rounding.py b/mac/utils/rounding.py new file mode 100644 index 0000000..5b5545a --- /dev/null +++ b/mac/utils/rounding.py @@ -0,0 +1,95 @@ +""" +Utilities for rounding solutions to convex optimization problems onto simple constraint sets. +""" + +import numpy as np + +def round_nearest(w, k, weights=None, break_ties_decimal_tol=None): + """ + Round a solution w to the relaxed problem, i.e. w \in [0,1]^m, |w| = k to a + solution to the original problem with w_i \in {0,1}. Ties between edges + are broken based on the original weight (all else being equal, we + prefer edges with larger weight). + + w: A solution in the feasible set for the relaxed problem + weights: The original weights of the edges to use as a tiebreaker + k: The number of edges to select + break_ties_decimal_tol: tolerance for determining floating point equality of weights w. If two selection weights are equal to this many decimal places, we break the tie based on the original weight. + + returns w': A solution in the feasible set for the original problem + """ + if weights is None or break_ties_decimal_tol is None: + # If there are no tiebreakers, just set the top k elements of w to 1, + # and the rest to 0 + idx = np.argpartition(w, -k)[-k:] + rounded = np.zeros(len(w)) + if k > 0: + rounded[idx] = 1.0 + return rounded + + # If there are tiebreakers, we truncate the selection weights to the + # specified number of decimal places, and then break ties based on the + # original weights + truncated_w = w.round(decimals=break_ties_decimal_tol) + zipped_vals = np.array( + [(truncated_w[i], weights[i]) for i in range(len(w))], + dtype=[("w", "float"), ("weight", "float")], + ) + idx = np.argpartition(zipped_vals, -k, order=["w", "weight"])[-k:] + rounded = np.zeros(len(w)) + if k > 0: + rounded[idx] = 1.0 + return rounded + +def round_random(w, k): + """ + Round a solution w to the relaxed problem, i.e. w \in [0,1]^m, |w| = k to + one with hard edge constraints and satisfying the constraint that the + expected number of selected edges is equal to k. + + w: A solution in the feasible set for the relaxed problem + k: The number of edges to select _in expectation_ + + returns w': A solution containing hard edge selections with an expected + number of selected edges equal to k. + """ + x = np.zeros(len(w)) + for i in range(len(w)): + r = np.random.rand() + if w[i] > r: + x[i] = 1.0 + return x + +def round_madow(w, k, seed=None, value_fn=None, max_iters=1): + if value_fn is None or max_iters == 1: + return round_madow_base(w, k, seed) + + best_x = None + best_val = -np.inf + for i in range(max_iters): + x = round_madow_base(w, k, seed) + val = value_fn(x) + if val > best_val: + best_val = val + best_x = x + return best_x + + +def round_madow_base(w, k, seed=None): + """ + Use Madow rounding + """ + if seed is None: + u = np.random.rand() + else: + u = seed.rand() + x = np.zeros(len(w)) + pi = np.zeros(len(w)) + sumw = np.cumsum(w) + pi[1:] = sumw[:-1] + for i in range(k): + total = u + i + x[np.where((pi <= total) & (total < sumw))] = 1.0 + + assert np.sum(x) == k, f"Error: {np.sum(x)} != {k}" + return x diff --git a/test/test_cholesky_utils.py b/test/test_cholesky_utils.py deleted file mode 100644 index d78a397..0000000 --- a/test/test_cholesky_utils.py +++ /dev/null @@ -1,231 +0,0 @@ -""" -Copyright 2023 MIT Marine Robotics Group - -Regression tests for Cholesky utilities - -Author: Kevin Doherty -""" -import unittest -import numpy as np - -np.set_printoptions(precision=4) -from mac.utils import ( - select_edges, - split_edges, - nx_to_mac, - mac_to_nx, - set_incidence_vector_for_edge_inplace, - weight_reduced_graph_lap_from_edge_list, - weight_graph_lap_from_edge_list - ) -from mac.cholesky_utils import ( - update_cholesky_factorization_inplace, - get_cholesky_forward_solve, - get_matrix_from_chol_factor_with_original_ordering, - find_fiedler_pair_cholesky, - CholeskyFiedlerSolver -) -from mac.fiedler import find_fiedler_pair -from scipy.sparse import spmatrix -from sksparse.cholmod import cholesky, Factor, analyze, cholesky_AAt -from .utils import get_split_petersen_graph, get_split_erdos_renyi_graph - - -ORDERING_METHODS = ["natural", "best", "amd", "metis", "nesdis", "colamd", "default"] - -class TestCholeskyUtils(unittest.TestCase): - def test_set_auv_inplace(self): - # Get the Erdos-Renyi graph - fixed_edges, candidate_edges, n = get_split_erdos_renyi_graph() - - auv_inplace_edge = np.zeros(n - 1) - auv_inplace_tuple = np.zeros(n - 1) - # iterate over all possible edges - for e in candidate_edges: - # set the AUV to be the edge - auv_on_demand = np.zeros_like(auv_inplace_edge) - set_incidence_vector_for_edge_inplace(auv_inplace_edge, e, n) - i, j, _ = e - edge_indices = (i, j) - set_incidence_vector_for_edge_inplace(auv_inplace_tuple, edge_indices, n) - - i -= 1 - j -= 1 - if i >= 0: - auv_on_demand[i] = 1 - if j >= 0: - auv_on_demand[j] = -1 - - # check that the two are the same - self.assertTrue(np.allclose(auv_inplace_edge, auv_on_demand)) - self.assertTrue(np.allclose(auv_inplace_tuple, auv_on_demand)) - - def test_update_cholesky_inplace_and_get_factored_matrix(self): - # TODO right now we're testing that these 2 functions agree, but we - # should try to make these tests independent of each other - - fixed_edges, candidate_edges, num_nodes = get_split_erdos_renyi_graph() - for ord_method in ORDERING_METHODS: - with self.subTest(ord_method=ord_method): - laplacian = weight_reduced_graph_lap_from_edge_list(fixed_edges, num_nodes) - laplacian_factorization = cholesky( - laplacian, beta=0, ordering_method=ord_method - ) - laplacian = laplacian.toarray() - - auv = np.zeros(num_nodes - 1) - for e in candidate_edges: - # set the AUV to be the edge - set_incidence_vector_for_edge_inplace(auv, e, num_nodes) - lap_update = np.outer(auv, auv) * e.weight - laplacian += lap_update - assert isinstance(laplacian, np.ndarray) - - # update the factorization - update_cholesky_factorization_inplace( - laplacian_factorization, e, num_nodes, reduced=True, subtract=False - ) - sksparse_matrix = get_matrix_from_chol_factor_with_original_ordering( - laplacian_factorization - ) - if isinstance(sksparse_matrix, spmatrix): - sksparse_matrix = sksparse_matrix.toarray() - - # check that the factorization is still valid - self.assertTrue( - np.allclose(sksparse_matrix, laplacian, atol=1e-6), - msg=f"""The factorization is not valid after updating with edge - {e}\n - The sksparse matrix is: \n {sksparse_matrix}\n - The laplacian is: \n {laplacian}\n - The difference is: \n {np.round(sksparse_matrix - laplacian,2)}""", - ) - - def test_lower_triangular_solve_norm(self): - """Test that the lower triangular solve returns a vector with the - correct norm - - We are solving the system Lx = b, where L is a lower triangular matrix - corresponding to the Cholesky factorization of a matrix and b is a - vector. - - We will test that the solution is correct by comparing it to the - solution from the numpy solve function. - """ - - def verify_solutions(np_laplacian, chol_factor, b): - np_lower = np.linalg.cholesky(np_laplacian) - np_solve = np.linalg.lstsq(np_lower, b, rcond=None)[0] - sksparse_solve = get_cholesky_forward_solve(chol_factor, b) - - # verify that the norms are the same - self.assertTrue( - np.allclose(np.linalg.norm(np_solve), np.linalg.norm(sksparse_solve)) - ) - - # double check the norm by using the inv function - chol_factor_inv = chol_factor.inv().toarray() - reff = b.T @ chol_factor_inv @ b - self.assertAlmostEqual(reff, np.linalg.norm(sksparse_solve) ** 2) - - - fixed_edges, candidate_edges, num_nodes = get_split_erdos_renyi_graph() - - for ord_method in ORDERING_METHODS: - with self.subTest(ord_method=ord_method): - laplacian = weight_reduced_graph_lap_from_edge_list(fixed_edges, num_nodes) - laplacian_factorization = cholesky( - laplacian, beta=0, ordering_method=ord_method - ) - laplacian = laplacian.toarray() - - b = np.random.rand(num_nodes - 1) - idx = -1 - verify_solutions(laplacian, laplacian_factorization, b) - - auv = np.zeros(num_nodes - 1) - for idx, e in enumerate(candidate_edges): - # get a random vector (b) - b = np.random.rand(num_nodes - 1) - - # set the AUV to be the edge - set_incidence_vector_for_edge_inplace(auv, e, num_nodes) - lap_update = np.outer(auv, auv) * e.weight - laplacian += lap_update - - # update the factorization - update_cholesky_factorization_inplace( - laplacian_factorization, e, num_nodes, reduced=True, subtract=False - ) - - # verify that the solutions are the same - verify_solutions(laplacian, laplacian_factorization, b) - pass - pass - pass - pass - - def test_cholesky_fiedler(self): - """Test that the Fiedler vector and value computed via Cholesky is correct. To - do this we will compare to the Fiedler vector and value computed - using LU decomposition. - """ - # Get the Petersen graph - fixed_edges, candidate_edges, n = get_split_petersen_graph() - L = weight_graph_lap_from_edge_list(fixed_edges, n) - l2, fied = find_fiedler_pair(L) - l2_test, fied_test = find_fiedler_pair_cholesky(L, x=None, normalized=False, tol=1e-8, seed=np.random.RandomState(7)) - # Check that the Fiedler vector and value are the same for the - # fixed graph. - self.assertAlmostEqual(l2, l2_test) - self.assertTrue(np.allclose(fied, fied_test)) - for e in candidate_edges: - with self.subTest(e=e): - L += e.weight * weight_graph_lap_from_edge_list([e], n) - l2, fied = find_fiedler_pair(L) - l2_test, fied_test = find_fiedler_pair_cholesky(L, x=None, normalized=False, tol=1e-8, seed=np.random.RandomState(7)) - # Check that the Fiedler vector and value are the same for the - # new graph. - self.assertAlmostEqual(l2, l2_test) - self.assertTrue(np.allclose(fied, fied_test)) - pass - pass - pass - - def test_incremental_cholesky_fiedler(self): - """Test that the Fiedler vector and value computed via Cholesky is correct after incremental updates. - - To do this we will compare to the Fiedler - vector and value computed using LU decomposition. - - """ - # Get the Petersen graph - fixed_edges, candidate_edges, n = get_split_petersen_graph() - L = weight_graph_lap_from_edge_list(fixed_edges, n) - l2, fied = find_fiedler_pair(L) - - # Make a Cholesky Eig Solver - chol_fiedler = CholeskyFiedlerSolver(L, normalized=False, tol=1e-8, seed=np.random.RandomState(7)) - l2_test, fied_test = chol_fiedler.find_fiedler_pair() - # Check that the Fiedler vector and value are the same for the - # fixed graph. - self.assertAlmostEqual(l2, l2_test) - self.assertTrue(np.allclose(fied, fied_test)) - for e in candidate_edges[0:1]: - with self.subTest(e=e): - L += e.weight * weight_graph_lap_from_edge_list([e], n) - l2, fied = find_fiedler_pair(L) - - chol_fiedler.add_edge(e) - l2_test, fied_test = chol_fiedler.find_fiedler_pair() - # Check that the Fiedler vector and value are the same for the - # new graph. - self.assertAlmostEqual(l2, l2_test) - self.assertTrue(np.allclose(fied, fied_test)) - pass - pass - pass - - -if __name__ == "__main__": - unittest.main() diff --git a/test/test_greedy_eig.py b/test/test_greedy_eig.py deleted file mode 100644 index 7a67a57..0000000 --- a/test/test_greedy_eig.py +++ /dev/null @@ -1,116 +0,0 @@ -""" -Copyright 2023 MIT Marine Robotics Group - -Tests for GreedyEig - -Author: Kevin Doherty -""" - -import unittest -import numpy as np -from scipy.sparse import spmatrix -import networkx as nx - -from mac.greedy_eig_minimal import MinimalGreedyEig -from mac.greedy_eig import GreedyEig -from mac.utils import select_edges, split_edges, nx_to_mac, mac_to_nx -from .utils import get_split_petersen_graph, get_split_erdos_renyi_graph - - -class TestGreedyEig(unittest.TestCase): - """ - Tests for the GreedyEig class - """ - - def test_petersen(self): - """ - Test the Petersen graph - """ - # Get the Petersen graph - fixed_edges, candidate_edges, n = get_split_petersen_graph() - - # Construct and solve a MinimalGreedyEig instance. This is an - # implementation of GreedyEig that naively computes the Fiedler value - # from scratch each time. It does not scale, but is useful for testing. - minimal = MinimalGreedyEig(fixed_edges, candidate_edges, n) - - # Construct and solve a GreedyEig instance. This is our fancy - # implementation using Cholesky magic. - greedy_eig = GreedyEig(fixed_edges, candidate_edges, n) - - percentages = [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9] - for pct_candidates in percentages: - with self.subTest(pct_candidates=pct_candidates): - num_candidates = int(pct_candidates * len(candidate_edges)) - - minimal_result, minimal_edges = minimal.subset(num_candidates) - eig_result, eig_edges = greedy_eig.subset(num_candidates) - - # Result from both methods must be identical - self.assertTrue( - np.allclose(minimal_result, eig_result), - msg=f"""GreedyEig result must match the result from MinimalGreedyEig. \n - Actual (GreedyEig): {eig_result} \n - Expected (Minimal): {minimal_result}""", - ) - self.assertTrue(minimal_edges == eig_edges, - msg=f"""GreedyEig edges must match the edges - from MinimalGreedyEig. \n - minimal_edges: {minimal_edges} \n - eig_edges: {eig_edges}""") - - @unittest.skip("This test is excessively slow") - def test_erdos_renyi(self): - """ - Test a few Erdos-Renyi graphs - """ - # TODO: maybe look into seeding so tests are deterministic - num_nodes_list = [20] - percent_connected_list = [0.1, 0.4] - for num_nodes in num_nodes_list: - for percent_connected in percent_connected_list: - with self.subTest(num_nodes=num_nodes, percent_connected=percent_connected): - # Get the Erdos-Renyi graph - fixed_edges, candidate_edges, n = get_split_erdos_renyi_graph( - num_nodes, percent_connected - ) - - # Construct and solve a GreedyEig instance. This is our fancy - # implementation using Cholesky magic. - greedy_eig = GreedyEig(fixed_edges, candidate_edges, n) - - # Construct and solve a MinimalGreedyESP instance. This is an - # implementation of GreedyESP that naively computes effective - # resistances using a dense pseudoinverse. It does not scale, - # but is useful for our tests. - minimal = MinimalGreedyEig(fixed_edges, candidate_edges, n) - - percentages = [0.1, 0.3, 0.6, 0.8] - for pct_candidates in percentages: - with self.subTest(pct_candidates=pct_candidates): - num_candidates = int(pct_candidates * len(candidate_edges)) - eig_result, eig_edges = greedy_eig.subset(num_candidates) - minimal_result, minimal_edges = minimal.subset(num_candidates) - - # Result from both methods must be identical - self.assertTrue( - np.allclose(eig_result, minimal_result), - msg=f"""GreedyEig result must match the result from MinimalGreedyEig. \n - Actual (GreedyEig): {eig_result} \n - Expected (Minimal): {minimal_result}""", - ) - - self.assertTrue(eig_edges == minimal_edges, - msg=f"""GreedyEig edges must match - the edges from MinimalGreedyEig. \n - minimal_edges: {minimal_edges} \n - eig_edges: {eig_edges}""") - pass - pass - pass - pass - pass - pass - -if __name__ == "__main__": - unittest.main() diff --git a/test/test_greedy_esp.py b/test/test_greedy_esp.py deleted file mode 100644 index a31f8c5..0000000 --- a/test/test_greedy_esp.py +++ /dev/null @@ -1,146 +0,0 @@ -""" -Copyright 2023 MIT Marine Robotics Group - -Regression tests for GreedyESP - -Author: Kevin Doherty -""" - -import unittest -import numpy as np -from scipy.sparse import spmatrix -import networkx as nx - -from mac.greedy_esp_minimal import MinimalGreedyESP -from mac.greedy_esp import GreedyESP -from mac.utils import select_edges, split_edges, nx_to_mac, mac_to_nx -from .utils import get_split_petersen_graph, get_split_erdos_renyi_graph - - -class TestGreedyESP(unittest.TestCase): - # Test for regressions in MAC performance - - def test_petersen(self): - """ - Test the Petersen graph - """ - # Get the Petersen graph - fixed_edges, candidate_edges, n = get_split_petersen_graph() - - # Construct and solve a MinimalGreedyESP instance. This is an - # implementation of GreedyESP that naively computes effective - # resistances using a dense pseudoinverse. It does not scale, - # but is useful for our tests. - minimal = MinimalGreedyESP(fixed_edges, candidate_edges, n, use_reduced_laplacian=True) - - # Construct and solve a GreedyESP instance. This is our fancy - # implementation using Cholesky magic. - esp = GreedyESP(fixed_edges, candidate_edges, n) - - percentages = [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9] - for pct_candidates in percentages: - with self.subTest(pct_candidates=pct_candidates): - num_candidates = int(pct_candidates * len(candidate_edges)) - - minimal_result, minimal_edges = minimal.subset(num_candidates) - esp_result, esp_edges = esp.subset(num_candidates) - - # Result from both methods must be identical - self.assertTrue( - np.allclose(minimal_result, esp_result), - msg=f"""GreedyESP result must match the result from MinimalGreedyESP. \n - Actual (GreedyESP): {esp_result} \n - Expected (Minimal): {minimal_result}""", - ) - self.assertTrue(minimal_edges == esp_edges, - msg=f"""GreedyESP edges must match the edges - from MinimalGreedyESP. \n - minimal_edges: {minimal_edges} \n - esp_edges: {esp_edges}""") - - @unittest.skip("This test is excessively slow") - def test_erdos_renyi(self): - """ - Test a few Erdos-Renyi graphs - """ - # TODO: maybe look into seeding so tests are deterministic - num_nodes_list = [20] - percent_connected_list = [0.1, 0.4] - for num_nodes in num_nodes_list: - for percent_connected in percent_connected_list: - with self.subTest(num_nodes=num_nodes, percent_connected=percent_connected): - # Get the Erdos-Renyi graph - fixed_edges, candidate_edges, n = get_split_erdos_renyi_graph( - num_nodes, percent_connected - ) - - # Construct and solve a GreedyESP instance. This is our fancy - # implementation using Cholesky magic. - esp = GreedyESP(fixed_edges, candidate_edges, n) - - # Construct and solve a MinimalGreedyESP instance. This is an - # implementation of GreedyESP that naively computes effective - # resistances using a dense pseudoinverse. It does not scale, - # but is useful for our tests. - minimal = MinimalGreedyESP(fixed_edges, candidate_edges, n, use_reduced_laplacian=True) - - percentages = [0.1, 0.3, 0.6, 0.8] - for pct_candidates in percentages: - with self.subTest(pct_candidates=pct_candidates): - num_candidates = int(pct_candidates * len(candidate_edges)) - esp_result, esp_edges = esp.subset(num_candidates) - minimal_result, minimal_edges = minimal.subset(num_candidates) - - # Result from both methods must be identical - self.assertTrue( - np.allclose(esp_result, minimal_result), - msg=f"""GreedyESP result must match the result from MinimalGreedyESP. \n - Actual (GreedyESP): {esp_result} \n - Expected (Minimal): {minimal_result}""", - ) - - self.assertTrue(esp_edges == minimal_edges, - msg=f"""GreedyESP edges must match - the edges from MinimalGreedyESP. \n - minimal_edges: {minimal_edges} \n - esp_edges: {esp_edges}""") - - - def test_greedy_esp_is_stateless(self): - # Get the Petersen graph - fixed_edges, candidate_edges, n = get_split_petersen_graph() - esp = GreedyESP(fixed_edges, candidate_edges, n) - - # ensure that after solving problems none of the class variables are - # changed - - # copy the class variables - original_class_variables = esp.__dict__.copy() - pct_candidates = 0.4 - num_candidates = int(pct_candidates * len(candidate_edges)) - esp.subset(num_candidates) - - # ensure that the values of the class variables are unchanged - for key, value in esp.__dict__.items(): - with self.subTest(key=key, value=value): - orig_value = original_class_variables[key] - - # ensure that the values are the same type - self.assertTrue( - type(value) == type(orig_value), - msg=f"""Class variable {key} has changed type. \n - Actual: {type(value)} \n - Expected: {type(orig_value)}""", - ) - - if isinstance(value, np.ndarray): - self.assertTrue(np.allclose(value, orig_value)) - elif isinstance(value, spmatrix): - self.assertTrue(isinstance(orig_value, spmatrix)) - self.assertTrue(np.allclose(value.toarray(), orig_value.toarray())) - else: - self.assertEqual(value, orig_value) - - -if __name__ == "__main__": - unittest.main() diff --git a/test/test_greedy_esp_minimal.py b/test/test_greedy_esp_minimal.py deleted file mode 100644 index ec971f7..0000000 --- a/test/test_greedy_esp_minimal.py +++ /dev/null @@ -1,100 +0,0 @@ -""" -Copyright 2023 MIT Marine Robotics Group - -Regression tests to ensure that GreedyESP returns the same result whether or not -we are using the reduced Laplacian - -Author: Alan Papalia -""" -import unittest -import numpy as np -import networkx as nx - -from mac.greedy_esp_minimal import MinimalGreedyESP -from mac.greedy_esp import GreedyESP -from mac.utils import select_edges, split_edges, nx_to_mac, mac_to_nx, get_incidence_vector -from .utils import get_split_petersen_graph - - -class TestGreedyEspMinimal(unittest.TestCase): - - def test_petersen(self): - """ - Test the Petersen graph - """ - # Get the Petersen graph - fixed_edges, candidate_edges, n = get_split_petersen_graph() - - greedy_minimal_reduced = MinimalGreedyESP( - fixed_edges, candidate_edges, n, use_reduced_laplacian=True - ) - greedy_minimal_not_reduced = MinimalGreedyESP( - fixed_edges, candidate_edges, n, use_reduced_laplacian=False - ) - - for pct_candidates in [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9]: - with self.subTest(pct_candidates=pct_candidates): - num_candidates = int(pct_candidates * len(candidate_edges)) - - # Construct and solve a MinimalGreedyESP instance. This is an - # implementation of GreedyESP that naively computes effective - # resistances using a dense pseudoinverse. It does not scale, - # but is useful for our tests. - minimal_reduced_result, reduced_edges_selected = greedy_minimal_reduced.subset(num_candidates) - minimal_not_reduced_result, not_reduced_edges_selected = greedy_minimal_not_reduced.subset( - num_candidates - ) - - # Result from both methods must be identical - self.assertTrue( - np.allclose(minimal_reduced_result, minimal_not_reduced_result), - msg=f"""GreedyMinimalReduced result must match the result from GreedyMinimalNotReduced. \n - (GreedyMinimalNotReduced): {minimal_not_reduced_result} \n - (GreedyMinimalReduced): {minimal_reduced_result}""", - ) - self.assertTrue( - reduced_edges_selected == not_reduced_edges_selected, - msg=f"""GreedyMinimalReduced edges selected must match the - result from GreedyMinimalNotReduced. \n - (GreedyMinimalNotReduced): {not_reduced_edges_selected} \n - (GreedyMinimalReduced): {reduced_edges_selected}""", - ) - - def test_reff_computation(self): - # Get the Petersen graph - fixed_edges, candidate_edges, n = get_split_petersen_graph() - - greedy_minimal_reduced = MinimalGreedyESP( - fixed_edges, candidate_edges, n, use_reduced_laplacian=True - ) - greedy_minimal_not_reduced = MinimalGreedyESP( - fixed_edges, candidate_edges, n, use_reduced_laplacian=False - ) - - # make sure they have the same edge list so we can iterate over just one - # of them - assert (greedy_minimal_reduced.edge_list == greedy_minimal_not_reduced.edge_list).all() - - # run the Reff computation for 10 random assignments of weights and all - # of the edges - num_trials = 10 - edge_list = greedy_minimal_reduced.edge_list - for _ in range(num_trials): - rand_weights = np.random.rand(len(greedy_minimal_reduced.edge_list)) - - for edge in edge_list: - auv = get_incidence_vector(edge, n) - reff_reduced = greedy_minimal_reduced.compute_reff(rand_weights, auv) - reff_not_reduced = greedy_minimal_not_reduced.compute_reff(rand_weights, auv) - - self.assertTrue(np.allclose(reff_reduced, reff_not_reduced), - msg=f"""Reff computation must be the same for both - methods.\n - Reff (Reduced): {reff_reduced} \n - Reff (Not Reduced): {reff_not_reduced}""") - - - - -if __name__ == "__main__": - unittest.main() \ No newline at end of file diff --git a/test/test_mac_regression.py b/test/test_mac_regression.py deleted file mode 100644 index 957c307..0000000 --- a/test/test_mac_regression.py +++ /dev/null @@ -1,108 +0,0 @@ -""" -Copyright 2023 MIT Marine Robotics Group - -Regression tests for MAC - -Author: Kevin Doherty -""" - -import unittest -import numpy as np -import networkx as nx - -from mac.mac import MAC -from mac.utils import select_edges, split_edges, nx_to_mac, mac_to_nx -from .utils import get_split_petersen_graph - -class TestMACRegression(unittest.TestCase): - # Test for regressions in MAC performance - # This test is based on the following: - # - Algebraic connectivity for standard graphs - # - - - def test_petersen(self): - """ - Test the Petersen graph - """ - # Get the Petersen graph - fixed_edges, candidate_edges, n = get_split_petersen_graph() - - for pct_candidates in [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9]: - with self.subTest(pct_candidates=pct_candidates): - num_candidates = int(pct_candidates * len(candidate_edges)) - - # Construct an initial guess - w_init = np.zeros(len(candidate_edges)) - w_init[:num_candidates] = 1.0 - - mac = MAC(fixed_edges, candidate_edges, n) - result, unrounded, upper = mac.fw_subset(w_init, num_candidates, max_iters=100) - init_selected = select_edges(candidate_edges, w_init) - selected = select_edges(candidate_edges, result) - - init_l2 = mac.evaluate_objective(w_init) - mac_unrounded_l2 = mac.evaluate_objective(unrounded) - mac_l2 = mac.evaluate_objective(result) - - print("Initial L2: {}".format(init_l2)) - print("MAC Unrounded L2: {}".format(mac_unrounded_l2)) - print("MAC Rounded L2: {}".format(mac_l2)) - - self.assertGreater(mac_unrounded_l2, init_l2, msg="""MAC unrounded connectivity - should be greater than initial guess""") - - # NOTE: Stubbed out test. This assertion is not necessarily - # (mathematically) going to be true. In order to guarantee - # this, we would need to add a "fallback" check in the MAC - # solver that ensures that rounding does not result in a - # solution worse than the initial In that case, we would simply - # select the initial guess. - - # self.assertGreater(mac_l2, init_l2, msg="""MAC rounded - # connectivity should be greater than initial guess""") - pass - - @unittest.skip("Skipping, since this test will not pass [slight differences in eigvecs can change output]") - def test_cache(self): - # Get the Petersen graph - fixed_edges, candidate_edges, n = get_split_petersen_graph() - - for pct_candidates in [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9]: - with self.subTest(pct_candidates=pct_candidates): - num_candidates = int(pct_candidates * len(candidate_edges)) - - # Construct an initial guess - w_init = np.zeros(len(candidate_edges)) - w_init[:num_candidates] = 1.0 - - mac = MAC(fixed_edges, candidate_edges, n) - mac_cache = MAC(fixed_edges, candidate_edges, n, use_cache=True) - - result, unrounded, upper = mac.fw_subset(w_init, num_candidates, max_iters=100) - result_cache, unrounded_cache, upper_cache = mac_cache.fw_subset(w_init, num_candidates, max_iters=100) - - self.assertTrue(np.allclose(unrounded, unrounded_cache), msg=f"""Cached MAC unrounded - result should be the same as non-cached result\n MAC: {unrounded} \n Cache: {unrounded_cache}""") - - self.assertTrue(np.allclose(result, result_cache), msg=f"""Cached MAC result - should be the same as non-cached result\n MAC: {result} \n Cache: {result_cache}""") - - self.assertTrue(np.allclose(upper, upper_cache), msg=f"""Cached MAC upper - result should be the same as non-cached result\n MAC: {upper} \n Cache: {upper_cache}""") - - mac_unrounded_l2 = mac.evaluate_objective(unrounded) - mac_l2 = mac.evaluate_objective(result) - - mac_cache_unrounded_l2 = mac.evaluate_objective(unrounded_cache) - mac_cache_l2 = mac.evaluate_objective(result_cache) - - print("MAC Unrounded L2: {}".format(mac_unrounded_l2)) - print("MAC Rounded L2: {}".format(mac_l2)) - - print("Cache MAC Unrounded L2: {}".format(mac_cache_unrounded_l2)) - print("Cache MAC Rounded L2: {}".format(mac_cache_l2)) - - pass - -if __name__ == '__main__': - unittest.main() diff --git a/test/utils.py b/test/utils.py deleted file mode 100644 index 8bb520a..0000000 --- a/test/utils.py +++ /dev/null @@ -1,46 +0,0 @@ -import numpy as np -import networkx as nx -from typing import List, Tuple - -from mac.utils import select_edges, split_edges, nx_to_mac, mac_to_nx, Edge - -def get_split_petersen_graph() -> Tuple[List[Edge], List[Edge], int]: - """ - Get a "split" Petersen graph for testing purposes. - Returns: - fixed_edges: a chain of edges in the graph. - candidate_edges: the remaining edges in the graph. - n: the number of nodes in the graph. - """ - G = nx.petersen_graph() - n = len(G.nodes()) - - # Add a chain - for i in range(n-1): - if G.has_edge(i+1, i): - G.remove_edge(i+1, i) - pass - if not G.has_edge(i, i+1): - G.add_edge(i, i+1) - pass - pass - - edges = nx_to_mac(G) - - # Split chain and non-chain parts - fixed_edges , candidate_edges = split_edges(edges) - - return fixed_edges, candidate_edges, n - -def get_split_erdos_renyi_graph(n: int=20, p: float = 0.30) -> Tuple[List[Edge], List[Edge], int]: - G = nx.erdos_renyi_graph(n, p) - - for i in range(n-1): - if G.has_edge(i+1, i): - G.remove_edge(i+1, i) - if not G.has_edge(i, i+1): - G.add_edge(i, i+1) - - edges = nx_to_mac(G) - fixed_edges, candidate_edges = split_edges(edges) - return fixed_edges, candidate_edges, n \ No newline at end of file diff --git a/test/__init__.py b/tests/README.md similarity index 100% rename from test/__init__.py rename to tests/README.md diff --git a/tests/__init__.py b/tests/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/tests/benchmarks/__init__.py b/tests/benchmarks/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/tests/benchmarks/test_cache_performance.py b/tests/benchmarks/test_cache_performance.py new file mode 100644 index 0000000..560664f --- /dev/null +++ b/tests/benchmarks/test_cache_performance.py @@ -0,0 +1,49 @@ +import numpy as np +import networkx as nx +import pytest + +# MAC requirements +from mac.solvers.mac import MAC +from mac.utils.conversions import nx_to_mac + +@pytest.mark.parametrize("k", [5]) +def test_benchmark_solver(benchmark, k): + # Load graph + G = nx.petersen_graph() + n = len(G.nodes()) + + # Split the graph into a tree part and a "loop" part + spanning_tree = nx.minimum_spanning_tree(G) + loop_graph = nx.difference(G, spanning_tree) + + # Create solver + solver = MAC(fixed_edges=nx_to_mac(spanning_tree), + candidate_edges=nx_to_mac(loop_graph), num_nodes=n) + + # Set up initial guess + x_init = np.zeros(len(loop_graph.edges())) + x_init[:k] = 1.0 + + result = benchmark.pedantic(solver.solve, args=[k], kwargs={"x_init": x_init, + "use_cache": False}, rounds=5) + +@pytest.mark.parametrize("k", [5]) +def test_benchmark_solver_cache(benchmark, k): + # Load graph + G = nx.petersen_graph() + n = len(G.nodes()) + + # Split the graph into a tree part and a "loop" part + spanning_tree = nx.minimum_spanning_tree(G) + loop_graph = nx.difference(G, spanning_tree) + + # Create solver + solver = MAC(fixed_edges=nx_to_mac(spanning_tree), + candidate_edges=nx_to_mac(loop_graph), num_nodes=n) + + # Set up initial guess + x_init = np.zeros(len(loop_graph.edges())) + x_init[:k] = 1.0 + + result = benchmark.pedantic(solver.solve, args=[k], kwargs={"x_init": x_init, + "use_cache": True}, rounds=5) diff --git a/tests/optimization/__init__.py b/tests/optimization/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/tests/optimization/test_frankwolfe.py b/tests/optimization/test_frankwolfe.py new file mode 100644 index 0000000..78a60ce --- /dev/null +++ b/tests/optimization/test_frankwolfe.py @@ -0,0 +1,76 @@ +""" +Copyright 2023 MIT Marine Robotics Group + +Tests for Frank-Wolfe algorithm + +Author: Kevin Doherty +""" + +import unittest +import numpy as np +import networkx as nx + +from mac.utils.conversions import nx_to_mac +from mac.utils.graphs import weight_graph_lap_from_edge_list + +from mac.optimization.constraints import * + +# Code under test. +from mac.optimization.frankwolfe import * + +from networkx.linalg.algebraicconnectivity import algebraic_connectivity + +class TestSimpleQuadratic(unittest.TestCase): + def test_solve_box_constraint(self): + """ + Solve min -x^2 s.t. x_i in [0,1]. Trivially solved by x = 0. + """ + # Create a simple concave objective function f(x) = -x^2 with gradf(x) = -2x. + problem = lambda x: (-np.inner(x,x), -2*x) + # Solve + solve_lp = solve_box_lp + N = 10 + initial = 0.5 * np.ones(N) + x, u = frank_wolfe(initial, problem, solve_lp) + self.assertTrue(np.allclose(x, np.zeros(N))) + + def test_solve_subset_box_constraint(self): + """ + Solve min -x^2 s.t. x_i in [0,1], sum_i x_i = 1 + """ + problem = lambda x: (-np.inner(x,x), -2*x) + k = 1 + solve_lp = lambda g: solve_subset_box_lp(g, k) + N = 2 + initial = np.random.rand(N) + initial = (k / np.sum(initial)) * initial + + x, u = frank_wolfe(initial, problem, solve_lp) + expected = (k / N) * np.ones(N) + + # We're willing to accept a pretty loose tolerance here. + self.assertTrue(np.allclose(x, expected, atol=0.01)) + + def test_convergence_around_zero(self): + """ + Ensure that we handle the case where f(x) \approx 0 correctly (e.g. avoid + division by zero). + """ + # Create a simple concave objective function f(x) = -x^2 + 0.25 with gradf(x) = -2x. + problem = lambda x: (-np.inner(x,x) + 0.25, -2*x) + + # Set up box constraints for the solver + solve_lp = solve_box_lp + + # Set up an initialization with f(initial) = 0 + N = 10 + initial = np.zeros(N) + initial[0] = 0.5 + + # Solve + x, u = frank_wolfe(initial, problem, solve_lp) + + self.assertTrue(np.allclose(x, np.zeros(N))) + +if __name__ == '__main__': + unittest.main() diff --git a/tests/solvers/__init__.py b/tests/solvers/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/tests/solvers/test_mac.py b/tests/solvers/test_mac.py new file mode 100644 index 0000000..51e7f5f --- /dev/null +++ b/tests/solvers/test_mac.py @@ -0,0 +1,75 @@ +""" +Copyright 2023 MIT Marine Robotics Group + +Regression tests for MAC + +Author: Kevin Doherty +""" + +import unittest +import numpy as np +import networkx as nx + +from mac.utils.conversions import nx_to_mac +from mac.utils.graphs import select_edges + +# Code under test +from mac.solvers.mac import MAC + +class TestPetersenGraphConnected(unittest.TestCase): + """ + Test sparsification on the Petersen graph with a connected "fixed" base graph. + """ + def setUp(self): + graph = nx.petersen_graph() + # Split the graph into a tree part and a "loop" part + spanning_tree = nx.minimum_spanning_tree(graph) + loop_graph = nx.difference(graph, spanning_tree) + self.fixed_edges = nx_to_mac(spanning_tree) + self.candidate_edges = nx_to_mac(loop_graph) + self.n = graph.number_of_nodes() + + # Test for regressions in MAC performance + # This test is based on the following: + # - Algebraic connectivity for standard graphs + def test_petersen(self): + """ + Test the Petersen graph + """ + for pct_candidates in [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9]: + with self.subTest(pct_candidates=pct_candidates): + num_candidates = int(pct_candidates * len(self.candidate_edges)) + + # Construct an initial guess + x_init = np.zeros(len(self.candidate_edges)) + x_init[:num_candidates] = 1.0 + + mac = MAC(self.fixed_edges, self.candidate_edges, self.n) + result, unrounded, upper = mac.solve(num_candidates, x_init, max_iters=100) + init_selected = select_edges(self.candidate_edges, x_init) + selected = select_edges(self.candidate_edges, result) + + init_l2 = mac.evaluate_objective(x_init) + mac_unrounded_l2 = mac.evaluate_objective(unrounded) + mac_l2 = mac.evaluate_objective(result) + + print("Initial L2: {}".format(init_l2)) + print("MAC Unrounded L2: {}".format(mac_unrounded_l2)) + print("MAC Rounded L2: {}".format(mac_l2)) + + self.assertGreaterEqual(mac_unrounded_l2, init_l2, msg="""MAC unrounded connectivity + should be greater than initial guess""") + + # NOTE: Stubbed out test. This assertion is not necessarily + # (mathematically) going to be true. In order to guarantee + # this, we would need to add a "fallback" check in the MAC + # solver that ensures that rounding does not result in a + # solution worse than the initial In that case, we would simply + # select the initial guess. + + # self.assertGreater(mac_l2, init_l2, msg="""MAC rounded + # connectivity should be greater than initial guess""") + pass + +if __name__ == '__main__': + unittest.main() diff --git a/tests/utils/__init__.py b/tests/utils/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/tests/utils/test_conversions.py b/tests/utils/test_conversions.py new file mode 100644 index 0000000..0af620a --- /dev/null +++ b/tests/utils/test_conversions.py @@ -0,0 +1,67 @@ +""" +Copyright 2023 MIT Marine Robotics Group + +Tests for graph type conversion utilities + +Author: Kevin Doherty +""" + +import unittest +import numpy as np +import networkx as nx + +# Code under test +from mac.utils.conversions import * + +def edges_equal(x: List[Edge], y: List[Edge]) -> bool: + """ + Check for equality of two edge lists in the "MAC" format. + """ + x.sort(key=lambda e: (e.i, e.j)) + y.sort(key=lambda e: (e.i, e.j)) + return x == y + +class TestConversionsUtils(unittest.TestCase): + def setUp(self): + self.petersen_nx = nx.petersen_graph() + # Manually enumerate the 15 edges of the Petersen graph + # (cf. https://networkx.org/documentation/stable/_modules/networkx/generators/small.html#petersen_graph) + self.petersen_mac = [Edge(0, 1, weight=1), + Edge(0, 4, weight=1), + Edge(0, 5, weight=1), + Edge(1, 2, weight=1), + Edge(1, 6, weight=1), + Edge(2, 3, weight=1), + Edge(2, 7, weight=1), + Edge(3, 4, weight=1), + Edge(3, 8, weight=1), + Edge(4, 9, weight=1), + Edge(5, 7, weight=1), + Edge(5, 8, weight=1), + Edge(6, 8, weight=1), + Edge(6, 9, weight=1), + Edge(7, 9, weight=1)] + + self.weighted_graph = nx.petersen_graph() + rng = np.random.default_rng(seed=7) + weights = rng.random(self.weighted_graph.number_of_edges()) + for i, e in enumerate(self.weighted_graph.edges()): + self.weighted_graph[e[0]][e[1]]['weight'] = weights[i] + + def test_mac_to_mac_via_nx(self): + mac_to_mac_via_nx = nx_to_mac(mac_to_nx(self.petersen_mac)) + self.assertTrue(edges_equal(mac_to_mac_via_nx, self.petersen_mac)) + + def test_nx_to_mac_equals_mac(self): + self.assertTrue(edges_equal(nx_to_mac(self.petersen_nx), self.petersen_mac)) + + def test_nx_to_nx_via_mac_weighted(self): + nx_to_nx_via_mac = mac_to_nx(nx_to_mac(self.weighted_graph.copy())) + + # We need to overwrite the name of the graph or the NetworkX graphs_equal utility + # will not recognize these graphs as equal. + nx_to_nx_via_mac.name = self.weighted_graph.name + self.assertTrue(nx.utils.graphs_equal(nx_to_nx_via_mac, self.weighted_graph)) + +if __name__ == '__main__': + unittest.main() diff --git a/tests/utils/test_fiedler.py b/tests/utils/test_fiedler.py new file mode 100644 index 0000000..9b2296c --- /dev/null +++ b/tests/utils/test_fiedler.py @@ -0,0 +1,53 @@ +""" +Copyright 2023 MIT Marine Robotics Group + +Tests for Fiedler utilities + +Author: Kevin Doherty +""" + +import unittest +import numpy as np +import networkx as nx + +from mac.utils.conversions import nx_to_mac +from mac.utils.graphs import weight_graph_lap_from_edge_list + +# Code under test. +from mac.utils.fiedler import * + +from networkx.linalg.algebraicconnectivity import algebraic_connectivity + +class TestConnectedGraphs(unittest.TestCase): + def setUp(self): + # Complete graph on 5 nodes + self.complete_graph = nx.complete_graph(5) + + def test_fiedler(self): + # The algebraic connectivity of the complete graph K(N) on N nodes is + # exactly equal to N. + edge_list = nx_to_mac(self.complete_graph) + N = self.complete_graph.number_of_nodes() + L = weight_graph_lap_from_edge_list(edge_list, N) + fiedler_value, fiedler_vec, _ = find_fiedler_pair(L) + self.assertTrue(np.isclose(fiedler_value, N)) + +class TestDisconnectedGraphs(unittest.TestCase): + def setUp(self): + # Construct a disconnected graph with two connected components + component_size = 3 + self.disconnected_graph = nx.complete_graph(component_size) + self.disconnected_graph.add_edges_from( + (u,v) for u in range(component_size, 2*component_size) for v in range(u+1, 2*component_size)) + + @unittest.skip("Feature not yet supported.") + def test_fiedler(self): + # Test multiple connected components + edge_list = nx_to_mac(self.disconnected_graph) + N = self.disconnected_graph.number_of_nodes() + L = weight_graph_lap_from_edge_list(edge_list, N) + fiedler_value, fiedler_vec, _ = find_fiedler_pair(L) + self.assertTrue(np.isclose(fiedler_value, 0)) + +if __name__ == '__main__': + unittest.main() diff --git a/tests/utils/test_graphs.py b/tests/utils/test_graphs.py new file mode 100644 index 0000000..9da1dc4 --- /dev/null +++ b/tests/utils/test_graphs.py @@ -0,0 +1,53 @@ +""" +Copyright 2023 MIT Marine Robotics Group + +Tests for graph utilities. + +Author: Kevin Doherty +""" + +import unittest +import networkx as nx +import numpy as np + +from mac.utils.conversions import nx_to_mac + +# Code under test +from mac.utils.graphs import * + +class TestGraphUtils(unittest.TestCase): + def setUp(self): + self.graph = nx.petersen_graph() + self.weighted_graph = nx.petersen_graph() + rng = np.random.default_rng(seed=7) + weights = rng.random(self.weighted_graph.number_of_edges()) + for i, e in enumerate(self.weighted_graph.edges()): + self.weighted_graph[e[0]][e[1]]['weight'] = weights[i] + + def test_weight_graph_lap_from_edge_list_unweighted(self): + edge_list = nx_to_mac(self.graph) + L_ours = weight_graph_lap_from_edge_list(edge_list, self.graph.number_of_nodes()) + L_nx = nx.laplacian_matrix(self.graph) + self.assertTrue(np.allclose(L_ours.todense(), L_nx.todense())) + + def test_weight_graph_lap_from_edge_list_weighted(self): + edge_list = nx_to_mac(self.weighted_graph) + L_ours = weight_graph_lap_from_edge_list(edge_list, self.weighted_graph.number_of_nodes()) + L_nx = nx.laplacian_matrix(self.weighted_graph) + self.assertTrue(np.allclose(L_ours.todense(), L_nx.todense())) + + def test_weight_graph_lap_from_edges(self): + edge_list = nx_to_mac(self.weighted_graph) + edges = [] + weights = [] + for e in edge_list: + edges.append((e.i, e.j)) + weights.append(e.weight) + edges = np.array(edges) + weights = np.array(weights) + L_ours = weight_graph_lap_from_edges(edges, weights, self.weighted_graph.number_of_nodes()) + L_nx = nx.laplacian_matrix(self.weighted_graph) + self.assertTrue(np.allclose(L_ours.todense(), L_nx.todense())) + +if __name__ == '__main__': + unittest.main()