From dfdd6d323ae0f77c19440ff87e30c94e768cc6d5 Mon Sep 17 00:00:00 2001 From: Kevin Date: Fri, 2 Aug 2024 11:19:49 -0400 Subject: [PATCH 01/36] Clean up tests. --- mac/greedy_eig_minimal.py | 145 ------------------------- mac/greedy_esp_minimal.py | 108 ------------------ test/test_greedy_esp_minimal.py | 100 ----------------- {test => tests}/__init__.py | 0 {test => tests}/test_cholesky_utils.py | 0 {test => tests}/test_greedy_eig.py | 0 {test => tests}/test_greedy_esp.py | 0 {test => tests}/test_mac_regression.py | 0 {test => tests}/utils.py | 0 9 files changed, 353 deletions(-) delete mode 100644 mac/greedy_eig_minimal.py delete mode 100644 mac/greedy_esp_minimal.py delete mode 100644 test/test_greedy_esp_minimal.py rename {test => tests}/__init__.py (100%) rename {test => tests}/test_cholesky_utils.py (100%) rename {test => tests}/test_greedy_eig.py (100%) rename {test => tests}/test_greedy_esp.py (100%) rename {test => tests}/test_mac_regression.py (100%) rename {test => tests}/utils.py (100%) 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/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/__init__.py b/tests/__init__.py similarity index 100% rename from test/__init__.py rename to tests/__init__.py diff --git a/test/test_cholesky_utils.py b/tests/test_cholesky_utils.py similarity index 100% rename from test/test_cholesky_utils.py rename to tests/test_cholesky_utils.py diff --git a/test/test_greedy_eig.py b/tests/test_greedy_eig.py similarity index 100% rename from test/test_greedy_eig.py rename to tests/test_greedy_eig.py diff --git a/test/test_greedy_esp.py b/tests/test_greedy_esp.py similarity index 100% rename from test/test_greedy_esp.py rename to tests/test_greedy_esp.py diff --git a/test/test_mac_regression.py b/tests/test_mac_regression.py similarity index 100% rename from test/test_mac_regression.py rename to tests/test_mac_regression.py diff --git a/test/utils.py b/tests/utils.py similarity index 100% rename from test/utils.py rename to tests/utils.py From b0dbd0328b17c8497051e49a55f2e5e5ff24aa0b Mon Sep 17 00:00:00 2001 From: Kevin Date: Fri, 2 Aug 2024 13:51:34 -0400 Subject: [PATCH 02/36] Tidy up Petersen graph sparsification. --- examples/petersen_graph_sparsification.py | 49 +++++++++++------------ 1 file changed, 23 insertions(+), 26 deletions(-) diff --git a/examples/petersen_graph_sparsification.py b/examples/petersen_graph_sparsification.py index 57c80c8..8331f6e 100644 --- a/examples/petersen_graph_sparsification.py +++ b/examples/petersen_graph_sparsification.py @@ -5,7 +5,6 @@ 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 plt.rcParams['text.usetex'] = True @@ -13,50 +12,48 @@ 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) +num_candidates = int(pct_candidates * len(candidate_edges)) +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 = np.zeros(len(candidate_edges)) w_init[:num_candidates] = 1.0 -result, unrounded, upper = mac.fw_subset(w_init, num_candidates, max_iters=100) +result, unrounded, upper = mac.fw_subset(w_init, num_candidates, 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)}") From 1872b6dd705b23d774f6c0df89d0fc2a18775e12 Mon Sep 17 00:00:00 2001 From: Kevin Date: Fri, 2 Aug 2024 14:15:18 -0400 Subject: [PATCH 03/36] WIP more tidying. --- mac/mac.py | 18 +++++++++--------- mac/utils.py | 10 +++++----- 2 files changed, 14 insertions(+), 14 deletions(-) diff --git a/mac/mac.py b/mac/mac.py index d8c819b..dd28bb4 100644 --- a/mac/mac.py +++ b/mac/mac.py @@ -62,22 +62,22 @@ 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): + def laplacian(self, x): """ - Construct the combined Laplacian (fixed edges plus candidate edges weighted by w). + Construct the combined Laplacian (fixed edges plus candidate edges weighted by x). - w: An element of [0,1]^m; this is the edge selection to use + x: 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 > 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 + 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 find_fiedler_pair(self, w): """ @@ -95,7 +95,7 @@ def find_fiedler_pair(self, w): value and corresponding vector. """ - L = self.combined_laplacian(w) + L = self.laplacian(w) if L.shape[0] == 0: # If the graph is empty, then the Fiedler value is 0 and the # Fiedler vector is empty. diff --git a/mac/utils.py b/mac/utils.py index 11b156c..db8b767 100644 --- a/mac/utils.py +++ b/mac/utils.py @@ -181,18 +181,18 @@ def weight_reduced_graph_lap_from_edge_list( def weight_graph_lap_from_edges( - edges: List[Tuple[int, int]], weights: List[int], num_poses: int + edges: List[Tuple[int, int]], weights: List[int], num_nodes: int ) -> csr_matrix: - """Returns the sparse rotational weighted graph Laplacian matrix from a list + """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_poses (int): the number of poses + num_nodes (int): the number of nodes Returns: - csr_matrix: the rotational weighted graph Laplacian matrix + csr_matrix: the weighted graph Laplacian matrix """ assert len(edges) == len(weights) # Preallocate triplets @@ -220,7 +220,7 @@ def weight_graph_lap_from_edges( cols.append(edges[i, 0]) data.append(-weights[i]) - return csr_matrix(coo_matrix((data, (rows, cols)), shape=[num_poses, num_poses])) + return csr_matrix(coo_matrix((data, (rows, cols)), shape=[num_nodes, num_nodes])) def split_edges(edges: List[Edge]) -> Tuple[List[Edge], List[Edge]]: From 1c098f63bf95cfcc2109f3b172a2d4606d6cf1ec Mon Sep 17 00:00:00 2001 From: Kevin Date: Fri, 2 Aug 2024 15:03:57 -0400 Subject: [PATCH 04/36] Reorganizing MAC and Fiedler utilities. --- mac/fiedler.py | 52 +-------------------------- mac/mac.py | 98 ++++++++++++++++++++++---------------------------- 2 files changed, 43 insertions(+), 107 deletions(-) diff --git a/mac/fiedler.py b/mac/fiedler.py index bd07d46..c81b8a1 100644 --- a/mac/fiedler.py +++ b/mac/fiedler.py @@ -6,36 +6,7 @@ # 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): +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`. @@ -72,24 +43,3 @@ def find_fiedler_pair_and_eigvecs(L, X=None, method='tracemin_lu', tol=1e-8, see 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/mac.py b/mac/mac.py index dd28bb4..a109ac3 100644 --- a/mac/mac.py +++ b/mac/mac.py @@ -2,6 +2,8 @@ import mac.fiedler as fiedler import mac.frankwolfe as fw from timeit import default_timer as timer +from dataclasses import dataclass +from typing import Optional import numpy as np import networkx as nx @@ -12,6 +14,11 @@ from timeit import default_timer as timer +@dataclass +class Cache: + """Problem data cache""" + Q: Optional[np.ndarray] = None + class MAC: def __init__(self, fixed_edges, candidate_edges, num_nodes, fiedler_method='tracemin_lu', use_cache=False, fiedler_tol=1e-8, @@ -79,74 +86,48 @@ def laplacian(self, x): L_x = self.L_fixed + L_candidate return L_x - def find_fiedler_pair(self, w): - """ - 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. - - """ - L = self.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): + def evaluate_objective(self, x): """ - 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) - - returns F(w) = lambda_2(L(w)). - """ - return self.find_fiedler_pair(w)[0] - - 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. + x: Weights for each candidate edge (does not include fixed edges) - returns grad F(w) from equation (8) of our paper: https://arxiv.org/pdf/2203.13897.pdf. + returns F(x) = lambda_2(L(x)). """ - return fiedler.grad_from_fiedler(fiedler_vec, self.edge_list, self.weights) + return fiedler.find_fiedler_pair(self.laplacian(x), method=self.fiedler_method, tol=self.fiedler_tol)[0] 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 + Compute the algebraic connectivity of L(x) and a (super)gradient of the algebraic connectivity with respect to x. - def fw_subset(self, w_init, k, rounding="nearest", fallback=False, + returns x, grad F(x). + """ + Q = None if self.cache is None else self.cache.Q + fiedler_value, fiedler_vec, Qnew = fiedler.find_fiedler_pair(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 self.cache is not None: + self.cache.Q = Q + return f, gradf, cache + + 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 the Frank-Wolfe method to solve the subset selection problem,. Parameters ---------- - w_init : Array-like + x_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 @@ -195,10 +176,15 @@ def fw_subset(self, w_init, k, rounding="nearest", fallback=False, # Solution for the direction-finding subproblem solve_lp = lambda g: fw.solve_subset_box_lp(g, k) + # handle case where x is none + + # Set up problem to use cache (or not) + # problem = lambda x: self.problem() + # 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, + w, u = fw.frank_wolfe_with_recycling(initial=x_init, problem=self.problem_with_recycling, solve_lp=solve_lp, maxiter=max_iters, @@ -206,7 +192,7 @@ def fw_subset(self, w_init, k, rounding="nearest", fallback=False, grad_norm_tol=grad_norm_tol, verbose=verbose) else: - w, u = fw.frank_wolfe(initial=w_init, problem=self.problem, + w, u = fw.frank_wolfe(initial=x_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, @@ -222,7 +208,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 From debcfc86efb3334de0d6fc7d79f32d175c7dc977 Mon Sep 17 00:00:00 2001 From: Kevin Date: Fri, 2 Aug 2024 15:07:29 -0400 Subject: [PATCH 05/36] Cache setup. --- mac/mac.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/mac/mac.py b/mac/mac.py index a109ac3..519a339 100644 --- a/mac/mac.py +++ b/mac/mac.py @@ -60,7 +60,11 @@ def __init__(self, fixed_edges, candidate_edges, num_nodes, 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 + + # Set up caching. + self.cache = None + if use_cache: + self.cache = Cache(Q=None) # Configuration for Fiedler vector computation self.fiedler_method = fiedler_method From 9e8fc4915c148ef0963f02d8bc9efd6b36e0b5f2 Mon Sep 17 00:00:00 2001 From: Kevin Date: Fri, 2 Aug 2024 16:09:28 -0400 Subject: [PATCH 06/36] Full project reorg. --- mac/frankwolfe.py | 183 ------------------- mac/optimization/constraints.py | 29 +++ mac/optimization/frankwolfe.py | 78 ++++++++ mac/{ => solvers}/baseline.py | 0 mac/{ => solvers}/greedy_eig.py | 0 mac/{ => solvers}/greedy_esp.py | 0 mac/{ => solvers}/mac.py | 30 ++- mac/{cholesky_utils.py => utils/cholesky.py} | 0 mac/utils/conversions.py | 35 ++++ mac/{ => utils}/fiedler.py | 0 mac/{utils.py => utils/graphs.py} | 136 +------------- mac/utils/rounding.py | 93 ++++++++++ 12 files changed, 251 insertions(+), 333 deletions(-) delete mode 100644 mac/frankwolfe.py create mode 100644 mac/optimization/constraints.py create mode 100644 mac/optimization/frankwolfe.py rename mac/{ => solvers}/baseline.py (100%) rename mac/{ => solvers}/greedy_eig.py (100%) rename mac/{ => solvers}/greedy_esp.py (100%) rename mac/{ => solvers}/mac.py (91%) rename mac/{cholesky_utils.py => utils/cholesky.py} (100%) create mode 100644 mac/utils/conversions.py rename mac/{ => utils}/fiedler.py (100%) rename mac/{utils.py => utils/graphs.py} (58%) create mode 100644 mac/utils/rounding.py 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/optimization/constraints.py b/mac/optimization/constraints.py new file mode 100644 index 0000000..4975738 --- /dev/null +++ b/mac/optimization/constraints.py @@ -0,0 +1,29 @@ +from mac.utils 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(w) + 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..3c6fcc0 --- /dev/null +++ b/mac/optimization/frankwolfe.py @@ -0,0 +1,78 @@ +""" +Utilities for maximizing concave functions over simple feasible sets with Frank-Wolfe +""" + +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 + 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 100% rename from mac/greedy_eig.py rename to mac/solvers/greedy_eig.py diff --git a/mac/greedy_esp.py b/mac/solvers/greedy_esp.py similarity index 100% rename from mac/greedy_esp.py rename to mac/solvers/greedy_esp.py diff --git a/mac/mac.py b/mac/solvers/mac.py similarity index 91% rename from mac/mac.py rename to mac/solvers/mac.py index 519a339..add09d2 100644 --- a/mac/mac.py +++ b/mac/solvers/mac.py @@ -21,7 +21,7 @@ class Cache: class MAC: 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 ---------- @@ -36,8 +36,6 @@ 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 @@ -47,25 +45,16 @@ def __init__(self, fixed_edges, candidate_edges, num_nodes, assert(len(fixed_edges) == len(candidate_edges) == 0) 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) - # Set up caching. - self.cache = None - if use_cache: - self.cache = Cache(Q=None) - # Configuration for Fiedler vector computation self.fiedler_method = fiedler_method self.fiedler_tol = fiedler_tol @@ -102,13 +91,13 @@ def evaluate_objective(self, x): """ return fiedler.find_fiedler_pair(self.laplacian(x), method=self.fiedler_method, tol=self.fiedler_tol)[0] - def problem(self, x): + 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. returns x, grad F(x). """ - Q = None if self.cache is None else self.cache.Q + Q = None if cache is None else cache.Q fiedler_value, fiedler_vec, Qnew = fiedler.find_fiedler_pair(x, X=Q) gradf = np.zeros(len(self.weights)) @@ -120,13 +109,13 @@ def problem(self, x): kdelta = weight_k * (v_i - v_j) gradf[k] = kdelta * (v_i - v_j) - if self.cache is not None: - self.cache.Q = Q + if cache is not None: + cache.Q = Q return f, gradf, cache 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): + 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 @@ -183,13 +172,16 @@ def solve(self, k, x_init=None, rounding="nearest", fallback=False, # handle case where x is none # Set up problem to use cache (or not) - # problem = lambda x: self.problem() + cache = None + if use_cache: + cache = 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=x_init, - problem=self.problem_with_recycling, + problem=problem, solve_lp=solve_lp, maxiter=max_iters, relative_duality_gap_tol=relative_duality_gap_tol, diff --git a/mac/cholesky_utils.py b/mac/utils/cholesky.py similarity index 100% rename from mac/cholesky_utils.py rename to mac/utils/cholesky.py diff --git a/mac/utils/conversions.py b/mac/utils/conversions.py new file mode 100644 index 0000000..1618d4a --- /dev/null +++ b/mac/utils/conversions.py @@ -0,0 +1,35 @@ +""" +Utilities for converting between graph types. +""" +import networkx as nx +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(): + 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 diff --git a/mac/fiedler.py b/mac/utils/fiedler.py similarity index 100% rename from mac/fiedler.py rename to mac/utils/fiedler.py diff --git a/mac/utils.py b/mac/utils/graphs.py similarity index 58% rename from mac/utils.py rename to mac/utils/graphs.py index db8b767..a1985a4 100644 --- a/mac/utils.py +++ b/mac/utils/graphs.py @@ -1,140 +1,15 @@ +""" +Types and utilities for working with graphs. +""" + import numpy as np +from typing import List, Union, Tuple 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 @@ -332,4 +207,3 @@ def get_edge_selection_as_binary_mask( 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..6053eab --- /dev/null +++ b/mac/utils/rounding.py @@ -0,0 +1,93 @@ +""" +Utilities for rounding solutions to convex optimization problems onto simple constraint sets. +""" + +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 From 764de706e5a91cdf67ce8566a09c208b40f65cdc Mon Sep 17 00:00:00 2001 From: Kevin Date: Fri, 2 Aug 2024 16:14:49 -0400 Subject: [PATCH 07/36] Test structure now mirrors project structure. --- tests/README.md | 0 tests/{ => solvers}/test_greedy_eig.py | 0 tests/{ => solvers}/test_greedy_esp.py | 0 .../test_mac.py} | 0 .../test_cholesky.py} | 0 tests/utils/test_fiedler.py | 24 +++++++++++++++++++ 6 files changed, 24 insertions(+) create mode 100644 tests/README.md rename tests/{ => solvers}/test_greedy_eig.py (100%) rename tests/{ => solvers}/test_greedy_esp.py (100%) rename tests/{test_mac_regression.py => solvers/test_mac.py} (100%) rename tests/{test_cholesky_utils.py => utils/test_cholesky.py} (100%) create mode 100644 tests/utils/test_fiedler.py diff --git a/tests/README.md b/tests/README.md new file mode 100644 index 0000000..e69de29 diff --git a/tests/test_greedy_eig.py b/tests/solvers/test_greedy_eig.py similarity index 100% rename from tests/test_greedy_eig.py rename to tests/solvers/test_greedy_eig.py diff --git a/tests/test_greedy_esp.py b/tests/solvers/test_greedy_esp.py similarity index 100% rename from tests/test_greedy_esp.py rename to tests/solvers/test_greedy_esp.py diff --git a/tests/test_mac_regression.py b/tests/solvers/test_mac.py similarity index 100% rename from tests/test_mac_regression.py rename to tests/solvers/test_mac.py diff --git a/tests/test_cholesky_utils.py b/tests/utils/test_cholesky.py similarity index 100% rename from tests/test_cholesky_utils.py rename to tests/utils/test_cholesky.py diff --git a/tests/utils/test_fiedler.py b/tests/utils/test_fiedler.py new file mode 100644 index 0000000..8ee2d8d --- /dev/null +++ b/tests/utils/test_fiedler.py @@ -0,0 +1,24 @@ +""" +Copyright 2023 MIT Marine Robotics Group + +Tests for Fiedler utilities + +Author: Kevin Doherty +""" + +import unittest +import numpy as np + +import mac.fiedler as fiedler +from .utils import get_split_petersen_graph + +class TestFiedler(unittest.TestCase): + def setUp(self): + return + + def test_cache(self): + return + + +if __name__ == '__main__': + unittest.main() From fb451839f8da9ea2abc801d1a1c7c8b65bc55e3e Mon Sep 17 00:00:00 2001 From: Kevin Date: Fri, 2 Aug 2024 17:07:20 -0400 Subject: [PATCH 08/36] Update tests. --- mac/__init__.py | 1 - mac/utils/conversions.py | 14 +++++++-- mac/utils/graphs.py | 30 ------------------- mac/utils/rounding.py | 1 + tests/utils/test_conversions.py | 53 +++++++++++++++++++++++++++++++++ tests/utils/test_graphs.py | 50 +++++++++++++++++++++++++++++++ 6 files changed, 116 insertions(+), 33 deletions(-) create mode 100644 tests/utils/test_conversions.py create mode 100644 tests/utils/test_graphs.py 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/utils/conversions.py b/mac/utils/conversions.py index 1618d4a..fce3742 100644 --- a/mac/utils/conversions.py +++ b/mac/utils/conversions.py @@ -2,6 +2,8 @@ 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]: @@ -15,7 +17,12 @@ def nx_to_mac(G: nx.Graph) -> List[Edge]: """ edges = [] for nxedge in G.edges(): - edge = Edge(nxedge[0], nxedge[1], 1.0) + i = nxedge[0] + j = nxedge[1] + if i < j: + edge = Edge(i, j, 1.0) + else: + edge = Edge(j, i, 1.0) edges.append(edge) return edges @@ -31,5 +38,8 @@ def mac_to_nx(edges: List[Edge]) -> nx.Graph: """ G = nx.Graph() for edge in edges: - G.add_edge(edge.i, edge.j, weight=edge.weight) + 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/graphs.py b/mac/utils/graphs.py index a1985a4..0fb9cb0 100644 --- a/mac/utils/graphs.py +++ b/mac/utils/graphs.py @@ -98,36 +98,6 @@ def weight_graph_lap_from_edges( return csr_matrix(coo_matrix((data, (rows, cols)), shape=[num_nodes, num_nodes])) -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 diff --git a/mac/utils/rounding.py b/mac/utils/rounding.py index 6053eab..5cd62b2 100644 --- a/mac/utils/rounding.py +++ b/mac/utils/rounding.py @@ -1,6 +1,7 @@ """ 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): """ diff --git a/tests/utils/test_conversions.py b/tests/utils/test_conversions.py new file mode 100644 index 0000000..e87e657 --- /dev/null +++ b/tests/utils/test_conversions.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 + +# 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)] + + 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)) + +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..ed34f32 --- /dev/null +++ b/tests/utils/test_graphs.py @@ -0,0 +1,50 @@ +""" +Copyright 2023 MIT Marine Robotics Group + +Tests for Fiedler utilities + +Author: Kevin Doherty +""" + +import unittest +import numpy as np + +import mac.fiedler as fiedler +from .utils import get_split_petersen_graph + +class TestGraphUtils(unittest.TestCase): + def setUp(self): + return + + def test_cache(self): + # 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])) + return + + +if __name__ == '__main__': + unittest.main() From 541e05a91f43ff21cf3488ad534cb4ea163e2f4d Mon Sep 17 00:00:00 2001 From: Kevin Date: Fri, 2 Aug 2024 17:07:37 -0400 Subject: [PATCH 09/36] Add future benchmark dir. --- benchmarks/README.md | 0 benchmarks/cache_performance.py | 2 ++ 2 files changed, 2 insertions(+) create mode 100644 benchmarks/README.md create mode 100644 benchmarks/cache_performance.py diff --git a/benchmarks/README.md b/benchmarks/README.md new file mode 100644 index 0000000..e69de29 diff --git a/benchmarks/cache_performance.py b/benchmarks/cache_performance.py new file mode 100644 index 0000000..0340bd1 --- /dev/null +++ b/benchmarks/cache_performance.py @@ -0,0 +1,2 @@ +# TODO implement me +# We should benchmark performance of MAC with vs. without caching on some datasets From 128351a55b9e4eba0554ce4de98918e4cc1548ff Mon Sep 17 00:00:00 2001 From: Kevin Date: Fri, 2 Aug 2024 17:07:57 -0400 Subject: [PATCH 10/36] Move split edges utility into pose graph utils. --- examples/pose_graph_utils.py | 30 ++++++++++++++++++++++++++++++ 1 file changed, 30 insertions(+) diff --git a/examples/pose_graph_utils.py b/examples/pose_graph_utils.py index 7aa8cb2..0c83e1a 100644 --- a/examples/pose_graph_utils.py +++ b/examples/pose_graph_utils.py @@ -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: From 85e89ddf9036d567892152539c51ecce9f276d70 Mon Sep 17 00:00:00 2001 From: Kevin Date: Fri, 2 Aug 2024 17:42:30 -0400 Subject: [PATCH 11/36] Update tests and conversion code. --- mac/utils/conversions.py | 8 +++-- tests/utils/test_conversions.py | 16 ++++++++- tests/utils/test_graphs.py | 57 ++++++++++++++------------------- 3 files changed, 45 insertions(+), 36 deletions(-) diff --git a/mac/utils/conversions.py b/mac/utils/conversions.py index fce3742..e134e68 100644 --- a/mac/utils/conversions.py +++ b/mac/utils/conversions.py @@ -19,10 +19,14 @@ def nx_to_mac(G: nx.Graph) -> List[Edge]: 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, 1.0) + edge = Edge(i, j, weight) else: - edge = Edge(j, i, 1.0) + edge = Edge(j, i, weight) edges.append(edge) return edges diff --git a/tests/utils/test_conversions.py b/tests/utils/test_conversions.py index e87e657..0af620a 100644 --- a/tests/utils/test_conversions.py +++ b/tests/utils/test_conversions.py @@ -1,7 +1,7 @@ """ Copyright 2023 MIT Marine Robotics Group -Tests for Fiedler utilities +Tests for graph type conversion utilities Author: Kevin Doherty """ @@ -42,6 +42,12 @@ def setUp(self): 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)) @@ -49,5 +55,13 @@ def test_mac_to_mac_via_nx(self): 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_graphs.py b/tests/utils/test_graphs.py index ed34f32..5c3a1e6 100644 --- a/tests/utils/test_graphs.py +++ b/tests/utils/test_graphs.py @@ -1,49 +1,40 @@ """ Copyright 2023 MIT Marine Robotics Group -Tests for Fiedler utilities +Tests for graph utilities. Author: Kevin Doherty """ import unittest +import networkx as nx import numpy as np -import mac.fiedler as fiedler -from .utils import get_split_petersen_graph +from mac.utils.conversions import nx_to_mac + +# Code under test +from mac.utils.graphs import * class TestGraphUtils(unittest.TestCase): def setUp(self): - return - - def test_cache(self): - # 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])) - return + 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_unweighted(self): + edge_list = nx_to_mac(self.weighted_graph) + L_ours = weight_graph_lap_from_edge_list(edge_list, self.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__': From 1a81b821ed2731898472981738c343ddde0a33f1 Mon Sep 17 00:00:00 2001 From: Kevin Date: Fri, 2 Aug 2024 17:43:26 -0400 Subject: [PATCH 12/36] Fix typo. --- tests/utils/test_graphs.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/utils/test_graphs.py b/tests/utils/test_graphs.py index 5c3a1e6..d4388f9 100644 --- a/tests/utils/test_graphs.py +++ b/tests/utils/test_graphs.py @@ -30,7 +30,7 @@ def test_weight_graph_lap_from_edge_list_unweighted(self): 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_unweighted(self): + 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.graph.number_of_nodes()) L_nx = nx.laplacian_matrix(self.weighted_graph) From faf1ffc5f4733e7661f7d5c07d942381152e0e8d Mon Sep 17 00:00:00 2001 From: Kevin Date: Fri, 2 Aug 2024 17:47:56 -0400 Subject: [PATCH 13/36] Test main Laplacian builder function. --- tests/utils/test_graphs.py | 14 +++++++++++++- 1 file changed, 13 insertions(+), 1 deletion(-) diff --git a/tests/utils/test_graphs.py b/tests/utils/test_graphs.py index d4388f9..9da1dc4 100644 --- a/tests/utils/test_graphs.py +++ b/tests/utils/test_graphs.py @@ -32,10 +32,22 @@ def test_weight_graph_lap_from_edge_list_unweighted(self): 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.graph.number_of_nodes()) + 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() From 18831c5eb5ee38ca2f35f62731b8aa314204e659 Mon Sep 17 00:00:00 2001 From: Kevin Date: Fri, 2 Aug 2024 17:51:55 -0400 Subject: [PATCH 14/36] WIP tidying Cholesky tests. --- tests/utils/test_cholesky.py | 17 +++++------------ tests/utils/test_fiedler.py | 5 ++--- 2 files changed, 7 insertions(+), 15 deletions(-) diff --git a/tests/utils/test_cholesky.py b/tests/utils/test_cholesky.py index d78a397..3eb61de 100644 --- a/tests/utils/test_cholesky.py +++ b/tests/utils/test_cholesky.py @@ -10,26 +10,19 @@ 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 mac.utils.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 +# Code under test +from mac.utils.cholesky import * + ORDERING_METHODS = ["natural", "best", "amd", "metis", "nesdis", "colamd", "default"] diff --git a/tests/utils/test_fiedler.py b/tests/utils/test_fiedler.py index 8ee2d8d..b72b2f0 100644 --- a/tests/utils/test_fiedler.py +++ b/tests/utils/test_fiedler.py @@ -9,10 +9,9 @@ import unittest import numpy as np -import mac.fiedler as fiedler -from .utils import get_split_petersen_graph +import mac.utils.fiedler import * -class TestFiedler(unittest.TestCase): +class TestFiedlerUtils(unittest.TestCase): def setUp(self): return From 61ebc95b88a22a6d86d23d562e37ae4f874442cc Mon Sep 17 00:00:00 2001 From: Kevin Date: Fri, 2 Aug 2024 18:15:55 -0400 Subject: [PATCH 15/36] Move Cache to inner class. --- mac/solvers/mac.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/mac/solvers/mac.py b/mac/solvers/mac.py index add09d2..7699813 100644 --- a/mac/solvers/mac.py +++ b/mac/solvers/mac.py @@ -14,12 +14,12 @@ from timeit import default_timer as timer -@dataclass -class Cache: - """Problem data cache""" - Q: Optional[np.ndarray] = None - 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', fiedler_tol=1e-8, min_selection_weight_tol=1e-10): From 49d3a99c282f92ed9a615facac6e6667c55a8e2e Mon Sep 17 00:00:00 2001 From: Kevin Date: Fri, 2 Aug 2024 18:16:09 -0400 Subject: [PATCH 16/36] Note on future tests to add. --- tests/utils/test_fiedler.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/tests/utils/test_fiedler.py b/tests/utils/test_fiedler.py index b72b2f0..3e2f38b 100644 --- a/tests/utils/test_fiedler.py +++ b/tests/utils/test_fiedler.py @@ -15,6 +15,10 @@ class TestFiedlerUtils(unittest.TestCase): def setUp(self): return + def test_disconnected(self): + # Test multiple connected components + return + def test_cache(self): return From 7e21f939d21594e86f6b1418619a5f689273001e Mon Sep 17 00:00:00 2001 From: Kevin Date: Fri, 2 Aug 2024 18:39:49 -0400 Subject: [PATCH 17/36] Add test for complete graph. --- tests/utils/test_fiedler.py | 21 +++++++++++++++++++-- 1 file changed, 19 insertions(+), 2 deletions(-) diff --git a/tests/utils/test_fiedler.py b/tests/utils/test_fiedler.py index 3e2f38b..e4bb018 100644 --- a/tests/utils/test_fiedler.py +++ b/tests/utils/test_fiedler.py @@ -8,12 +8,29 @@ import unittest import numpy as np +import networkx as nx -import mac.utils.fiedler import * +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 TestFiedlerUtils(unittest.TestCase): def setUp(self): - return + # Complete graph on 5 nodes + self.complete_graph = nx.complete_graph(5) + + def test_complete_graph(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)) def test_disconnected(self): # Test multiple connected components From 3b154e3e4f283bba99cfe425f6238bc1d1bb601e Mon Sep 17 00:00:00 2001 From: Kevin Date: Fri, 2 Aug 2024 18:50:01 -0400 Subject: [PATCH 18/36] Add failing test for disconnected component case. --- tests/utils/test_fiedler.py | 24 ++++++++++++++++++------ 1 file changed, 18 insertions(+), 6 deletions(-) diff --git a/tests/utils/test_fiedler.py b/tests/utils/test_fiedler.py index e4bb018..06fdd4b 100644 --- a/tests/utils/test_fiedler.py +++ b/tests/utils/test_fiedler.py @@ -12,18 +12,19 @@ from mac.utils.conversions import nx_to_mac from mac.utils.graphs import weight_graph_lap_from_edge_list +import matplotlib.pyplot as plt # Code under test. from mac.utils.fiedler import * from networkx.linalg.algebraicconnectivity import algebraic_connectivity -class TestFiedlerUtils(unittest.TestCase): +class TestConnectedGraphs(unittest.TestCase): def setUp(self): # Complete graph on 5 nodes self.complete_graph = nx.complete_graph(5) - def test_complete_graph(self): + 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) @@ -32,13 +33,24 @@ def test_complete_graph(self): fiedler_value, fiedler_vec, _ = find_fiedler_pair(L) self.assertTrue(np.isclose(fiedler_value, N)) - def test_disconnected(self): - # Test multiple connected components - return - def test_cache(self): return +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)) + + 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() From e69b378811bcdbee33e9cd984376d3bf97f54e5e Mon Sep 17 00:00:00 2001 From: Kevin Date: Fri, 2 Aug 2024 19:43:49 -0400 Subject: [PATCH 19/36] Fixes for optimization code. --- mac/optimization/constraints.py | 4 ++-- mac/optimization/frankwolfe.py | 5 ++++- 2 files changed, 6 insertions(+), 3 deletions(-) diff --git a/mac/optimization/constraints.py b/mac/optimization/constraints.py index 4975738..63f7f5d 100644 --- a/mac/optimization/constraints.py +++ b/mac/optimization/constraints.py @@ -1,4 +1,4 @@ -from mac.utils import round_nearest +from mac.utils.rounding import round_nearest import numpy as np def solve_subset_box_lp(g, k): @@ -24,6 +24,6 @@ def solve_box_lp(g): the positions of all nonnegative elements of g. """ - solution = np.zeros_like(w) + 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 index 3c6fcc0..5a50831 100644 --- a/mac/optimization/frankwolfe.py +++ b/mac/optimization/frankwolfe.py @@ -4,6 +4,9 @@ import numpy as np +def naive_stepsize(k): + return 2.0 / (k + 2.0) + def frank_wolfe(initial, problem, solve_lp, @@ -65,7 +68,7 @@ def frank_wolfe(initial, return x, u # If the *relative* duality gap is sufficiently small, we are done. - if (u - f) / f < relative_duality_gap_tol: + if (u - f) / abs(f) < relative_duality_gap_tol: if verbose: print("Duality gap tolerance reached, found optimal solution") return x, u From b2ab8ac92abac5c439caa4e9e9d1ccc7ccd7ab2a Mon Sep 17 00:00:00 2001 From: Kevin Date: Fri, 2 Aug 2024 19:44:04 -0400 Subject: [PATCH 20/36] Remove unused import. --- tests/utils/test_fiedler.py | 1 - 1 file changed, 1 deletion(-) diff --git a/tests/utils/test_fiedler.py b/tests/utils/test_fiedler.py index 06fdd4b..d5dfed4 100644 --- a/tests/utils/test_fiedler.py +++ b/tests/utils/test_fiedler.py @@ -12,7 +12,6 @@ from mac.utils.conversions import nx_to_mac from mac.utils.graphs import weight_graph_lap_from_edge_list -import matplotlib.pyplot as plt # Code under test. from mac.utils.fiedler import * From 7a2691ab3fdb04e226f90e2d2500cb593f830671 Mon Sep 17 00:00:00 2001 From: Kevin Date: Fri, 2 Aug 2024 19:44:11 -0400 Subject: [PATCH 21/36] Start tests for optimization code. --- tests/optimization/test_frankwolfe.py | 55 +++++++++++++++++++++++++++ 1 file changed, 55 insertions(+) create mode 100644 tests/optimization/test_frankwolfe.py diff --git a/tests/optimization/test_frankwolfe.py b/tests/optimization/test_frankwolfe.py new file mode 100644 index 0000000..f36cd1d --- /dev/null +++ b/tests/optimization/test_frankwolfe.py @@ -0,0 +1,55 @@ +""" +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)) + +if __name__ == '__main__': + unittest.main() From 7b135db2936eeba7a792babe9c02c0bc32b564c9 Mon Sep 17 00:00:00 2001 From: Kevin Date: Fri, 2 Aug 2024 21:08:23 -0400 Subject: [PATCH 22/36] Fix relative duality gap check, update constraints docs, and add test. --- mac/optimization/constraints.py | 10 +++++++++- mac/optimization/frankwolfe.py | 3 +-- tests/optimization/test_frankwolfe.py | 21 +++++++++++++++++++++ 3 files changed, 31 insertions(+), 3 deletions(-) diff --git a/mac/optimization/constraints.py b/mac/optimization/constraints.py index 63f7f5d..41232cf 100644 --- a/mac/optimization/constraints.py +++ b/mac/optimization/constraints.py @@ -1,3 +1,11 @@ +""" +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 @@ -25,5 +33,5 @@ def solve_box_lp(g): """ solution = np.zeros_like(g) - solution[g >= 0.0] = 1.0 + solution[g > 0.0] = 1.0 return solution diff --git a/mac/optimization/frankwolfe.py b/mac/optimization/frankwolfe.py index 5a50831..db34673 100644 --- a/mac/optimization/frankwolfe.py +++ b/mac/optimization/frankwolfe.py @@ -68,13 +68,12 @@ def frank_wolfe(initial, return x, u # If the *relative* duality gap is sufficiently small, we are done. - if (u - f) / abs(f) < relative_duality_gap_tol: + 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) - pass if verbose: print("Reached maximum number of iterations, returning best solution") return x, u diff --git a/tests/optimization/test_frankwolfe.py b/tests/optimization/test_frankwolfe.py index f36cd1d..78a60ce 100644 --- a/tests/optimization/test_frankwolfe.py +++ b/tests/optimization/test_frankwolfe.py @@ -51,5 +51,26 @@ def test_solve_subset_box_constraint(self): # 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() From c75dcab16e690fca5e338272cd67248ebb0d3f77 Mon Sep 17 00:00:00 2001 From: Kevin Date: Sat, 3 Aug 2024 14:40:40 -0400 Subject: [PATCH 23/36] MAC now working. --- mac/solvers/mac.py | 47 ++++++++++++++++----------------------- mac/utils/rounding.py | 1 + tests/solvers/test_mac.py | 47 ++++++++++++++++++++++----------------- 3 files changed, 47 insertions(+), 48 deletions(-) diff --git a/mac/solvers/mac.py b/mac/solvers/mac.py index 7699813..edc1636 100644 --- a/mac/solvers/mac.py +++ b/mac/solvers/mac.py @@ -1,19 +1,19 @@ -from mac.utils import * -import mac.fiedler as fiedler -import mac.frankwolfe as fw -from timeit import default_timer as timer -from dataclasses import dataclass -from typing import Optional - import numpy as np import networkx as nx import networkx.linalg as la -import matplotlib.pyplot as plt +from dataclasses import dataclass +from typing import Optional from scipy.sparse import csc_matrix, csr_matrix - 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: @@ -89,7 +89,7 @@ def evaluate_objective(self, x): returns F(x) = lambda_2(L(x)). """ - return fiedler.find_fiedler_pair(self.laplacian(x), method=self.fiedler_method, tol=self.fiedler_tol)[0] + return fiedler.find_fiedler_pair(L=self.laplacian(x), method=self.fiedler_method, tol=self.fiedler_tol)[0] def problem(self, x, cache=None): """ @@ -98,7 +98,7 @@ def problem(self, x, cache=None): returns x, grad F(x). """ Q = None if cache is None else cache.Q - fiedler_value, fiedler_vec, Qnew = fiedler.find_fiedler_pair(x, X=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)): @@ -111,7 +111,7 @@ def problem(self, x, cache=None): if cache is not None: cache.Q = Q - return f, gradf, cache + return f, gradf def solve(self, k, x_init=None, rounding="nearest", fallback=False, max_iters=5, relative_duality_gap_tol=1e-4, @@ -164,10 +164,10 @@ def solve(self, k, x_init=None, rounding="nearest", fallback=False, return result, result, self.evaluate_objective(np.ones(len(self.weights))) - assert(len(w_init) == len(self.weights)) + 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) # handle case where x is none @@ -179,20 +179,11 @@ def solve(self, k, x_init=None, rounding="nearest", fallback=False, # 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=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) - else: - w, u = fw.frank_wolfe(initial=x_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": diff --git a/mac/utils/rounding.py b/mac/utils/rounding.py index 5cd62b2..5b5545a 100644 --- a/mac/utils/rounding.py +++ b/mac/utils/rounding.py @@ -1,6 +1,7 @@ """ 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): diff --git a/tests/solvers/test_mac.py b/tests/solvers/test_mac.py index 957c307..c6b0ebc 100644 --- a/tests/solvers/test_mac.py +++ b/tests/solvers/test_mac.py @@ -10,37 +10,46 @@ 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 +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() -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)) + num_candidates = int(pct_candidates * len(self.candidate_edges)) # Construct an initial guess - w_init = np.zeros(len(candidate_edges)) - w_init[:num_candidates] = 1.0 + x_init = np.zeros(len(self.candidate_edges)) + x_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) + 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(w_init) + init_l2 = mac.evaluate_objective(x_init) mac_unrounded_l2 = mac.evaluate_objective(unrounded) mac_l2 = mac.evaluate_objective(result) @@ -48,7 +57,7 @@ def test_petersen(self): 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 + 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 @@ -102,7 +111,5 @@ def test_cache(self): 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() From 0b1bcf46935bb81a96c64ae5308c22d4712651ff Mon Sep 17 00:00:00 2001 From: Kevin Date: Sun, 4 Aug 2024 16:52:50 -0400 Subject: [PATCH 24/36] Some doc updates and minor cleanup. --- mac/solvers/mac.py | 53 +++++++++++++++++++++++++++++----------------- 1 file changed, 34 insertions(+), 19 deletions(-) diff --git a/mac/solvers/mac.py b/mac/solvers/mac.py index edc1636..add826a 100644 --- a/mac/solvers/mac.py +++ b/mac/solvers/mac.py @@ -4,7 +4,6 @@ from dataclasses import dataclass from typing import Optional -from scipy.sparse import csc_matrix, csr_matrix from timeit import default_timer as timer from mac.utils.graphs import * @@ -41,13 +40,23 @@ def __init__(self, fixed_edges, candidate_edges, num_nodes, 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.weights = [] self.edge_list = [] - for edge in candidate_edges: self.weights.append(edge.weight) self.edge_list.append((edge.i, edge.j)) @@ -63,15 +72,15 @@ def __init__(self, fixed_edges, candidate_edges, num_nodes, self.min_selection_weight_tol = min_selection_weight_tol def laplacian(self, x): - """ - Construct the combined Laplacian (fixed edges plus candidate edges weighted by 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 - 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) + 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. + + returns the matrix L(x) """ idx = np.where(x > self.min_selection_weight_tol) prod = x[idx]*self.weights[idx] @@ -89,11 +98,16 @@ def evaluate_objective(self, x): returns F(x) = lambda_2(L(x)). """ - return fiedler.find_fiedler_pair(L=self.laplacian(x), method=self.fiedler_method, tol=self.fiedler_tol)[0] + return fiedler.find_fiedler_pair(L=self.laplacian(x), + method=self.fiedler_method, tol=self.fiedler_tol)[0] 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. + """Compute the algebraic connectivity of L(x) and a (super)gradient of the + algebraic connectivity with respect to x. + + 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 x, grad F(x). """ @@ -114,17 +128,18 @@ def problem(self, x, cache=None): 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): + 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 ---------- - x_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). From 587feef496eb09e38f10e95434f62da78832f1a2 Mon Sep 17 00:00:00 2001 From: Kevin Date: Sat, 10 Aug 2024 14:22:43 -0400 Subject: [PATCH 25/36] Add __init__ to support test discovery. --- tests/optimization/__init__.py | 0 tests/solvers/__init__.py | 0 tests/utils.py | 46 ---------------------------------- tests/utils/__init__.py | 0 4 files changed, 46 deletions(-) create mode 100644 tests/optimization/__init__.py create mode 100644 tests/solvers/__init__.py delete mode 100644 tests/utils.py create mode 100644 tests/utils/__init__.py diff --git a/tests/optimization/__init__.py b/tests/optimization/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/tests/solvers/__init__.py b/tests/solvers/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/tests/utils.py b/tests/utils.py deleted file mode 100644 index 8bb520a..0000000 --- a/tests/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/tests/utils/__init__.py b/tests/utils/__init__.py new file mode 100644 index 0000000..e69de29 From 628be78291d677dc3dedf83a18cc12e220596f40 Mon Sep 17 00:00:00 2001 From: Kevin Date: Wed, 14 Aug 2024 19:36:48 -0400 Subject: [PATCH 26/36] Tidying unused names. --- mac/solvers/greedy_esp.py | 2 +- tests/solvers/test_mac.py | 40 --------------------------------------- 2 files changed, 1 insertion(+), 41 deletions(-) diff --git a/mac/solvers/greedy_esp.py b/mac/solvers/greedy_esp.py index eb6ed16..36384b5 100644 --- a/mac/solvers/greedy_esp.py +++ b/mac/solvers/greedy_esp.py @@ -272,7 +272,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/tests/solvers/test_mac.py b/tests/solvers/test_mac.py index c6b0ebc..51e7f5f 100644 --- a/tests/solvers/test_mac.py +++ b/tests/solvers/test_mac.py @@ -71,45 +71,5 @@ def test_petersen(self): # 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)) - if __name__ == '__main__': unittest.main() From 46c27ab8b9d078b06c5ded0376976ee6a1eb22b0 Mon Sep 17 00:00:00 2001 From: Kevin Date: Tue, 8 Oct 2024 17:29:56 -0400 Subject: [PATCH 27/36] More cleanup. --- examples/greedy_eig_test.py | 61 --------- tests/solvers/test_greedy_eig.py | 116 ---------------- tests/solvers/test_greedy_esp.py | 146 -------------------- tests/utils/test_cholesky.py | 224 ------------------------------- 4 files changed, 547 deletions(-) delete mode 100644 examples/greedy_eig_test.py delete mode 100644 tests/solvers/test_greedy_eig.py delete mode 100644 tests/solvers/test_greedy_esp.py delete mode 100644 tests/utils/test_cholesky.py 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/tests/solvers/test_greedy_eig.py b/tests/solvers/test_greedy_eig.py deleted file mode 100644 index 7a67a57..0000000 --- a/tests/solvers/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/tests/solvers/test_greedy_esp.py b/tests/solvers/test_greedy_esp.py deleted file mode 100644 index a31f8c5..0000000 --- a/tests/solvers/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/tests/utils/test_cholesky.py b/tests/utils/test_cholesky.py deleted file mode 100644 index 3eb61de..0000000 --- a/tests/utils/test_cholesky.py +++ /dev/null @@ -1,224 +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 ( - set_incidence_vector_for_edge_inplace, - weight_reduced_graph_lap_from_edge_list, - weight_graph_lap_from_edge_list - ) - -from mac.utils.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 - -# Code under test -from mac.utils.cholesky import * - - -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() From 6613c4bfe8c50aa1241c090c4185c24852671a2e Mon Sep 17 00:00:00 2001 From: Kevin Date: Wed, 23 Oct 2024 20:30:17 -0400 Subject: [PATCH 28/36] Fix comment. --- mac/utils/cholesky.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/mac/utils/cholesky.py b/mac/utils/cholesky.py index 223646a..9bcfe2b 100644 --- a/mac/utils/cholesky.py +++ b/mac/utils/cholesky.py @@ -107,7 +107,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 From 96f175eb0b376749290b454d37a06fae74dc9bc5 Mon Sep 17 00:00:00 2001 From: Kevin Date: Wed, 23 Oct 2024 20:30:40 -0400 Subject: [PATCH 29/36] Skip test for unimplemented feature. --- tests/utils/test_fiedler.py | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/tests/utils/test_fiedler.py b/tests/utils/test_fiedler.py index d5dfed4..9b2296c 100644 --- a/tests/utils/test_fiedler.py +++ b/tests/utils/test_fiedler.py @@ -32,9 +32,6 @@ def test_fiedler(self): fiedler_value, fiedler_vec, _ = find_fiedler_pair(L) self.assertTrue(np.isclose(fiedler_value, N)) - def test_cache(self): - return - class TestDisconnectedGraphs(unittest.TestCase): def setUp(self): # Construct a disconnected graph with two connected components @@ -43,6 +40,7 @@ def setUp(self): 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) From bf5e3b33f9d99e4af4ec7f549fa315fe9971a752 Mon Sep 17 00:00:00 2001 From: Kevin Date: Sat, 2 Nov 2024 13:26:14 -0400 Subject: [PATCH 30/36] Rm old benchmarks dir. --- benchmarks/README.md | 0 benchmarks/cache_performance.py | 2 -- 2 files changed, 2 deletions(-) delete mode 100644 benchmarks/README.md delete mode 100644 benchmarks/cache_performance.py diff --git a/benchmarks/README.md b/benchmarks/README.md deleted file mode 100644 index e69de29..0000000 diff --git a/benchmarks/cache_performance.py b/benchmarks/cache_performance.py deleted file mode 100644 index 0340bd1..0000000 --- a/benchmarks/cache_performance.py +++ /dev/null @@ -1,2 +0,0 @@ -# TODO implement me -# We should benchmark performance of MAC with vs. without caching on some datasets From 70c7bf4b83aac49c32f8971e4709a380b5ea6720 Mon Sep 17 00:00:00 2001 From: Kevin Date: Sat, 2 Nov 2024 13:26:27 -0400 Subject: [PATCH 31/36] Fix bug in cache usage. --- mac/solvers/mac.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/mac/solvers/mac.py b/mac/solvers/mac.py index add826a..4f786f9 100644 --- a/mac/solvers/mac.py +++ b/mac/solvers/mac.py @@ -179,17 +179,16 @@ def solve(self, k, x_init=None, rounding="nearest", fallback=False, return result, result, self.evaluate_objective(np.ones(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: constraints.solve_subset_box_lp(g, k) - # handle case where x is none - # Set up problem to use cache (or not) cache = None if use_cache: - cache = Cache() + cache = MAC.Cache() problem = lambda x: self.problem(x, cache=cache) # Run Frank-Wolfe to solve the relaxation of subset constrained From 49c1959e84623c4d1b09cfdd8299a1d58aae7c84 Mon Sep 17 00:00:00 2001 From: Kevin Date: Sat, 2 Nov 2024 13:26:41 -0400 Subject: [PATCH 32/36] Add new benchmarks. --- tests/benchmarks/__init__.py | 0 tests/benchmarks/test_cache_performance.py | 49 ++++++++++++++++++++++ 2 files changed, 49 insertions(+) create mode 100644 tests/benchmarks/__init__.py create mode 100644 tests/benchmarks/test_cache_performance.py 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) From 576c66a8ca3fbc4d54876abe67ba6ccd817e334b Mon Sep 17 00:00:00 2001 From: Kevin Date: Sat, 2 Nov 2024 13:28:32 -0400 Subject: [PATCH 33/36] Add pytest-benchmark to deps. --- .github/workflows/test.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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: | From 55be6adf5d7a1df52e259126aace1d010e3fae32 Mon Sep 17 00:00:00 2001 From: Kevin Date: Sat, 2 Nov 2024 13:49:17 -0400 Subject: [PATCH 34/36] Clean up for Petersen graph example. --- examples/petersen_graph_sparsification.py | 21 +++++++++++++-------- mac/solvers/greedy_eig.py | 4 ++-- mac/solvers/greedy_esp.py | 5 ++--- mac/utils/cholesky.py | 3 ++- 4 files changed, 19 insertions(+), 14 deletions(-) diff --git a/examples/petersen_graph_sparsification.py b/examples/petersen_graph_sparsification.py index 8331f6e..cc56d94 100644 --- a/examples/petersen_graph_sparsification.py +++ b/examples/petersen_graph_sparsification.py @@ -1,11 +1,11 @@ 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.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 @@ -30,16 +30,21 @@ fixed_edges = nx_to_mac(spanning_tree) candidate_edges = nx_to_mac(loop_graph) -pct_candidates = 0.4 +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_edges)) -w_init[:num_candidates] = 1.0 +# 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, rounding="nearest") +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) diff --git a/mac/solvers/greedy_eig.py b/mac/solvers/greedy_eig.py index aa3cdbe..4108871 100644 --- a/mac/solvers/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/solvers/greedy_esp.py b/mac/solvers/greedy_esp.py index 36384b5..7565e33 100644 --- a/mac/solvers/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( diff --git a/mac/utils/cholesky.py b/mac/utils/cholesky.py index 9bcfe2b..34cecb6 100644 --- a/mac/utils/cholesky.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 From 312ebaf8cfafb04bb5077cfad23a367fdd2aa735 Mon Sep 17 00:00:00 2001 From: Kevin Date: Thu, 7 Nov 2024 21:15:33 -0500 Subject: [PATCH 35/36] Everything should be functional now. --- examples/g2o_experiment.py | 15 +-- examples/g2o_mac_cache.py | 148 ------------------------ examples/greedy_esp_example.py | 50 -------- examples/pose_graph_utils.py | 2 +- examples/random_graph_sparsification.py | 28 ++--- mac/solvers/__init__.py | 3 + mac/utils/fiedler.py | 2 +- 7 files changed, 26 insertions(+), 222 deletions(-) delete mode 100644 examples/g2o_mac_cache.py delete mode 100644 examples/greedy_esp_example.py create mode 100644 mac/solvers/__init__.py diff --git a/examples/g2o_experiment.py b/examples/g2o_experiment.py index 81c65b8..df28510 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, GreedyESP +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,13 @@ 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) greedy_esp = GreedyESP(odom_edges, lc_edges, num_poses, lazy=True) ############################# @@ -318,7 +315,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) 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_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/pose_graph_utils.py b/examples/pose_graph_utils.py index 0c83e1a..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 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/solvers/__init__.py b/mac/solvers/__init__.py new file mode 100644 index 0000000..a8bad2c --- /dev/null +++ b/mac/solvers/__init__.py @@ -0,0 +1,3 @@ +from .mac import MAC +from .baseline import NaiveGreedy +from .greedy_esp import GreedyESP diff --git a/mac/utils/fiedler.py b/mac/utils/fiedler.py index c81b8a1..98539ce 100644 --- a/mac/utils/fiedler.py +++ b/mac/utils/fiedler.py @@ -36,7 +36,7 @@ def find_fiedler_pair(L, X=None, method='tracemin_lu', tol=1e-8, seed=None): assert X.shape[1] == q if method == 'tracemin_cholesky': - from mac.cholesky_utils import tracemin_fiedler_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) From 2977f0dd0c27dd3d9639a39289b68917b4ae1d64 Mon Sep 17 00:00:00 2001 From: Kevin Date: Thu, 7 Nov 2024 21:20:27 -0500 Subject: [PATCH 36/36] Rm imported GreedyESP depending on cholesky utils. --- examples/g2o_experiment.py | 11 +++-------- mac/solvers/__init__.py | 1 - 2 files changed, 3 insertions(+), 9 deletions(-) diff --git a/examples/g2o_experiment.py b/examples/g2o_experiment.py index df28510..9b9c240 100644 --- a/examples/g2o_experiment.py +++ b/examples/g2o_experiment.py @@ -6,7 +6,7 @@ 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.solvers import MAC, NaiveGreedy, GreedyESP +from mac.solvers import MAC, NaiveGreedy from mac.utils.graphs import Edge from mac.utils.rounding import round_madow @@ -272,6 +272,7 @@ def to_sesync_format(measurements): # Make a GreedyEig Solver if run_greedy: + from mac.solvers.greedy_esp import GreedyESP greedy_esp = GreedyESP(odom_edges, lc_edges, num_poses, lazy=True) ############################# @@ -334,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/mac/solvers/__init__.py b/mac/solvers/__init__.py index a8bad2c..b54ea33 100644 --- a/mac/solvers/__init__.py +++ b/mac/solvers/__init__.py @@ -1,3 +1,2 @@ from .mac import MAC from .baseline import NaiveGreedy -from .greedy_esp import GreedyESP