#!/usr/bin/env python # coding: utf-8 # For better rendering of the DataFrames [view this notebook on nbviewer.org](https://nbviewer.org/github/klezm/QuantumAnnealingPlayground/blob/main/choosing_boxes_problem/choosing_boxes.ipynb) # or [open it on the D-Wave Leap IDE](https://ide.dwavesys.io/#https://github.com/klezm/QuantumAnnealingPlayground) # In[ ]: import dimod try: dimod.ExactCQMSolver except: get_ipython().system('pip install -U dwave-ocean-sdk dimod') # !pip install pandas numpy sympy matplotlib tabulate import os os._exit(0) # restarts the ipython kernel # # Choosing Boxes Problem # as shown in the [problem formulation guide (chapt. 4.4)](https://www.dwavesys.com/media/bu0lh5ee/problem-formulation-guide-2022-01-10.pdf). # # # - https://github.com/dwave-training/choosing-boxes # ## Problem # # Pick 2 Boxes and minimize the summed weight. # # | Box | Weight | # | ----- | ------ | # | $x_1$ | 15 | # | $x_2$ | 20 | # | $x_3$ | 25 | # # QUBO: # # $$ # \begin{align*} # \text{Objective:}\quad & min(15x_1 + 20x_2 + 25x_3) \\ # \text{Constraint:}\quad & x_1 + x_2 + x_3 = 2 # \end{align*} # $$ # # ## Mathematical formulation # # For $\gamma = 1$ (Lagrange Multiplier) # # $$ # \begin{align} # & min # ( # (15x_1 + 20x_2 + 25x_3) # # +\gamma # { # \underbrace{ # (x_1+x_2+x_3-2)^2 # }_{ # \underbrace{ # x_1^2 # - 4 x_1 # + x_2^2 # - 4 x_2 # + x_3^2 # - 4 x_3 # # + 2 x_1 x_2 # + 2 x_1 x_3 # + 2 x_2 x_3 # # + 4 # }_{ # - 3x_1 - 3x_2 - 3x_3 # + 2x_1x_2 + 2x_1x_3 + 2x_2x_3 # + 4 # \quad # (x_i^2 = x_i \ \text{with} \ x_i \in \{0,1\}) # } # } # } # ) \\ # # & = min # \left( # 15x_1 # - 3x_1 # + 20x_2 # - 3x_2 # + 25x_3 # - 3x_3 # + 2x_1x_2 # + 2x_1x_3 # + 2x_2x_3 # + 4 # \right) \\ # # & = min # \left( # \underbrace{ # 12x_1 # + 17x_2 # + 22x_3 # }_{\text{linear terms}} # + # \underbrace{ # 2x_1x_2 # + 2x_1x_3 # + 2x_2x_3 # }_{\text{quadratic terms}} # + 4 # \right) \\ # # & \Rightarrow # Q = # \left( # \begin{matrix} # 12 & 2 & 2 \\ # & 17 & 2 \\ # & & 22 \\ # \end{matrix} # \right) # \end{align} # $$ # ## Content # # 1. [BQM & CQM](#BQM-%26-CQM) - Just passing the QUBO to a solver and let it do the rest # 1. [Lagrange Multiplier (γ) Comparison](#Lagrange-Multiplier-%28%CE%B3%29-Comparison) - How does γ effect the QUBO? # 2. [QUBO Matrix](#QUBO-Matrix) - Calculating the QUBO matrix from scratch and then passing it to a solver. # ## BQM & CQM # # We can use equation (1) to build a BQM or the objective & constraint functions to build a CQM. # Both can be used to solve the problem on a D-Wave QA. # # [D-Wave Docs: Quadratic Models](https://docs.dwavesys.com/docs/latest/c_gs_workflow.html#quadratic-models) # [D-Wave Docs: Problem Parameters](https://docs.dwavesys.com/docs/latest/c_solver_problems.html#bqm) # In[17]: from IPython.display import Markdown from pprint import pprint from tabulate import tabulate # from collections import defaultdict import numpy as np import pandas as pd import sympy import dwave.system from dwave.system import DWaveSampler, EmbeddingComposite, LeapHybridCQMSampler import dimod # import dwave.inspector import networkx as nx import matplotlib.pyplot as plt import seaborn as sns # In[2]: import os # print("DWAVE_API_TOKEN" in os.environ) # print(os.environ.get("DWAVE_API_TOKEN")) # # os.environ["DWAVE_API_TOKEN"] = "YOUR_API_TOKEN" # print(os.environ.get("DWAVE_API_TOKEN")) if "DWAVE_API_TOKEN" not in locals() or "DWAVE_API_TOKEN" not in os.environ: DWAVE_API_TOKEN = input("paste your D-Wave API token here:") os.environ["DWAVE_API_TOKEN"] = DWAVE_API_TOKEN # In[3]: def bqm2matrix(bqm: dimod.BinaryQuadraticModel) -> np.ndarray: bqm_np = bqm.to_numpy_vectors() # bqm.variables.to_serializable() bqm_m = np.diag(bqm_np.linear_biases) for ir, ic, b in zip(bqm_np.quadratic.row_indices, bqm_np.quadratic.col_indices, bqm_np.quadratic.biases): bqm_m[ir, ic] = b return bqm_m # bqm2matrix(bqm) # ### BQM # In[283]: γ = 35 x1 = dimod.Binary("x1") x2 = dimod.Binary("x2") x3 = dimod.Binary("x3") bqm = 15 * x1 + 20 * x2 + 25 * x3 + γ * (x1 + x2 + x3 - 2)**2 print(bqm, "\n") # In[81]: # G = bqm.to_networkx_graph() # nx.draw_networkx(G) # G = nx.from_numpy_matrix(bqm.to_numpy_matrix(bqm.variables.to_serializable())) G = nx.from_numpy_matrix(bqm2matrix(bqm)) pos = nx.circular_layout(G) pos_labels = {k: v * 1.1 for k, v in pos.items()} nx.draw(G, pos, with_labels = True) nx.draw_networkx_edge_labels( G, pos_labels, # edge_labels = {e: G.edges[e]["bias"] for e in G.edges}, edge_labels = {e: G.edges[e]["weight"] for e in G.edges}, ) plt.title("QUBO Matrix as Graph") plt.show() nx.to_numpy_array(G) # In[289]: sampler_bqm = EmbeddingComposite(DWaveSampler()) # token = DWAVE_API_TOKEN sampleset_bqm = sampler_bqm.sample(bqm, num_reads = 100, label = f'Choosing Boxes BQM (γ = {γ})') sampleset_bqm_df = sampleset_bqm.to_pandas_dataframe() # Reintroduce the offset to get the "real" energies (not the scaled ones) sampleset_bqm_df["energy_no_offset"] = sampleset_bqm_df.energy - bqm.offset display(sampleset_bqm_df.sort_values(by = "energy")) print(str(sampleset_bqm).splitlines()[-1]) # #### Exact Solver # In[ ]: γ = 60 x1, x2, x3 = [dimod.Binary(f'x{i}') for i in range(1, 4)] bqm = 15 * x1 + 20 * x2 + 25 * x3 + γ * (x1 + x2 + x3 - 2)**2 qvars = bqm.variables.to_serializable() sampleset_bqm_exact = dimod.ExactSolver().sample(bqm) # sampleset_bqm_exact = dimod.ExactDQMSolver().sample_dqm(bqm) # sampleset_bqm_exact = dimod.ExactPolySolver().sample(bqm) sampleset_bqm_exact_df = sampleset_bqm_exact.to_pandas_dataframe().assign(constraint_satisfied = lambda r: r.loc[:, qvars].sum(axis = 1) == 2) display(sampleset_bqm_exact_df.sort_values(by = "energy")) display(sampleset_bqm_exact_df.sort_values(by = qvars)) print(str(sampleset_bqm_exact).splitlines()[-1]) # #### Eigenvalues Visualization # In[ ]: def scale_cmap_log(cmap = plt.cm.gist_ncar, plot = False) -> plt.matplotlib.colors.LinearSegmentedColormap: """ Scales a colormap logarithmically. """ n = 10000 # number of samples from original colormap logs = np.geomspace(1, 256 + 1, n, dtype = int) - 1 log_colors = [cmap(i) for i in range(n)] log_colors = [log_colors[i] for i in logs] cmap_log = cmap.from_list(cmap.name + "_log", log_colors, len(log_colors)) if plot: plt.imshow([np.linspace(0, 1, 256)] , aspect = "auto", cmap = cmap_log) plt.gcf().set_size_inches(8.5, .5) plt.title("Logarithmic Colormap") plt.gca().yaxis.set_visible(False) plt.show() return cmap_log # In[556]: dfs = [] for γ in range(26, 35): x1, x2, x3 = [dimod.Binary(f'x{i}') for i in range(1, 4)] bqm = 15 * x1 + 20 * x2 + 25 * x3 + γ * (x1 + x2 + x3 - 2)**2 qvars = bqm.variables.to_serializable() sampleset_bqm_exact = dimod.ExactSolver().sample(bqm) sampleset_bqm_exact_df = sampleset_bqm_exact.to_pandas_dataframe().assign(constraint_satisfied = lambda r: r.loc[:, qvars].sum(axis = 1) == 2) if True: # Reintroduce the offset to get the "real" energies (not the scaled ones) sampleset_bqm_exact_df.energy -= bqm.offset # print(bqm.energies((dimod.reference.samplers.exact_solver._graycode(bqm), list(bqm.variables))) - bqm.offset) dfs.append([γ, sampleset_bqm_exact_df.sort_values(by = "energy")]) # display(pd.concat([x[1] for x in dfs], axis = 1, ).T.reset_index().set_index("index", append = True).T.style.background_gradient(subset = [(i, "energy") for i in [3, 9, 15]], cmap = "RdBu_r")) # pprint(df.columns.get_locs([slice(None), "energy"])) # pprint(df.columns.get_loc_level("energy", 1, drop_level = False)[1]) display( pd.concat( [x[1] for x in dfs], keys = [x[0] for x in dfs], axis = 1 ).T.drop_duplicates( ).T.astype( int ).rename_axis( columns = ["γ", ""], ).style.background_gradient( subset = pd.IndexSlice[:, pd.IndexSlice[:, "energy"]], cmap = scale_cmap_log(plt.cm.gist_ncar, True).reversed(), ) ) display(Markdown("Here we can see ... TODO ")) display(pd.concat( [x.energy.rename(index = f'γ = {γ}') for γ, x in dfs], axis = 1 ).rename_axis( index = "solution set", columns = "energy @", ).astype( int ).style.background_gradient( axis = None, cmap = scale_cmap_log(plt.cm.gist_ncar).reversed(), ) ) # In[558]: # Here: with added offset # Reintroduce the offset to get the "real" energies (not the scaled ones) # sampleset_bqm_exact_df.energy -= bqm.offset # ### CQM # # CQM can be used with the `LeapHybridCQMSampler()`. # With the `ExactCQMSolver()` the problem can be solved without a QPU. # In[158]: x1, x2, x3 = [dimod.Binary(f'x{i}') for i in range(1, 4)] cqm = dimod.ConstrainedQuadraticModel() cqm.set_objective(15 * x1 + 20 * x2 + 25 * x3) cqm.add_constraint(x1 + x2 + x3 == 2, label = "pick exactly 2 boxes") # pprint(vars(cqm)); print() # pprint(vars(dimod.cqm_to_bqm(cqm)[0])); print() sampler_cqm = LeapHybridCQMSampler() # print(sampler_cqm.parameters, "\n") sampleset_cqm = sampler_cqm.sample_cqm(cqm, label = f'Choosing Boxes CQM (γ = {γ})') display(sampleset_cqm.to_pandas_dataframe().sort_values(by = "energy").convert_dtypes().style.background_gradient(subset = "energy")) print(str(sampleset_cqm).splitlines()[-1]) # In[279]: # pprint(vars(cqm)) # pprint(cqm.objective); print(cqm.constraints); pprint(cqm.discrete); pprint(cqm.variables) # print(f'{cqm.check_feasible(sampleset_cqm.first.sample) = }') # # for sample in sampleset_cqm.samples(): print(f'{cqm.check_feasible(sample) = }') # pprint(list(cqm.iter_constraint_data(sampleset_cqm.first.sample))) # # for sample in sampleset_cqm.samples(): pprint(list(cqm.iter_constraint_data(sample))) # pprint(list(cqm.iter_violations(sampleset_cqm.first.sample))) # # for sample in sampleset_cqm.samples(): pprint(list(cqm.iter_violations(sample))) # print(f'{cqm.num_biases() = }') # print(f'{cqm.num_quadratic_variables() = }') # pprint(cqm.violations(sampleset_cqm.first.sample)) # # for sample in sampleset_cqm.samples(): pprint(cqm.violations(sample)) # pprint(sampleset_cqm.info) # sampleset_cqm.record # list(sampleset_cqm.samples()) # dimod.as_samples(sampleset_cqm) # pprint(dimod.cqm_to_bqm(cqm)) # #### Exact CQM Solver # In[43]: x1, x2, x3 = [dimod.Binary(f'x{i}') for i in range(1, 4)] # symbols = [x1, x2, x3] cqm = dimod.ConstrainedQuadraticModel() cqm.set_objective(15 * x1 + 20 * x2 + 25 * x3) cqm.add_constraint(x1 + x2 + x3 == 2, label = "pick exactly 2 boxes") qvars = cqm.variables.to_serializable() # pprint(vars(cqm)); print() pprint(vars(dimod.cqm_to_bqm(cqm)[0])); print() sampleset_cqm = dimod.ExactCQMSolver().sample_cqm(cqm) # print(sampleset_cqm) # print(sampleset_cqm.to_pandas_dataframe().to_string()) display(sampleset_cqm.to_pandas_dataframe().sort_values(by = "energy")) display(sampleset_cqm.to_pandas_dataframe().sort_values(by = qvars)) print(str(sampleset_cqm).splitlines()[-1]) # In[ ]: # ### Lagrange Multiplier (γ) Comparison # # Comparison of different values for the Lagrange Multiplier. # # The [D-Wave "Problem Formulation Guide"](https://www.dwavesys.com/media/bu0lh5ee/problem-formulation-guide-2022-01-10.pdf) suggests in chapter 4.3.4: # > Note that you may need to try a few different values to identify the best Lagrange parameter value for your specific BQM. # > **A good starting value is to set γ equal to your best estimate of the value of your objective function.** # > If you find that your constraints are not satisfied in the solutions returned, you may need to increase γ. # > On the other hand, if your constraints are all satisfied but your solutions are not close to optimal, you may need to decrease γ. # # In our problem we see that the optimal solution would be to pick box 1 & 2 with a weight of 35, which we can pick as the "best estimate" for γ. # In[66]: samplesets = {} for γ in [1, 10, 20, 30, 35, 40, 60, 80, 100, 120, 150, 200]: sampler = EmbeddingComposite(DWaveSampler()) x1, x2, x3 = [dimod.Binary(f'x{i}') for i in range(1, 4)] bqm = 15 * x1 + 20 * x2 + 25 * x3 + γ * (x1 + x2 + x3 - 2)**2 # Q = make_Q_dict(make_qubo(γ = γ)) sampleset = sampler.sample(bqm, num_reads = 100, label = f'Choosing Boxes (γ = {γ})') # sampleset = sampler.sample_qubo(Q, num_reads = 100) print(f'{sampleset}\t\twith\tγ = {γ}\n') samplesets[γ] = sampleset # In the following DataFrame the results of the sampleset with differen $\gamma$ are shown. # # We see that for $\gamma = 1$ and $\gamma = 10$ no satisfying solution were found. # We should have picked 2 boxes (not just 1 or even 0). # Therefore we should pick a higher Lagrange Multiplier to give the constraint more weight. # # For $\gamma = 35$ we find combinations for picking 1 or 2 Boxes. # First all solutions where 2 boxes were picked with the weight in descending order and then the same for cases where only one box was picked. # # For very high $\gamma$ we see that the constraint has "enough" weight, so no solutions violating the constraint have been found. # In[72]: df = pd.concat({γ: s.to_pandas_dataframe() for γ, s in samplesets.items()}).rename_axis(index = ["γ", "#"]) df = df.convert_dtypes() df = df.rename(columns = {0: "x1", 1: "x2", 2: "x3"}) # df["result"] = (df[0].astype(str) + df[1].astype(str) + df[2].astype(str)).apply(lambda x: int(f'0b{x}', 2)) df["result"] = df["x1"] * 15 + df["x2"] * 20 + df["x3"] * 25 # summed weight df["#picked"] = df["x1"] + df["x2"] + df["x3"] # number of picked boxes (constraint) # display(df.style.bar(subset = ["energy", "num_occurrences"]).background_gradient(subset = "result", cmap = "Blues_r").background_gradient(subset = "#picked")) display( df.style.bar( subset = ["chain_break_fraction", "energy", "num_occurrences"], color = "#ff92a585", ).set_properties( **{"background-color": "#adadad69"}, # #004ba752 subset = pd.IndexSlice[pd.IndexSlice[df.index.unique(level = 0)[::2], :], :] ).background_gradient( subset = "result", cmap = "Blues_r" ).set_properties( subset = "#picked", **{"font-weight": "bold"} ).text_gradient( subset = "#picked", cmap = "tab10" ) ) # In[ ]: # Some not very helpful plots :D df2 = df.T.stack() # df2.drop(columns = [0, 1, 2, "#"]).groupby(by = "γ").plot() # df2.loc[["energy", "num_occurrences"]].plot() df2.loc[["num_occurrences"]].plot(); df2.loc[["energy", ]].plot(); # df2.loc[["chain_break_fraction", ]].plot() # ## QUBO Matrix # # Here we repeat the calculation from above, but we use `sympy` to create the QUBO matrix ourselves. # In[256]: # σx = np.array([[0, 1], # [1, 0]]) # σy = np.array([[0, -1j], # [1j, 0]]) # σz = np.array([[1, 0], # [0, -1]]) # In[82]: # In[83]: def make_Q(qubo: sympy.Expr) -> np.ndarray: """ takes a QUBO sympy expression and returns the QUBO matrix. """ qubo = sympy.cancel(qubo) # simplified symbols = qubo.free_symbols qubo = qubo.subs({x ** 2: x for x in symbols}) # x_i^2 = x_i with x_i ∈ {0,1} qubo = qubo.as_poly() qubo_dict = qubo.as_dict() qubo_dict = sorted(tuple(qubo_dict.items()), reverse = True) num_vars = len(symbols) Q = np.zeros((num_vars, num_vars), dtype = int) def group_rows(qubo: dict): rows = {} for i in range(num_vars): row = [(coord, coef) for ix, (coord, coef) in enumerate(qubo) if coord[i] == 1] rows[i] = row for elem in row: qubo.remove(elem) rows = {k: {"linear": [e for e in v if sum(e[0]) == 1], "quadratic": [e for e in v if sum(e[0]) > 1]} for k, v in rows.items()} return rows rows = group_rows(qubo_dict) def fill_Q(rows, Q): for ir, row in rows.items(): for (ic, col) in row["quadratic"]: Q[(ir, ic.index(1, ir + 1))] = col for (ic, col) in row["linear"]: Q[(ir, ir)] = col return Q Q = fill_Q(rows, Q) return Q # make_Q(sympy.cancel(15 * sympy.Symbol("x1") + 20 * sympy.Symbol("x2") + 25 * sympy.Symbol("x3") + 1 * (sympy.Symbol("x1") + sympy.Symbol("x2") + sympy.Symbol("x3") - 2)**2)) # In[95]: def make_qubo(γ = 35, verbose = False) -> np.ndarray: """ make QUBO for the choosing boxes problem: min(15a + 20b + 25c + γ * (a + b + c - 2)^2) Args: γ: int or list of ints ( in the latter case the outputs for all will be printed and the last one returned) """ γs = γ if isinstance(γ, list) else [γ] γ = γs[-1] symbols = sympy.symbols("x1 x2 x3") x1, x2, x3 = symbols if verbose: print("QUBO:") γ_sym = sympy.Symbol("γ") display(15 * x1 + 20 * x2 + 25 * x3 + γ_sym * (x1 + x2 + x3 - 2)**2) for γ in γs: qubo = sympy.cancel(15 * x1 + 20 * x2 + 25 * x3 + γ * (x1 + x2 + x3 - 2)**2) # simplified qubo = qubo.subs({x ** 2: x for x in symbols}) # x_i^2 = x_i with x_i ∈ {0,1} Q = make_Q(qubo) if verbose: display(Markdown(f'#### Simplified QUBO with $γ = {γ}$')) display(qubo) print(tabulate(pd.DataFrame([np.array2string(m) for m in [Q, *np.linalg.eig(Q)]]).T, ["QUBO Matrix", "Eigenvalues", "Eigenvectors"], showindex = False)) if verbose > 1: plt.imshow(Q, cmap = "turbo", interpolation = "none") #, vmin = Q[np.triu_indices_from(Q)].min()) plt.axis("off") plt.colorbar() plt.gcf().set_size_inches(2, 1) plt.show() return Q # Q = make_qubo(γ = [1, 35], verbose = 2) # In[88]: Q = make_qubo(γ = [1, 35, 80, 200], verbose = 1) # for i in [1,4,40]: make_qubo(γ = i, verbose = True) # In[ ]: Q = make_qubo(γ = [23, 24, 25, 26, 27], verbose = 1) # Calculate all possible solutions for the choosing boxes problem with 3 variables for different $γ$ : # In[ ]: # Calcuate the QUBO for various γ # For high γ (200) the dynamic of the hardware is not utilized for γ in [1, 35, 50, 200]: # for γ in [23, 24, 25, 26, 27]: display(Markdown(f'All possible solution values for $γ = {γ}$')) print(f' x1 x2 x3 γ {" " * len(str(γ))} x1 x2 x3') for x1 in [0, 1]: for x2 in [0, 1]: for x3 in [0, 1]: print(f'15 * {x1} + 20 * {x2} + 25 * {x3} + {γ} * ({x1} + {x2} + {x3} - 2)**2 \t= ', 15 * x1 + 20 * x2 + 25 * x3 + γ * (x1 + x2 + x3 - 2)**2) # Here again all possible solutions are shown, including the duplicate counts of the results per $γ$ : # In[96]: # Calcuate the QUBO for various γ # Note: for γ in [5, 10, 15, 20, 25, 30] there are duplicate solutions num_vars = 3 γs = [22, 23, 24, 25, 26, 27] γs = [1, 10, 20, 21, 22, 23, 24, 25, 26, 27, 35, 40, 80, 200][:-2] # γs = [5, 10, 15, 20, 25, 30] # duplicates # γs = range(2100) # vars_combs = np.array(sorted(np.array(np.meshgrid(*([[0, 1]] * num_vars))).reshape(-1, num_vars).tolist())) vars_combs = np.array([[int(b) for b in ("0" * num_vars + bin(x)[2:])[-num_vars:]] for x in range(2 ** num_vars)]) vars_combs = {f'x{k}': v for k, v in zip(range(1, num_vars + 1), vars_combs.T)} df = pd.DataFrame(vars_combs) df2 = pd.Series([γs] * len(df), name = "γ", index = df.index) df = df.join(df2.explode()).set_index("γ", append = True).swaplevel(0, 1, axis = 0).sort_index() df["|"] = "|" df["15"] = "" df["* x1"] = df.x1 df["+ 20"] = "" df["* x2"] = df.x2 df["+ 25"] = "" df["* x3"] = df.x3 df["+ γ"] = df.index.get_level_values(0) df["(x1"] = df.x1 df["+ x2"] = df.x2 df["+ x3"] = df.x3 df["- 2"] = "" #2 df[")^2"] = "" df["="] = "" df["result"] = 15 * df.x1 + 20 * df.x2 + 25 * df.x3 + df.index.get_level_values(0) * (df.x1 + df.x2 + df.x3 - 2)**2 # df["Equation"] = "15 * " + df.x1.astype(str) + " + 20 * {df.x2} + 25 * {df.x3} + {γ} * ({df.x1} + {df.x2} + {df.x3} - 2)**2" df["duplicates"] = df.groupby(level = "γ", axis = 0).apply(lambda x: x.duplicated(subset = "result", keep = False)).values # df_dups = df.duplicates.groupby(level = 0).sum() df_dups = df.duplicates.groupby(level = "γ").sum() display(df_dups[df_dups != 0].to_frame().rename(columns = {"duplicates": "Duplicate_count"})) # display(df.result.groupby(level = "γ").agg(["min", "max", "var", "std", "mean", "median", "sem", "sum", "skew", np.ptp]).style.bar()) display( df.style.bar( subset = ["result"], color = "#ff92a585", ).set_properties( **{"background-color": "#adadad69"}, subset = pd.IndexSlice[pd.IndexSlice[df.index.unique(level = 0)[::2], :], :] ).text_gradient( subset = [c for c in df.columns if any([x in c for x in ["x1", "x2", "x3"]]) and c not in ["x1", "x2", "x3"]], cmap = "bwr_r", # tab10 PiYG ).text_gradient( subset = "duplicates", cmap = "bwr", ) ) # In[90]: def make_Q_dict(Q: np.ndarray, verbose = False) -> dict: """ Prepare the Q matrix for DWaveSampler() """ Q = Q.astype(float) # for nan masking if verbose: pprint(Q) triu_indices = np.triu_indices(Q.shape[0]) # triu_indices = np.dstack(triu_indices)[0] mask = np.ones_like(Q).astype(bool) mask[triu_indices] = False assert np.isnan(Q[mask]).all() or (Q[mask] == 0).all(), "Lower triangle should not contain any values" # set lower triangle to np.nan Q[mask] = np.nan Q_dict = {k: v for (k, v) in np.ndenumerate(Q) if not np.isnan(v)} if verbose: pprint(Q_dict) return Q_dict Q = make_Q_dict(make_qubo(γ = 40, verbose = False), verbose = True) # for i in [1,4,40]: make_Q_dict(make_qubo(γ = i, verbose = False), verbose = True) # ### Run Problem on D-Wave QA # In[91]: Q = make_Q_dict(make_qubo(γ = 35)) sampler = EmbeddingComposite(DWaveSampler()) # Create a bqm object... bqm = dimod.BQM.from_qubo(Q) sampleset = sampler.sample(bqm, num_reads = 10, label = f'Choosing Boxes QUBO (γ = {γ})') # ...or use sample_qubo() directly. # sampleset = sampler.sample_qubo(Q, num_reads = 10) # Note that the energies returned by the QA are negative here, unlike in the solutions above. # In[7]: print(sampleset, end = "\n\n") pprint(vars(sampleset)); print() # sampleset # ### Lagrange Multiplier (γ) Comparison # # Note that the energies are not scaled here. # In[ ]: samplesets = {} # for γ in [23, 24, 25, 26, 27]: # for γ in reversed([23, 24, 25, 26, 27]): for γ in [1, 10, 20, 30, 35, 40, 60, 80, 100, 120, 150, 200]: sampler = EmbeddingComposite(DWaveSampler()) # x1, x2, x3 = [dimod.Binary(f'x{i}') for i in range(1, 4)] # bqm = 15 * x1 + 20 * x2 + 25 * x3 + γ * (x1 + x2 + x3 - 2)**2 # sampleset = sampler.sample(bqm, num_reads = 100) Q = make_Q_dict(make_qubo(γ = γ)) sampleset = sampler.sample_qubo(Q, num_reads = 100, label = f'Choosing Boxes (γ = {γ})') # print(f'>>> γ = {γ} <<<\n{sampleset}\n') samplesets[γ] = sampleset # In[ ]: df = pd.concat({γ: s.to_pandas_dataframe() for γ, s in samplesets.items()}).rename_axis(index = ["γ", "#"]) df = df.convert_dtypes() df = df.rename(columns = {0: "x1", 1: "x2", 2: "x3"}) df["result"] = df["x1"] * 15 + df["x2"] * 20 + df["x3"] * 25 # summed weight df["#picked"] = df["x1"] + df["x2"] + df["x3"] # number of picked boxes (constraint) display( df.style.bar( subset = ["chain_break_fraction", "energy", "num_occurrences"], color = "#ff92a585", ).set_properties( **{"background-color": "#adadad69"}, # #004ba752 subset = pd.IndexSlice[pd.IndexSlice[df.index.unique(level = 0)[::2], :], :] ).background_gradient( subset = "result", cmap = "Blues_r" ).set_properties( subset = "#picked", **{"font-weight": "bold"} ).text_gradient( subset = "#picked", cmap = "tab10" ) ) # Note: Changing the order of the γ's or inserting a pause in between seems to have quite an impact on the outcome. Maybe `num_reads = 100` is also not sufficient. # In[ ]: # # Single Qubit Problem # Here it should be shown if the Quantum Annealer by D-Wave exclusively finds the correct solution for all trials. # In[4]: γ = 1 x1 = dimod.Binary("x1") bqm = x1 #+ γ * (x1 - 1)**2 bqm # In[5]: sampler_bqm = EmbeddingComposite(DWaveSampler()) sampleset_bqm = sampler_bqm.sample(bqm, num_reads = 10_000, label = "Trivial 1-Qubit Problem") display(sampleset_bqm.to_pandas_dataframe()) print(str(sampleset_bqm).splitlines()[-1]) # In[32]: pd.DataFrame.from_dict(sampleset_bqm.info["timing"], orient = "index", columns = ["time [μs]"]) # The Exact Solver shows all possible solutions: # In[108]: sampleset_bqm_exact = dimod.ExactSolver().sample(bqm) # sampleset_bqm_exact = dimod.ExactDQMSolver().sample_dqm(bqm) # sampleset_bqm_exact = dimod.ExactPolySolver().sample(bqm) qvars = bqm.variables.to_serializable() sampleset_bqm_exact_df = sampleset_bqm_exact.to_pandas_dataframe() # sampleset_bqm_exact_df["constraint_satisfied"] = sampleset_bqm_exact_df.loc[:, qvars].sum(axis = 1) == 0 display(sampleset_bqm_exact_df.sort_values(by = "energy")) # display(sampleset_bqm_exact_df.sort_values(by = qvars)) print(str(sampleset_bqm_exact).splitlines()[-1])