From d56ee52ddd8e97d7a409437a7b6b3d1da9cb0403 Mon Sep 17 00:00:00 2001 From: Stefan Hannie Date: Tue, 31 Mar 2020 19:16:24 -0700 Subject: [PATCH] Layout p norm (#121) * P_norm layout A Kamada-Kawai-like algorithm with a variable p-norm distance function. Now also passing dim, center, scale as kwargs to layout callables. * Set p_norm as find_embedding default * p_norm unit test test_dimension also changed as the behavior is slightly different now. * Minor changes Layout.d changed to Layout.dim, documentation changes, split a test. * Timeout changed to perf_counter * Fixed typo added TODO #122 --- minorminer/layout/__init__.py | 15 +- minorminer/layout/layout.py | 191 ++++++++++++++++-- minorminer/layout/placement.py | 10 +- minorminer/layout/tests/common.py | 6 +- .../layout/tests/test_find_embedding.py | 35 ++-- minorminer/layout/tests/test_layout.py | 80 ++++++-- minorminer/layout/tests/test_placement.py | 4 - 7 files changed, 272 insertions(+), 69 deletions(-) diff --git a/minorminer/layout/__init__.py b/minorminer/layout/__init__.py index 50483565..4dfd262d 100644 --- a/minorminer/layout/__init__.py +++ b/minorminer/layout/__init__.py @@ -1,16 +1,17 @@ import time -import minorminer as mm import networkx as nx -from .layout import Layout +import minorminer as mm + +from .layout import Layout, p_norm from .placement import Placement, closest def find_embedding( S, T, - layout=None, + layout=p_norm, placement=closest, mm_hint_type="initial_chains", return_layouts=False, @@ -51,7 +52,7 @@ def find_embedding( Output is dependent upon kwargs passed to minonminer, but more or less emb is a mapping from vertices of S (keys) to chains in T (values). """ - start = time.process_time() + start = time.perf_counter() # Parse kwargs layout_kwargs, placement_kwargs = _parse_kwargs(kwargs) @@ -63,7 +64,7 @@ def find_embedding( S_T_placement = Placement( S_layout, T_layout, placement, **placement_kwargs) - end = time.process_time() + end = time.perf_counter() timeout = kwargs.get("timeout") if timeout is not None: time_remaining = timeout - (end - start) @@ -94,8 +95,8 @@ def _parse_kwargs(kwargs): """ layout_kwargs = {} # For the layout object - if "d" in kwargs: - layout_kwargs["d"] = kwargs.pop("d") + if "dim" in kwargs: + layout_kwargs["dim"] = kwargs.pop("dim") if "center" in kwargs: layout_kwargs["center"] = kwargs.pop("center") if "scale" in kwargs: diff --git a/minorminer/layout/layout.py b/minorminer/layout/layout.py index 6668f0b8..37113f1a 100644 --- a/minorminer/layout/layout.py +++ b/minorminer/layout/layout.py @@ -2,6 +2,158 @@ import networkx as nx import numpy as np +from scipy import optimize + + +def p_norm(G, p=2, starting_layout=None, G_distances=None, dim=None, center=None, scale=None, **kwargs): + """ + Embeds a graph in R^d with the p-norm and minimizes a Kamada-Kawai-esque objective function to achieve + an embedding with low distortion. This computes a layout where the graph distance and the p-distance are + very close to each other. + + Parameters + ---------- + G : NetworkX graph + The graph you want to compute the layout for. + p : int (default 2) + The order of the p-norm to use as a metric. + starting_layout : dict + A mapping from the vertices of G to points in R^d. + G_distances : dict (default None) + A dictionary of dictionaries representing distances from every vertex in G to every other vertex in G. If None, + it is computed. + + Returns + ------- + layout : dict + A mapping from vertices of G (keys) to points in R^d (values). + """ + # Use the user provided starting_layout, or a spectral_layout if the dimension is low enough. + # If neither, use a random_layout. + if starting_layout: + pass + elif dim: + if dim >= len(G): + starting_layout = nx.random_layout + else: + starting_layout = nx.spectral_layout + else: + starting_layout = nx.spectral_layout + + # Make a layout object + layout = Layout(G, starting_layout, dim=dim, + center=center, scale=scale) + + # Update dim to the starting layout's + dim = layout.dim + + # Save on distance calculations by passing them in + G_distances = _graph_distance_matrix(G, G_distances) + + # Solve the Kamada-Kawai-esque minimization function + X = optimize.minimize( + _p_norm_objective, + layout.layout_array.ravel(), + method='L-BFGS-B', + args=(G_distances, dim, p), + jac=True, + ) + + # Read out the solution to the minimization problem and save layouts + layout.layout_array = X.x.reshape(len(G), dim) + + return layout.layout + + +def _graph_distance_matrix(G, all_pairs_shortest_path_length=None): + """ + Compute the distance matrix of G. + + Parameters + ---------- + G : NetworkX Graph + The graph to find the distance matrix of. + all_pairs_shortest_path_length : dict (default None) + If None, it is computed by calling nx.all_pairs_shortest_path_length. + + Returns + ------- + G_distances : Numpy 2d array + An array indexed by sorted vertices of G whose i,j value is d_G(i,j). + """ + if all_pairs_shortest_path_length is None: + all_pairs_shortest_path_length = nx.all_pairs_shortest_path_length(G) + + return np.array( + [ + [V[v] for v in sorted(G)] for u, V in all_pairs_shortest_path_length + ] + ) + + +def _p_norm_objective(layout_vector, G_distances, dim, p): + """ + Compute the sum of differences squared between the p-norm and the graph distance as well as the gradient. + + Parameters + ---------- + layout : Numpy Array + A vector indexed by sorted vertices of G whose values are points in some metric space. + G_distances : Numpy 2d array + An array indexed by sorted vertices of G whose i,j value is d_G(i,j). + dim : int + The dimension of the metric space. This will reshape the flattened array passed in to the cost function. + p : int + The order of the p-norm to use. + + Returns + ------- + cost : float + The sum of differences squared between the metric distance and the graph distance. + """ + # Reconstitute the flattened array that scipy.optimize.minimize passed in + n = len(G_distances) + layout = layout_vector.reshape(n, dim) + + # Difference between pairs of points in a 3d matrix + diff = layout[:, np.newaxis, :] - layout[np.newaxis, :, :] + + # A 2d matrix of the distances between points + dist = np.linalg.norm(diff, ord=p, axis=-1) + + # TODO: Compare this division-by-zero strategy to adding epsilon. + # A vectorized version of the gradient function + with np.errstate(divide='ignore', invalid='ignore'): # handle division by 0 + if p == 1: + grad = np.einsum( + 'ijk,ij,ijk->ik', + 2*diff, + dist - G_distances, + np.nan_to_num(1/np.abs(diff)) + ) + elif p == float("inf"): + # Note: It may not be faster to do this outside of einsum + abs_diff = np.abs(diff) + x_bigger = abs_diff[:, :, 0] > abs_diff[:, :, 1] + + grad = np.einsum( + 'ijk,ij,ijk,ijk->ik', + 2*diff, + dist - G_distances, + np.nan_to_num(1/abs_diff), + np.dstack((x_bigger, np.logical_not(x_bigger))) + ) + else: + grad = np.einsum( + 'ijk,ijk,ij,ij->ik', + 2*diff, + np.abs(diff)**(p-2), + dist - G_distances, + np.nan_to_num((1/dist)**(p-1)) + ) + + # Return the cost and the gradient + return np.sum((G_distances - dist)**2), grad.ravel() class Layout(abc.MutableMapping): @@ -9,13 +161,14 @@ def __init__( self, G, layout=None, - d=None, + dim=None, center=None, scale=None, **kwargs ): """ Compute a layout for G, i.e., a map from G to R^d. + Parameters ---------- G : NetworkX graph or NetworkX supported edges data structure (dict, list, ...) @@ -23,8 +176,8 @@ def __init__( layout : dict or function (default None) If a dict, this specifies a pre-computed layout for G. If a function, the function is called on G `layout(G)` and should return a layout of G. If None, nx.spectral_layout is called. - d : int (default None) - The desired dimension of the layout, R^d. If None, set d to be the dimension of layout. + dim : int (default None) + The desired dimension of the layout, R^dim. If None, set dim to be the dimension of layout. center : tuple (default None) The desired center point of the layout. If None, set center to be the center of layout. scale : float (default None) @@ -36,9 +189,17 @@ def __init__( # Ensure G is a graph object self.G = _parse_graph(G) + # Add dim, center, and scale to kwargs if not None + if dim is not None: + kwargs["dim"] = dim + if center is not None: + kwargs["center"] = center + if scale is not None: + kwargs["scale"] = scale + # If passed in, save or compute the layout if layout is None: - self.layout = nx.spectral_layout(self.G) + self.layout = nx.spectral_layout(self.G, **kwargs) elif callable(layout): self.layout = layout(self.G, **kwargs) else: @@ -46,7 +207,7 @@ def __init__( self.layout = layout # Set specs in the order of (user input, precomputed layout) - self.d = d or self._d + self.dim = dim or self._dim self.scale = scale or self._scale if center is not None: self.center = center @@ -87,17 +248,17 @@ def layout_array(self, value): self.layout = {v: p for v, p in zip(sorted(self.G), value)} @property - def d(self): - return self._d + def dim(self): + return self._dim - @d.setter - def d(self, value): + @dim.setter + def dim(self, value): """ If the dimension is changed, change the dimension of the layout, if possible. """ if value: self.layout_array = _dimension_layout( - self.layout_array, value, self._d) + self.layout_array, value, self._dim) @property def center(self): @@ -170,11 +331,11 @@ def _set_layout_specs(self, empty=False): Set the dimension, center, and scale of the layout_array currently in the Layout object. """ if self.layout_array.size == 0: - self._d = 0 + self._dim = 0 self._center = np.array([]) self._scale = 0 else: - self._d = self.layout_array.shape[1] + self._dim = self.layout_array.shape[1] self._center = np.mean(self.layout_array, axis=0) self._scale = np.max( np.linalg.norm( @@ -186,6 +347,7 @@ def _set_layout_specs(self, empty=False): def _dimension_layout(layout_array, new_d, old_d=None): """ This helper function transforms a layout from R^old_d to R^new_d by padding extra dimensions with 0's. + Parameters ---------- layout_array : numpy array @@ -194,6 +356,7 @@ def _dimension_layout(layout_array, new_d, old_d=None): The new dimension to convert the layout to, must be larger than old_dim. old_d : int (default None) The current dimension of the laoyut. If None, the dimension is looked up via layout_array. + Returns ------- layout_array : numpy array @@ -221,6 +384,7 @@ def _center_layout(layout_array, new_center, old_center=None): """ This helper function transforms a layout from [old_center - scale, old_center + scale]^d to [new_center - scale, new_center + scale]^d. + Parameters ---------- layout_array : numpy array @@ -230,6 +394,7 @@ def _center_layout(layout_array, new_center, old_center=None): old_center : tuple or numpy array (default None) A point in R^d representing the center of the layout. If None, the approximate center of layout is computed by calculating the center of mass (or centroid). + Returns ------- layout_array : numpy array @@ -247,6 +412,7 @@ def _scale_layout(layout_array, new_scale, old_scale=None, center=None): """ This helper function transforms a layout from [center - old_scale, center + old_scale]^d to [center - new_scale, center + new_scale]^d. + Parameters ---------- layout_array : numpy array (default None) @@ -259,6 +425,7 @@ def _scale_layout(layout_array, new_scale, old_scale=None, center=None): center : tuple or numpy array (default None) A point in R^d representing the center of the layout. If None, the approximate center of layout is computed by calculating the center of mass (centroid). + Returns ------- layout_array : numpy array diff --git a/minorminer/layout/placement.py b/minorminer/layout/placement.py index be3ef867..ae5cb4ca 100644 --- a/minorminer/layout/placement.py +++ b/minorminer/layout/placement.py @@ -15,6 +15,10 @@ def closest(S_layout, T_layout, subset_size=(1, 1), num_neighbors=1, **kwargs): Parameters ---------- + S_layout : layout.Layout + A layout for S; i.e. a map from S to R^d. + T_layout : layout.Layout + A layout for T; i.e. a map from T to R^d. subset_size : tuple (default (1, 1)) A lower (subset_size[0]) and upper (subset_size[1]) bound on the size of subets of T that will be considered when mapping vertices of S. @@ -141,7 +145,7 @@ def __init__( ---------- S_layout : layout.Layout A layout for S; i.e. a map from S to R^d. - T_layout : layout.Layout or networkx.Graph + T_layout : layout.Layout A layout for T; i.e. a map from T to R^d. placement : dict or function (default None) If a dict, this specifies a pre-computed placement for S in T. If a function, the function is called on @@ -156,10 +160,10 @@ def __init__( self.T_layout = _parse_layout(T_layout) # Layout dimensions should match - if self.S_layout.d != self.T_layout.d: + if self.S_layout.dim != self.T_layout.dim: raise ValueError( "S_layout has dimension {} but T_layout has dimension {}. These must match.".format( - self.S_layout.d, self.T_layout.d) + self.S_layout.dim, self.T_layout.dim) ) # Scale S if S_layout is bigger than T_layout diff --git a/minorminer/layout/tests/common.py b/minorminer/layout/tests/common.py index 936a4293..a269450a 100644 --- a/minorminer/layout/tests/common.py +++ b/minorminer/layout/tests/common.py @@ -13,15 +13,17 @@ def __init__(self, *args, **kwargs): # Graphs for testing self.S = nx.random_regular_graph(3, 50) + self.S_small = nx.random_regular_graph(3, 10) self.G = nx.Graph() self.H = nx.complete_graph(1) self.C = dnx.chimera_graph(4) # Compute some layouts self.S_layout = mml.Layout(self.S) + self.S_small_layout = mml.Layout(self.S_small) self.G_layout = mml.Layout(self.G) self.C_layout = mml.Layout(self.C) - self.C_layout_3 = mml.Layout(self.C, d=3) + self.C_layout_3 = mml.Layout(self.C, dim=3) def assertArrayEqual(self, a, b): """ @@ -41,7 +43,7 @@ def assertIsLayout(self, S, layout): Tests that layout is a mapping from S to R^d """ for u in S: - self.assertEqual(len(layout[u]), layout.d) + self.assertEqual(len(layout[u]), layout.dim) def assertIsPlacement(self, S, T, placement): """ diff --git a/minorminer/layout/tests/test_find_embedding.py b/minorminer/layout/tests/test_find_embedding.py index 708508be..2a3ba094 100644 --- a/minorminer/layout/tests/test_find_embedding.py +++ b/minorminer/layout/tests/test_find_embedding.py @@ -11,19 +11,6 @@ class TestFindEmb(TestLayoutPlacement): - def assertIsLayout(self, S, layout): - """ - Tests that layout is a mapping from S to R^d - """ - for u in S: - self.assertEqual(len(layout[u]), layout.d) - - def assertArrayEqual(self, a, b): - """ - Tests that two arrays are equal via numpy. - """ - np.testing.assert_almost_equal(a, b) - def test_default(self): """ Minimal find_embedding call @@ -65,18 +52,18 @@ def test_layout_kwargs(self): Pass in layout kwargs. """ # Pick some values to pass in - d = 3 + dim = 3 center = (0, 0, 0) scale = 2 _, (S_layout, C_layout) = mml.find_embedding(self.S, self.C, - d=d, center=center, scale=scale, return_layouts=True) + dim=dim, center=center, scale=scale, return_layouts=True) # Test that S_layout matches - self.assertEqual(S_layout.d, d) + self.assertEqual(S_layout.dim, dim) self.assertArrayEqual(S_layout.center, center) self.assertAlmostEqual(S_layout.scale, scale) # Test that C_layout matches - self.assertEqual(C_layout.d, d) + self.assertEqual(C_layout.dim, dim) self.assertArrayEqual(C_layout.center, center) self.assertAlmostEqual(C_layout.scale, scale) @@ -86,10 +73,18 @@ def test_placement_kwargs(self): """ # Pick some values to pass in scale_ratio = .8 + + mml.find_embedding(self.S, self.C, scale_ratio=scale_ratio) + + def test_placement_closest(self): + """ + Test the closest placement strategy + """ + # Pick some values to pass in subset_size = (1, 2) num_neighbors = 5 - mml.find_embedding(self.S, self.C, scale_ratio=scale_ratio, + mml.find_embedding(self.S, self.C, placement=mml.closest, subset_size=subset_size, num_neighbors=num_neighbors) def test_layout_parameter(self): @@ -117,7 +112,3 @@ def test_layout_parameter(self): # Too many things in layout self.assertRaises(ValueError, mml.find_embedding, self.S, self.C, layout=(1, 2, 3)) - - -if __name__ == '__main__': - unittest.main() diff --git a/minorminer/layout/tests/test_layout.py b/minorminer/layout/tests/test_layout.py index 6d5dfceb..e411b436 100644 --- a/minorminer/layout/tests/test_layout.py +++ b/minorminer/layout/tests/test_layout.py @@ -13,6 +13,50 @@ class TestLayout(TestLayoutPlacement): + def test_pnorm(self): + """ + Test the p_norm layout strategy. + """ + # Some specs to test with + low_dim = random.randint(3, 9) + high_dim = len(self.S_small) + center = (1, 1) + scale = random.random()*random.randint(1, 10) + + # Default behavior + mml.p_norm(self.S_small) + + # Using a starting_layout + mml.p_norm(self.S_small, starting_layout=nx.random_layout(self.S_small)) + + # Passing in G_distances + mml.p_norm( + self.S_small, G_distances=nx.all_pairs_shortest_path_length(self.S_small)) + + # Passing in dim + mml.p_norm(self.S_small, dim=low_dim) + mml.p_norm(self.S_small, dim=high_dim) + + # Passing in center + mml.p_norm(self.S_small, center=center) + + # Passing in scale + mml.p_norm(self.S_small, scale=scale) + + # Different p-norms + mml.p_norm(self.S_small, p=1) + mml.p_norm(self.S_small, p=3) + mml.p_norm(self.S_small, p=float("inf")) + + # # Playing with a starting_layout and a dimension + # starting_layout = nx.spectral_layout(self.S, dim=dim) + # layout_in = mml.Layout( + # self.S, mml.p_norm, starting_layout=starting_layout) + # layout_out = mml.Layout( + # self.S, mml.p_norm, starting_layout=nx.spectral_layout, dim=dim) + # self.assertLayoutEqual(self.S, layout_in, layout_out) + # self.assertIsLayout(self.S, layout_in) + def test_precomputed_layout(self): """ Pass in a precomputed layout to the Layout class. @@ -30,29 +74,31 @@ def test_dimension(self): """ Change the dimension of a layout. """ - d = random.randint(3, 10) + dim = random.randint(3, 10) - # Pass in d as an argument - layout_pre = mml.Layout(self.S, d=d) - self.assertEqual(layout_pre.d, d) + # Pass in dim as an argument + layout_pre = mml.Layout(self.S, dim=dim) + self.assertEqual(layout_pre.dim, dim) + self.assertIsLayout(self.S, layout_pre) - # Change the layout to have the d + # Change the layout to have the dim layout_post = mml.Layout(self.S) - layout_post.d = d - self.assertEqual(layout_post.d, d) + layout_post.dim = dim + self.assertEqual(layout_post.dim, dim) + self.assertIsLayout(self.S, layout_post) # Change the dimension without changing the object, - layout = mml.Layout(self.S, d=2) - new_layout_array = _dimension_layout(layout.layout_array, d) - self.assertEqual(layout.d, 2) + layout = mml.Layout(self.S, dim=2) + new_layout_array = _dimension_layout(layout.layout_array, dim) + self.assertEqual(layout.dim, 2) - # The layouts should match each other - self.assertLayoutEqual(self.S, layout_pre, layout_post) - self.assertIsLayout(self.S, layout_pre) - self.assertArrayEqual(new_layout_array, layout_pre.layout_array) + # The layout_arrays changed after the fact should match each other + self.assertArrayEqual(new_layout_array, layout_post.layout_array) # Test dimension too small - self.assertRaises(ValueError, mml.Layout, self.S, d=1) + layout = mml.Layout(self.S, dim=2) + self.assertRaises(ValueError, _dimension_layout, + layout.layout_array, 1) def test_center(self): """ @@ -157,7 +203,3 @@ def test_layout_class(self): # Test __repr__ self.assertEqual(repr(L), "{}") - - -if __name__ == '__main__': - unittest.main() diff --git a/minorminer/layout/tests/test_placement.py b/minorminer/layout/tests/test_placement.py index 96c6f502..840af54f 100644 --- a/minorminer/layout/tests/test_placement.py +++ b/minorminer/layout/tests/test_placement.py @@ -106,7 +106,3 @@ def test_failures(self): def rando_placer(S_layout, T_layout): T_vertices = list(T_layout) return {v: [random.choice(T_vertices)] for v in S_layout} - - -if __name__ == '__main__': - unittest.main()