Merge pull request #10 from vgoliber/add-numpy-version

Add numpy version

README.md

@@ -111,15 +111,18 @@ problems of a wide-variety of sizes. For more information on setting this
 parameter, see D-Wave's [Problem Formulation
 Guide](https://www.dwavesys.com/practical-quantum-computing-developers).
 
-## Ocean Features
-
-This code example utilizes Ocean's ```AdjVectorBQM``` functionality. For
-smaller problems we can use Python dictionaries to store a BQM. However, for
-large, real-world sized problems, using dictionaries to store the BQM biases
-can become quite slow. Using NumPy arrays instead allows Python to run quickly,
-and is much more efficient on large problems. The Ocean ```AdjVectorBQM```
-functions allow the user to store biases as numpy arrays and load them quickly
-to build a BQM object, suitable for both quantum and hybrid solvers.
+## Faster BQM Construction
+
+An alternative demo file, `demo_numpy.py`, shows how the BQM for this problem
+can be constructed using NumPy arrays and vectors. Utilizing NumPy and matrix
+operations allows for a much faster construction of the BQM than building it
+with for-loops. As problem instances become larger and larger, it becomes more
+and more important to efficiently build the BQM to save time in the
+initialization and setup of the model. The chart below demonstrates the savings
+in classical compute time when setting up the BQM for this problem using
+for-loops versus using efficient NumPy operations in the Leap IDE.
+
+![Classical comparison](readme_imgs/runtimes.png "Classical Runtime Comparison")
 
 ## References
 

demo.py

@@ -28,7 +28,7 @@
     import matplotlib.pyplot as plt
 
 def read_in_args():
-    """ Read in user specified parameters or use defaults."""
+    """Read in user specified parameters or use defaults."""
 
     # Set up user-specified optional arguments
     parser = argparse.ArgumentParser()
@@ -64,7 +64,22 @@ def read_in_args():
     return args
 
 def set_up_scenario(w, h, num_poi, num_cs):
-    """ Build scenario set up with specified parameters. """
+    """Build scenario set up with specified parameters.
+    
+    Args:
+        w (int): Width of grid
+        h (int): Height of grid
+        num_poi (int): Number of points of interest
+        num_cs (int): Number of existing charging stations
+    
+    Returns:
+        G (networkx graph): Grid graph of size w by h
+        pois (list of tuples of ints): A fixed set of points of interest
+        charging_stations (list of tuples of ints): 
+            Set of current charging locations
+        potential_new_cs_nodes (list of tuples of ints): 
+            Potential new charging locations
+    """
 
     G = nx.grid_2d_graph(w, h)
     nodes = list(G.nodes)
@@ -84,7 +99,21 @@ def distance(a, b):
     return (a[0]**2 - 2*a[0]*b[0] + b[0]**2) + (a[1]**2 - 2*a[1]*b[1] + b[1]**2)
 
 def build_bqm(potential_new_cs_nodes, num_poi, pois, num_cs, charging_stations, num_new_cs):
-    """ Build bqm that models our problem scenario for the hybrid sampler. """
+    """Build bqm that models our problem scenario for the hybrid sampler. 
+
+    Args:
+        potential_new_cs_nodes (list of tuples of ints):
+            Potential new charging locations
+        num_poi (int): Number of points of interest
+        pois (list of tuples of ints): A fixed set of points of interest
+        num_cs (int): Number of existing charging stations
+        charging_stations (list of tuples of ints): 
+            Set of current charging locations
+        num_new_cs (int): Number of new charging stations desired
+    
+    Returns:
+        bqm_np (BinaryQuadraticModel): QUBO model for the input scenario
+    """
 
     # Tunable parameters
     gamma1 = len(potential_new_cs_nodes) * 4
@@ -94,7 +123,7 @@ def build_bqm(potential_new_cs_nodes, num_poi, pois, num_cs, charging_stations,
 
     # Build BQM using adjVectors to find best new charging location s.t. min
     # distance to POIs and max distance to existing charging locations
-    bqm = dimod.AdjVectorBQM(len(potential_new_cs_nodes), 'BINARY')
+    bqm = dimod.BinaryQuadraticModel(len(potential_new_cs_nodes), 'BINARY')
 
     # Constraint 1: Min average distance to POIs
     if num_poi > 0:
@@ -128,7 +157,19 @@ def build_bqm(potential_new_cs_nodes, num_poi, pois, num_cs, charging_stations,
     return bqm
 
 def run_bqm_and_collect_solutions(bqm, sampler, potential_new_cs_nodes, **kwargs):
-    """ Solve the bqm with the provided sampler to find new charger locations. """
+    """Solve the bqm with the provided sampler to find new charger locations. 
+    
+    Args:
+        bqm (BinaryQuadraticModel): The QUBO model for the problem instance
+        sampler: Sampler or solver to be used
+        potential_new_cs_nodes (list of tuples of ints): 
+            Potential new charging locations
+        **kwargs: Sampler-specific parameters to be used
+    
+    Returns:
+        new_charging_nodes (list of tuples of ints): 
+            Locations of new charging stations
+    """
 
     sampleset = sampler.sample(bqm,
                                label='Example - EV Charger Placement',
@@ -140,7 +181,21 @@ def run_bqm_and_collect_solutions(bqm, sampler, potential_new_cs_nodes, **kwargs
     return new_charging_nodes
 
 def printout_solution_to_cmdline(pois, num_poi, charging_stations, num_cs, new_charging_nodes, num_new_cs):
-    """ Print solution statistics to command line. """
+    """Print solution statistics to command line.
+    
+    Args:
+        pois (list of tuples of ints): A fixed set of points of interest
+        num_poi (int): Number of points of interest
+        charging_stations (list of tuples of ints): 
+            A fixed set of current charging locations
+        num_cs (int): Number of existing charging stations
+        new_charging_nodes (list of tuples of ints): 
+            Locations of new charging stations
+        num_new_cs (int): Number of new charging stations desired
+    
+    Returns:
+        None.
+    """
 
     print("\nSolution returned: \n------------------")
 
@@ -170,6 +225,17 @@ def save_output_image(G, pois, charging_stations, new_charging_nodes):
             - Red nodes: current charger location
             - Nodes marked 'P': POI locations
             - Blue nodes: new charger locations
+
+    Args:
+        G (networkx graph): Grid graph of size w by h
+        pois (list of tuples of ints): A fixed set of points of interest
+        charging_stations (list of tuples of ints): 
+            A fixed set of current charging locations
+        new_charging_nodes (list of tuples of ints): 
+            Locations of new charging stations
+    
+    Returns:
+        None. Output saved to file "map.png".
     """
 
     fig, (ax1, ax2) = plt.subplots(1, 2)

demo_numpy.py

@@ -0,0 +1,131 @@
+# Copyright 2021 D-Wave Systems Inc.
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+import numpy as np
+import dimod
+from dwave.system import LeapHybridSampler
+
+import demo
+
+def build_bqm(potential_new_cs_nodes, num_poi, pois, num_cs, charging_stations, num_new_cs):
+    """Build bqm that models our problem scenario using NumPy. 
+
+    Args:
+        potential_new_cs_nodes (list of tuples of ints): 
+            Potential new charging locations
+        num_poi (int): Number of points of interest
+        pois (list of tuples of ints): A fixed set of points of interest
+        num_cs (int): Number of existing charging stations
+        charging_stations (list of tuples of ints): 
+            A fixed set of current charging locations
+        num_new_cs (int): Number of new charging stations desired
+    
+    Returns:
+        bqm_np (BinaryQuadraticModel): QUBO model for the input scenario
+    """
+
+    # Tunable parameters
+    gamma1 = len(potential_new_cs_nodes) * 4.
+    gamma2 = len(potential_new_cs_nodes) / 3.
+    gamma3 = len(potential_new_cs_nodes) * 1.7
+    gamma4 = len(potential_new_cs_nodes) ** 3
+
+    # Build BQM using adjVectors to find best new charging location s.t. min
+    # distance to POIs and max distance to existing charging locations
+    linear = np.zeros(len(potential_new_cs_nodes))
+
+    nodes_array = np.asarray(potential_new_cs_nodes)
+    pois_array = np.asarray(pois)
+    cs_array = np.asarray(charging_stations)
+
+    # Constraint 1: Min average distance to POIs
+    if num_poi > 0:
+
+        ct_matrix = (np.matmul(nodes_array, pois_array.T)*(-2.) 
+                    + np.sum(np.square(pois_array), axis=1).astype(float) 
+                    + np.sum(np.square(nodes_array), axis=1).reshape(-1,1).astype(float))
+
+        linear += np.sum(ct_matrix, axis=1) / num_poi * gamma1
+
+    # Constraint 2: Max distance to existing chargers
+    if num_cs > 0:    
+
+        dist_mat = (np.matmul(nodes_array, cs_array.T)*(-2.) 
+                    + np.sum(np.square(cs_array), axis=1).astype(float) 
+                    + np.sum(np.square(nodes_array), axis=1).reshape(-1,1).astype(float))
+
+        linear += -1 * np.sum(dist_mat, axis=1) / num_cs * gamma2 
+
+    # Constraint 3: Max distance to other new charging locations
+    if num_new_cs > 1:
+
+        dist_mat = -gamma3*((np.matmul(nodes_array, nodes_array.T)*(-2.) 
+                    + np.sum(np.square(nodes_array), axis=1)).astype(float) 
+                    + np.sum(np.square(nodes_array), axis=1).reshape(-1,1).astype(float))
+
+    # Constraint 4: Choose exactly num_new_cs new charging locations
+    linear += (1-2*num_new_cs)*gamma4
+    dist_mat += 2*gamma4
+    dist_mat = np.triu(dist_mat, k=1).flatten()
+
+    quad_col = np.tile(np.arange(len(potential_new_cs_nodes)), len(potential_new_cs_nodes))
+    quad_row = np.tile(np.arange(len(potential_new_cs_nodes)), 
+                (len(potential_new_cs_nodes),1)).flatten('F')
+
+    q2 = quad_col[dist_mat != 0]
+    q1 = quad_row[dist_mat != 0]
+    q3 = dist_mat[dist_mat != 0]
+
+    bqm_np = dimod.BinaryQuadraticModel.from_numpy_vectors(linear=linear, 
+                                                            quadratic=(q1, q2, q3), 
+                                                            offset=0, 
+                                                            vartype=dimod.BINARY)
+
+    return bqm_np
+
+if __name__ == '__main__':
+
+    # Collect user inputs
+    args = demo.read_in_args()
+
+    # Build large grid graph for city
+    G, pois, charging_stations, potential_new_cs_nodes = demo.set_up_scenario(args.width, 
+                                                                            args.height, 
+                                                                            args.poi, 
+                                                                            args.chargers)
+
+    # Build BQM
+    bqm = build_bqm(potential_new_cs_nodes, 
+                    args.poi, 
+                    pois, 
+                    args.chargers, 
+                    charging_stations, 
+                    args.new_chargers)
+
+    # Run BQM on HSS
+    sampler = LeapHybridSampler()
+    print("\nRunning scenario on", sampler.solver.id, "solver...")
+
+    new_charging_nodes = demo.run_bqm_and_collect_solutions(bqm, sampler, potential_new_cs_nodes)
+
+    # Print results to commnand-line for user
+    demo.printout_solution_to_cmdline(pois, 
+                                    args.poi, 
+                                    charging_stations, 
+                                    args.chargers, 
+                                    new_charging_nodes, 
+                                    args.new_chargers)
+
+    # Create scenario output image
+    demo.save_output_image(G, pois, charging_stations, new_charging_nodes)

tests/test_integration.py

@@ -19,9 +19,11 @@
 import random
 
 import neal
+import dimod
 import numpy as np
 
 import demo
+import demo_numpy
 
 project_dir = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))
 
@@ -105,5 +107,19 @@ def test_solution_quality(self):
 
         self.assertGreater(new_cs_dist, 10)
 
+    def test_same_bqm(self):
+        """Run demo.py and demo_numpy.py with same inputs to check same BQM created."""
+
+        w, h = (random.randint(10,20), random.randint(10,20))
+        num_poi, num_cs, num_new_cs = (random.randint(1,4), random.randint(1,4), random.randint(1,4))
+
+        G, pois, charging_stations, potential_new_cs_nodes = demo.set_up_scenario(w, h, num_poi, num_cs)
+
+        bqm = demo.build_bqm(potential_new_cs_nodes, num_poi, pois, num_cs, charging_stations, num_new_cs)
+        bqm_np = demo_numpy.build_bqm(potential_new_cs_nodes, num_poi, pois, num_cs, charging_stations, num_new_cs)
+        bqm_np.add_offset(bqm.offset)
+
+        dimod.testing.asserts.assert_bqm_almost_equal(bqm, bqm_np)
+
 if __name__ == '__main__':
     unittest.main()