Merge pull request #24 from jpmorganchase/fix-default-behavior-overla…

…p-objective

fix default behavior of overlap objective

qokit/maxcut.py

@@ -35,14 +35,14 @@ def get_maxcut_terms(G: nx.Graph) -> TermsType:
         terms to be used in the simulation
     """
     if nx.is_weighted(G):
-        terms = [(-float(G[u][v]["weight"]) / 2, (int(u), int(v))) for u, v, *_ in G.edges()]
+        terms = [(float(G[u][v]["weight"]) / 2, (int(u), int(v))) for u, v, *_ in G.edges()]
         total_w = sum([float(G[u][v]["weight"]) for u, v, *_ in G.edges()])
 
     else:
-        terms = [(-1 / 2, (int(e[0]), int(e[1]))) for e in G.edges()]
+        terms = [(1 / 2, (int(e[0]), int(e[1]))) for e in G.edges()]
         total_w = int(G.number_of_edges())
     N = G.number_of_nodes()
-    terms.append((+total_w / 2, tuple()))
+    terms.append((-total_w / 2, tuple()))
     return terms
 
 

qokit/parameter_utils.py

@@ -11,6 +11,8 @@
 from importlib_resources import files
 from enum import Enum
 from typing import Callable
+from functools import lru_cache, cached_property, cache
+from datetime import datetime
 
 
 def from_fourier_basis(u, v):
@@ -206,9 +208,25 @@ def get_sk_gamma_beta(p, parameterization: QAOAParameterization | str = "gamma b
         raise ValueError(f"p={p} not supported, try lower p")
     parameterization = QAOAParameterization(parameterization)
     if parameterization == QAOAParameterization.THETA:
-        return np.concatenate((4 * gamma, beta), axis=0)
+        return np.concatenate((-1 * (4 * gamma), beta), axis=0)
     elif parameterization == QAOAParameterization.GAMMA_BETA:
-        return 4 * gamma, beta
+        return -1 * (4 * gamma), beta
+
+
+@cache
+def _get_gamma_beta_from_file():
+    """
+    Caches the dataframe after the first call to load JSon, subsequent calls will get from cache and save I/O
+
+    Parameters
+    ----------
+    None
+
+    Returns
+    -------
+    df: Pandas Dataframe
+    """
+    return pd.read_json(str(files("qokit.assets.maxcut_datasets").joinpath("fixed_angles_for_regular_graphs.json")), orient="index")
 
 
 def get_fixed_gamma_beta(d, p, return_AR=False):
@@ -231,16 +249,15 @@ def get_fixed_gamma_beta(d, p, return_AR=False):
     AR : float
         Only returned is flag return_AR is raised
     """
-    df = pd.read_json(str(files("qokit.assets.maxcut_datasets").joinpath("fixed_angles_for_regular_graphs.json")), orient="index")
-
+    df = _get_gamma_beta_from_file()
     row = df[(df["d"] == d) & (df["p"] == p)]
     if len(row) != 1:
         raise ValueError(f"Failed to retrieve fixed angles for d={d}, p={p}")
     row = row.squeeze()
     if return_AR:
-        return row["gamma"], row["beta"], row["AR"]
+        return -1 * np.asarray(row["gamma"]), row["beta"], row["AR"]
     else:
-        return row["gamma"], row["beta"]
+        return -1 * np.asarray(row["gamma"]), row["beta"]
 
 
 def get_best_known_parameters_for_LABS_wrt_overlap(N: int) -> pd.DataFrame:

qokit/portfolio_optimization.py

@@ -275,7 +275,7 @@ def get_sk_ini(p: int):
     """
     scaled the sk look-up table for the application of portfolio optimziation
     """
-    gamma_scale, beta_scale = -0.5, 1
+    gamma_scale, beta_scale = 0.5, 1
     gamma, beta = get_sk_gamma_beta(p, parameterization="gamma beta")
     scaled_gamma, scaled_beta = gamma_scale * gamma, beta_scale * beta
     X0 = np.concatenate((scaled_gamma, scaled_beta), axis=0)

qokit/qaoa_circuit_maxcut.py

@@ -11,7 +11,7 @@
 
 
 def append_zz_term(qc, q1, q2, gamma):
-    qc.rzz(-gamma / 2, q1, q2)
+    qc.rzz(gamma / 2, q1, q2)
 
 
 def append_maxcut_cost_operator_circuit(qc, G, gamma):

qokit/qaoa_objective.py

@@ -117,7 +117,7 @@ def get_qaoa_objective(
     N: int,
     p: int | None = None,
     precomputed_diagonal_hamiltonian=None,
-    precomputed_objectives=None,
+    precomputed_costs=None,
     terms=None,
     precomputed_optimal_bitstrings=None,
     parameterization: str | QAOAParameterization = "theta",
@@ -166,7 +166,7 @@ def get_qaoa_objective(
 
     # -- Qiskit edge case
     if simulator == "qiskit":
-        g = _get_qiskit_objective(parameterized_circuit, precomputed_objectives, precomputed_optimal_bitstrings, objective, terms, parameterization, mixer)
+        g = _get_qiskit_objective(parameterized_circuit, precomputed_costs, precomputed_optimal_bitstrings, objective, terms, parameterization, mixer)
 
         def fq(*args):
             gamma, beta = qokit.parameter_utils.convert_to_gamma_beta(*args, parameterization=parameterization)
@@ -183,8 +183,8 @@ def fq(*args):
         raise ValueError(f"Unknown mixer type passed to get_qaoa_objective: {mixer}, allowed ['x', 'xy']")
 
     sim = simulator_cls(N, terms=terms, costs=precomputed_diagonal_hamiltonian)
-    if precomputed_objectives is None:
-        precomputed_objectives = sim.get_cost_diagonal()
+    if precomputed_costs is None:
+        precomputed_costs = sim.get_cost_diagonal()
 
     bitstring_loc = None
     if precomputed_optimal_bitstrings is not None and objective != "expectation":
@@ -195,13 +195,13 @@ def f(*args):
         gamma, beta = qokit.parameter_utils.convert_to_gamma_beta(*args, parameterization=parameterization)
         result = sim.simulate_qaoa(gamma, beta, initial_state, n_trotters=n_trotters)
         if objective == "expectation":
-            return sim.get_expectation(result, costs=precomputed_objectives, preserve_state=False)
+            return sim.get_expectation(result, costs=precomputed_costs, preserve_state=False)
         elif objective == "overlap":
-            overlap = sim.get_overlap(result, costs=precomputed_objectives, indices=bitstring_loc, preserve_state=False)
+            overlap = sim.get_overlap(result, costs=precomputed_costs, indices=bitstring_loc, preserve_state=False)
             return 1 - overlap
         elif objective == "expectation and overlap":
-            overlap = sim.get_overlap(result, costs=precomputed_objectives, indices=bitstring_loc, preserve_state=True)
-            expectation = sim.get_expectation(result, costs=precomputed_objectives)
+            overlap = sim.get_overlap(result, costs=precomputed_costs, indices=bitstring_loc, preserve_state=True)
+            expectation = sim.get_expectation(result, costs=precomputed_costs)
             return expectation, 1 - overlap
 
     return f

qokit/qaoa_objective_labs.py

@@ -142,7 +142,7 @@ def get_random_guess_merit_factor(N: int) -> float:
 def get_qaoa_labs_objective(
     N: int,
     p: int,
-    precomputed_merit_factors: np.ndarray | None = None,
+    precomputed_negative_merit_factors: np.ndarray | None = None,
     parameterization: str = "theta",
     objective: str = "expectation",
     precomputed_optimal_bitstrings: np.ndarray | None = None,
@@ -156,7 +156,7 @@ def get_qaoa_labs_objective(
         Number of qubits
     p : int
         Number of QAOA layers (number of parameters will be 2*p)
-    precomputed_merit_factors : np.array
+    precomputed_negative_merit_factors : np.array
         precomputed merit factors to compute the QAOA expectation
     parameterization : str
         If parameterization == 'theta', then f takes one parameter (gamma and beta concatenated)
@@ -186,19 +186,19 @@ def get_qaoa_labs_objective(
 
     terms_ix, offset = get_energy_term_indices(N)
 
-    if precomputed_merit_factors is None:
-        precomputed_merit_factors = get_precomputed_labs_merit_factors(N)
+    if precomputed_negative_merit_factors is None:
+        precomputed_negative_merit_factors = get_precomputed_labs_merit_factors(N)
 
     if objective in ["overlap", "expectation and overlap"] and precomputed_optimal_bitstrings is None:
         precomputed_optimal_bitstrings = get_precomputed_optimal_bitstrings(N)
 
-    precomputed_diagonal_hamiltonian = -(N**2) / (2 * precomputed_merit_factors) - offset
+    precomputed_diagonal_hamiltonian = -(N**2) / (2 * precomputed_negative_merit_factors) - offset
 
     return get_qaoa_objective(
         N=N,
         p=p,
         precomputed_diagonal_hamiltonian=precomputed_diagonal_hamiltonian,
-        precomputed_objectives=precomputed_merit_factors,
+        precomputed_costs=precomputed_negative_merit_factors,
         precomputed_optimal_bitstrings=precomputed_optimal_bitstrings,
         parameterization=parameterization,
         objective=objective,

qokit/qaoa_objective_maxcut.py

@@ -35,7 +35,8 @@ def get_qaoa_maxcut_objective(
     G : nx.Graph
         graph on which MaxCut will be solved
     precomputed_cuts : np.array
-        precomputed cuts to compute the QAOA expectation
+        precomputed cuts to compute the QAOA expectation, for maximization problem
+        send the precomputed cuts/energies as negative
     parameterization : str
         If parameterization == 'theta', then f takes one parameter (gamma and beta concatenated)
         If parameterization == 'gamma beta', then f takes two parameters (gamma and beta)
@@ -68,12 +69,13 @@ def get_qaoa_maxcut_objective(
         parameterized_circuit = get_parameterized_qaoa_circuit(G, p)
     else:
         parameterized_circuit = None
-
+    # Reverse the sign as get_qaoa_objective assumes that the problem is minimization
+    precomputed_costs = precomputed_cuts * -1 if precomputed_cuts is not None else None
     return get_qaoa_objective(
         N=N,
         p=p,
         precomputed_diagonal_hamiltonian=precomputed_cuts,
-        precomputed_objectives=precomputed_cuts,
+        precomputed_costs=precomputed_costs,
         terms=terms,
         precomputed_optimal_bitstrings=precomputed_optimal_bitstrings,
         parameterized_circuit=parameterized_circuit,

qokit/qaoa_objective_portfolio.py

@@ -78,7 +78,7 @@ def get_qaoa_portfolio_objective(
         N=N,
         p=p,
         precomputed_diagonal_hamiltonian=po_problem["scale"] * precomputed_energies,
-        precomputed_objectives=precomputed_energies,
+        precomputed_costs=precomputed_energies,
         precomputed_optimal_bitstrings=precomputed_optimal_bitstrings,
         parameterized_circuit=parameterized_circuit,
         parameterization=parameterization,

tests/test_maxcut.py

@@ -50,12 +50,11 @@ def test_maxcut_qaoa_obj_fixed_angles():
 
         obj = partial(maxcut_obj, w=get_adjacency_matrix(G))
         optimal_cut, x = brute_force(obj, N, function_takes="bits")
-
         for p in range(1, max_p + 1):
             gamma, beta, AR = get_fixed_gamma_beta(d, p, return_AR=True)
             for simulator in ["auto", "qiskit"]:
                 f = get_qaoa_maxcut_objective(N, p, G=G, parameterization="gamma beta", simulator=simulator)
-                assert f(gamma, beta) / optimal_cut > AR
+                assert -f(gamma, beta) / optimal_cut > AR
 
 
 def test_maxcut_weighted_qaoa_obj():
@@ -66,17 +65,18 @@ def test_maxcut_weighted_qaoa_obj():
         lambda row: nx.node_link_graph(row["G_json"]),
         axis=1,
     )
+    # changing the sign of gamma got from file
+    df["gamma"] = df["gamma"].apply(lambda x: [-y for y in x])
 
     for _, row in df.iterrows():
         for simulator in ["auto", "qiskit"]:
             f = get_qaoa_maxcut_objective(row["G"].number_of_nodes(), row["p"], G=row["G"], parameterization="gamma beta", simulator=simulator)
-            assert np.isclose(f(row["gamma"], row["beta"]), row["Expected cut of QAOA"])
+            assert np.isclose(-f(row["gamma"], row["beta"]), row["Expected cut of QAOA"])
 
         # Qiskit non-parameterized circuit must be tested separately
         precomputed_cuts = precompute_energies(maxcut_obj, row["G"].number_of_nodes(), w=get_adjacency_matrix(row["G"]))
         qc = get_qaoa_circuit(row["G"], row["beta"], row["gamma"])
         qc_param = get_parameterized_qaoa_circuit(row["G"], row["p"]).bind_parameters(np.hstack([row["beta"], row["gamma"]]))
-
         sv = np.asarray(qiskit_backend.run(qc).result().get_statevector())
         sv_param = np.asarray(qiskit_backend.run(qc_param).result().get_statevector())
 
@@ -92,6 +92,7 @@ def test_maxcut_precompute(simclass):
     for u, v, w in G.edges(data=True):
         w["weight"] = np.random.rand()
     precomputed_cuts = precompute_energies(maxcut_obj, N, w=get_adjacency_matrix(G))
+    precomputed_cuts = precomputed_cuts * -1
     terms = get_maxcut_terms(G)
     sim = simclass(N, terms=terms)
     cuts = sim.get_cost_diagonal()
@@ -105,14 +106,15 @@ def test_sk_ini_maxcut():
         obj = partial(maxcut_obj, w=get_adjacency_matrix(G))
         optimal_cut, x = brute_force(obj, N, function_takes="bits")
         precomputed_energies = precompute_energies(obj, N)
+        precomputed_energies = precomputed_energies * -1
         last_ar = 0
         for p in range(1, max_p + 1):
             gamma, beta = get_sk_gamma_beta(p)
             for simulator in ["auto"]:
                 f = get_qaoa_maxcut_objective(N, p, G=G, parameterization="gamma beta", simulator=simulator)
                 cur_ar = f(gamma / np.sqrt(d), beta) / optimal_cut
             if p == 1:
-                assert cur_ar > np.mean(precomputed_energies) / optimal_cut
+                assert cur_ar < np.mean(precomputed_energies) / optimal_cut
             else:
-                assert cur_ar > last_ar
+                assert cur_ar < last_ar
                 last_ar = cur_ar

tests/test_qaoa_objective_maxcut.py

@@ -0,0 +1,35 @@
+###############################################################################
+# // SPDX-License-Identifier: Apache-2.0
+# // Copyright : JP Morgan Chase & Co
+###############################################################################
+import networkx as nx
+import numpy as np
+from qokit.qaoa_objective_maxcut import get_qaoa_maxcut_objective
+from qokit.utils import precompute_energies
+from qokit.parameter_utils import get_fixed_gamma_beta
+from qokit.maxcut import maxcut_obj, get_adjacency_matrix
+
+
+def test_validate_energy_for_terms_and_precomputedcuts_are_same():
+    G = nx.random_regular_graph(3, 16)
+    # without precomputed cuts
+    f_terms = get_qaoa_maxcut_objective(G.number_of_nodes(), 1, G=G, parameterization="gamma beta")
+    # with precomputed cuts
+    precomputed_cuts = precompute_energies(maxcut_obj, G.number_of_nodes(), w=get_adjacency_matrix(G))
+    f_precomputedcuts = get_qaoa_maxcut_objective(G.number_of_nodes(), 1, precomputed_cuts=precomputed_cuts, parameterization="gamma beta")
+    p = 1
+    gamma, beta = get_fixed_gamma_beta(3, p)
+    energy_terms = f_terms(-1 * np.asarray(gamma), beta)
+    energy_precomputedcuts = f_precomputedcuts(gamma, beta)
+    assert energy_terms == energy_precomputedcuts
+
+
+def test_validate_energy_for_terms_with_simulators_are_same():
+    G = nx.random_regular_graph(3, 16)
+    f = f = get_qaoa_maxcut_objective(G.number_of_nodes(), 1, G=G, parameterization="gamma beta", simulator="auto")
+    g = get_qaoa_maxcut_objective(G.number_of_nodes(), 1, G=G, parameterization="gamma beta", simulator="qiskit")
+    p = 1
+    gamma, beta = get_fixed_gamma_beta(3, p)
+    auto = f(gamma, beta)
+    qiskit = g(gamma, beta)
+    assert np.isclose(auto, qiskit)