diff --git a/README.md b/README.md index 29ffbc0..a6a43d4 100644 --- a/README.md +++ b/README.md @@ -13,6 +13,7 @@ A flexible, modular Python library for the [Quantum Approximate Optimization Alg - [Installation](#installation) - [Requirements](#requirements) - [Quick Example](#quick-example) +- [Objective direction and sign convention](#objective-direction-and-sign-convention) - [Background](#background) - [Custom Ansatz](#custom-ansatz) - [Running Optimization](#running-optimization-at-depth-p) @@ -81,7 +82,8 @@ qaoa.sample_cost_landscape() qaoa.optimize(depth=3) # Extract results -print("Optimal expectation value:", qaoa.get_Exp(depth=3)) +print("Optimal energy:", qaoa.get_energy(depth=3)) +print("Optimal objective:", qaoa.get_objective(depth=3)) print("Optimal parameters (gamma):", qaoa.get_gamma(depth=3)) print("Optimal parameters (beta):", qaoa.get_beta(depth=3)) ``` @@ -90,14 +92,41 @@ See [examples/](examples/) for more complete worked examples. --- +## Objective direction and sign convention + +Each problem declares an explicit objective direction via `problem.objective_sense` (`"minimize"` or `"maximize"`). + +- `problem.objective_value(x)` is the natural mathematical objective. +- `problem.energy(x)` is the canonical quantity minimized by QAOA. + +Conversion is centralized: + +- `MINIMIZE`: `energy(x) = objective_value(x)` +- `MAXIMIZE`: `energy(x) = -objective_value(x)` + +Phase separators follow: + +$$U_P(\gamma)|x\rangle = e^{-i\gamma\,\mathrm{energy}(x)}|x\rangle.$$ + +Examples: + +- **MaxCut**: natural objective is positive cut value; energy is its negative. +- **QUBO**: natural objective is the un-negated polynomial + $x^\top Q x + c^\top x + b$. + +Legacy `cost(x)` is kept as a deprecated compatibility score with old behavior: +`cost(x) = -energy(x)`. + +--- + ## Background -Given a **cost function** -$$c: \lbrace 0, 1\rbrace^n \rightarrow \mathbb{R}$$ +Given an **energy function** +$$E: \lbrace 0, 1\rbrace^n \rightarrow \mathbb{R}$$ one defines a **problem Hamiltonian** $H_P$ through the action on computational basis states via -$$ H_P |x\rangle = c(x) |x\rangle,$$ +$$ H_P |x\rangle = E(x) |x\rangle,$$ -which means that ground states minimize the cost function $c$. +which means that ground states minimize the canonical energy. Given a parametrized ansatz $| \gamma, \beta \rangle$, a classical optimizer is used to minimize the energy $$ \langle \gamma, \beta | H_P | \gamma, \beta \rangle.$$ @@ -119,7 +148,7 @@ $U_M(\beta_l)=e^{-i\beta_l X^{\otimes n}}$, $U_P(\gamma_l)=e^{-i\gamma_l H_P}$, ## Custom Ansatz -To create a custom QAOA ansatz, specify a [problem](qaoa/problems/base_problem.py), a [mixer](qaoa/mixers/base_mixer.py), and an [initial state](qaoa/initialstates/base_initialstate.py). These base classes each have an abstract method `def create_circuit:` that must be implemented. The problem base class additionally requires `def cost:`. +To create a custom QAOA ansatz, specify a [problem](qaoa/problems/base_problem.py), a [mixer](qaoa/mixers/base_mixer.py), and an [initial state](qaoa/initialstates/base_initialstate.py). These base classes each have an abstract method `def create_circuit:` that must be implemented. A custom problem should define `objective_sense` and `objective_value()`. The phase separator must encode `energy(x)` as above. This library already contains several standard implementations. @@ -177,12 +206,12 @@ qaoa.sample_cost_landscape() Sampling high-dimensional target functions quickly becomes intractable for depth $p>1$. The library therefore **iteratively increases the depth**. At each depth a **local optimization** algorithm (e.g. COBYLA) finds a local minimum, using the following **initial guess**: -- At depth $p=1$: parameters $(\gamma, \beta)$ are taken from the minimum of the sampled cost landscape. +- At depth $p=1$: parameters $(\gamma, \beta)$ are taken from the minimum of the sampled energy landscape. - At depth $p>1$: two strategies are available, controlled by the `interpolate` parameter: * **Interpolation** (`interpolate=True`, default): uses the [INTERP heuristic](https://arxiv.org/pdf/1812.01041.pdf) to produce a smooth initial guess by interpolating the optimal angles from depth $p-1$. Works well for vanilla QAOA. - * **Layer-by-layer grid scan** (`interpolate=False`): the best angles from depth $p-1$ are *locked* and a 2-D grid search is performed over the new layer's parameters. Because the grid includes $(γ=0, β=0)$ — which adds an identity layer reproducing the depth-$(p-1)$ result — the initial cost at depth $p$ is guaranteed to be ≤ cost at depth $p-1$, ensuring a monotonically increasing approximation ratio. Recommended for multi-angle and orbit ansätze. + * **Layer-by-layer grid scan** (`interpolate=False`): the best angles from depth $p-1$ are *locked* and a 2-D grid search is performed over the new layer's parameters. Because the grid includes $(γ=0, β=0)$ — which adds an identity layer reproducing the depth-$(p-1)$ result — the initial energy at depth $p$ is guaranteed to be ≤ energy at depth $p-1$, ensuring a monotonically increasing approximation ratio. Recommended for multi-angle and orbit ansätze. ```python # Interpolation (default) diff --git a/qaoa/problems/__init__.py b/qaoa/problems/__init__.py index 8c5680b..454a2b0 100644 --- a/qaoa/problems/__init__.py +++ b/qaoa/problems/__init__.py @@ -1,4 +1,4 @@ -from .base_problem import Problem +from .base_problem import Problem, ObjectiveSense from .qubo_problem import QUBO from .graph_problem import GraphProblem from .exactcover_problem import ExactCover diff --git a/qaoa/problems/base_problem.py b/qaoa/problems/base_problem.py index dc6d8e0..d6582ba 100644 --- a/qaoa/problems/base_problem.py +++ b/qaoa/problems/base_problem.py @@ -1,8 +1,16 @@ from abc import ABC, abstractmethod +from enum import Enum +import itertools +import warnings from qaoa.utils import validation +class ObjectiveSense(str, Enum): + MINIMIZE = "minimize" + MAXIMIZE = "maximize" + + class BaseProblem(ABC): """ Base class for defining optimization problems. @@ -65,21 +73,59 @@ def create_circuit(self): ``` """ - @abstractmethod - def cost(self, string): - """ - Abstract method to calculate the cost of a solution. - - Subclasses must implement this method to define how the cost of a - solution is calculated for the specific optimization problem. - - Args: - string (str): A solution string or configuration to evaluate. + def __init__(self, objective_sense: ObjectiveSense = ObjectiveSense.MINIMIZE) -> None: + super().__init__() + if not isinstance(objective_sense, ObjectiveSense): + try: + objective_sense = ObjectiveSense(objective_sense) + except Exception as exc: + raise ValueError( + "objective_sense must be one of " + f"{[s.value for s in ObjectiveSense]}" + ) from exc + self.objective_sense = objective_sense + + def _has_objective_value_override(self) -> bool: + return self.__class__.objective_value is not Problem.objective_value + + def _has_legacy_cost_override(self) -> bool: + return self.__class__.cost is not Problem.cost + + def objective_value(self, string): + if self._has_legacy_cost_override(): + warnings.warn( + "Implement objective_value() and objective_sense for custom problems. " + "Legacy cost() fallback is deprecated.", + DeprecationWarning, + stacklevel=2, + ) + return self.cost(string) + raise NotImplementedError("Subclasses must implement objective_value().") + + def energy(self, string): + if not self._has_objective_value_override() and self._has_legacy_cost_override(): + return -self.cost(string) + value = self.objective_value(string) + if self.objective_sense is ObjectiveSense.MINIMIZE: + return value + return -value + + def objective_from_energy(self, energy): + if self.objective_sense is ObjectiveSense.MINIMIZE: + return energy + return -energy + + def score(self, string): + return -self.energy(string) - Returns: - float: The cost of the given solution. - """ - pass + def cost(self, string): + warnings.warn( + "cost() is deprecated and kept for backward compatibility. " + "Use objective_value() or energy(). cost(x) == -energy(x).", + DeprecationWarning, + stacklevel=2, + ) + return self.score(string) @abstractmethod def create_circuit(self): @@ -120,25 +166,47 @@ def get_num_parameters(self): def computeMinMaxCosts(self): """ - Brute force method to compute min and max cost of feasible solution + Deprecated wrapper for objective_bounds(). Kept for backward compatibility. """ - import itertools - - max_cost = float("-inf") - min_cost = float("inf") - for s in ["".join(i) for i in itertools.product("01", repeat=self.N_qubits)]: + warnings.warn( + "computeMinMaxCosts() is deprecated. Use objective_bounds().", + DeprecationWarning, + stacklevel=2, + ) + return self.objective_bounds() + + def objective_bounds(self): + min_objective = float("inf") + max_objective = float("-inf") + for s in map("".join, itertools.product("01", repeat=self.N_qubits)): + if self.isFeasible(s): + value = self.objective_value(s) + min_objective = min(min_objective, value) + max_objective = max(max_objective, value) + return min_objective, max_objective + + def optimal_objective(self): + min_objective, max_objective = self.objective_bounds() + if self.objective_sense is ObjectiveSense.MINIMIZE: + return min_objective + return max_objective + + def energy_bounds(self): + min_energy = float("inf") + max_energy = float("-inf") + for s in map("".join, itertools.product("01", repeat=self.N_qubits)): if self.isFeasible(s): - cost = -self.cost(s) - max_cost = max(max_cost, cost) - min_cost = min(min_cost, cost) - return min_cost, max_cost + value = self.energy(s) + min_energy = min(min_energy, value) + max_energy = max(max_energy, value) + return min_energy, max_energy def validate_circuit(self, t=1, flip=True, atol=1e-8, rtol=1e-8): """ Exact check that the problem's circuit represents the problem's cost function. This tests checks that the unitary operator represented by the quantum circuit is - equal to the expected matrix with diagonal elements - exp(-j*t*cost(e)), + equal to the expected matrix with diagonal elements + exp(-j*t*energy(e)), where e is the corresponding binary state, up to a global phase. Suitable for <= 10 qubits as this check uses the full unitary matrix of size 2^n x 2^n). diff --git a/qaoa/problems/exactcover_problem.py b/qaoa/problems/exactcover_problem.py index e501a64..12af67f 100644 --- a/qaoa/problems/exactcover_problem.py +++ b/qaoa/problems/exactcover_problem.py @@ -86,9 +86,9 @@ def __init__( self.N_qubits = numColumns - def cost(self, string): + def objective_value(self, string): """ - Calculates the cost so that states where an element is not covered, or covered more than once, will be penalized, whereas + Calculates the natural objective so that states where an element is not covered, or covered more than once, will be penalized, whereas sets that contain elements that are covered exactly once are favored. Args: @@ -97,7 +97,7 @@ def cost(self, string): x = np.array(list(map(int, string))) c_e = self.__exactCover(x) - return -(self.weights @ x + self.penalty_factor * c_e) + return self.weights @ x + self.penalty_factor * c_e def isFeasible(self, string): @@ -161,14 +161,14 @@ def bitstrings_all_generator(n, k): bitstrings_generator = bitstrings_all_generator if self.hamming_weight is not None: bitstrings_generator = bitstrings_hamming_weight_generator - opt_val = -np.inf + opt_val = np.inf opt_sol = None num_feasible = 0 for bs in bitstrings_generator(self.N_qubits, self.hamming_weight): - cost = self.cost(bs) - if cost > opt_val: - opt_val = cost + energy = self.energy(bs) + if energy < opt_val: + opt_val = energy opt_sol = bs num_feasible += self.isFeasible(bs) diff --git a/qaoa/problems/graph_problem.py b/qaoa/problems/graph_problem.py index 64fc8f4..55026d6 100644 --- a/qaoa/problems/graph_problem.py +++ b/qaoa/problems/graph_problem.py @@ -2,7 +2,7 @@ from qiskit.circuit import Parameter from abc import abstractmethod -from .base_problem import Problem +from .base_problem import ObjectiveSense, Problem from qaoa.utils import * @@ -43,6 +43,7 @@ def __init__( entangling gates eliminated (proportional to the degree of the fixed node). """ super().__init__() + self.objective_sense = ObjectiveSense.MAXIMIZE # fixes the highest-degree node (node n-1 after relabeling) to "color1" self.fix_one_node = fix_one_node @@ -149,9 +150,9 @@ def slice_string(self, string: str) -> list: labels.append(self.colors["color1"][0]) return labels - def cost(self, string: str) -> float | int: + def objective_value(self, string: str) -> float | int: """ - Compute the cost for a given solution. + Compute the natural objective value for a given solution. Args: string (str): Binary string. @@ -160,7 +161,7 @@ def cost(self, string: str) -> float | int: ValueError: If the length of the string does not match the number of qubits. Returns: - float | int: The cost of the given solution. + float | int: The weighted cut value of the given solution. """ if len(string) != self.N_qubits: raise ValueError( diff --git a/qaoa/problems/maxkcut_one_hot_problem.py b/qaoa/problems/maxkcut_one_hot_problem.py index 91037fb..f17324e 100644 --- a/qaoa/problems/maxkcut_one_hot_problem.py +++ b/qaoa/problems/maxkcut_one_hot_problem.py @@ -2,7 +2,7 @@ from qiskit.circuit import Parameter import networkx as nx -from .base_problem import Problem +from .base_problem import ObjectiveSense, Problem class MaxKCutOneHot(Problem): @@ -32,7 +32,7 @@ def __init__(self, G: nx.Graph, k_cuts: int) -> None: Raises: ValueError: If k_cuts is less than 2 or greater than 8. """ - super().__init__() + super().__init__(objective_sense=ObjectiveSense.MAXIMIZE) if (k_cuts < 2) or (k_cuts > 8): raise ValueError( "k_cuts must be 2 or more, and is not implemented for k_cuts > 8" @@ -68,9 +68,9 @@ def binstringToLabels(self, string: str) -> str: labels += str(idx) return labels - def cost(self, string: str) -> float | int: + def objective_value(self, string: str) -> float | int: """ - Computes the Max k-Cut cost for a given binary string representing a coloring. + Computes the Max k-Cut objective for a given binary string representing a coloring. Args: string (str): The binary string representing the one-hot encoding of node colors. diff --git a/qaoa/problems/portfolio_problem.py b/qaoa/problems/portfolio_problem.py index 9c2fc07..fd8acf6 100644 --- a/qaoa/problems/portfolio_problem.py +++ b/qaoa/problems/portfolio_problem.py @@ -53,22 +53,21 @@ def __init__(self, risk, budget, cov_matrix, exp_return, penalty=0) -> None: super().__init__(Q=Q, c=c, b=b) - def cost(self, string): + def objective_value(self, string): """ - Computes the portfolio cost of a given bitstring. This overrides the QUBO base class - cost to use the problem-specific formula directly. + Computes the portfolio natural objective of a given bitstring. Args: string (str): Bitstring representing the selected assets (portfolio). Returns: - cost (float): The negative of the portfolio objective value (including penalty). + float: The portfolio objective value (including penalty). """ x = np.array(list(map(int, string))) cost = self.risk * (x.T @ self.cov_matrix @ x) - self.exp_return.T @ x cost += self.penalty * (x.sum() - self.budget) ** 2 - return -cost + return cost def isFeasible(self, string): """ @@ -97,12 +96,12 @@ def bitstrings_hamming_weight_generator(n, k): s[i] = '1' yield ''.join(s) - opt_val = -np.inf + opt_val = np.inf opt_sol = None for bs in bitstrings_hamming_weight_generator(self.N_qubits, self.budget): - cost = self.cost(bs) - if cost > opt_val: - opt_val = cost + energy = self.energy(bs) + if energy < opt_val: + opt_val = energy opt_sol = bs return opt_sol diff --git a/qaoa/problems/qubo_problem.py b/qaoa/problems/qubo_problem.py index 5fb4be6..c48ba52 100644 --- a/qaoa/problems/qubo_problem.py +++ b/qaoa/problems/qubo_problem.py @@ -1,10 +1,11 @@ import math import numpy as np +import warnings from qiskit import QuantumCircuit, QuantumRegister from qiskit.circuit import Parameter -from .base_problem import Problem +from .base_problem import ObjectiveSense, Problem import structlog LOG = structlog.get_logger(file=__name__) @@ -32,7 +33,7 @@ class QUBO(Problem): create_circuit(): Creates a parametrized quantum circuit corresponding to the cost function of the QUBO problem. createParameterizedCostCircuitTril(): Creates a parameterized circuit of the triangularized QUBO problem. """ - def __init__(self, Q=None, c=None, b=None) -> None: + def __init__(self, Q=None, c=None, b=None, objective_sense=ObjectiveSense.MINIMIZE) -> None: """ Implements the mapping from the parameters in params to the QUBO problem. Is expected to be called by the child class. @@ -49,7 +50,7 @@ def __init__(self, Q=None, c=None, b=None) -> None: AssertionError: If c is not a 1D numpy ndarray of compatible size. AssertionError: If b is not a scalar. """ - super().__init__() + super().__init__(objective_sense=objective_sense) assert type(Q) is np.ndarray, "Q needs to be a numpy ndarray, but is " + str( type(Q) ) @@ -86,17 +87,22 @@ def __init__(self, Q=None, c=None, b=None) -> None: assert np.isscalar(b), "b is expected to be scalar, but is " + str(b) self.QUBO_b = b - def cost(self, string): + def objective_value(self, string): """ - Computes the cost of a given binary string according to the QUBO formulation. + Computes the natural objective value of a given binary string according to the QUBO formulation. Args: string (str): Binary string representing a candidate solution to the QUBO problem. Returns: - float: The cost of the solution. + float: The natural objective value of the solution. """ - return self.qubo_cost(string) + x = np.array(list(map(int, string))) + return x.T @ self.QUBO_Q @ x + self.QUBO_c.T @ x + self.QUBO_b + + def qubo_objective(self, string): + x = np.array(list(map(int, string))) + return x.T @ self.QUBO_Q @ x + self.QUBO_c.T @ x + self.QUBO_b def qubo_cost(self, string): """ @@ -104,8 +110,12 @@ def qubo_cost(self, string): original cost function is equivalent to the qubo-transformed cost function. This wrapper enables that validation check """ - x = np.array(list(map(int, string))) - return -(x.T @ self.QUBO_Q @ x + self.QUBO_c.T @ x + self.QUBO_b) + warnings.warn( + "qubo_cost() is deprecated. Use objective_value() or energy().", + DeprecationWarning, + stacklevel=2, + ) + return self.objective_value(string) def create_circuit(self): @@ -118,8 +128,9 @@ def create_circuit(self): # To simplify notation: N = self.N_qubits - Q = self.QUBO_Q - c = self.QUBO_c + sign = 1.0 if self.objective_sense is ObjectiveSense.MINIMIZE else -1.0 + Q = sign * self.QUBO_Q + c = sign * self.QUBO_c gamma = Parameter("x_gamma") # Ensure that Q is symmetric and add c to its diagonal @@ -153,13 +164,13 @@ def validate_circuit(self, t=1, flip=True, atol=1e-8, rtol=1e-8): """ Validates two elements: - 1) That the QUBO cost function (self.qubo_cost) is equivalent to the problem-specific - cost function (self.cost) + 1) That the natural QUBO polynomial is equivalent to the problem-specific + natural objective function (self.objective_value) - 2) Exact check that the problem's circuit represents the problem's cost function. + 2) Exact check that the problem's circuit represents the problem's canonical energy. This tests checks that the unitary operator represented by the quantum circuit is equal to the expected matrix with diagonal elements - exp(-j*t*cost(e)), + exp(-j*t*energy(e)), where e is the corresponding binary state, up to a global phase. Suitable for <= 10 qubits as this check uses the full unitary matrix of size 2^n x 2^n). @@ -173,17 +184,17 @@ def validate_circuit(self, t=1, flip=True, atol=1e-8, rtol=1e-8): n = self.N_qubits for i in range(2**n): bitstring = format(i, f'0{n}b') - cost = self.cost(bitstring) - qubo_cost = self.qubo_cost(bitstring) - abs_error = np.abs(cost - qubo_cost) + objective = self.objective_value(bitstring) + qubo_objective = self.qubo_objective(bitstring) + abs_error = np.abs(objective - qubo_objective) if (abs_error > atol): qubo_mapping_errors += 1 max_abs_error = max(max_abs_error, abs_error) if qubo_mapping_errors < 9: mismatches.append({ "bitstring": list(bitstring), - "cost": cost, - "qubo_cost": qubo_cost, + "objective": objective, + "qubo_objective": qubo_objective, "abs_error": abs_error }) @@ -191,12 +202,12 @@ def validate_circuit(self, t=1, flip=True, atol=1e-8, rtol=1e-8): if qubo_mapping_errors > 0: report = { "n_qubits": self.N_qubits, - "max_error": abs_error, + "max_error": max_abs_error, "examples": mismatches } return False, report - # Validate mapping from cost function to quantum circuit: + # Validate mapping from energy function to quantum circuit: circ_ok, circ_report = super().validate_circuit(t=t, flip=flip, atol=atol, rtol=rtol) - circ_report["max_qubo_cost_vs_cost_error"] = max_abs_error + circ_report["max_qubo_objective_error"] = max_abs_error return circ_ok, circ_report \ No newline at end of file diff --git a/qaoa/qaoa.py b/qaoa/qaoa.py index da28bde..bb56762 100644 --- a/qaoa/qaoa.py +++ b/qaoa/qaoa.py @@ -1,4 +1,5 @@ import structlog +import warnings LOG = structlog.get_logger(file=__name__) @@ -54,7 +55,7 @@ class OptResult: get_best_solution(): Returns the best solutions and their corresponding cost. """ - def __init__(self, depth): + def __init__(self, depth, problem): """ Initializes the OptResult object with the given depth. @@ -62,12 +63,19 @@ def __init__(self, depth): depth (int): The depth p of the optimization. """ self.depth = depth + self.problem = problem self.angles = [] self.Exp = [] + self.energy_history = self.Exp + self.objective_history = [] self.Var = [] self.WorstCost = [] self.BestCost = [] + self.best_energy = self.BestCost + self.worst_energy = self.WorstCost + self.best_objective = [] + self.worst_objective = [] self.BestSols = [] self.shots = [] @@ -84,11 +92,17 @@ def add_iteration(self, angles, stat, shots): shots (int): Number of shots taken in the iteration. """ self.angles.append(angles) - self.Exp.append(-stat.get_CVaR()) + cvar_energy = stat.get_CVaR() + self.Exp.append(cvar_energy) + self.objective_history.append(self.problem.objective_from_energy(cvar_energy)) self.Var.append(stat.get_Variance()) - self.BestCost.append(-stat.get_max()) - self.WorstCost.append(-stat.get_min()) - self.BestSols.append(stat.get_max_sols()) + best_energy = stat.get_min() + worst_energy = stat.get_max() + self.BestCost.append(best_energy) + self.WorstCost.append(worst_energy) + self.best_objective.append(self.problem.objective_from_energy(best_energy)) + self.worst_objective.append(self.problem.objective_from_energy(worst_energy)) + self.BestSols.append(stat.get_min_sols()) self.shots.append(shots) def compute_best_index(self): @@ -105,6 +119,12 @@ def get_best_Exp(self): """ return self.Exp[self.index_Exp_min] + def get_best_energy(self): + return self.get_best_Exp() + + def get_best_objective(self): + return self.objective_history[self.index_Exp_min] + def get_best_Var(self): """ Returns: @@ -142,8 +162,8 @@ def get_best_solution(self): - list: The best solutions (bit-strings) that yield the best cost. - float: The best cost found. """ - best_cost = np.min(self.BestCost) - iterations_with_best_cost = np.where(self.BestCost == best_cost)[0] + best_energy = np.min(self.BestCost) + iterations_with_best_cost = np.where(self.BestCost == best_energy)[0] all_best_sols = [] for i in iterations_with_best_cost: @@ -151,7 +171,7 @@ def get_best_solution(self): # flatten the list: all_best_sols = [item for sublist in all_best_sols for item in sublist] best_sols = np.unique(all_best_sols) - return best_sols, best_cost + return best_sols, best_energy class QAOA: @@ -295,6 +315,9 @@ def __init__( self.n_init = 0 self.Exp_sampled_p1 = None + self.Energy_sampled_p1 = None + self.MinEnergy_sampled_p1 = None + self.MaxEnergy_sampled_p1 = None self.landscape_p1_angles = {} self.Var_sampled_p1 = None self.MaxCost_sampled_p1 = None @@ -309,6 +332,7 @@ def __init__( self.post = post self.Exp_post_processed = None + self.Energy_post_processed = None self.Var_post_processed = None self.samplecount_hists = {} self.last_hist = {} @@ -317,10 +341,10 @@ def __init__( def exp_landscape(self): """ Returns: - float: The expected value of the cost landscape at depth p = 1. + float: The expected energy landscape at depth p = 1. """ ### at depth p = 1 - return self.Exp_sampled_p1 + return self.Energy_sampled_p1 def var_landscape(self): """ @@ -332,6 +356,8 @@ def var_landscape(self): def get_Exp(self, depth=None): """ + Deprecated alias for get_energy(). + Args: depth (int, optional): The depth at which to retrieve the expected value. @@ -340,14 +366,32 @@ def get_Exp(self, depth=None): If depth is None, returns a list of best expected values for all depths up to the current depth. If depth is specified, returns the best expected value at that depth. """ + warnings.warn( + "get_Exp() is deprecated. Use get_energy().", + DeprecationWarning, + stacklevel=2, + ) + return self.get_energy(depth=depth) + + def get_energy(self, depth=None): + if not depth: + ret = [] + for i in range(1, self.current_depth + 1): + ret.append(self.optimization_results[i].get_best_energy()) + return ret + if depth > self.current_depth + 1: + raise ValueError + return self.optimization_results[depth].get_best_energy() + + def get_objective(self, depth=None): if not depth: ret = [] for i in range(1, self.current_depth + 1): - ret.append(self.optimization_results[i].get_best_Exp()) + ret.append(self.optimization_results[i].get_best_objective()) return ret if depth > self.current_depth + 1: raise ValueError - return self.optimization_results[depth].get_best_Exp() + return self.optimization_results[depth].get_best_objective() def get_Var(self, depth): """ @@ -565,8 +609,8 @@ def sample_cost_landscape( if self.sequential: expectations = [] variances = [] - maxcosts = [] - mincosts = [] + max_energies = [] + min_energies = [] self.createParameterizedCircuit(depth) logger.info("Executing sample_cost_landscape") @@ -602,30 +646,33 @@ def sample_cost_landscape( self.stat.reset() for string in counts_list: - # qiskit binary strings use little endian encoding, but our cost function expects big endian encoding. Therefore, we reverse the order - cost = self.problem.cost(string[::-1]) + # qiskit binary strings use little endian encoding, but our energy function expects big endian encoding. Therefore, we reverse the order + energy = self.problem.energy(string[::-1]) self.stat.add_sample( - cost, counts_list[string], string[::-1] + energy, counts_list[string], string[::-1] ) expectations.append(self.stat.get_CVaR()) variances.append(self.stat.get_Variance()) - maxcosts.append(self.stat.get_max()) - mincosts.append(self.stat.get_min()) + max_energies.append(self.stat.get_max()) + min_energies.append(self.stat.get_min()) angles = self.landscape_p1_angles - self.Exp_sampled_p1 = -np.array(expectations).reshape( + self.Energy_sampled_p1 = np.array(expectations).reshape( angles["beta"][2], angles["gamma"][2] ) + self.Exp_sampled_p1 = self.Energy_sampled_p1 self.Var_sampled_p1 = np.array(variances).reshape( angles["beta"][2], angles["gamma"][2] ) - self.MaxCost_sampled_p1 = -np.array(maxcosts).reshape( + self.MaxEnergy_sampled_p1 = np.array(max_energies).reshape( angles["beta"][2], angles["gamma"][2] ) - self.MinCost_sampled_p1 = -np.array(mincosts).reshape( + self.MinEnergy_sampled_p1 = np.array(min_energies).reshape( angles["beta"][2], angles["gamma"][2] ) + self.MaxCost_sampled_p1 = self.MinEnergy_sampled_p1 + self.MinCost_sampled_p1 = self.MaxEnergy_sampled_p1 logger.info("Done measurement") else: self.createParameterizedCircuit(depth) @@ -697,7 +744,7 @@ def measurementStatistics(self, job): if self.memorysize > 0: for measurement in memory_list: self.memory_lists.append( - [measurement, self.problem.cost(measurement[::-1])] + [measurement, self.problem.energy(measurement[::-1])] ) self.memorysize -= 1 if self.memorysize < 1: @@ -708,36 +755,39 @@ def measurementStatistics(self, job): if isinstance(counts_list, list): expectations = [] variances = [] - maxcosts = [] - mincosts = [] + max_energies = [] + min_energies = [] for i, counts in enumerate(counts_list): self.stat.reset() for string in counts: - # qiskit binary strings use little endian encoding, but our cost function expects big endian encoding. Therefore, we reverse the order - cost = self.problem.cost(string[::-1]) - self.stat.add_sample(cost, counts[string], string[::-1]) + # qiskit binary strings use little endian encoding, but our energy function expects big endian encoding. Therefore, we reverse the order + energy = self.problem.energy(string[::-1]) + self.stat.add_sample(energy, counts[string], string[::-1]) expectations.append(self.stat.get_CVaR()) variances.append(self.stat.get_Variance()) - maxcosts.append(self.stat.get_max()) - mincosts.append(self.stat.get_min()) + max_energies.append(self.stat.get_max()) + min_energies.append(self.stat.get_min()) angles = self.landscape_p1_angles - self.Exp_sampled_p1 = -np.array(expectations).reshape( + self.Energy_sampled_p1 = np.array(expectations).reshape( angles["beta"][2], angles["gamma"][2] ) + self.Exp_sampled_p1 = self.Energy_sampled_p1 self.Var_sampled_p1 = np.array(variances).reshape( angles["beta"][2], angles["gamma"][2] ) - self.MaxCost_sampled_p1 = -np.array(maxcosts).reshape( + self.MaxEnergy_sampled_p1 = np.array(max_energies).reshape( angles["beta"][2], angles["gamma"][2] ) - self.MinCost_sampled_p1 = -np.array(mincosts).reshape( + self.MinEnergy_sampled_p1 = np.array(min_energies).reshape( angles["beta"][2], angles["gamma"][2] ) + self.MaxCost_sampled_p1 = self.MinEnergy_sampled_p1 + self.MinCost_sampled_p1 = self.MaxEnergy_sampled_p1 else: for string in counts_list: - # qiskit binary strings use little endian encoding, but our cost function expects big endian encoding. Therefore, we reverse the order - cost = self.problem.cost(string[::-1]) - self.stat.add_sample(cost, counts_list[string], string[::-1]) + # qiskit binary strings use little endian encoding, but our energy function expects big endian encoding. Therefore, we reverse the order + energy = self.problem.energy(string[::-1]) + self.stat.add_sample(energy, counts_list[string], string[::-1]) def optimize( self, @@ -819,7 +869,7 @@ def optimize( angles0 = self._grid_search_layer(best_angles, angles) self.optimization_results[self.current_depth + 1] = OptResult( - self.current_depth + 1 + self.current_depth + 1, self.problem ) # Create parameterized circuit at new depth new_depth = int((len(angles0) - n_init) / n_per_layer) @@ -833,7 +883,7 @@ def optimize( self.optimization_results[self.current_depth + 1].opt_time = time.perf_counter() - start_time LOG.info( - f"cost(depth { self.current_depth + 1} = {res.fun}", + f"energy(depth { self.current_depth + 1}) = {res.fun}", func=self.optimize.__name__, ) @@ -852,7 +902,8 @@ def optimize( if self.post: samples = self.samplecount_hists[self.current_depth] post_processing(self, samples=samples, K=self.post) - self.Exp_post_processed = -self.stat.get_CVaR() + self.Energy_post_processed = self.stat.get_CVaR() + self.Exp_post_processed = self.Energy_post_processed self.Var_post_processed = self.stat.get_Variance() def local_opt(self, angles0): @@ -890,7 +941,7 @@ def loss(self, angles): angles (list): List of angles (gamma and beta). Returns: - float: The negative expected value of the cost function (CVaR) for the given angles. + float: The expected energy CVaR for the given angles. This is used as the objective function to be minimized during optimization. Raises: @@ -932,7 +983,7 @@ def loss(self, angles): angles.copy(), self.stat, shots_taken ) - return -self.stat.get_CVaR() + return self.stat.get_CVaR() def getParametersToBind(self, angles, depth, asList=False): """ @@ -1028,7 +1079,7 @@ def interp(self, angles): def _eval_cost(self, angle_array): """ - Evaluate the expected cost (CVaR) for a specific angle array without + Evaluate the expected energy (CVaR) for a specific angle array without recording the result in ``optimization_results``. Intended for use during grid searches where many candidate points are @@ -1039,7 +1090,7 @@ def _eval_cost(self, angle_array): be consistent with the current ``parametrized_circuit_depth``. Returns: - float: Negative expected cost (CVaR), i.e. the value to minimise. + float: Expected energy (CVaR), i.e. the value to minimise. Raises: NotImplementedError: If the backend is not local. @@ -1061,9 +1112,9 @@ def _eval_cost(self, angle_array): counts = jres.get_counts() self.stat.reset() for string in counts: - cost = self.problem.cost(string[::-1]) - self.stat.add_sample(cost, counts[string], string[::-1]) - return -self.stat.get_CVaR() + energy = self.problem.energy(string[::-1]) + self.stat.add_sample(energy, counts[string], string[::-1]) + return self.stat.get_CVaR() def _grid_search_layer(self, prev_angles, angles): """ @@ -1128,7 +1179,7 @@ def _grid_search_layer(self, prev_angles, angles): best_cost = cost best_angles = candidate.copy() - logger.info(f"Layer grid search done, best cost: {-best_cost:.6f}") + logger.info(f"Layer grid search done, best energy: {best_cost:.6f}") return best_angles def hist(self, angles, shots): diff --git a/qaoa/utils/flip.py b/qaoa/utils/flip.py index dd95fbd..9827f88 100644 --- a/qaoa/utils/flip.py +++ b/qaoa/utils/flip.py @@ -4,7 +4,7 @@ class BitFlip: """ - BitFlip class for performing random bit flips on a string to increase cost. + BitFlip class for performing random bit flips on a string to reduce energy. Attributes: circuit (QuantumCircuit): Quantum circuit for bit flips. @@ -24,7 +24,7 @@ def __init__(self, n): def boost_samples(self, problem, string, K=5): """ - Random bitflips on string/list of strings to increase cost. + Random bitflips on string/list of strings to reduce energy. Args: problem: BaseType Problem. @@ -36,7 +36,7 @@ def boost_samples(self, problem, string, K=5): """ string_arr = np.array([int(bit) for bit in string]) old_string = string - cost = problem.cost(string[::-1]) + energy = problem.energy(string[::-1]) for _ in range(K): shuffled_indices = np.arange(self.N_qubits) @@ -46,10 +46,10 @@ def boost_samples(self, problem, string, K=5): string_arr_altered = np.copy(string_arr) string_arr_altered[i] = not (string_arr[i]) string_altered = "".join(map(str, string_arr_altered)) - new_cost = problem.cost(string_altered[::-1]) + new_energy = problem.energy(string_altered[::-1]) - if new_cost > cost: - cost = new_cost + if new_energy < energy: + energy = new_energy string_arr = string_arr_altered string = string_altered diff --git a/qaoa/utils/plotroutines.py b/qaoa/utils/plotroutines.py index cd55054..f6382f1 100644 --- a/qaoa/utils/plotroutines.py +++ b/qaoa/utils/plotroutines.py @@ -120,7 +120,7 @@ def plot_ApproximationRatio( tuple: ``(fig, ax)``. """ if not shots: - exp = np.array(qaoa_instance.get_Exp()) + exp = np.array(qaoa_instance.get_objective()) else: exp = [] for p in range(1, qaoa_instance.current_depth + 1): @@ -130,13 +130,16 @@ def plot_ApproximationRatio( fig, ax = _get_fig_ax(fig) ax.hlines(1, 1, maxdepth, linestyles="solid", colors="black") - # Normalized approximation ratio for a minimization objective. - # Here mincost is the optimal (most negative) value and maxcost the worst. - # This maps exp = maxcost → 0 (worst) and exp = mincost → 1 (optimal). - # Hence we use (maxcost - exp) / (maxcost - mincost). + if np.isclose(maxcost, mincost): + appr_ratio = np.ones_like(exp) + elif qaoa_instance.problem.objective_sense.value == "maximize": + appr_ratio = (exp - mincost) / (maxcost - mincost) + else: + appr_ratio = (maxcost - exp) / (maxcost - mincost) + ax.plot( np.arange(1, maxdepth + 1), - (maxcost - exp) / (maxcost - mincost), + appr_ratio, style, label=label, ) @@ -197,11 +200,11 @@ def _apprrat_successprob(qaoa_instance, depth, shots=10**4): for string in hist: if qaoa_instance.problem.isFeasible(string): - cost = qaoa_instance.problem.cost(string) + cost = qaoa_instance.problem.objective_value(string) counts += hist[string] stat.add_sample(cost, hist[string], string) - return -stat.get_CVaR(), counts / shots + return stat.get_CVaR(), counts / shots # Keep the old private name as an alias for internal backward compatibility. __apprrat_successprob = _apprrat_successprob @@ -441,7 +444,7 @@ def printBestHistogramEntries(qaoa, classical_solution=None, num_solutions=10, s best_classical_sol = None if classical_solution is not None: best_classical_sol = _np2str(classical_solution) - print("Classical best result: ", (best_classical_sol, qaoa.problem.cost(best_classical_sol))) + print("Classical best result: ", (best_classical_sol, qaoa.problem.objective_value(best_classical_sol))) print(" --> points to the classical solution ") print(" * marks feasible solutions ") for p in range(1, qaoa.current_depth + 1): @@ -454,13 +457,13 @@ def printBestHistogramEntries(qaoa, classical_solution=None, num_solutions=10, s best_classical_sol_i = None print("Results for depth " + str(p) + " using best angles:") for s, freq in sorted_hist.items(): - cost = qaoa.problem.cost(s) + cost = qaoa.problem.objective_value(s) if i == 1: best_sol = s best_cost = cost best_freq = freq best_i = i - elif cost > best_cost: + elif qaoa.problem.energy(s) < qaoa.problem.energy(best_sol): best_sol = s best_cost = cost best_freq = freq @@ -562,4 +565,3 @@ def prob_hit_ones(p, n): ax.grid(True, which="both", ls="--") return fig, ax - diff --git a/qaoa/utils/post.py b/qaoa/utils/post.py index ba2c6b0..52a76ea 100644 --- a/qaoa/utils/post.py +++ b/qaoa/utils/post.py @@ -34,7 +34,7 @@ def post_processing(instance, samples, K=5): count = 1 instance.stat.add_sample( - instance.problem.cost(boosted[::-1]), count, boosted[::-1] + instance.problem.energy(boosted[::-1]), count, boosted[::-1] ) hist_post[boosted] = hist_post.get(boosted, 0) + count return hist_post @@ -65,7 +65,7 @@ def post_process_all_depths(instance, K=5): samples=hist, K=K, ) - exp_in_layers[d] = exp_in_layers.get(d, []) + [-instance.stat.get_CVaR()] + exp_in_layers[d] = exp_in_layers.get(d, []) + [instance.stat.get_CVaR()] exp.append(stat.mean(exp_in_layers[d])) var.append(stat.variance(exp_in_layers[d])) return (np.array(exp), np.array(var)) diff --git a/qaoa/utils/qaoaIO.py b/qaoa/utils/qaoaIO.py index dc0d366..cc9dc7e 100644 --- a/qaoa/utils/qaoaIO.py +++ b/qaoa/utils/qaoaIO.py @@ -83,15 +83,17 @@ class ExactCoverProblemData(ProblemData): weights: np.ndarray = None solution: np.ndarray = None hamming_weight: int = None + objective_sense: str = "minimize" problem_type: str = "ExactCover" @dataclass class PortfolioOptimizationProblemData(ProblemData): risk: float = 0.0 - exp_returns: np.ndarray = None + exp_return: np.ndarray = None cov_matrix: np.ndarray = None budget: int = 0 + objective_sense: str = "minimize" problem_type: str = "PortfolioOptimization" @@ -121,6 +123,8 @@ class DepthResult: optimal_angles: List[float] histogram: Dict[str, int] opt_time: float # runtime in seconds + best_energy: float | None = None + best_objective: float | None = None @dataclass @@ -139,6 +143,7 @@ class QAOAResult: problem: ProblemData qaoa_params: QAOAParameters metadata: Dict[str, str] = field(default_factory=dict) + schema_version: int = 2 def __post_init__(self): """Automatically populate metadata if not provided.""" @@ -150,6 +155,7 @@ def __post_init__(self): def save(self, filename: str): """Save result (including problem type info) to JSON file.""" data = { + "schema_version": self.schema_version, "problem": self.problem.to_dict(), "qaoa_params": _numpy_to_list(asdict(self.qaoa_params)), "metadata": self.metadata @@ -166,7 +172,7 @@ def load(cls, filename: str) -> "QAOAResult": # Rebuild problem instance from its dict problem = ProblemData.from_dict(data["problem"]) - #depths = [DepthResult(**d) for d in data["qaoa_params"]["depths"]] + schema_version = data.get("schema_version", 1) depths_data = data["qaoa_params"]["depths"] depths = {int(k): DepthResult(**v) for k, v in depths_data.items()} @@ -180,7 +186,12 @@ def load(cls, filename: str) -> "QAOAResult": depths=depths, ) - return cls(problem=problem, qaoa_params=qaoa_params, metadata=data.get("metadata", {})) + return cls( + problem=problem, + qaoa_params=qaoa_params, + metadata=data.get("metadata", {}), + schema_version=schema_version, + ) def _generate_metadata(self) -> dict: @@ -224,14 +235,17 @@ def from_qaoa(cls, qaoa: QAOA, depths[k] = DepthResult( optimal_angles = qaoa.optimization_results[k].get_best_angles(), histogram = qaoa.hist(qaoa.optimization_results[k].get_best_angles(), hist_shots), - opt_time = qaoa.optimization_results[k].opt_time + opt_time = qaoa.optimization_results[k].opt_time, + best_energy = qaoa.get_energy(k), + best_objective = qaoa.get_objective(k), ) problem_data = ExactCoverProblemData( columns = qaoa.problem.columns, weights = qaoa.problem.weights, solution = solution, - hamming_weight = qaoa.problem.hamming_weight + hamming_weight = qaoa.problem.hamming_weight, + objective_sense = qaoa.problem.objective_sense.value, ) init_method = InitMethod(str(qaoa.initialstate).split(" ")[0].split(".")[-1].upper()) @@ -258,7 +272,7 @@ def from_qaoa(cls, qaoa: QAOA, depths = depths ) - return cls(problem=problem_data, qaoa_params=qaoa_params) + return cls(problem=problem_data, qaoa_params=qaoa_params, schema_version=2) # TODO: Implement # def generate_qaoa_object(self) -> "QAOA" diff --git a/qaoa/utils/statistic.py b/qaoa/utils/statistic.py index 88b4f75..0a09193 100644 --- a/qaoa/utils/statistic.py +++ b/qaoa/utils/statistic.py @@ -135,8 +135,8 @@ def get_CVaR(self): float: The CVaR based on the samples. """ if self.cvar < 1: - cvarK = int(np.round(self.cvar * len(self.all_values))) - cvar = np.sum(self.all_values[-cvarK:]) / cvarK + cvarK = max(1, int(np.round(self.cvar * len(self.all_values)))) + cvar = np.sum(self.all_values[:cvarK]) / cvarK return cvar else: return self.get_E() diff --git a/qaoa/utils/validation.py b/qaoa/utils/validation.py index 37a1a4d..f8ac5b1 100644 --- a/qaoa/utils/validation.py +++ b/qaoa/utils/validation.py @@ -15,10 +15,10 @@ def check_phase_separator_exact_qaoa(qaoa, *arg, **kwarg): def check_phase_separator_exact_problem(problem, t=1, flip=True, atol=1e-8, rtol=1e-8): """ - Exact check that the problem's circuit represents the problem's cost function. + Exact check that the problem's circuit represents the problem's energy function. This tests checks that the unitary operator represented by the quantum circuit is equal to the expected matrix with diagonal elements - exp(-j*t*cost(e)), + exp(-j*t*energy(e)), where e is the corresponding binary state, up to a global phase. Suitable for <= 10 qubits as this check uses the full unitary matrix of size 2^n x 2^n). @@ -30,17 +30,17 @@ def check_phase_separator_exact_problem(problem, t=1, flip=True, atol=1e-8, rtol {problem.circuit.parameters[0]: t}, inplace = False ) - cost_fn = problem.cost + energy_fn = problem.energy U = Operator(circ).data # complex ndarray n = circ.num_qubits d = 2**n # Compare diagonal phases to expected, modulo a global phase # expected diag entries - costs = [] + energies = [] for i in range(d): - costs.append(cost_fn(_bitstring(i, n, flip=flip))) - expected = np.exp(1j * t * np.asarray(costs, dtype=float)) + energies.append(energy_fn(_bitstring(i, n, flip=flip))) + expected = np.exp(-1j * t * np.asarray(energies, dtype=float)) diag = np.diag(U) diff --git a/unittests/test_objective_sense_migration.py b/unittests/test_objective_sense_migration.py new file mode 100644 index 0000000..4f0bcba --- /dev/null +++ b/unittests/test_objective_sense_migration.py @@ -0,0 +1,209 @@ +import tempfile +import unittest + +import networkx as nx +import numpy as np +from qiskit import QuantumCircuit, QuantumRegister +from qiskit.circuit import Parameter + +from qaoa import QAOA, initialstates, mixers, problems +from qaoa.problems.base_problem import ObjectiveSense, Problem +from qaoa.utils import BitFlip, Statistic, qaoaIO + + +class MinToyProblem(Problem): + def __init__(self): + super().__init__(objective_sense=ObjectiveSense.MINIMIZE) + self.N_qubits = 2 + + def objective_value(self, string): + return float(sum(int(b) for b in string)) + + def create_circuit(self): + q = QuantumRegister(2) + self.circuit = QuantumCircuit(q) + gamma = Parameter("x_gamma") + self.circuit.p(-gamma, q[0]) + self.circuit.p(-gamma, q[1]) + + +class MaxToyProblem(Problem): + def __init__(self): + super().__init__(objective_sense=ObjectiveSense.MAXIMIZE) + self.N_qubits = 2 + + def objective_value(self, string): + return float(sum(int(b) for b in string)) + + def create_circuit(self): + q = QuantumRegister(2) + self.circuit = QuantumCircuit(q) + gamma = Parameter("x_gamma") + self.circuit.p(gamma, q[0]) + self.circuit.p(gamma, q[1]) + + +class WrongSignMinProblem(MinToyProblem): + def create_circuit(self): + q = QuantumRegister(2) + self.circuit = QuantumCircuit(q) + gamma = Parameter("x_gamma") + self.circuit.p(gamma, q[0]) + self.circuit.p(gamma, q[1]) + + +class TestObjectiveSenseMigration(unittest.TestCase): + def test_objective_sense_enum(self): + self.assertEqual(ObjectiveSense.MINIMIZE.value, "minimize") + self.assertEqual(ObjectiveSense.MAXIMIZE.value, "maximize") + with self.assertRaises(ValueError): + MinToyProblem().objective_sense = ObjectiveSense("invalid") + + def test_energy_objective_invariants(self): + pmin = MinToyProblem() + pmax = MaxToyProblem() + self.assertEqual(pmin.energy("10"), pmin.objective_value("10")) + self.assertEqual(pmax.energy("10"), -pmax.objective_value("10")) + self.assertEqual(pmin.cost("10"), -pmin.energy("10")) + self.assertEqual(pmax.cost("10"), -pmax.energy("10")) + + def test_maxcut_objective_and_energy(self): + G = nx.Graph() + G.add_edge(0, 1, weight=2.0) + problem = problems.MaxCut(G) + self.assertEqual(problem.objective_value("01"), 2.0) + self.assertEqual(problem.energy("01"), -2.0) + + def test_qubo_default_minimize_and_objective(self): + Q = np.array([[1.0, 0.0], [0.0, 2.0]]) + problem = problems.QUBO(Q) + self.assertEqual(problem.objective_sense, ObjectiveSense.MINIMIZE) + self.assertEqual(problem.objective_value("11"), 3.0) + self.assertEqual(problem.energy("11"), 3.0) + + def test_exact_cover_objective(self): + columns = np.array([[1, 0], [0, 1]]) + weights = np.array([2.0, 3.0]) + problem = problems.ExactCover(columns, weights=weights, penalty_factor=5.0) + self.assertEqual(problem.objective_value("11"), 5.0) + + def test_portfolio_objective(self): + cov = np.array([[1.0, 0.2], [0.2, 1.0]]) + exp_ret = np.array([0.1, 0.3]) + problem = problems.PortfolioOptimization( + risk=0.5, budget=1, cov_matrix=cov, exp_return=exp_ret, penalty=2.0 + ) + # x = 01 => 0.5*1 - 0.3 + 0 penalty = 0.2 + self.assertAlmostEqual(problem.objective_value("01"), 0.2) + + def test_lower_tail_cvar_and_maximize_mapping(self): + stat = Statistic(cvar=0.5) + for v in [1.0, 2.0, 3.0, 4.0]: + stat.add_sample(v, 1.0, str(v)) + self.assertAlmostEqual(stat.get_CVaR(), 1.5) + + max_problem = MaxToyProblem() + stat2 = Statistic(cvar=0.5) + for v in [-1.0, -2.0, -3.0, -4.0]: + stat2.add_sample(v, 1.0, str(v)) + energy_cvar = stat2.get_CVaR() + self.assertAlmostEqual(energy_cvar, -3.5) + self.assertAlmostEqual(max_problem.objective_from_energy(energy_cvar), 3.5) + + def test_deterministic_selection_by_energy(self): + pmin = MinToyProblem() + pmax = MaxToyProblem() + self.assertLess(pmin.energy("00"), pmin.energy("11")) + self.assertLess(pmax.energy("11"), pmax.energy("00")) + self.assertLess(pmin.objective_value("00"), pmin.objective_value("11")) + self.assertGreater(pmax.objective_value("11"), pmax.objective_value("00")) + + def test_phase_validation_pass_and_fail(self): + ok_min, _ = MinToyProblem().validate_circuit() + ok_max, _ = MaxToyProblem().validate_circuit() + ok_bad, _ = WrongSignMinProblem().validate_circuit() + self.assertTrue(ok_min) + self.assertTrue(ok_max) + self.assertFalse(ok_bad) + + def test_flip_boosting_both_senses(self): + flipper = BitFlip(2) + np.random.seed(0) + pmin = MinToyProblem() + s0 = "11" + s1 = flipper.boost_samples(problem=pmin, string=s0, K=10) + self.assertLessEqual(pmin.energy(s1[::-1]), pmin.energy(s0[::-1])) + + np.random.seed(0) + pmax = MaxToyProblem() + s0 = "00" + s1 = flipper.boost_samples(problem=pmax, string=s0, K=10) + self.assertLessEqual(pmax.energy(s1[::-1]), pmax.energy(s0[::-1])) + + def test_objective_and_energy_bounds(self): + pmin = MinToyProblem() + pmax = MaxToyProblem() + self.assertEqual(pmin.objective_bounds(), (0.0, 2.0)) + self.assertEqual(pmax.objective_bounds(), (0.0, 2.0)) + self.assertEqual(pmin.energy_bounds(), (0.0, 2.0)) + self.assertEqual(pmax.energy_bounds(), (-2.0, -0.0)) + self.assertEqual(pmin.optimal_objective(), 0.0) + self.assertEqual(pmax.optimal_objective(), 2.0) + + def test_qaoa_get_exp_alias_and_energy(self): + from qiskit_aer import AerSimulator + + G = nx.path_graph(3) + for u, v in G.edges(): + G[u][v]["weight"] = 1.0 + q = QAOA( + problems.MaxCut(G), + mixers.X(), + initialstates.Plus(), + backend=AerSimulator(), + shots=128, + ) + q.optimize(depth=1, angles={"gamma": [0, np.pi, 3], "beta": [0, np.pi, 3]}) + self.assertAlmostEqual(q.get_Exp(depth=1), q.get_energy(depth=1)) + self.assertAlmostEqual( + q.get_objective(depth=1), -q.get_energy(depth=1) + ) + + def test_serialization_roundtrip_preserves_objective_sense(self): + pd = qaoaIO.ExactCoverProblemData( + columns=np.array([[1, 0], [0, 1]]), + weights=np.array([1.0, 2.0]), + solution=np.array([1, 1]), + hamming_weight=1, + objective_sense="minimize", + ) + params = qaoaIO.QAOAParameters( + cvar=0.5, + init_method=qaoaIO.InitMethod.PLUS, + mixer_method=qaoaIO.MixerMethod.X, + backend="sim", + optimizer="COBYLA", + N_qubits=2, + depths={1: qaoaIO.DepthResult([0.1, 0.2], {"00": 1}, 0.01, 0.0, 0.0)}, + ) + result = qaoaIO.QAOAResult(problem=pd, qaoa_params=params) + with tempfile.NamedTemporaryFile(suffix=".json") as fp: + result.save(fp.name) + loaded = qaoaIO.QAOAResult.load(fp.name) + self.assertEqual(loaded.schema_version, 2) + self.assertEqual(loaded.problem.objective_sense, "minimize") + + def test_approximation_ratio_mapping_formula(self): + best, worst = 10.0, 2.0 + value_best, value_worst = 10.0, 2.0 + ratio_max_best = (value_best - worst) / (best - worst) + ratio_max_worst = (value_worst - worst) / (best - worst) + self.assertAlmostEqual(ratio_max_best, 1.0) + self.assertAlmostEqual(ratio_max_worst, 0.0) + + best, worst = 2.0, 10.0 + value_best, value_worst = 2.0, 10.0 + ratio_min_best = (worst - value_best) / (worst - best) + ratio_min_worst = (worst - value_worst) / (worst - best) + self.assertAlmostEqual(ratio_min_best, 1.0) + self.assertAlmostEqual(ratio_min_worst, 0.0) diff --git a/unittests/test_statistic_utility.py b/unittests/test_statistic_utility.py index 273b3ea..2298d96 100644 --- a/unittests/test_statistic_utility.py +++ b/unittests/test_statistic_utility.py @@ -63,13 +63,13 @@ def test_variance(self): self.assertAlmostEqual(stat.get_Variance(), 2.0) def test_cvar_below_one(self): - """CVaR with alpha=0.5 should average over top 50% of values.""" + """CVaR with alpha=0.5 should average over lower 50% of values.""" stat = self.Statistic(cvar=0.5) # add 4 samples: 1, 2, 3, 4 for v in [1.0, 2.0, 3.0, 4.0]: stat.add_sample(v, 1.0, str(v)) - # Top 50% = [3, 4], CVaR = 3.5 - self.assertAlmostEqual(stat.get_CVaR(), 3.5) + # Lower 50% = [1, 2], CVaR = 1.5 + self.assertAlmostEqual(stat.get_CVaR(), 1.5) def test_cvar_equal_one_is_expectation(self): stat = self.Statistic(cvar=1)