diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 1b2d69b5..d374792b 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -23,8 +23,10 @@ repos: hooks: - id: mypy exclude: examples|docs + args: [--always-false=USE_DOUBLE_PRECISION] additional_dependencies: - "qoolqit[extras]>=1.1" # The version in pyproject.toml must match + - "jaxtyping" # Cleanup jupyter notebooks - repo: https://github.com/srstevenson/nb-clean diff --git a/docs/content/classical/cplex_solver/cplex_solving.md b/docs/content/classical/cplex_solver/cplex_solving.md index 61a3b155..ec48ceab 100644 --- a/docs/content/classical/cplex_solver/cplex_solving.md +++ b/docs/content/classical/cplex_solver/cplex_solving.md @@ -11,13 +11,13 @@ class CplexSolver(BaseClassicalSolver): def __init__( self, instance: QUBOInstance, - config: Optional[Dict[str, Any]] = None + config: Optional[dict[str, Any]] = None ) ``` - **Parameters**: - `instance` (`QUBOInstance`): QUBO problem containing a square coefficient matrix (`torch.Tensor`). - - `config` (`Optional[Dict[str, Any]]`): Dictionary supporting: + - `config` (`Optional[dict[str, Any]]`): Dictionary supporting: - `cplex_maxtime` (`float`, default `600.0`): Maximum solve time in seconds. - `cplex_log_path` (`str`, default `"solver.log"`): Path for CPLEX log output. diff --git a/docs/tutorial/greedy_demo.gif b/docs/tutorial/greedy_demo.gif index 957cf75e..8d58ad85 100644 Binary files a/docs/tutorial/greedy_demo.gif and b/docs/tutorial/greedy_demo.gif differ diff --git a/pyproject.toml b/pyproject.toml index bf011ba3..a33675bd 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -64,6 +64,8 @@ dev = [ "jupyter", "tqdm", "mypy", + "jaxtyping", + "beartype", ] doc = [ "mkdocs", diff --git a/pyrightconfig.json b/pyrightconfig.json new file mode 100644 index 00000000..2eae3371 --- /dev/null +++ b/pyrightconfig.json @@ -0,0 +1,5 @@ +{ + "defineConstant": { + "USE_DOUBLE_PRECISION": false + } +} diff --git a/qubosolver/__init__.py b/qubosolver/__init__.py index f6eab914..1e627116 100644 --- a/qubosolver/__init__.py +++ b/qubosolver/__init__.py @@ -1,16 +1,65 @@ from __future__ import annotations -from importlib import import_module +from qubosolver.types._checks import _RUNTIME_CHECKS -from .data import * # noqa: F403 -from .qubo_instance import * # noqa: F403 -from .qubo_types import * # noqa: F403 -from .config import * # noqa: F403 -from .utils import * # noqa: F403 +from qubosolver.types import ( + bitstring, + bitstrings, + matrix, + tensor, + vector, + vectori, + Bitstring, + Bitstrings, + Matrix, + Tensor, + Vector, + Vectori, + QUBOSolution, + QUBOAnalyzer, + QUBOInstance, + EmbedderType, + DriveType, + LayoutType, + SolutionStatusType, + DensityType, +) +from qubosolver.config import SolverConfig, EmbeddingConfig, DriveShapingConfig -list_of_submodules = [".data", ".utils", ".qubo_instance", ".config"] +__all__ = [ + # Submodules + "bitstring", + "bitstrings", + "matrix", + "tensor", + "vector", + "vectori", + # Type Aliases + "Bitstring", + "Bitstrings", + "Matrix", + "Tensor", + "Vector", + "Vectori", + # Classes + "QUBOSolution", + "QUBOAnalyzer", + "QUBOInstance", + # Enums + "EmbedderType", + "DriveType", + "LayoutType", + "SolutionStatusType", + "DensityType", + # Configs + "SolverConfig", + "EmbeddingConfig", + "DriveShapingConfig", + # Functions +] -__all__ = [] -for submodule in list_of_submodules: - __all_submodule__ = getattr(import_module(submodule, package="qubosolver"), "__all__") - __all__ += __all_submodule__ +if _RUNTIME_CHECKS: + from beartype import BeartypeConf + from beartype.claw import beartype_this_package + + beartype_this_package(conf=BeartypeConf(warning_cls_on_decorator_exception=None)) diff --git a/qubosolver/_drive_shaping/__init__.py b/qubosolver/_drive_shaping/__init__.py new file mode 100644 index 00000000..5c1a0a28 --- /dev/null +++ b/qubosolver/_drive_shaping/__init__.py @@ -0,0 +1,4 @@ +from __future__ import annotations + +from qubosolver.drive_shaping import heuristic as heuristic +from qubosolver.drive_shaping import optimized as optimized diff --git a/qubosolver/_drive_shaping/device_specs.py b/qubosolver/_drive_shaping/device_specs.py new file mode 100644 index 00000000..0c852fd9 --- /dev/null +++ b/qubosolver/_drive_shaping/device_specs.py @@ -0,0 +1,103 @@ +from __future__ import annotations + +import math +import warnings + +import qoolqit + + +def compare_specs(a: dict, b: dict) -> None: + keys_a = set(a.keys()) + keys_b = set(b.keys()) + + if keys_a != keys_b: + only_in_a = keys_a - keys_b + only_in_b = keys_b - keys_a + if only_in_a: + warnings.warn(f"Keys present in a but missing in b: {only_in_a}") + if only_in_b: + warnings.warn(f"Keys present in b but missing in a: {only_in_b}") + + for key in keys_a & keys_b: + v1 = a[key] + v2 = b[key] + if v1 is None or v2 is None: + continue + if not math.isclose(float(v1), float(v2)): + warnings.warn(f"Value mismatch for key '{key}': a={v1}, b={v2}") + + +def pulser_specs( + device: qoolqit.Device, normalize: bool = False, check_against_qoolqit: bool = False +) -> dict[str, float | None]: + pulser_device = device._device + channel = pulser_device.channels["rydberg_global"] + specs: dict[str, float | None] = {} + specs["max_duration"] = ( + float(pulser_device.max_sequence_duration) + if pulser_device.max_sequence_duration is not None + else None + ) + specs["max_amplitude"] = channel.max_amp + specs["max_abs_detuning"] = channel.max_abs_detuning + specs["min_distance"] = pulser_device.min_atom_distance + specs["max_radial_distance"] = pulser_device.max_radial_distance + specs["min_avg_amp"] = channel.min_avg_amp + specs["dmm_bottom_detuning"] = None + + dmm_channels = list(getattr(pulser_device, "dmm_channels", {}).values()) + if dmm_channels: + specs["dmm_bottom_detuning"] = getattr(dmm_channels[0], "bottom_detuning", None) + + if not normalize: + return specs + + def _normalize(name: str, scale: float) -> None: + if specs[name] is not None: + specs[name] /= scale # type: ignore[operator] + + r0 = specs["min_distance"] + assert r0 is not None + _normalize("min_distance", r0) + _normalize("max_radial_distance", r0) + + C6 = pulser_device.interaction_coeff + J0 = C6 / (r0**6) + + _normalize("max_amplitude", J0) + _normalize("max_abs_detuning", J0) + _normalize("min_avg_amp", J0) + _normalize("dmm_bottom_detuning", J0) + + # J0 is in rad/us and t in ns, hence the factor 1000 + _normalize("max_duration", 1000.0 / J0) + + if check_against_qoolqit: + qq_specs = qoolqit_specs(device) + compare_specs(specs, qq_specs) + + return specs + + +def qoolqit_specs( + device: qoolqit.Device, complete_with_pulser: bool = False, check_against_pulser: bool = False +) -> dict[str, float | None]: + specs = device.specs + _pulser_specs = pulser_specs(device, normalize=True) + + def import_from_pulser_or_set(name: str, fallback: float | None = None) -> None: + if name in specs.keys(): + return + if complete_with_pulser: + specs[name] = _pulser_specs.get(name, fallback) + else: + specs[name] = fallback + + # Don't merge, update specific keys that are known not be in Qoolqit specs + import_from_pulser_or_set("min_avg_amp") + import_from_pulser_or_set("dmm_bottom_detuning") + + if check_against_pulser: + compare_specs(specs, _pulser_specs) + + return specs diff --git a/qubosolver/_drive_shaping/drive_shaper.py b/qubosolver/_drive_shaping/drive_shaper.py new file mode 100644 index 00000000..cefb4aee --- /dev/null +++ b/qubosolver/_drive_shaping/drive_shaper.py @@ -0,0 +1,217 @@ +from __future__ import annotations + +from abc import ABC, abstractmethod +from typing import cast + +import torch + +from qoolqit import Register, Drive + +from qubosolver import drive_shaping, QUBOInstance, SolverConfig, QUBOSolution, DriveType +from qubosolver.types import protocols + + +class BaseDriveShaper(ABC): + """ + Abstract base class for generating Qoolqit drives based on a QUBO problem. + + This class transforms the structure of a QUBOInstance into a quantum + waveform sequence or drive that can be applied to a physical register. The register + is passed at the time of drive generation, not during initialization. + + Attributes: + instance (QUBOInstance): The QUBO problem instance. + config (SolverConfig): The solver configuration. + drive (Drive, optional): A saved current drive obtained by `generate`. + backend (Backend): Backend to use. + device (Device): Device from backend. + + """ + + def __init__(self, instance: QUBOInstance, config: SolverConfig, backend: protocols.Backend): + """ + Initialize the drive shaping module with a QUBO instance. + + Args: + instance (QUBOInstance): The QUBO problem instance. + config (SolverConfig): The solver configuration. + backend (Backend): Backend to use. + """ + self.instance: QUBOInstance = instance + self.config: SolverConfig = config + self.drive: Drive | None = None + self.backend = backend + self.device = self.config.device + + # check if device allow DMM + self.dmm = self.config.drive_shaping.dmm and ( + len(list(self.config.device._device.dmm_channels.keys())) > 0 + ) + + @property + def qubo_coefficients(self) -> torch.Tensor: + return self.instance.coefficients + + @property + def qubo_normalized_coefficients(self) -> torch.Tensor: + return self.instance.normalized_coefficients + + @abstractmethod + def generate( + self, + register: Register, + ) -> tuple[Drive, QUBOSolution]: + """ + Generate a drive based on the problem and the provided register. + + Args: + register (Register): The physical register layout. + + Returns: + Drive: A generated Drive. + QUBOSolution: An instance of the qubo solution + """ + pass + + +class HeuristicDriveShaper(BaseDriveShaper): + """ + Heuristic schedule drive shaper. + + With DMM: + Final target encoding: + d_i = -alpha * Q_ii + + DMM convention in this stack: + WeightedDetuning waveform must be <= 0 + + Hence we encode the local final detuning as: + delta_i(T) = delta_g(T) + delta_dmm(T) * w_i + + with: + delta_g(T) = d_max + delta_dmm(T) = -(d_max - d_min) <= 0 + w_i = (d_max - d_i) / (d_max - d_min) in [0, 1] + + so that: + delta_i(T) = d_i + + Without DMM: + Only a global detuning is available, so the final detuning is chosen as: + delta_g(T) = mean(d_i) + and no weighted detunings are declared. + """ + + def generate(self, register: Register) -> tuple[Drive, QUBOSolution]: + device = self.config.device + dmm = self.config.drive_shaping.dmm + # Heuristic coefficient for omega + kappa = self.config.drive_shaping.heuristic_kappa + return drive_shaping.heuristic.generate( + register, self.instance, device=device, dmm=dmm, kappa=kappa + ), QUBOSolution(torch.Tensor(), torch.Tensor()) + + +class OptimizedDriveShaper(BaseDriveShaper): + """ + Drive shaper that uses optimization to find the best drive parameters for solving QUBOs. + Returns an optimized drive, the bitstrings, their counts, probabilities, and costs. + + Attributes: + drive (Drive): current drive. + best_cost (float): Current best cost. + best_bitstring (Tensor | list): Current best bitstring. + bitstrings (Tensor | list): List of current bitstrings obtained. + counts (Tensor | list): Frequencies of bitstrings. + probabilities (Tensor | list): Probabilities of bitstrings. + costs (Tensor | list): Qubo cost. + optimized_custom_qubo_cost (Callable[[str, torch.Tensor], float], optional): + Apply a different qubo cost evaluation during optimization. + Must be defined as: + `def optimized_custom_qubo_cost(bitstring: str, QUBO: torch.Tensor) -> float`. + Defaults to None, meaning we use the default QUBO evaluation. + optimized_custom_objective (Callable[[list, list, list, list, float, str], float], optional): + For bayesian optimization, one can change the output of + `self.run_simulation` to optimize differently. Instead of using the best cost + out of the samples, one can change the objective for an average, + or any function out of the form + `cost_eval = optimized_custom_objective(bitstrings, + counts, probabilities, costs, best_cost, best_bitstring)` + Defaults to None, which means we optimize using the best cost + out of the samples. + optimized_callback_objective (Callable[..., None], optional): Apply a callback + during bayesian optimization. Only accepts one input dictionary + created during optimization `d = {"x": x, "cost_eval": cost_eval}` + hence should be defined as: + `def callback_fn(d: dict) -> None:` + Defaults to None, which means no callback is applied. + """ + + def __init__( + self, + instance: QUBOInstance, + config: SolverConfig, + backend: protocols.Backend, + ): + """Instantiate an `OptimizedDriveShaper`. + + Args: + instance (QUBOInstance): Qubo instance. + config (SolverConfig): Configuration for solving. + backend (Backend): Backend to use during optimization. + + """ + super().__init__(instance, config, backend) + + def generate( + self, + register: Register, + ) -> tuple[Drive, QUBOSolution]: + """ + Generate a drive via optimization. + + Args: + register (Register): The physical register layout. + + Returns: + Drive: A generated Drive. + QUBOSolution: An instance of the qubo solution + """ + + config = drive_shaping.optimized.Config.from_drive_shaping_config(self.config.drive_shaping) + + return drive_shaping.optimized.generate( + register, self.instance, self.device, dmm=self.dmm, backend=self.backend, config=config + ) + + +def get_drive_shaper( + instance: QUBOInstance, + config: SolverConfig, + backend: protocols.Backend, +) -> BaseDriveShaper: + """ + Method that returns the correct DriveShaper based on configuration. + The correct drive shaping method can be identified using the config, and an + object of this driveshaper can be returned using this function. + + Args: + instance (QUBOInstance): The QUBO problem to embed. + config (SolverConfig): The solver configuration used. + backend (Backend): Backend to extract device from or to use + during drive shaping. + + Returns: + (BaseDriveShaper): The representative Drive Shaper object. + """ + if config.drive_shaping.drive_shaping_method == DriveType.HEURISTIC: + return HeuristicDriveShaper(instance, config, backend) + elif config.drive_shaping.drive_shaping_method == DriveType.OPTIMIZED: + return OptimizedDriveShaper(instance, config, backend) + elif issubclass(config.drive_shaping.drive_shaping_method, BaseDriveShaper): + return cast( + BaseDriveShaper, + config.drive_shaping.drive_shaping_method(instance, config, backend), + ) + else: + raise NotImplementedError diff --git a/qubosolver/_drive_shaping/heuristic.py b/qubosolver/_drive_shaping/heuristic.py new file mode 100644 index 00000000..b8e31d86 --- /dev/null +++ b/qubosolver/_drive_shaping/heuristic.py @@ -0,0 +1,97 @@ +from __future__ import annotations + +import numpy as np +import torch +import warnings + +import qoolqit +from qubosolver import QUBOInstance + +from .device_specs import pulser_specs as _pulser_specs +from .waveforms import weighted_detunings + + +def generate( + register: qoolqit.Register, + Q: QUBOInstance, + device: qoolqit.Device, + dmm: bool, + kappa: float = 0.25, +) -> qoolqit.Drive: + + # Hardware bounds + specs = device.specs + max_seq_duration: float = specs["max_duration"] or 1000.0 + pulser_specs = _pulser_specs(device) + use_dmm = dmm and (pulser_specs["dmm_bottom_detuning"] is not None) + + max_amplitude = specs["max_amplitude"] or 1.0 + max_abs_detuning = specs["max_abs_detuning"] or 9.0 + + det_amp_ratio = max_amplitude / max_abs_detuning + if kappa < det_amp_ratio: + warnings.warn( + f"heuristic_kappa is too small ({kappa}), you're likely to get a qoolqit CompilationError. Set it above {det_amp_ratio}." + ) + + n = Q.size + + # Target local final detunings + d = (-0.5 * torch.diag(Q.normalized_coefficients)).cpu().numpy() + d_min = np.min(d) + d_max = np.max(d) + + if use_dmm: + # Final global detuning is the top value, DMM pulls down locally + delta_g_T = d_max + + # DMM convention required by WeightedDetuning: + # waveform must be <= 0 + spread = max(0.0, d_max - d_min) + # if spread > 1e-15 and delta_dmm_max > 0.0: + if spread > 1e-15: + delta_dmm_T = -spread # must be <= 0 + denom = d_max - d_min + weights = ((d_max - d) / denom).clip(0.0, 1.0).tolist() + else: + use_dmm = False + delta_dmm_T = 0.0 + weights = [0.0] * n + else: + # No DMM: use a single global final detuning + delta_g_T = np.mean(d) + delta_dmm_T = 0.0 + weights = [0.0] * n + + omega_max = kappa * np.max(np.abs(d)) + # How to get max detuning ? + delta_0 = -np.max(np.abs(d)) + + # Amplitude waveform: 0 -> plateau -> 0 + eps = 1e-9 + amp_wave = qoolqit.waveforms.Interpolated( + max_seq_duration, + [eps, omega_max, omega_max, eps], + ) + + # Global detuning waveform: initial negative -> final target + det_wave = qoolqit.waveforms.Interpolated( + max_seq_duration, + [delta_0, delta_0, delta_g_T, delta_g_T], + ) + + # DMM weighted detunings + wdetunings = None + if use_dmm: + wdetunings = weighted_detunings( + register, + max_seq_duration, + weights, + final_detuning=delta_dmm_T, + ) + + return qoolqit.Drive( + amplitude=amp_wave, + detuning=det_wave, + dmm=wdetunings, + ) diff --git a/qubosolver/_drive_shaping/optimized.py b/qubosolver/_drive_shaping/optimized.py new file mode 100644 index 00000000..5db4e7d3 --- /dev/null +++ b/qubosolver/_drive_shaping/optimized.py @@ -0,0 +1,224 @@ +from __future__ import annotations + +from dataclasses import dataclass, field +import numpy as np +from skopt import gp_minimize +from collections.abc import Callable, Sequence +from typing import TypedDict +import torch + +import qoolqit +from qoolqit.execution.compilation_functions import CompilerProfile + +from qubosolver import utils, QUBOSolution, DriveShapingConfig, QUBOInstance, solvers +from .waveforms import weighted_detunings +from qubosolver.types import protocols, Bitstring, Matrix + + +def _default_objective(solution: QUBOSolution) -> float: + return solution.costs[0].item() if solution else float("inf") + + +class _CallbackObjectiveInput(TypedDict): + x: Sequence[float] + cost_eval: float + + +@dataclass +class Config: + + x0: list[float] = field( + default_factory=lambda: [ + 0.5, + 0.9, + 0.5, + -0.8, + 0.0, + 0.8, + ] + ) + n_calls: int = 20 + seed: int | None = None + qubo_cost: Callable[[Bitstring, Matrix], float] = utils.qubo_bitstring_cost + objective: Callable[[QUBOSolution], float] = _default_objective + callback_objective: Callable[[_CallbackObjectiveInput], None] = lambda data: None + + @staticmethod + def from_drive_shaping_config(config: DriveShapingConfig) -> Config: + cfg = Config() + cfg.x0 = ( + config.optimized_initial_omega_parameters + config.optimized_initial_detuning_parameters + ) + cfg.n_calls = config.optimized_n_calls + cfg.seed = config.optimized_seed + if config.optimized_custom_qubo_cost is not None: + cfg.qubo_cost = config.optimized_custom_qubo_cost + if config.optimized_custom_objective is not None: + cfg.objective = config.optimized_custom_objective + if config.optimized_callback_objective is not None: + cfg.callback_objective = config.optimized_callback_objective + + return cfg + + +def _compute_norm_weights(Q: QUBOInstance) -> list[float]: + """Compute normalization weights. + + Returns: + list[float]: normalization weights. + """ + weights_list = torch.abs(torch.diag(Q.coefficients)).tolist() + max_node_weight = max(weights_list) if weights_list else 1.0 + norm_weights_list = [ + (1 - (w / max_node_weight)) if max_node_weight != 0 else 0.0 for w in weights_list + ] + return norm_weights_list + + +def _build_drive( + params: Sequence[float], + Q: QUBOInstance, + register: qoolqit.Register, + device_specs: dict[str, float | None], + dmm: bool, +) -> qoolqit.Drive: + """Build the drive from a list of parameters for the objective. + + Args: + params (list): List of parameters. + + Returns: + Drive: Drive sequence. + """ + max_seq_duration: float = device_specs["max_duration"] or 1e3 + max_amplitude: float = device_specs["max_amplitude"] or 1e4 + max_detuning: float = device_specs["max_abs_detuning"] or 1e4 + + amp_params = [1e-9] + list(params[:3]) + [1e-9] + # FIXME: det_params of length 4 ? with last param as final det for dmm? + det_params = list(params[3:]) + amp_params = [p * max_amplitude for p in amp_params] + det_params = [p * max_detuning for p in det_params] + + amp_wave = qoolqit.waveforms.Interpolated(max_seq_duration, amp_params) + det_wave = qoolqit.waveforms.Interpolated(max_seq_duration, det_params) + + wdetunings = None + final_detuning = det_params[-1] + if dmm and final_detuning > 0: + wdetunings = weighted_detunings( + register, + max_seq_duration, + _compute_norm_weights(Q), + final_detuning=-final_detuning, + ) + + shaped_drive = qoolqit.Drive(amplitude=amp_wave, detuning=det_wave, dmm=wdetunings) + + return shaped_drive + + +def _run_simulation( + register: qoolqit.Register, + drive: qoolqit.Drive, + Q: torch.Tensor, + device: qoolqit.Device, + backend: protocols.Backend, + config: Config, +) -> QUBOSolution: + """Run a quantum program using backend and returns + a tuple of (bitstrings, counts, probabilities, costs, best cost, best bitstring). + + Args: + register (Register): register of quantum program. + drive (Drive): drive to run on backend. + QUBO (torch.Tensor): Qubo coefficients. + convert_to_tensor (bool, optional): Convert tuple components to tensors. + Defaults to True. + + Returns: + tuple: tuple of (bitstrings, counts, probabilities, costs, best cost, best bitstring) + """ + try: + job = solvers.analog_quantum_sample( + backend, device, drive, register, CompilerProfile.WORKING_POINT + ) + solution = QUBOSolution.from_results(job.results()) + costs = [config.qubo_cost(b, Q) for b in solution.bitstrings] + solution.costs = torch.Tensor(costs) + solution.probabilities = solution.compute_probabilities() + solution.sort_by_cost() + return solution + except Exception as e: + print(f"Simulation failed: {e}") + return QUBOSolution() + + +def generate( + register: qoolqit.Register, + Q: QUBOInstance, + device: qoolqit.Device, + dmm: bool, + backend: protocols.Backend, + config: Config, +) -> tuple[qoolqit.Drive, QUBOSolution]: + n_amp = 3 + n_det = 3 + + eps = 0.0001 + zero = eps + one = 1.0 - eps + + bounds = [(zero, one)] * n_amp + [(-one, -zero)] + [(-one, one)] * (n_det - 2) + [(zero, one)] + + def run(x: list[float], eval: bool = True) -> tuple[float, QUBOSolution, qoolqit.Drive]: + + solution = QUBOSolution() + drive = _build_drive( + x, + Q, + register, + device.specs, + dmm, + ) + + try: + solution = _run_simulation( + register, + drive, + Q.coefficients, + device, + backend, + config, + ) + if eval: + cost_eval = config.objective(solution) + if not np.isfinite(cost_eval): + print(f"[Warning] Non-finite cost encountered: {cost_eval} at x={x}") + cost_eval = 1e4 + else: + cost_eval = float("nan") + + except Exception as e: + print(f"[Exception] Error during simulation at x={x}: {e}") + cost_eval = 1e4 + return cost_eval, solution, drive + + def objective(x: list[float]) -> float: + cost_eval, _, _ = run(x) + config.callback_objective({"x": x, "cost_eval": cost_eval}) + + return cost_eval + + opt_result = gp_minimize( + objective, + bounds, + x0=config.x0, + n_calls=config.n_calls, + random_state=config.seed, + ) + + best_params = opt_result.x if opt_result else config.x0 + _, solution, drive = run(best_params, eval=False) + + return drive, solution diff --git a/qubosolver/pipeline/waveforms.py b/qubosolver/_drive_shaping/waveforms.py similarity index 71% rename from qubosolver/pipeline/waveforms.py rename to qubosolver/_drive_shaping/waveforms.py index a8b89182..f5b3176b 100644 --- a/qubosolver/pipeline/waveforms.py +++ b/qubosolver/_drive_shaping/waveforms.py @@ -1,16 +1,14 @@ from __future__ import annotations -from qoolqit.register import Register -from qoolqit.drive import DetuningMapModulator -from qoolqit.waveforms import Constant +import qoolqit def weighted_detunings( - embedding: Register, + embedding: qoolqit.Register, duration: float, norm_weights: list[float], final_detuning: float, -) -> DetuningMapModulator: +) -> qoolqit.drive.DetuningMapModulator: """Create the list of weighted detuning for a drive. Args: @@ -23,8 +21,8 @@ def weighted_detunings( DetuningMapModulator: DetuningMapModulator with a constant waveform for QUBO solving. """ - waveform = Constant(duration, final_detuning) - return DetuningMapModulator( + waveform = qoolqit.waveforms.Constant(duration, final_detuning) + return qoolqit.drive.DetuningMapModulator( weights={embedding.qubits_ids[i]: w for i, w in enumerate(norm_weights)}, waveform=waveform, ) diff --git a/qubosolver/_embedding/__init__.py b/qubosolver/_embedding/__init__.py new file mode 100644 index 00000000..9199c867 --- /dev/null +++ b/qubosolver/_embedding/__init__.py @@ -0,0 +1,8 @@ +from __future__ import annotations + +from qubosolver.embedding import blade, greedy + +__all__ = [ + "blade", + "greedy", +] diff --git a/qubosolver/algorithms/greedy/__init__.py b/qubosolver/_embedding/algorithms/greedy/__init__.py similarity index 100% rename from qubosolver/algorithms/greedy/__init__.py rename to qubosolver/_embedding/algorithms/greedy/__init__.py diff --git a/qubosolver/algorithms/greedy/greedy.py b/qubosolver/_embedding/algorithms/greedy/greedy.py similarity index 96% rename from qubosolver/algorithms/greedy/greedy.py rename to qubosolver/_embedding/algorithms/greedy/greedy.py index 6d427758..6348c972 100644 --- a/qubosolver/algorithms/greedy/greedy.py +++ b/qubosolver/_embedding/algorithms/greedy/greedy.py @@ -1,14 +1,15 @@ from __future__ import annotations import copy +from collections.abc import Callable import typing -from typing import Any, Callable, Dict, List, Optional, Tuple, cast +from typing import Any, Optional, cast import torch from pulser.register.register_layout import RegisterLayout from qoolqit.devices.device import DigitalAnalogDevice -from qubosolver.qubo_types import LayoutType +from qubosolver.types.enums import LayoutType # Optional imports for animation; guarded so library usage stays safe in non-notebook envs. try: # pragma: no cover @@ -113,7 +114,7 @@ def get_best(self, Q: torch.Tensor, positioned: set, all_vertices: set) -> Any: to the already-positioned set. """ all_vertices = all_vertices.difference(positioned) - node_contributes: List[Tuple[int, float]] = [] + node_contributes: list[tuple[int, float]] = [] for u in all_vertices: s: float = 0.0 for j in positioned: @@ -134,7 +135,7 @@ def optimize_position( all_traps: set, used_traps: set, return_candidates: bool = False, - ) -> tuple[Any, Any, Any] | tuple[Any, Any, Any, List[Tuple[int, float]]]: + ) -> tuple[Any, Any, Any] | tuple[Any, Any, Any, list[tuple[int, float]]]: """ Evaluate all available traps p for node u and pick the one that minimizes: s(p) = sum_{j in positioned} Z[u, j, p, trap(j)]. @@ -150,7 +151,7 @@ def optimize_position( choice_p: int = -1 choice_coordinates: tuple = (None, None) min_val: float = float("inf") - candidates: List[Tuple[int, float]] = [] + candidates: list[tuple[int, float]] = [] for p in available_traps: s = 0.0 @@ -180,7 +181,7 @@ def greedy_algorithm( v: int, results: dict, params: dict, - on_step: Optional[Callable[[Dict[str, Any]], None]] = None, + on_step: Optional[Callable[[dict[str, Any]], None]] = None, ) -> dict: """ Greedy loop starting from node v. If `on_step` is provided, emit a @@ -207,8 +208,8 @@ def greedy_algorithm( n_extra_traps = n_traps - n # helpers for instrumentation - def _trap_of_from_coords() -> Dict[int, int]: # pragma: no cover - out: Dict[int, int] = {} + def _trap_of_from_coords() -> dict[int, int]: # pragma: no cover + out: dict[int, int] = {} for node_id, coord in positioned_coords.items(): out[node_id] = int(self.MAPPING_COORDS_POSITIONS[coord]) return out @@ -254,7 +255,7 @@ def _trap_of_from_coords() -> Dict[int, int]: # pragma: no cover ) # Help mypy: explicitly cast 4-tuple _, u_coordinates, _, candidates = typing.cast( - Tuple[Any, Any, Any, List[Tuple[int, float]]], res4 + tuple[Any, Any, Any, list[tuple[int, float]]], res4 ) candidates.sort(key=lambda t: t[1]) # ascending by mismatch else: @@ -268,7 +269,7 @@ def _trap_of_from_coords() -> Dict[int, int]: # pragma: no cover return_candidates=False, ) # Help mypy: explicitly cast 3-tuple - _, u_coordinates, _ = typing.cast(Tuple[Any, Any, Any], res3) + _, u_coordinates, _ = typing.cast(tuple[Any, Any, Any], res3) candidates = [] distance = torch.tensor(u_coordinates).norm().item() @@ -372,7 +373,7 @@ def _trap_of_from_coords() -> Dict[int, int]: # pragma: no cover # ---------------------------- def _render_animation( # pragma: no cover self, - frames: List[Dict[str, Any]], + frames: list[dict[str, Any]], all_coords_np: "np.ndarray", spacing: float, layout_name: str, @@ -416,7 +417,7 @@ def _render_animation( # pragma: no cover [], [], s=130, color="tab:green", edgecolor="k", linewidths=0.8, zorder=3 ) # We manage labels manually to remove/recreate them each frame - labels: List[Any] = [] + labels: list[Any] = [] # ---- Info panel ax_info.axis("off") @@ -629,7 +630,7 @@ def launch_greedy( self, Q: torch.Tensor, params: dict, - on_step: Optional[Callable[[Dict[str, Any]], None]] = None, + on_step: Optional[Callable[[dict[str, Any]], None]] = None, ) -> Any: """ Run greedy from each start node and keep the best result. @@ -661,11 +662,11 @@ def launch_greedy( anim_flag = bool(params.get("animation", False) or params.get("draw_steps", False)) instrument = bool(params.get("draw_steps", False) or on_step is not None or anim_flag) - frames: List[Dict[str, Any]] = [] + frames: list[dict[str, Any]] = [] if instrument: # pragma: no cover - def _collector(state: Dict[str, Any]) -> None: + def _collector(state: dict[str, Any]) -> None: if on_step is not None: try: on_step(state) diff --git a/qubosolver/_embedding/blade.py b/qubosolver/_embedding/blade.py new file mode 100644 index 00000000..58b2909a --- /dev/null +++ b/qubosolver/_embedding/blade.py @@ -0,0 +1,20 @@ +from __future__ import annotations + +from typing import TypeAlias + +from qubosolver import QUBOInstance +from qoolqit import Register +from qoolqit.embedding import Blade, BladeConfig + +#  TODO: Using `type` statement when Python >= 3.12 +Config: TypeAlias = BladeConfig + + +def embed(instance: QUBOInstance, config: Config, normalize: bool = True) -> Register: + _blade = Blade(config) + graph = _blade.embed(instance.coefficients.numpy()) + if normalize: + graph.rescale_coords(spacing=1.0001) + register = Register.from_graph(graph) + + return register diff --git a/qubosolver/_embedding/embedder.py b/qubosolver/_embedding/embedder.py new file mode 100644 index 00000000..56d298e1 --- /dev/null +++ b/qubosolver/_embedding/embedder.py @@ -0,0 +1,128 @@ +from __future__ import annotations + +import typing +from abc import ABC, abstractmethod +import warnings + +from qoolqit import Register + + +from qubosolver import QUBOInstance, embedding, EmbedderType, SolverConfig +from qubosolver.types import protocols + +warnings.filterwarnings("ignore", module="pulser") + + +class BaseEmbedder(ABC): + """ + Abstract base class for all embedders. + + Prepares the geometry (register) of atoms based on the QUBO instance. + Returns a Register compatible with Pasqal/Pulser devices. + """ + + def __init__(self, instance: QUBOInstance, config: SolverConfig, backend: protocols.Backend): + """ + Args: + instance (QUBOInstance): The QUBO problem to embed. + config (SolverConfig): The Solver Configuration. + """ + self.instance: QUBOInstance = instance + self.config: SolverConfig = config + self.register: Register | None = None + self.backend = backend + + # TODO: remove when bumping to qoolqit v1 + # for converting to qoolqit + self._distance_conversion = self.config.device.converter.factors[2] + + @abstractmethod + def embed(self) -> Register: + """ + Creates a layout of atoms as the register. + + Returns: + Register: The register. + """ + ... + + +class BLaDEmbedder(BaseEmbedder): + + def embed(self) -> Register: + + embed_config = self.config.embedding + default = embedding.blade.Config() + step_per_round = embed_config.blade_steps_per_round + if step_per_round is None: + step_per_round = default.steps_per_round + if embed_config.blade_starting_positions is not None: + starting_positions = embed_config.blade_starting_positions.numpy() + else: + starting_positions = None + + min_distance = self.config.embedding.min_distance + max_radial_distance = self.config.device.specs["max_radial_distance"] + if min_distance is None or max_radial_distance is None: + device = self.config.device + max_min_dist_ratio = None + else: + device = None + max_min_dist_ratio = max_radial_distance / min_distance + + config = embedding.blade.Config( + steps_per_round=step_per_round, + starting_positions=starting_positions, + dimensions=tuple(embed_config.blade_dimensions), + max_min_dist_ratio=max_min_dist_ratio, + device=device, + ) + return embedding.blade.embed(self.instance, config, normalize=(min_distance is not None)) + + +class GreedyEmbedder(BaseEmbedder): + """Create an embedding in a greedy fashion. + + At each step, place one logical node onto one trap to minimize the + incremental mismatch between the logical QUBO matrix Q and the physical + interaction matrix U (approx. C / ||r_i - r_j||^6). + """ + + def embed(self) -> Register: + """ + Creates a layout of atoms as the register. + + Returns: + Register: The register. + """ + config = embedding.greedy.Config.from_embedding_config(self.config.embedding) + normalize = self.config.embedding.min_distance is not None + return embedding.greedy.embed(self.instance, self.config.device, config, normalize) + + +def get_embedder( + instance: QUBOInstance, config: SolverConfig, backend: protocols.Backend +) -> BaseEmbedder: + """ + Method that returns the correct embedder based on configuration. + The correct embedding method can be identified using the config, and an + object of this embedding can be returned using this function. + + Args: + instance (QUBOInstance): The QUBO problem to embed. + config (Device): The quantum device to target. + + Returns: + (BaseEmbedder): The representative embedder object. + """ + + if config.embedding.embedding_method == EmbedderType.BLADE: + return BLaDEmbedder(instance, config, backend) + elif config.embedding.embedding_method == EmbedderType.GREEDY: + return GreedyEmbedder(instance, config, backend) + elif issubclass(config.embedding.embedding_method, BaseEmbedder): + return typing.cast( + BaseEmbedder, config.embedding.embedding_method(instance, config, backend) + ) + else: + raise NotImplementedError diff --git a/qubosolver/_embedding/greedy.py b/qubosolver/_embedding/greedy.py new file mode 100644 index 00000000..10b851f3 --- /dev/null +++ b/qubosolver/_embedding/greedy.py @@ -0,0 +1,117 @@ +from __future__ import annotations + +from dataclasses import dataclass +import numpy as np +import pathlib +import torch + +import qoolqit + +from qubosolver import QUBOInstance, LayoutType, EmbeddingConfig +from qubosolver.embedding.algorithms import greedy + + +@dataclass +class Config: + traps: int = -1 + spacing: float = 7.0 + layout: LayoutType = LayoutType.TRIANGULAR + draw_steps: bool = False + animation_save_path: pathlib.Path | None = None + + def update_from_device(self, device: qoolqit.Device) -> None: + if self.traps == -1: + self.traps = _number_of_traps_from_device(device) + + _device = device._device + if hasattr(_device, "min_atom_distance"): + self.spacing = max(self.spacing, float(_device.min_atom_distance)) + + @staticmethod + def from_embedding_config(config: EmbeddingConfig) -> Config: + cfg = Config() + cfg.traps = config.greedy_traps + cfg.spacing = config.greedy_spacing + cfg.layout = EmbeddingConfig._normalize_layout(config.greedy_layout) + cfg.draw_steps = config.draw_steps + path = config.animation_save_path + cfg.animation_save_path = pathlib.Path(path) if path else None + + return cfg + + +def _number_of_traps_from_device(device: qoolqit.Device) -> int: + """Determine the number of traps to use based on the device constraints. + + Inspects the device's layout and atom number limits to derive an + appropriate trap count. The resolution order is: + + 1. ``max_layout_traps`` – if the device exposes a hard trap limit, use it directly. + 2. ``max_atom_num`` / ``max_layout_filling`` – if only an atom-number limit is + available, derive the minimum number of traps needed to accommodate that + many atoms at the device's maximum filling ratio. + 3. Fallback – return ``200`` when neither property is set. + + Args: + device (Device): The quantum device whose constraints are inspected. + + Returns: + int: The number of traps to allocate for the embedding. + """ + + if device._device.max_layout_traps: + return device._device.max_layout_traps + + if device._device.max_atom_num: + return int(np.ceil(device._device.max_atom_num / device._device.max_layout_filling)) + + return 200 + + +def embed( + instance: QUBOInstance, + device: qoolqit.Device, + config: Config = Config(), + normalize: bool = True, +) -> qoolqit.Register: + """ + Creates a layout of atoms as the register. + + Returns: + Register: The register. + """ + config.update_from_device(device) + + if config.traps < instance.size: + raise ValueError( + "Number of traps must be at least equal to the number of atoms on the register." + ) + + # build params for the Greedy algorithm + params = { + "device": device._device, + "layout": config.layout, + "traps": config.traps, + "spacing": config.spacing, + # animation controls (all read by Greedy) + "draw_steps": config.draw_steps, # collect per-step data + "animation": config.draw_steps, # render animation after run + "animation_save_path": config.animation_save_path, # optional export + } + + # --- Call Greedy (unchanged public signature) + _, _, coords, _, _ = greedy.Greedy().launch_greedy( + Q=instance.coefficients, + params=params, + ) + if normalize: + min_reg_distance = torch.cdist(coords, coords).fill_diagonal_(float("inf")).min() + coords *= 1.0001 / min_reg_distance + else: + distance_conversion = device.converter.factors[2] + coords /= distance_conversion + + # build the register (unchanged) + qubits = {f"q{i}": coord for i, coord in enumerate(coords)} + register = qoolqit.Register(qubits) + return register diff --git a/qubosolver/io/utils.py b/qubosolver/_io/utils.py similarity index 99% rename from qubosolver/io/utils.py rename to qubosolver/_io/utils.py index 0bf3c342..f4ddce5c 100644 --- a/qubosolver/io/utils.py +++ b/qubosolver/_io/utils.py @@ -7,8 +7,9 @@ from contextlib import nullcontext from typing import TYPE_CHECKING, overload, Sized +from qubosolver.types._checks import _RUNTIME_CHECKS -if TYPE_CHECKING: +if TYPE_CHECKING or _RUNTIME_CHECKS: from typing import Any, Union, IO, Literal, TypeVar from typing_extensions import Buffer from contextlib import AbstractContextManager diff --git a/qubosolver/_solvers/__init__.py b/qubosolver/_solvers/__init__.py new file mode 100644 index 00000000..fcb653e2 --- /dev/null +++ b/qubosolver/_solvers/__init__.py @@ -0,0 +1,17 @@ +from __future__ import annotations + +from .bitflip import iterative_bitflip_local_search +from .quantum import analog_quantum_sample +from .trivial import trivial_solution_search +from .cplex import cplex +from .tabu_search import tabu_search +from .simulated_annealing import simulated_annealing + +__all__ = [ + "iterative_bitflip_local_search", + "analog_quantum_sample", + "trivial_solution_search", + "cplex", + "tabu_search", + "simulated_annealing", +] diff --git a/qubosolver/pipeline/basesolver.py b/qubosolver/_solvers/basesolver.py similarity index 58% rename from qubosolver/pipeline/basesolver.py rename to qubosolver/_solvers/basesolver.py index 6759b41b..8d469a29 100644 --- a/qubosolver/pipeline/basesolver.py +++ b/qubosolver/_solvers/basesolver.py @@ -2,28 +2,25 @@ from abc import ABC, abstractmethod import inspect -import json -import torch -from qoolqit import QuantumProgram from qoolqit.execution.job import Job from qubosolver import QUBOInstance -from qubosolver.config import SolverConfig, compiler_profile, max_duration_ratio +from qubosolver.config import SolverConfig, compiler_profile from qubosolver.data import QUBOSolution -from qubosolver.qubo_types import SolutionStatusType -from qubosolver.pipeline.fixtures import Fixtures +from qubosolver.types.enums import SolutionStatusType import qubosolver.io.utils as io_utils +from qubosolver import transforms, solvers from typing import TYPE_CHECKING if TYPE_CHECKING: - from typing import Optional, Callable, Any + from collections.abc import Callable + from typing import Optional, Any from typing_extensions import Self from qoolqit import Register, Drive - from pulser.backend import Results class BaseSolver(ABC): @@ -54,8 +51,12 @@ def __init__(self, instance: QUBOInstance, config: SolverConfig | None = None): self.backend = self.config.backend self.device = self.config.device - self.fixtures = Fixtures(self.instance, self.config) - self.n_fixed_variables_preprocessing = 0 + @property + def n_fixed_variables_preprocessing(self) -> int: + if isinstance(self.instance, transforms.variable_fixing.QUBOInstance): + return self.instance.n_fixed_indices + else: + return 0 @abstractmethod def solve(self) -> QUBOSolution: @@ -121,52 +122,11 @@ def submit( RuntimeError: If wait=False is specified for local backends, as async execution is not supported on local backends. """ - program = QuantumProgram( - register=embedding, - drive=drive, - ) - program.compile_to( - self.device, - profile=compiler_profile(self.config), - device_max_duration_ratio=max_duration_ratio(self.config), + return solvers.analog_quantum_sample( + self.backend, self.device, drive, embedding, compiler_profile(self.config) ) - return self.backend.run(program) - - @staticmethod - def parse_results( - results: Results, - ) -> tuple[torch.Tensor, torch.Tensor]: - """ - Parse execution results from quantum backends into standardized tensor format. - - Extracts bitstring measurements and their corresponding counts from either - remote or local backend execution results. Handles different result formats - and converts them into PyTorch tensors for further processing. - - Args: - results (RemoteResults | Sequence[Results]): Execution results from the - quantum backend. Can be either RemoteResults from a remote backend - or a sequence of Results from a local emulator. - - Returns: - tuple[torch.Tensor, torch.Tensor]: A tuple containing: - - bitstrings (torch.Tensor): Tensor of shape (n_samples, n_qubits) - containing the measured bitstrings as integers (0 or 1). Returns - empty tensor (0, 0) if no measurements were obtained. - - counts (torch.Tensor): Tensor of shape (n_samples,) containing - the number of times each corresponding bitstring was measured. - """ - counter = results.final_bitstrings - bitstrings = torch.tensor( - [list(map(int, list(b))) for b in list(counter.keys())], dtype=torch.int64 - ) - if bitstrings.numel() == 0: - bitstrings = torch.empty((0, 0), dtype=torch.int64) - counts = torch.tensor(list(map(int, list(counter.values()))), dtype=torch.int64) - return bitstrings, counts - - def execute(self, drive: Drive, embedding: Register) -> tuple[torch.Tensor, torch.Tensor]: + def execute(self, drive: Drive, embedding: Register) -> QUBOSolution: """ Execute the drive schedule on the backend and retrieve the solution. # TODO: We do not currently execute using the async run. @@ -181,15 +141,7 @@ def execute(self, drive: Drive, embedding: Register) -> tuple[torch.Tensor, torc tuple: A tuple of (bitstrings, counts) from the execution. """ job = self.submit(drive, embedding) - bitstrings, counts = self.parse_results(job.results()) - - if self.config.drive_shaping.optimized_re_execute_opt_drive and ( - bitstrings.numel() == 0 or counts.numel() == 0 - ): - job = self.submit(drive, embedding) - bitstrings, counts = self.parse_results(job.results()) - - return bitstrings, counts + return QUBOSolution.from_results(job.results()) def draw_sequence(self, drive: Drive, embedding: Register) -> None: """Draw sequence of the `QuantumProgram` submitted. @@ -199,82 +151,18 @@ def draw_sequence(self, drive: Drive, embedding: Register) -> None: embedding (Register): embedding program is defined over. """ if self.config.use_quantum: - program = QuantumProgram( - register=embedding, - drive=drive, - ) - program.compile_to( - self.device, - profile=compiler_profile(self.config), - device_max_duration_ratio=max_duration_ratio(self.config), + program = solvers.quantum._quantum_program( + self.device, drive, embedding, compiler_profile(self.config) ) program.draw(compiled=True) def _trivial_solution(self) -> Optional[QUBOSolution]: - """ - Check for the two trivial QUBO cases: - 1) all coefficients >= 0 → solution = 0^n - 2) all coefficients <= 0 → solution = 1^n - 3) diagonal qubo, negative coeffs gets 1, positive gets 0 - - Returns: - QUBOSolution if a trivial case applies, else None. - """ - coeffs = self.instance.coefficients # torch.Tensor (n, n) - n = self.instance.size - device, dtype = coeffs.device, coeffs.dtype - - # Case 1: all coeffs >= 0 → x = [0,...,0] - if torch.all(coeffs >= 0): - raw = torch.zeros(n, dtype=torch.int64, device=device) - # always make a batch of one: shape (1, n) - batch = raw.unsqueeze(0) - cost = self.instance.evaluate_solution(raw) - return QUBOSolution( - bitstrings=batch, - costs=torch.tensor([cost], dtype=dtype, device=device), - solution_status=SolutionStatusType.TRIVIALZERO, - ) - - # Case 2: all coeffs <= 0 → x = [1,...,1] - if torch.all(coeffs <= 0): - raw = torch.ones(n, dtype=torch.int64, device=device) - # always make a batch of one: shape (1, n) - batch = raw.unsqueeze(0) - cost = self.instance.evaluate_solution(raw) - return QUBOSolution( - bitstrings=batch, - costs=torch.tensor([cost], dtype=dtype, device=device), - solution_status=SolutionStatusType.TRIVIALONE, - ) - - # Case 3: diagonal cases - # negative coeffs gets 1, positive gets 0 - diagonal = torch.diag(coeffs) - if (torch.diag(diagonal) == coeffs).all(): - raw = (diagonal < 0).long() - cost = self.instance.evaluate_solution(raw) - batch = raw.unsqueeze(0) - return QUBOSolution( - bitstrings=batch, - costs=torch.tensor([cost], dtype=dtype, device=device), - solution_status=SolutionStatusType.TRIVIALDIAGONAL, - ) - return None + return solvers.trivial_solution_search(self.instance) def preprocess(self) -> None: """Apply preprocessing on instance to reduce its size.""" if self.config.do_preprocessing: - # Apply preprocessing and change the solved QUBO by the reduced one - self.fixtures.preprocess() - if ( - self.fixtures.reduced_qubo.coefficients is not None - and len(self.fixtures.reduced_qubo.coefficients) > 0 - and self.fixtures.n_fixed_variables < self.instance.size - ): - - self.instance = self.fixtures.reduced_qubo - self.n_fixed_variables_preprocessing = self.fixtures.n_fixed_variables + self.instance = transforms.variable_fixing.apply_recursively(self.instance) def post_process_fixation(self, solution: QUBOSolution) -> QUBOSolution: """Post-process fixations of the preprocessing and restore the original QUBO. @@ -285,10 +173,17 @@ def post_process_fixation(self, solution: QUBOSolution) -> QUBOSolution: Returns: QUBOSolution: New restored solution if preprocessing was applied. """ - if self.config.do_preprocessing: - solution = self.fixtures.post_process_fixation(solution) - self.instance = self.fixtures.instance - return solution + # Means that preprocessing was not applied + if not self.config.do_preprocessing: + return solution + + assert isinstance(self.instance, transforms.variable_fixing.QUBOInstance) + self.instance.solution = solution + new_solution = transforms.variable_fixing.unapply(self.instance) + new_solution.solution_status = SolutionStatusType.PREPROCESSED + self.instance = self.instance._parent_instance + + return new_solution def post_process(self, solution: QUBOSolution) -> QUBOSolution: """Apply post-processing. @@ -299,10 +194,18 @@ def post_process(self, solution: QUBOSolution) -> QUBOSolution: Returns: QUBOSolution: New postprocessed solution. """ + if not self.config.do_postprocessing: + return solution + + new_solution = solvers.iterative_bitflip_local_search(solution, self.instance) + + # FIXME: Use bitfield ? + if new_solution.solution_status == SolutionStatusType.PREPROCESSED: + new_solution.solution_status = SolutionStatusType.PREPOSTPROCESSED + else: + new_solution.solution_status = SolutionStatusType.POSTPROCESSED - if self.config.do_postprocessing: - solution = self.fixtures.postprocess(solution) - return solution + return new_solution @classmethod def save(cls, file_like: io_utils.FileLike[bytes], solver: Self) -> None: @@ -327,15 +230,13 @@ def save(cls, file_like: io_utils.FileLike[bytes], solver: Self) -> None: directly to the provided file-like object. """ with io_utils.open(file_like, "wb") as f: + io_utils.save(f, "?", solver.config.do_preprocessing) + io_utils.save(f, "?", solver.config.do_postprocessing) if solver.config.do_preprocessing: - QUBOInstance.save(f, solver.fixtures.instance) + assert isinstance(solver.instance, transforms.variable_fixing.QUBOInstance) + transforms.variable_fixing.QUBOInstance.save(f, solver.instance) else: QUBOInstance.save(f, solver.instance) - io_utils.save(f, "?", solver.config.do_preprocessing) - io_utils.save(f, "?", solver.config.do_postprocessing) - - fixed_var_json = json.dumps(solver.fixtures.fixed_var_dict_list) - io_utils.save_string(f, fixed_var_json) @classmethod def load(cls, file_like: io_utils.FileLike[bytes]) -> Self: @@ -364,23 +265,19 @@ def load(cls, file_like: io_utils.FileLike[bytes]) -> Self: methods are disabled to prevent incorrect usage of an incompletely loaded solver. """ with io_utils.open(file_like, "rb") as f: - instance = QUBOInstance.load(f) do_preprocessing: bool = io_utils.load(f, "?") do_postprocessing: bool = io_utils.load(f, "?") - fixed_var_json = io_utils.load_string(f) + instance = ( + transforms.variable_fixing.QUBOInstance.load(f) + if do_preprocessing + else QUBOInstance.load(f) + ) config = SolverConfig( do_preprocessing=do_preprocessing, do_postprocessing=do_postprocessing ) solver = cls(instance, config) - def decode_int_keys(obj: dict) -> dict: - return {int(k): v for k, v in obj.items()} - - solver.fixtures.fixed_var_dict_list = json.loads( - fixed_var_json, object_hook=decode_int_keys - ) - # Solver is incompletely loaded, most functions are unvailable for name, _ in inspect.getmembers(solver, predicate=inspect.ismethod): if not name.startswith("__") and name not in ( diff --git a/qubosolver/_solvers/bitflip.py b/qubosolver/_solvers/bitflip.py new file mode 100644 index 00000000..1a2bac90 --- /dev/null +++ b/qubosolver/_solvers/bitflip.py @@ -0,0 +1,87 @@ +from __future__ import annotations + +import numpy as np +import torch +from collections.abc import Callable + + +from qubosolver import QUBOInstance, QUBOSolution +from qubosolver.types import make_vector, make_vectori, make_bitstrings + + +def _bit_flip_local_search( + qubo_func: Callable[[np.ndarray], float], s: np.ndarray, shuffle: bool = True +) -> tuple[np.ndarray, float]: + """ + Performs a local search by flipping bits to improve the objective value. + + Args: + qubo_func: Function that computes the objective value for a solution array. + s (np.ndarray): Binary array representing a candidate solution. + shuffle (bool, optional): Shuffle to diversify + + Returns: + tuple[np.ndarray, float]: The improved solution and its objective value. + """ + s_current = s.copy() + current_objective = qubo_func(s_current) + while True: + best_idx = None + best_obj = current_objective + indices = list(range(len(s_current))) + if shuffle: + np.random.shuffle(indices) # option to diversify + # Evaluate all possible flips, keep best + for i in indices: + s_new = s_current.copy() + s_new[i] = 1 - s_new[i] + new_objective = qubo_func(s_new) + if new_objective < best_obj: + best_obj = new_objective + best_idx = i + # If no improvements, stop + if best_idx is None: + break + # Apply best flip + s_current[best_idx] = 1 - s_current[best_idx] + current_objective = best_obj + return s_current, current_objective + + +def iterative_bitflip_local_search(solution: QUBOSolution, Q: QUBOInstance) -> QUBOSolution: + + # If there are no bitstrings, return the solution unchanged. + if solution.bitstrings.numel() == 0: + return solution + + # Define an objective function that uses the existing evaluate_solution method. + def qubo_objective(s_arr: np.ndarray) -> float: + # Convert the solution array to a list of integers + return Q.evaluate_solution(s_arr.tolist()) + + num_solutions = solution.bitstrings.shape[0] + improved_bitstrings = make_bitstrings(num_solutions, Q.size) + improved_costs = make_vector(num_solutions) + + for idx in range(num_solutions): + # Get the current solution (row) as a numpy array of integers. + s_orig = solution.bitstrings[idx].detach().cpu().numpy().astype(int) + # Apply bit-flip local search to improve the solution. + s_improved, new_cost = _bit_flip_local_search(qubo_objective, s_orig) + improved_bitstrings[idx, :] = torch.from_numpy(s_improved).to(torch.int8) + improved_costs[idx] = new_cost + + unique_bitstrings, inverse = torch.unique(improved_bitstrings, dim=0, return_inverse=True) + n = unique_bitstrings.shape[0] + + # Update the solution object. + solution.bitstrings = unique_bitstrings + solution.costs = make_vector(n).scatter_reduce( + dim=0, index=inverse, src=improved_costs, reduce="amin", include_self=False + ) + solution.counts = make_vectori(n).scatter_reduce( + dim=0, index=inverse, src=solution.counts, reduce="sum", include_self=False + ) + solution.probabilities = solution.compute_probabilities() + + return solution diff --git a/qubosolver/classical_solver/classical_solver.py b/qubosolver/_solvers/classical_solver.py similarity index 67% rename from qubosolver/classical_solver/classical_solver.py rename to qubosolver/_solvers/classical_solver.py index 10b20517..1a8359ec 100644 --- a/qubosolver/classical_solver/classical_solver.py +++ b/qubosolver/_solvers/classical_solver.py @@ -12,20 +12,14 @@ from __future__ import annotations from abc import ABC, abstractmethod -from typing import List -import cplex import torch from qubosolver.config import ClassicalConfig from qubosolver import QUBOInstance, QUBOSolution -from qubosolver.classical_solver.simulated_annealing import qubo_simulated_annealing -from qubosolver.classical_solver.tabu_search import qubo_tabu_search -from qubosolver.classical_solver.classical_solver_conversion_tools import ( - qubo_instance_to_sparsepairs, -) -from qubosolver.qubo_types import ClassicalSolverType -from qubosolver.utils.qubo_eval import qubo_cost +from qubosolver.types.enums import ClassicalSolverType +from qubosolver.utils.qubo_eval import batched_qubo_cost +from qubosolver import solvers class BaseClassicalSolver(ABC): @@ -40,7 +34,7 @@ def __init__(self, instance: QUBOInstance, config: ClassicalConfig): Args: instance (QUBOInstance): The QUBO problem instance to solve. - config (Optional[Dict[str, Any]]): Solver configuration + config (Optional[dict[str, Any]]): Solver configuration (e.g., cplex_maxtime, cplex_log_path, classical_solver_type). """ self.instance = instance @@ -64,56 +58,12 @@ class CplexSolver(BaseClassicalSolver): """ def solve(self) -> QUBOSolution: - # Extract configuration parameters using new keys. + from qubosolver.solvers import cplex + log_path: str = self.config.cplex_log_path maxtime: float = self.config.cplex_maxtime - if self.instance.coefficients is None: - raise ValueError("The QUBO instance does not contain coefficients.") - - # Determine the number of variables. - N: int = self.instance.coefficients.shape[0] - # If there are no variables, return an empty solution. - if N == 0: - bitstring_tensor = torch.empty((0, 0), dtype=torch.float32) - cost_tensor = torch.empty((0,), dtype=torch.float32) - return QUBOSolution(bitstrings=bitstring_tensor, costs=cost_tensor) - - # Convert the coefficient matrix into CPLEX sparse pairs format using the conversion tool. - sparsepairs: List[cplex.SparsePair] = qubo_instance_to_sparsepairs(self.instance) - - # Open a log file. - log_file = open(log_path, "w") - problem = cplex.Cplex() - - # Redirect logging streams. - problem.set_log_stream(log_file) - problem.set_error_stream(log_file) - problem.set_warning_stream(log_file) - problem.set_results_stream(log_file) - - problem.parameters.timelimit.set(maxtime) - problem.objective.set_sense(problem.objective.sense.minimize) - - # Add binary variables. - problem.variables.add(types="B" * N) - - # Set the quadratic objective. - problem.objective.set_quadratic(sparsepairs) - - problem.solve() - - # Retrieve solution. - solution_values = problem.solution.get_values() - solution_cost = problem.solution.get_objective_value() - - log_file.close() - - # Convert the solution into a QUBOSolution. - bitstring_tensor = torch.tensor([[int(b) for b in solution_values]], dtype=torch.float32) - cost_tensor = torch.tensor([solution_cost], dtype=torch.float32) - - return QUBOSolution(bitstrings=bitstring_tensor, costs=cost_tensor) + return cplex(self.instance, maxtime=maxtime, log_path=log_path) class SimulatedAnnealingSolver(BaseClassicalSolver): @@ -122,7 +72,7 @@ class SimulatedAnnealingSolver(BaseClassicalSolver): """ def solve(self) -> QUBOSolution: - simulated_annealing_solution = qubo_simulated_annealing( + simulated_annealing_solution = solvers.simulated_annealing( qubo=self.instance, top_k=self.config.max_bitstrings, max_iter=self.config.max_iter, @@ -147,7 +97,7 @@ def solve(self) -> QUBOSolution: x0 = torch.randint(0, 2, size=(self.instance.size,)) else: x0 = self.config.tabu_x0 - tabu_search_solution = qubo_tabu_search( + tabu_search_solution = solvers.tabu_search( qubo=self.instance, x0=x0, max_iter=self.config.max_iter, @@ -195,7 +145,7 @@ def solve(self) -> QUBOSolution: torch.float32 ) unique_bits, counts = torch.unique(bitstrings, dim=0, return_counts=True) - costs = qubo_cost(unique_bits, self.instance.coefficients) + costs = batched_qubo_cost(unique_bits, self.instance.coefficients) return QUBOSolution( bitstrings=unique_bits, costs=costs, @@ -210,7 +160,7 @@ def get_classical_solver(instance: QUBOInstance, config: ClassicalConfig) -> Bas Args: instance (QUBOInstance): The QUBO problem instance. - config (Optional[Dict[str, Any]]): Configuration, + config (Optional[dict[str, Any]]): Configuration, possibly including 'classical_solver_type'. Returns: diff --git a/qubosolver/_solvers/cplex.py b/qubosolver/_solvers/cplex.py new file mode 100644 index 00000000..96f28ebb --- /dev/null +++ b/qubosolver/_solvers/cplex.py @@ -0,0 +1,73 @@ +from __future__ import annotations + +import cplex as CPLEX +import torch + +from qubosolver import QUBOInstance, QUBOSolution + + +def _qubo_instance_to_sparsepairs( + instance: QUBOInstance, tol: float = 1e-8 +) -> list[CPLEX.SparsePair]: + + matrix = instance.coefficients.cpu().numpy() + size = matrix.shape[0] + sparsepairs: list[CPLEX.SparsePair] = [] + + for i in range(size): + indices: list[int] = [] + values: list[float] = [] + for j in range(size): + coeff = matrix[i, j] * 2 + if abs(coeff) > tol: + indices.append(j) + values.append(float(coeff)) # <<< conversion ici + sparsepairs.append(CPLEX.SparsePair(ind=indices, val=values)) + + return sparsepairs + + +def cplex(Q: QUBOInstance, maxtime: float = 600.0, log_path: str = "solver.log") -> QUBOSolution: + + # Determine the number of variables. + N: int = Q.size + # If there are no variables, return an empty solution. + if N == 0: + return QUBOSolution() + + # Convert the coefficient matrix into CPLEX sparse pairs format using the conversion tool. + sparsepairs: list[CPLEX.SparsePair] = _qubo_instance_to_sparsepairs(Q) + + # Open a log file. + log_file = open(log_path, "w") + problem = CPLEX.Cplex() + + # Redirect logging streams. + problem.set_log_stream(log_file) + problem.set_error_stream(log_file) + problem.set_warning_stream(log_file) + problem.set_results_stream(log_file) + + problem.parameters.timelimit.set(maxtime) + problem.objective.set_sense(problem.objective.sense.minimize) + + # Add binary variables. + problem.variables.add(types="B" * N) + + # Set the quadratic objective. + problem.objective.set_quadratic(sparsepairs) + + problem.solve() + + # Retrieve solution. + solution_values = problem.solution.get_values() + solution_cost = problem.solution.get_objective_value() + + log_file.close() + + # Convert the solution into a QUBOSolution. + bitstring_tensor = torch.tensor([[int(b) for b in solution_values]], dtype=torch.int8) + counts = torch.tensor([1], dtype=torch.int64) + cost_tensor = torch.tensor([solution_cost], dtype=torch.float32) + + return QUBOSolution(bitstrings=bitstring_tensor, costs=cost_tensor, counts=counts) diff --git a/qubosolver/_solvers/quantum.py b/qubosolver/_solvers/quantum.py new file mode 100644 index 00000000..53410593 --- /dev/null +++ b/qubosolver/_solvers/quantum.py @@ -0,0 +1,41 @@ +from __future__ import annotations + +from qubosolver.types import protocols +import qoolqit +from qoolqit.execution.compilation_functions import CompilerProfile +from qoolqit.execution import job + + +def _quantum_program( + device: qoolqit.Device, + drive: qoolqit.Drive, + register: qoolqit.Register, + compiler_profile: CompilerProfile = CompilerProfile.MAX_ENERGY, +) -> qoolqit.QuantumProgram: + program = qoolqit.QuantumProgram( + register=register, + drive=drive, + ) + max_duration_ratio = 0.99 if device.specs["max_duration"] is not None else None + + program.compile_to( + device, + profile=compiler_profile, + device_max_duration_ratio=max_duration_ratio, + ) + return program + + +def analog_quantum_sample( + backend: protocols.Backend, + device: qoolqit.Device, + drive: qoolqit.Drive, + register: qoolqit.Register, + compiler_profile: CompilerProfile = CompilerProfile.MAX_ENERGY, +) -> job.Job: + """ + Submit a quantum program for execution on the configured backend.xecution is not supported on local backends. + """ + program = _quantum_program(device, drive, register, compiler_profile) + + return backend.run(program) diff --git a/qubosolver/classical_solver/simulated_annealing.py b/qubosolver/_solvers/simulated_annealing.py similarity index 80% rename from qubosolver/classical_solver/simulated_annealing.py rename to qubosolver/_solvers/simulated_annealing.py index 333b0e48..9c0d53bd 100644 --- a/qubosolver/classical_solver/simulated_annealing.py +++ b/qubosolver/_solvers/simulated_annealing.py @@ -5,57 +5,9 @@ from qubosolver import QUBOInstance, QUBOSolution -def qubo_simulated_annealing( - qubo: QUBOInstance, - top_k: int = 5, - max_iter: int = 1000, - initial_temp: float = 10.0, - final_temp: float = 1e-3, - cooling_rate: float | None = None, - seed: int | None = None, - start: torch.Tensor | None = None, - energy_tol: float = 0.0, -) -> QUBOSolution: - """ - Solve a QUBO instance using the Simulated Annealing metaheuristic. - - This function wraps the low-level `simulated_annealing()` routine - and converts its output into a standardized `QUBOSolution` object. - - The algorithm gradually lowers the system temperature to reduce - the probability of accepting worse solutions, balancing exploration - and exploitation. - - Returns: - A `QUBOSolution` object containing: - - `bitstrings`: The best binary solution(s) found. - - `costs`: Corresponding objective values as tensors. - - `counts`: Corresponding frequencies. - - `probabilities`: Frequencies divided by `top_k`. - - Example: - >>> solution = qubo_simulated_annealing(qubo) - >>> print(solution.bitstrings, solution.costs) - """ - bitstrings, costs, counts = simulated_annealing( - Q=qubo.coefficients, - top_k=top_k, - max_iter=max_iter, - initial_temp=initial_temp, - final_temp=final_temp, - cooling_rate=cooling_rate, - seed=seed, - start=start, - energy_tol=energy_tol, - ) - return QUBOSolution( - bitstrings=bitstrings, costs=costs, counts=counts, probabilities=counts.float() / top_k - ) - - @torch.no_grad() def simulated_annealing( - Q: torch.Tensor, + qubo: QUBOInstance, top_k: int = 5, max_iter: int = 1000, initial_temp: float = 5.0, @@ -64,7 +16,7 @@ def simulated_annealing( seed: int | None = None, start: torch.Tensor | None = None, energy_tol: float = 0.0, -) -> tuple[torch.Tensor, torch.Tensor, torch.Tensor]: +) -> QUBOSolution: """ Perform Simulated Annealing (SA) for a Quadratic Unconstrained Binary Optimization (QUBO) problem using PyTorch, returning the top-K unique low-energy solutions found. @@ -144,6 +96,7 @@ def simulated_annealing( if seed is not None: rng.manual_seed(seed) + Q = qubo.coefficients Q = 0.5 * (Q + Q.T) n = int(Q.shape[0]) @@ -224,4 +177,10 @@ def maybe_insert(b_u8: torch.Tensor, e: float) -> None: bitstrings, dim=0, return_inverse=True, return_counts=True ) energies = energies[inverse_indices] - return unique_bits, energies, counts + + solution = QUBOSolution( + bitstrings=unique_bits.to(torch.int8), counts=counts, costs=energies.to(torch.float32) + ) + solution.probabilities = solution.compute_probabilities() + + return solution diff --git a/qubosolver/solver.py b/qubosolver/_solvers/solver.py similarity index 63% rename from qubosolver/solver.py rename to qubosolver/_solvers/solver.py index 83cc951e..05492fb7 100644 --- a/qubosolver/solver.py +++ b/qubosolver/_solvers/solver.py @@ -1,21 +1,15 @@ from __future__ import annotations -from typing import Callable -from collections import Counter +from collections.abc import Callable import torch -import random -# Import the classical solver factory from our classical_solver module. from qoolqit import Register, Drive -from qubosolver.qubo_instance import QUBOInstance -from qubosolver.data import QUBOSolution +from qubosolver import QUBOInstance, QUBOSolution, SolverConfig, DecompositionConfig, transforms from qubosolver.classical_solver import get_classical_solver -from qubosolver.config import SolverConfig, DecompositionConfig from qubosolver.pipeline import ( BaseSolver, - Fixtures, get_embedder, get_drive_shaper, ) @@ -139,7 +133,7 @@ def solve(self) -> QUBOSolution: # 1) try trivial else delegate to quantum solver if self.config.activate_trivial_solutions: trivial = self._trivial_solution() - if trivial is not None: + if trivial: return trivial self._check_size_limit() @@ -148,37 +142,10 @@ def solve(self) -> QUBOSolution: embedding = self.embedding() - drive, qubo_solution = self.drive(embedding) + drive, solution = self.drive(embedding) - bitstrings, counts, _, _ = ( - qubo_solution.bitstrings, - qubo_solution.counts, - qubo_solution.probabilities, - qubo_solution.costs, - ) - if ( - len(bitstrings) == 0 and qubo_solution.counts is None - ) or self.config.drive_shaping.optimized_re_execute_opt_drive: - bitstrings, counts = self.execute(drive, embedding) - - bitstring_strs = bitstrings - bitstrings_tensor = torch.tensor( - [list(map(int, bs)) for bs in bitstring_strs], dtype=torch.float32 - ) - if counts is None: - counts_tensor = torch.empty((0,), dtype=torch.int32) - elif isinstance(counts, dict) or isinstance(counts, Counter): - count_values = [counts.get(bs, 0) for bs in bitstring_strs] - counts_tensor = torch.tensor(count_values, dtype=torch.int32) - else: - counts_tensor = counts - - solution = QUBOSolution( - bitstrings=bitstrings_tensor, - counts=counts_tensor, - costs=torch.Tensor(), - probabilities=None, - ) + if not solution or self.config.drive_shaping.optimized_re_execute_opt_drive: + solution = self.execute(drive, embedding) # Post-process fixations of the preprocessing and restore the original QUBO solution = self.post_process_fixation(solution) @@ -205,8 +172,6 @@ class QuboSolverClassical(BaseSolver): def __init__(self, instance: QUBOInstance, config: SolverConfig | None = None): super().__init__(instance, config) - # Optionally, you could instantiate Fixtures here for postprocessing: - self.fixtures = Fixtures(self.instance, self.config) def embedding(self) -> Register: # Classical solvers do not require an embedding. @@ -219,7 +184,7 @@ def solve(self) -> QUBOSolution: # 1) try trivial if self.config.activate_trivial_solutions: trivial = self._trivial_solution() - if trivial is not None: + if trivial: return trivial self.preprocess() @@ -339,109 +304,47 @@ def solve(self) -> QUBOSolution: return solver.solve() else: - from qubosolver.algorithms.decompose import ( - compute_distance_interaction_matrix, - geometric_search, - interaction_matrix_from_placed, - last_target_matrix, - transfer_edge_values, - update_global_solution, - vertices_to_place, - positive_vertices_update, + config = transforms.decompositions.Config.from_decomposition_config( + self.decomposition_config ) - - self._decomposition = [] - - global_solution = torch.full((self.instance.size,), -1) - qubo_mat = self.instance.coefficients.clone() - dist_matrix = compute_distance_interaction_matrix( - self.device._pulser_device, - qubo_mat, - neglecting_inter_distance=self.decomposition_config.neglecting_inter_distance, - neglecting_max_coefficient=self.decomposition_config.neglecting_max_coefficient, + decomposed_qubo = transforms.decompositions.QUBOInstance( + self.instance, self.device, config ) - # The following dictionary contain vertices to apply the decomposition search - # where each vertex key maps to other blocking, separated and neighbors vertices - # and gets updated as we iterate in the decomposition - dict_vertices_to_place = vertices_to_place( - dist_matrix, - qubo_mat, - separation_threshold=self.decomposition_config.neglecting_inter_distance, - ) + while len(decomposed_qubo._vertices_to_place) > config.decompose_stop_number: - transfer_edge_values(dict_vertices_to_place, dict(), global_solution, qubo_mat) - positive_vertices = positive_vertices_update(dict_vertices_to_place, global_solution) - for v in positive_vertices: - self._decomposition.append([v]) - while len(dict_vertices_to_place) > self.decomposition_config.decompose_stop_number: - # find a first vertex to start the geometric search - # random works better according to some performed numerics - # sort to have reproducibility when setting the seed - first_vertex_search = random.choice(sorted(dict_vertices_to_place.keys())) - placed_vertices = geometric_search( - qubo_mat, - dict_vertices_to_place, - first_vertex_search, - self.decomposition_config.decompose_threshold, - self.device._pulser_device, + subqubo = transforms.decompositions.extract_subqubo( + decomposed_qubo, self.device, config ) - if len(placed_vertices) <= self.decomposition_config.decompose_break_placement: + + if subqubo.size == 0: break self.number_iterations += 1 - matrix_to_solve, map_index_vertices = interaction_matrix_from_placed( - placed_vertices, self.device._pulser_device - ) - qubo = QUBOInstance(matrix_to_solve) - subsolver = self._solver_factory(qubo, self._config_subproblems) + subsolver = self._solver_factory(subqubo, self._config_subproblems) # only one bitstring is kept as per design choice of the # decomposition algorithm - sub_solution = subsolver.solve().bitstrings[0] - self._decomposition.append(list(map_index_vertices.keys())) - update_global_solution( - global_solution=global_solution, - sub_solution=sub_solution, - mapping=map_index_vertices, - ) - - transfer_edge_values( - dict_vertices_to_place, - placed_vertices, - global_solution, - qubo_mat, - ) - positive_vertices = positive_vertices_update( - dict_vertices_to_place, global_solution - ) - for v in positive_vertices: - self._decomposition.append([v]) + subqubo.solution = subsolver.solve() + subqubo.solution.costs = subqubo.solution.compute_costs(subqubo) + subqubo.solution.sort_by_cost() + transforms.decompositions.update(decomposed_qubo, subqubo) # classical resolution of last matrix - matrix_to_solve, mapping_target_vertices = last_target_matrix( - list(dict_vertices_to_place.keys()), qubo_mat - ) - qubo = QUBOInstance(matrix_to_solve) - lastsolver = QuboSolverClassical( - qubo, - SolverConfig(use_quantum=False, decompose=None), - ) - sub_solution = lastsolver.solve().bitstrings[0] - if mapping_target_vertices: - self._decomposition.append(list(mapping_target_vertices.keys())) - update_global_solution( - global_solution=global_solution, - sub_solution=sub_solution, - mapping=mapping_target_vertices, + subqubo = transforms.decompositions.extract_subqubo( + decomposed_qubo, self.device, config, last=True ) + if subqubo.size != 0: + lastsolver = QuboSolverClassical( + subqubo, + SolverConfig(use_quantum=False, decompose=None), + ) + subqubo.solution = lastsolver.solve() + subqubo.solution.costs = subqubo.solution.compute_costs(subqubo) + subqubo.solution.sort_by_cost() + transforms.decompositions.update(decomposed_qubo, subqubo) + + self._decomposition = decomposed_qubo._decomposition + # Probabilities and counts are ignored as we return one solution - qubosol = QUBOSolution( - bitstrings=global_solution.unsqueeze(0).to(dtype=torch.float32), - counts=torch.tensor([1], dtype=torch.int), - costs=torch.Tensor(), - probabilities=None, - ) - qubosol.costs = qubosol.compute_costs(self.instance) - qubosol.bitstrings = qubosol.bitstrings.int() - return qubosol + return decomposed_qubo.solution diff --git a/qubosolver/classical_solver/tabu_search.py b/qubosolver/_solvers/tabu_search.py similarity index 65% rename from qubosolver/classical_solver/tabu_search.py rename to qubosolver/_solvers/tabu_search.py index 0c08d7a5..68a4fd8d 100644 --- a/qubosolver/classical_solver/tabu_search.py +++ b/qubosolver/_solvers/tabu_search.py @@ -3,68 +3,18 @@ import torch from qubosolver import QUBOInstance, QUBOSolution -from qubosolver.utils.qubo_eval import qubo_cost - - -def qubo_tabu_search( - qubo: QUBOInstance, - x0: torch.Tensor, - max_iter: int = 100, - tabu_tenure: int = 7, - max_no_improve: int = 20, - max_bitstrings: int = 1, -) -> QUBOSolution: - """ - Solve a QUBO problem using a simple Tabu Search heuristic. - - This function wraps the core `tabu_search()` routine and converts - its output into a standardized `QUBOSolution` object, including - the bitstrings and their evaluated costs. - - Args: - qubo: The QUBO instance to optimize, providing the cost matrix - and an evaluation method. - x0: The initial solution as a binary tensor of shape (n,). - max_iter: Maximum number of iterations to perform. - tabu_tenure: Number of iterations a flipped variable remains tabu. - max_no_improve: Stop criterion based on consecutive iterations - without improvement. - - Returns: - A `QUBOSolution` object containing: - - `bitstrings`: The best solution found as a tensor. - - `costs`: The corresponding objective value tensor. - - `counts`: The frequencies of each bitstring. - - `probabilities`: Frequencies divided by max_bitstrings. - - Example: - >>> solution = qubo_tabu_search(qubo, x0=torch.randint(0, 2, (10,))) - >>> print(solution) - """ - best_solutions, costs, counts = tabu_search( - qubo=qubo, - x0=x0, - max_iter=max_iter, - tabu_tenure=tabu_tenure, - max_no_improve=max_no_improve, - max_bitstrings=max_bitstrings, - ) - return QUBOSolution( - bitstrings=best_solutions, - costs=costs, - counts=counts, - probabilities=counts.float() / max_bitstrings, - ) +from qubosolver.utils.qubo_eval import batched_qubo_cost +from qubosolver.types import Bitstring def tabu_search( qubo: QUBOInstance, - x0: torch.Tensor, + x0: Bitstring, max_iter: int = 100, tabu_tenure: int = 7, max_no_improve: int = 20, max_bitstrings: int = 1, -) -> tuple[torch.Tensor, torch.Tensor, torch.Tensor]: +) -> QUBOSolution: """ Perform Tabu Search on a QUBO instance to find low-cost bitstrings. @@ -101,7 +51,7 @@ def tabu_search( # Repeat x0 for each parallel run x_current = x0.clone().to(torch.int64).unsqueeze(0).repeat(max_bitstrings, 1) - f_current = qubo_cost(x_current, Q) + f_current = batched_qubo_cost(x_current, Q) x_best = x_current.clone() f_best = f_current.clone() @@ -113,7 +63,7 @@ def tabu_search( flips = torch.eye(n, dtype=torch.int64, device=device).unsqueeze(0) x_neighbors = x_current.unsqueeze(1).clone() x_neighbors = (x_neighbors + flips) % 2 # each bit flipped - f_candidates = qubo_cost(x_neighbors.view(-1, n), Q).view(max_bitstrings, n) + f_candidates = batched_qubo_cost(x_neighbors.view(-1, n), Q).view(max_bitstrings, n) # Tabu and aspiration tabu_mask = tabu_list > iteration @@ -144,6 +94,13 @@ def tabu_search( # Get unique final solutions uniq, counts = torch.unique(x_best, dim=0, return_counts=True) - costs = qubo_cost(uniq, Q) + costs = batched_qubo_cost(uniq, Q) + + solution = QUBOSolution( + bitstrings=uniq.to(torch.int8), + costs=costs, + counts=counts, + ) + solution.probabilities = solution.compute_probabilities() - return uniq, costs, counts + return solution diff --git a/qubosolver/_solvers/trivial.py b/qubosolver/_solvers/trivial.py new file mode 100644 index 00000000..c65edad0 --- /dev/null +++ b/qubosolver/_solvers/trivial.py @@ -0,0 +1,62 @@ +from __future__ import annotations + +import torch + +from qubosolver import QUBOInstance, QUBOSolution, SolutionStatusType + + +def trivial_solution_search(Q: QUBOInstance) -> QUBOSolution: + """ + Check for the two trivial QUBO cases: + 1) all coefficients >= 0 → solution = 0^n + 2) all coefficients <= 0 → solution = 1^n + 3) diagonal qubo, negative coeffs gets 1, positive gets 0 + + Returns: + QUBOSolution if a trivial case applies, else None. + """ + coeffs = Q.coefficients + n = Q.size + device, dtype = Q.device, Q.dtype + + # Case 1: all coeffs >= 0 → x = [0,...,0] + if torch.all(coeffs >= 0): + raw = torch.zeros(n, dtype=torch.int64, device=device) + # always make a batch of one: shape (1, n) + batch = raw.unsqueeze(0) + cost = Q.evaluate_solution(raw) + return QUBOSolution( + bitstrings=batch, + counts=torch.tensor([1], dtype=torch.int64, device=device), + costs=torch.tensor([cost], dtype=dtype, device=device), + solution_status=SolutionStatusType.TRIVIALZERO, + ) + + # Case 2: all coeffs <= 0 → x = [1,...,1] + if torch.all(coeffs <= 0): + raw = torch.ones(n, dtype=torch.int64, device=device) + # always make a batch of one: shape (1, n) + batch = raw.unsqueeze(0) + cost = Q.evaluate_solution(raw) + return QUBOSolution( + bitstrings=batch, + counts=torch.tensor([1], dtype=torch.int64, device=device), + costs=torch.tensor([cost], dtype=dtype, device=device), + solution_status=SolutionStatusType.TRIVIALONE, + ) + + # Case 3: diagonal cases + # negative coeffs gets 1, positive gets 0 + diagonal = torch.diag(coeffs) + if (torch.diag(diagonal) == coeffs).all(): + raw = (diagonal < 0).long() + cost = Q.evaluate_solution(raw) + batch = raw.unsqueeze(0) + return QUBOSolution( + bitstrings=batch, + counts=torch.tensor([1], dtype=torch.int64, device=device), + costs=torch.tensor([cost], dtype=dtype, device=device), + solution_status=SolutionStatusType.TRIVIALDIAGONAL, + ) + + return QUBOSolution() diff --git a/qubosolver/_transforms/__init__.py b/qubosolver/_transforms/__init__.py new file mode 100644 index 00000000..ba9c591a --- /dev/null +++ b/qubosolver/_transforms/__init__.py @@ -0,0 +1,4 @@ +from __future__ import annotations + +from qubosolver.transforms import variable_fixing as variable_fixing +from qubosolver.transforms import decompositions as decompositions diff --git a/qubosolver/algorithms/__init__.py b/qubosolver/_transforms/algorithms/__init__.py similarity index 100% rename from qubosolver/algorithms/__init__.py rename to qubosolver/_transforms/algorithms/__init__.py diff --git a/qubosolver/algorithms/decompose.py b/qubosolver/_transforms/algorithms/decompose.py similarity index 99% rename from qubosolver/algorithms/decompose.py rename to qubosolver/_transforms/algorithms/decompose.py index 473a6960..d6e0b658 100644 --- a/qubosolver/algorithms/decompose.py +++ b/qubosolver/_transforms/algorithms/decompose.py @@ -2,7 +2,7 @@ import math import random -from typing import TypedDict, List +from typing import TypedDict import numpy as np from pulser.devices._device_datacls import BaseDevice as PulserBaseDevice @@ -241,7 +241,7 @@ def update_vertex_info_from_placed( def positive_vertices_update( dict_vertices_to_place: dict[int, VertexToPlace], global_solution: torch.Tensor -) -> List[int]: +) -> list[int]: """Remove vertices whose weight became positive during decomposition. Args: diff --git a/qubosolver/_transforms/decompositions.py b/qubosolver/_transforms/decompositions.py new file mode 100644 index 00000000..4f1c3042 --- /dev/null +++ b/qubosolver/_transforms/decompositions.py @@ -0,0 +1,158 @@ +from __future__ import annotations + +import copy +from dataclasses import dataclass +import random +import torch + +import qoolqit + +from qubosolver import QUBOSolution, DecompositionConfig +from qubosolver import QUBOInstance as QUBOInstanceBase +from qubosolver.transforms.algorithms.decompose import ( + compute_distance_interaction_matrix, + geometric_search, + interaction_matrix_from_placed, + last_target_matrix, + transfer_edge_values, + update_global_solution, + vertices_to_place, + positive_vertices_update, + VertexToPlace, + WeightedZone, +) +from qubosolver.types import Matrix, make_matrix + + +@dataclass +class Config: + neglecting_inter_distance: float = 15.0 + neglecting_max_coefficient: float = 1.0 + decompose_stop_number: int = 15 + decompose_threshold: float = 25.0 + decompose_break_placement: int = 3 + + @staticmethod + def from_decomposition_config(config: DecompositionConfig) -> Config: + return Config( + neglecting_inter_distance=config.neglecting_inter_distance, + neglecting_max_coefficient=config.neglecting_max_coefficient, + decompose_stop_number=config.decompose_stop_number, + decompose_threshold=config.decompose_threshold, + decompose_break_placement=config.decompose_break_placement, + ) + + +class QUBOInstance(QUBOInstanceBase): + + def __init__( + self, + parent_instance: QUBOInstanceBase, + device: qoolqit.Device, + config: Config = Config(), + ): + + super().__init__(parent_instance.coefficients) + self._parent_instance = copy.deepcopy(parent_instance) + + self._global_solution: torch.Tensor = torch.full( + (parent_instance.size,), -1, dtype=torch.int64 + ) + self._qubo_matrix: torch.Tensor = parent_instance.coefficients.clone() + self._vertices_to_place: dict[int, VertexToPlace] = {} + self._placed_vertices: dict[int, WeightedZone] = {} + self._decomposition: list[list[int]] = [] + + dist_matrix = compute_distance_interaction_matrix( + device._pulser_device, + self._qubo_matrix, + neglecting_inter_distance=config.neglecting_inter_distance, + neglecting_max_coefficient=config.neglecting_max_coefficient, + ) + # The following dictionary contain vertices to apply the decomposition search + # where each vertex key maps to other blocking, separated and neighbors vertices + # and gets updated as we iterate in the decomposition + self._vertices_to_place = vertices_to_place( + dist_matrix, + self._qubo_matrix, + separation_threshold=config.neglecting_inter_distance, + ) + + update(self, SubQUBOInstance()) + + +class SubQUBOInstance(QUBOInstanceBase): + def __init__( + self, + coefficients: Matrix = make_matrix(0), + map_index_vertices: dict[int, int] = {}, + ): + super().__init__(coefficients) + self._map_index_vertices = map_index_vertices + + +def extract_subqubo( + qubo: QUBOInstance, + device: qoolqit.Device, + config: Config, + last: bool = False, +) -> SubQUBOInstance: + + if last: + matrix_to_solve, map_index_vertices = last_target_matrix( + list(qubo._vertices_to_place.keys()), + qubo._qubo_matrix, + ) + # qubo._vertices_to_place = {} + return SubQUBOInstance(matrix_to_solve, map_index_vertices) + + # find a first vertex to start the geometric search + # random works better according to some performed numerics + # sort to have reproducibility when setting the seed + first_vertex_search = random.choice(sorted(qubo._vertices_to_place.keys())) + + qubo._placed_vertices = geometric_search( + qubo._qubo_matrix, + qubo._vertices_to_place, + first_vertex_search, + config.decompose_threshold, + device._pulser_device, + ) + if len(qubo._placed_vertices) <= config.decompose_break_placement: + return SubQUBOInstance() + + matrix_to_solve, map_index_vertices = interaction_matrix_from_placed( + qubo._placed_vertices, device._pulser_device + ) + return SubQUBOInstance(matrix_to_solve, map_index_vertices) + + +def update(qubo: QUBOInstance, subqubo: SubQUBOInstance) -> None: + if subqubo._map_index_vertices: + qubo._decomposition.append(list(subqubo._map_index_vertices.keys())) + if subqubo.solution: + update_global_solution( + global_solution=qubo._global_solution, + sub_solution=subqubo.solution.bitstrings[0], + mapping=subqubo._map_index_vertices, + ) + + transfer_edge_values( + qubo._vertices_to_place, + qubo._placed_vertices, + qubo._global_solution, + qubo._qubo_matrix, + ) + + positive_vertices = positive_vertices_update(qubo._vertices_to_place, qubo._global_solution) + + for v in positive_vertices: + qubo._decomposition.append([v]) + + if (qubo._global_solution != -1).all(): + # Probabilities and counts are ignored as we return one solution + qubo.solution = QUBOSolution( + bitstrings=qubo._global_solution.unsqueeze(0), + counts=torch.tensor([1], dtype=torch.int), + ) + qubo.solution.costs = qubo.solution.compute_costs(qubo) diff --git a/qubosolver/_transforms/variable_fixing.py b/qubosolver/_transforms/variable_fixing.py new file mode 100644 index 00000000..3a70b7ce --- /dev/null +++ b/qubosolver/_transforms/variable_fixing.py @@ -0,0 +1,279 @@ +from __future__ import annotations + +from collections.abc import Callable, Sequence +from typing import cast, TypeAlias + +import copy +import json +import torch + +from qubosolver import QUBOSolution +from qubosolver.io import utils as io_utils +from qubosolver import QUBOInstance as QUBOInstanceBase + +#  TODO: Using `type` statement when Python >= 3.12 +Rule: TypeAlias = Callable[[QUBOInstanceBase], dict[int, int]] + + +def hansen_fixing(qubo: QUBOInstanceBase) -> dict[int, int]: + """ + Identifies and fixes variables in a QUBO instance based on threshold conditions. + + This method determines whether a variable should be fixed to 0 or 1 by computing + lower and upper bounds from the diagonal and off-diagonal elements of the QUBO matrix. + + Args: + qubo (QUBOInstance): The QUBO instance containing the coefficients matrix. + + Returns: + dict[int, int]: A dictionary mapping variable indices to fixed values (0 or 1). + """ + if qubo.coefficients is None: + raise ValueError("QUBO coefficients are not initialized.") + + fixed_dict: dict[int, int] = {} + size_raw = qubo.size + size: int = cast(int, size_raw) + epsilon: float = 1e-8 # Tolerance to avoid floating-point precision issues + + for i in range(size): + ci = qubo.coefficients[i, i].item() # Diagonal element + + q_minus = sum(min(0, qubo.coefficients[i, j].item()) for j in range(size) if j != i) + q_plus = sum(max(0, qubo.coefficients[i, j].item()) for j in range(size) if j != i) + + if ci + q_minus * 2 >= -epsilon: + fixed_dict[i] = 0 + elif ci + q_plus * 2 <= epsilon: + fixed_dict[i] = 1 + + return fixed_dict + + +class QUBOInstance(QUBOInstanceBase): + + def __init__(self, parent_instance: QUBOInstanceBase): + super().__init__( + parent_instance.coefficients, + ) + self._parent_instance = copy.deepcopy(parent_instance) + self._fixed_indices: list[dict[int, int]] = [] + + @property + def n_fixed_indices(self) -> int: + """Returns the number of fixed variables. + + Returns: + int: The number of fixed variables. + """ + return sum([len(fixed) for fixed in self._fixed_indices]) + + @staticmethod + def save(file_like: io_utils.FileLike[bytes], instance: QUBOInstanceBase) -> None: + """ + Saves a QUBOInstance to a file-like object. + + Args: + file_like (io_utils.FileLike[bytes]): + File-like object opened in binary write mode where the instance will be saved. + instance (QUBOInstance): + The QUBOInstance object to be saved. + + Returns: + None + """ + _check_QUBOInstance(instance) + assert isinstance(instance, QUBOInstance) + + with io_utils.open(file_like, "wb") as f: + QUBOInstanceBase.save(f, instance) + QUBOInstanceBase.save(f, instance._parent_instance) + fixed_var_json = json.dumps(instance._fixed_indices) + io_utils.save_string(f, fixed_var_json) + + @staticmethod + def load(file_like: io_utils.FileLike[bytes]) -> QUBOInstance: + """ + Loads a QUBOInstance from a file-like object. + + Args: + file_like (io_utils.FileLike[bytes]): + File-like object opened in binary read mode containing the saved QUBOInstance data. + + Returns: + QUBOInstance: + A new QUBOInstance object reconstructed from the saved data. + """ + + def decode_int_keys(obj: dict) -> dict: + return {int(k): v for k, v in obj.items()} + + with io_utils.open(file_like, "rb") as f: + instance = QUBOInstance(QUBOInstanceBase.load(f)) + instance._parent_instance = QUBOInstanceBase.load(f) + fixed_var_json = io_utils.load_string(f) + instance._fixed_indices = json.loads(fixed_var_json, object_hook=decode_int_keys) + + return instance + + +def _check_QUBOInstance(qubo: QUBOInstanceBase) -> None: + if not isinstance(qubo, QUBOInstance): + raise TypeError("Input must be an instance of _QUBOInstance.") + + +def default_rules() -> tuple[Rule]: + return (hansen_fixing,) + + +def reduce_qubo( + qubo: QUBOInstanceBase, fixed_indices: dict[int, int], *, inplace: bool = False +) -> QUBOInstance: + """ + Applies variable fixation to reduce the size of the QUBO problem. + + This function modifies the QUBO coefficient matrix by: + - Removing rows and columns corresponding to fixed variables. + - Adjusting diagonal elements to account for fixed variables. + + Args: + fixed_dict (dict[int, int]): A dictionary of fixed variable assignments. + - Keys are variable indices. + - Values are fixed binary values (0 or 1). + + Returns: + None: Modifies `self.reduced_qubo` in place. + """ + if not inplace: + qubo = QUBOInstance(qubo) + + _check_QUBOInstance(qubo) + assert isinstance(qubo, QUBOInstance) + + if not fixed_indices: + return qubo + + Q = qubo.coefficients.clone() + + fixed_to_0 = {i for i, v in fixed_indices.items() if v == 0} + fixed_to_1 = {i for i, v in fixed_indices.items() if v == 1} + fixed_vars = sorted(fixed_to_0 | fixed_to_1, reverse=True) + + for i in fixed_vars: + if i >= Q.shape[0]: + continue + + if i in fixed_to_1: + for j in range(Q.shape[0]): + if j != i: + Q[j, j] += Q[i, j] * 2 + + Q = torch.cat((Q[:i, :], Q[i + 1 :, :]), dim=0) + Q = torch.cat((Q[:, :i], Q[:, i + 1 :]), dim=1) + + qubo.coefficients = Q + qubo.update_metrics() + + qubo._fixed_indices.append(fixed_indices) + return qubo + + +def apply( + qubo: QUBOInstanceBase, + fixation_rules: Sequence[Rule] = default_rules(), + *, + inplace: bool = False, +) -> QUBOInstance: + """ + Applies a sequence of variable fixation rules to the QUBO instance. + + Args: + qubo (QUBOInstance): The QUBO instance to apply rules to. + fixation_rules (Sequence[FixationRule]): A sequence of functions that + return dictionaries mapping variable indices to fixed values. + + Returns: + list[dict[int, int]]: A list of fixation dictionaries, one per rule that fixed variables. + """ + if not inplace: + qubo = QUBOInstance(qubo) + + _check_QUBOInstance(qubo) + assert isinstance(qubo, QUBOInstance) + + for rule in fixation_rules: + fixed = rule(qubo) + reduce_qubo(qubo, fixed, inplace=True) + + return qubo + + +def apply_recursively( + qubo: QUBOInstanceBase, + fixation_rules: Sequence[Rule] = default_rules(), + *, + inplace: bool = False, +) -> QUBOInstance: + """ + Iteratively applies all fixation rules until no more variables can be fixed. + + This function repeatedly applies all rules in `self.fixation_rule_list` + until no further reduction is possible. + """ + if not inplace: + qubo = QUBOInstance(qubo) + + _check_QUBOInstance(qubo) + assert isinstance(qubo, QUBOInstance) + + while True: + prev_n_fixations = len(qubo._fixed_indices) + apply(qubo, fixation_rules, inplace=True) + n_fixations = len(qubo._fixed_indices) + assert n_fixations >= prev_n_fixations + if n_fixations == prev_n_fixations: + return qubo + + +def unapply(reduced_qubo: QUBOInstance) -> QUBOSolution: + """ + Restores fixed variables in the solution bitstrings after QUBO reduction. + + This method reconstructs the full-length bitstrings by reinserting the fixed + variables at their original positions. + + Args: + solution (QUBOSolution): The solution object from the reduced QUBO problem. + + Returns: + QUBOSolution: A solution object with bitstrings restored to their original size. + """ + # FIXME: raise if empty solution ? + bitstring_list = reduced_qubo.solution.bitstrings.tolist() or [[]] + + def reinsert_fixed_variables(bitstring: list[int]) -> list[int]: + for fixation_dict in reversed(reduced_qubo._fixed_indices): + for position, bit_value in sorted(fixation_dict.items()): + bitstring.insert(position, bit_value) + return bitstring + + bits_to_reinsert = sum(len(fixation_dict) for fixation_dict in reduced_qubo._fixed_indices) + assert (bits_to_reinsert + len(bitstring_list[0])) == reduced_qubo._parent_instance.size + + if bits_to_reinsert == 0: + return reduced_qubo.solution + + bitstring_list = [reinsert_fixed_variables(bitstring) for bitstring in bitstring_list] + bitstrings = torch.tensor(bitstring_list, dtype=torch.int8) + + costs = torch.tensor( + [reduced_qubo._parent_instance.evaluate_solution(sample.tolist()) for sample in bitstrings], + dtype=torch.float32, + ) + + return QUBOSolution( + bitstrings=bitstrings, + costs=costs, + counts=reduced_qubo.solution.counts, + probabilities=reduced_qubo.solution.probabilities, + ) diff --git a/qubosolver/_utils/__init__.py b/qubosolver/_utils/__init__.py new file mode 100644 index 00000000..80b57420 --- /dev/null +++ b/qubosolver/_utils/__init__.py @@ -0,0 +1,8 @@ +from __future__ import annotations + + +from qubosolver._utils import costs + +__all__ = [ + "costs", +] diff --git a/qubosolver/utils/qubo_eval.py b/qubosolver/_utils/costs.py similarity index 77% rename from qubosolver/utils/qubo_eval.py rename to qubosolver/_utils/costs.py index 8e694c73..a65aaaa5 100644 --- a/qubosolver/utils/qubo_eval.py +++ b/qubosolver/_utils/costs.py @@ -1,9 +1,10 @@ from __future__ import annotations import torch +from qubosolver.types import Vector, Matrix, Bitstring -def qubo_cost(x: torch.Tensor, Q: torch.Tensor) -> torch.Tensor: +def batched_qubo_cost(x: Matrix, Q: Matrix) -> Matrix: """ Compute the quadratic cost of a given binary vector under a QUBO matrix. @@ -29,6 +30,15 @@ def qubo_cost(x: torch.Tensor, Q: torch.Tensor) -> torch.Tensor: return torch.einsum("bi,ij,bj->b", x, Q, x) +def qubo_cost(x: Vector, Q: Matrix) -> float: + return batched_qubo_cost(x, Q).item() + + +def qubo_bitstring_cost(b: Bitstring, Q: Matrix) -> float: + dtype = Q.dtype + return qubo_cost(b.to(dtype), Q) + + def calculate_qubo_cost(bitstring: str, QUBO: torch.Tensor) -> float: """Apply the default qubo evaluation b Q b^T. diff --git a/qubosolver/classical_solver/__init__.py b/qubosolver/classical_solver/__init__.py deleted file mode 100644 index 0ce10c9f..00000000 --- a/qubosolver/classical_solver/__init__.py +++ /dev/null @@ -1,9 +0,0 @@ -# qubosolver/classical_solver/__init__.py - -from __future__ import annotations - -from .classical_solver import get_classical_solver - -__all__ = [ - "get_classical_solver", -] diff --git a/qubosolver/classical_solver/classical_solver_conversion_tools.py b/qubosolver/classical_solver/classical_solver_conversion_tools.py deleted file mode 100644 index 04d3b7a0..00000000 --- a/qubosolver/classical_solver/classical_solver_conversion_tools.py +++ /dev/null @@ -1,29 +0,0 @@ -from __future__ import annotations -from typing import List - -import cplex - -from qubosolver import QUBOInstance - - -def qubo_instance_to_sparsepairs( - instance: QUBOInstance, tol: float = 1e-8 -) -> List[cplex.SparsePair]: - if instance.coefficients is None: - raise ValueError("The QUBO instance does not have coefficients.") - - matrix = instance.coefficients.cpu().numpy() - size = matrix.shape[0] - sparsepairs: List[cplex.SparsePair] = [] - - for i in range(size): - indices: List[int] = [] - values: List[float] = [] - for j in range(size): - coeff = matrix[i, j] * 2 - if abs(coeff) > tol: - indices.append(j) - values.append(float(coeff)) # <<< conversion ici - sparsepairs.append(cplex.SparsePair(ind=indices, val=values)) - - return sparsepairs diff --git a/qubosolver/concepts/backend.py b/qubosolver/concepts/backend.py deleted file mode 100644 index 6fe1ef57..00000000 --- a/qubosolver/concepts/backend.py +++ /dev/null @@ -1,15 +0,0 @@ -from __future__ import annotations - -from typing import Protocol, TYPE_CHECKING - -if TYPE_CHECKING: - from typing import Self - -from qoolqit import QuantumProgram -from qoolqit.execution import job -from pulser.backend import Results - - -class Backend(Protocol): - - def run(self: Self, program: QuantumProgram) -> job.Job[Results]: ... diff --git a/qubosolver/config.py b/qubosolver/config.py index 015b1ecc..191cef6b 100644 --- a/qubosolver/config.py +++ b/qubosolver/config.py @@ -2,8 +2,9 @@ import inspect from abc import ABC +from collections.abc import Callable from dataclasses import field -from typing import Any, Callable +from typing import Any import torch from pydantic import BaseModel, ConfigDict, field_validator, model_validator, model_serializer @@ -14,12 +15,13 @@ from qoolqit.execution.compilation_functions import CompilerProfile from pulser_pasqal import PasqalCloud -from qubosolver.qubo_types import ( +from .types._enums import ( EmbedderType, LayoutType, DriveType, ClassicalSolverType, ) +from .types import Bitstring, Matrix, QUBOSolution # to handle torch Tensor BaseModel.model_config["arbitrary_types_allowed"] = True @@ -266,9 +268,9 @@ class DriveShapingConfig(Config): after optimization. Defaults to False. optimized_n_calls (int, optional): Number of calls for the optimization process. Defaults to 20. Note the optimizer accepts a minimal value of 12. - optimized_initial_omega_parameters (List[float], optional): Default initial omega parameters + optimized_initial_omega_parameters (list[float], optional): Default initial omega parameters for the drive. Defaults to Omega = (1, 2, 1). - optimized_initial_detuning_parameters (List[float], optional): Default initial detuning parameters + optimized_initial_detuning_parameters (list[float], optional): Default initial detuning parameters for the drive. Defaults to delta = (-2, 0, 2). optimized_custom_qubo_cost (Callable[[str, torch.Tensor], float], optional): Apply a different qubo cost evaluation @@ -277,13 +279,13 @@ class DriveShapingConfig(Config): Must be defined as: `def optimized_custom_qubo_cost(bitstring: str, QUBO: torch.Tensor) -> float`. Defaults to None, meaning we use the default QUBO evaluation. - optimized_custom_objective_fn (Callable[[list, list, list, list, float, str], float], optional): + optimized_custom_objective (Callable[[list, list, list, list, float, str], float], optional): For bayesian optimization, one can change the output of `qubosolver/pipeline/drive.py:OptimizedDriveShaper.run_simulation` to optimize differently. Instead of using the best cost out of the samples, one can change the objective for an average, or any function out of the form - `cost_eval = optimized_custom_objective_fn(bitstrings, + `cost_eval = optimized_custom_objective(bitstrings, counts, probabilities, costs, best_cost, best_bitstring)` Defaults to None, which means we optimize using the best cost out of the samples. @@ -297,7 +299,6 @@ class DriveShapingConfig(Config): drive_shaping_method: Any = DriveType.HEURISTIC dmm: bool = True - optimized_re_execute_opt_drive: bool = False optimized_n_calls: int = 20 optimized_initial_omega_parameters: list[float] = field( default_factory=lambda: [ @@ -313,10 +314,11 @@ class DriveShapingConfig(Config): 0.8, ] ) # ---> default initial drive parameters: delta = (-2, 0, 2) - optimized_custom_qubo_cost: Callable[[str, torch.Tensor], float] | None = None - optimized_custom_objective: Callable[[list, list, list, list, float, str], float] | None = None + optimized_custom_qubo_cost: Callable[[Bitstring, Matrix], float] | None = None + optimized_custom_objective: Callable[[QUBOSolution], float] | None = None optimized_callback_objective: Callable[..., None] | None = None optimized_seed: int | None = None + optimized_re_execute_opt_drive: bool = False # Heuristic coefficient for omega heuristic_kappa: float = 0.5 @@ -526,9 +528,3 @@ def compiler_profile(config: SolverConfig) -> CompilerProfile: if config.drive_shaping.drive_shaping_method == DriveType.OPTIMIZED: return CompilerProfile.WORKING_POINT return CompilerProfile.MAX_ENERGY - - -def max_duration_ratio(config: SolverConfig) -> float | None: - if config.device.specs["max_duration"] is None: - return None - return 0.99 diff --git a/qubosolver/data_utils.py b/qubosolver/data_utils.py index 790df91a..25a5bc53 100644 --- a/qubosolver/data_utils.py +++ b/qubosolver/data_utils.py @@ -41,7 +41,7 @@ def from_dict(data: dict, device: str = "cpu", dtype: torch.dtype = torch.float3 return tensor raise TypeError( "Unsupported dictionary format: Expected dict[int, float] for 1D tensors " - "or dict[Tuple[int, int], float] for 2D tensors." + "or dict[tuple[int, int], float] for 2D tensors." ) @@ -71,51 +71,3 @@ def from_tensor( data: torch.Tensor, device: str = "cpu", dtype: torch.dtype = torch.float32 ) -> torch.Tensor: return data.to(dtype=dtype, device=device) - - -def generate_symmetric_mask( - size: int, target: int, device: str, generator: torch.Generator -) -> torch.Tensor: - """Generate a symmetric boolean mask with an exact number of True elements - to match a certain density of QUBO. - Used in the `from_random` method of `QUBODataset`. - - Args: - size (int): Size of problem. - target (int): Target number of elements. - device (str): Torch device. - generator (torch.Generator): generator for randomness. - - Returns: - torch.Tensor: Mask. - """ - possible_x = [] - for x in range(1, min(size, target) + 1): - if (target - x) % 2 == 0: - y = (target - x) // 2 - if y <= (size * (size - 1)) // 2: - possible_x.append(x) - if not possible_x: - x, y = 1, 0 - else: - x = possible_x[ - int(torch.randint(0, len(possible_x), (1,), device=device, generator=generator).item()) - ] - y = (target - x) // 2 - - mask = torch.zeros((size, size), dtype=torch.bool, device=device) - diag_indices = torch.randperm(size, device=device, generator=generator)[:x] - for i in diag_indices.tolist(): - mask[i, i] = True - - upper_indices = torch.tensor( - [(i, j) for i in range(size) for j in range(i + 1, size)], - device=device, - ) - if upper_indices.size(0) > 0 and y > 0: - perm = torch.randperm(upper_indices.size(0), device=device, generator=generator)[:y] - chosen_upper = upper_indices[perm] - for i, j in chosen_upper.tolist(): - mask[i, j] = True - mask[j, i] = True - return mask diff --git a/qubosolver/main.py b/qubosolver/main.py deleted file mode 100644 index c72d893e..00000000 --- a/qubosolver/main.py +++ /dev/null @@ -1,15 +0,0 @@ -from __future__ import annotations - -from typing import Optional - - -def main(str_to_add: Optional[str] = None) -> str: - msg = "Welcome to Qubo Solver!" - if str_to_add is not None: - msg += str_to_add - return msg - - -if __name__ == "__main__": - msg = main(str_to_add="\n(Executed from main.py)") - print(msg) diff --git a/qubosolver/pipeline/__init__.py b/qubosolver/pipeline/__init__.py deleted file mode 100644 index 8a8f01b3..00000000 --- a/qubosolver/pipeline/__init__.py +++ /dev/null @@ -1,13 +0,0 @@ -from __future__ import annotations - -from .basesolver import BaseSolver -from .embedder import get_embedder -from .fixtures import Fixtures -from .drive import get_drive_shaper - -__all__ = [ - "get_drive_shaper", - "get_embedder", - "BaseSolver", - "Fixtures", -] diff --git a/qubosolver/pipeline/drive.py b/qubosolver/pipeline/drive.py deleted file mode 100644 index 904bb394..00000000 --- a/qubosolver/pipeline/drive.py +++ /dev/null @@ -1,638 +0,0 @@ -from __future__ import annotations - -from abc import ABC, abstractmethod -from typing import cast - -import numpy as np -import torch -from skopt import gp_minimize -import math -import warnings - -from qoolqit import Register, QuantumProgram, Drive, Device -from qoolqit.waveforms import Interpolated as InterpolatedWaveform -from qubosolver import concepts - - -from qubosolver import QUBOInstance -from qubosolver.config import SolverConfig, compiler_profile, max_duration_ratio -from qubosolver.data import QUBOSolution -from qubosolver.qubo_types import DriveType -from qubosolver.utils import calculate_qubo_cost -from qubosolver.pipeline.waveforms import weighted_detunings - - -class BaseDriveShaper(ABC): - """ - Abstract base class for generating Qoolqit drives based on a QUBO problem. - - This class transforms the structure of a QUBOInstance into a quantum - waveform sequence or drive that can be applied to a physical register. The register - is passed at the time of drive generation, not during initialization. - - Attributes: - instance (QUBOInstance): The QUBO problem instance. - config (SolverConfig): The solver configuration. - drive (Drive, optional): A saved current drive obtained by `generate`. - backend (Backend): Backend to use. - device (Device): Device from backend. - - """ - - def __init__(self, instance: QUBOInstance, config: SolverConfig, backend: concepts.Backend): - """ - Initialize the drive shaping module with a QUBO instance. - - Args: - instance (QUBOInstance): The QUBO problem instance. - config (SolverConfig): The solver configuration. - backend (Backend): Backend to use. - """ - self.instance: QUBOInstance = instance - self.config: SolverConfig = config - self.drive: Drive | None = None - self.backend = backend - self.device = self.config.device - - # check if device allow DMM - self.dmm = self.config.drive_shaping.dmm and ( - len(list(self.config.device._device.dmm_channels.keys())) > 0 - ) - - @property - def qubo_coefficients(self) -> torch.Tensor: - return self.instance.coefficients - - @property - def qubo_normalized_coefficients(self) -> torch.Tensor: - return self.instance.normalized_coefficients - - def _compute_norm_weights(self) -> list[float]: - """Compute normalization weights. - - Returns: - list[float]: normalization weights. - """ - TIME, _, _ = self.device.converter.factors - weights_list = torch.abs(torch.diag(self.qubo_coefficients)).tolist() - max_node_weight = max(weights_list) if weights_list else 1.0 - norm_weights_list = [ - (1 - (w / max_node_weight)) if max_node_weight != 0 else 0.0 for w in weights_list - ] - return norm_weights_list - - @abstractmethod - def generate( - self, - register: Register, - ) -> tuple[Drive, QUBOSolution]: - """ - Generate a drive based on the problem and the provided register. - - Args: - register (Register): The physical register layout. - - Returns: - Drive: A generated Drive. - QUBOSolution: An instance of the qubo solution - """ - pass - - -def _compare_specs(a: dict, b: dict) -> None: - keys_a = set(a.keys()) - keys_b = set(b.keys()) - - if keys_a != keys_b: - only_in_a = keys_a - keys_b - only_in_b = keys_b - keys_a - if only_in_a: - warnings.warn(f"Keys present in a but missing in b: {only_in_a}") - if only_in_b: - warnings.warn(f"Keys present in b but missing in a: {only_in_b}") - - for key in keys_a & keys_b: - v1 = a[key] - v2 = b[key] - if v1 is None or v2 is None: - continue - if not math.isclose(float(v1), float(v2)): - warnings.warn(f"Value mismatch for key '{key}': a={v1}, b={v2}") - - -def _pulser_specs( - device: Device, normalize: bool = False, check_against_qoolqit: bool = False -) -> dict[str, float | None]: - pulser_device = device._device - channel = pulser_device.channels["rydberg_global"] - specs: dict[str, float | None] = {} - specs["max_duration"] = pulser_device.max_sequence_duration - specs["max_amplitude"] = channel.max_amp - specs["max_abs_detuning"] = channel.max_abs_detuning - specs["min_distance"] = pulser_device.min_atom_distance - specs["max_radial_distance"] = pulser_device.max_radial_distance - specs["min_avg_amp"] = channel.min_avg_amp - specs["dmm_bottom_detuning"] = None - - dmm_channels = list(getattr(pulser_device, "dmm_channels", {}).values()) - if dmm_channels: - specs["dmm_bottom_detuning"] = getattr(dmm_channels[0], "bottom_detuning", None) - - if not normalize: - return specs - - def _normalize(name: str, scale: float) -> None: - if specs[name] is not None: - specs[name] /= scale # type: ignore[operator] - - r0 = specs["min_distance"] - assert r0 is not None - _normalize("min_distance", r0) - _normalize("max_radial_distance", r0) - - C6 = pulser_device.interaction_coeff - J0 = C6 / (r0**6) - - _normalize("max_amplitude", J0) - _normalize("max_abs_detuning", J0) - _normalize("min_avg_amp", J0) - _normalize("dmm_bottom_detuning", J0) - - # J0 is in rad/us and t in ns, hence the factor 1000 - _normalize("max_duration", 1000.0 / J0) - - if check_against_qoolqit: - qq_specs = _qoolqit_specs(device) - _compare_specs(specs, qq_specs) - - return specs - - -def _qoolqit_specs( - device: Device, complete_with_pulser: bool = False, check_against_pulser: bool = False -) -> dict[str, float | None]: - specs = device.specs - pulser_specs = _pulser_specs(device, normalize=True) - - def import_from_pulser_or_set(name: str, fallback: float | None = None) -> None: - if name in specs.keys(): - return - if complete_with_pulser: - specs[name] = pulser_specs.get(name, fallback) - else: - specs[name] = fallback - - # Don't merge, update specific keys that are known not be in Qoolqit specs - import_from_pulser_or_set("min_avg_amp") - import_from_pulser_or_set("dmm_bottom_detuning") - - if check_against_pulser: - _compare_specs(specs, pulser_specs) - - return specs - - -class HeuristicDriveShaper(BaseDriveShaper): - """ - Heuristic schedule drive shaper. - - With DMM: - Final target encoding: - d_i = -alpha * Q_ii - - DMM convention in this stack: - WeightedDetuning waveform must be <= 0 - - Hence we encode the local final detuning as: - delta_i(T) = delta_g(T) + delta_dmm(T) * w_i - - with: - delta_g(T) = d_max - delta_dmm(T) = -(d_max - d_min) <= 0 - w_i = (d_max - d_i) / (d_max - d_min) in [0, 1] - - so that: - delta_i(T) = d_i - - Without DMM: - Only a global detuning is available, so the final detuning is chosen as: - delta_g(T) = mean(d_i) - and no weighted detunings are declared. - """ - - def generate(self, register: Register) -> tuple[Drive, QUBOSolution]: - - # Heuristic coefficient for omega - kappa = self.config.drive_shaping.heuristic_kappa - # Hardware bounds - specs = self.device.specs - max_seq_duration: float = specs["max_duration"] or 1000.0 - pulser_specs = _pulser_specs(self.device) - use_dmm = self.dmm and (pulser_specs["dmm_bottom_detuning"] is not None) - - max_amplitude = specs["max_amplitude"] or 1.0 - max_abs_detuning = specs["max_abs_detuning"] or 9.0 - - det_amp_ratio = max_amplitude / max_abs_detuning - if kappa < det_amp_ratio: - warnings.warn( - f"heuristic_kappa is too small ({kappa}), you're likely to get a qoolqit CompilationError. Set it above {det_amp_ratio}." - ) - - Q = self.qubo_normalized_coefficients - n = self.instance.size - - # Target local final detunings - d = (-0.5 * torch.diag(Q)).cpu().numpy() - d_min = np.min(d) - d_max = np.max(d) - - if use_dmm: - # Final global detuning is the top value, DMM pulls down locally - delta_g_T = d_max - - # DMM convention required by WeightedDetuning: - # waveform must be <= 0 - spread = max(0.0, d_max - d_min) - # if spread > 1e-15 and delta_dmm_max > 0.0: - if spread > 1e-15: - delta_dmm_T = -spread # must be <= 0 - denom = d_max - d_min - weights = ((d_max - d) / denom).clip(0.0, 1.0).tolist() - else: - use_dmm = False - delta_dmm_T = 0.0 - weights = [0.0] * n - else: - # No DMM: use a single global final detuning - delta_g_T = np.mean(d) - delta_dmm_T = 0.0 - weights = [0.0] * n - - omega_max = kappa * np.max(np.abs(d)) - # How to get max detuning ? - delta_0 = -np.max(np.abs(d)) - - # Amplitude waveform: 0 -> plateau -> 0 - eps = 1e-9 - amp_wave = InterpolatedWaveform( - max_seq_duration, - [eps, omega_max, omega_max, eps], - ) - - # Global detuning waveform: initial negative -> final target - det_wave = InterpolatedWaveform( - max_seq_duration, - [delta_0, delta_0, delta_g_T, delta_g_T], - ) - - # DMM weighted detunings - wdetunings = None - if use_dmm: - wdetunings = weighted_detunings( - register, - max_seq_duration, - weights, - final_detuning=delta_dmm_T, - ) - - shaped_drive = Drive( - amplitude=amp_wave, - detuning=det_wave, - dmm=wdetunings, - ) - solution = QUBOSolution(torch.Tensor(), torch.Tensor()) - return shaped_drive, solution - - -class OptimizedDriveShaper(BaseDriveShaper): - """ - Drive shaper that uses optimization to find the best drive parameters for solving QUBOs. - Returns an optimized drive, the bitstrings, their counts, probabilities, and costs. - - Attributes: - drive (Drive): current drive. - best_cost (float): Current best cost. - best_bitstring (Tensor | list): Current best bitstring. - bitstrings (Tensor | list): List of current bitstrings obtained. - counts (Tensor | list): Frequencies of bitstrings. - probabilities (Tensor | list): Probabilities of bitstrings. - costs (Tensor | list): Qubo cost. - optimized_custom_qubo_cost (Callable[[str, torch.Tensor], float], optional): - Apply a different qubo cost evaluation during optimization. - Must be defined as: - `def optimized_custom_qubo_cost(bitstring: str, QUBO: torch.Tensor) -> float`. - Defaults to None, meaning we use the default QUBO evaluation. - optimized_custom_objective_fn (Callable[[list, list, list, list, float, str], float], optional): - For bayesian optimization, one can change the output of - `self.run_simulation` to optimize differently. Instead of using the best cost - out of the samples, one can change the objective for an average, - or any function out of the form - `cost_eval = optimized_custom_objective_fn(bitstrings, - counts, probabilities, costs, best_cost, best_bitstring)` - Defaults to None, which means we optimize using the best cost - out of the samples. - optimized_callback_objective (Callable[..., None], optional): Apply a callback - during bayesian optimization. Only accepts one input dictionary - created during optimization `d = {"x": x, "cost_eval": cost_eval}` - hence should be defined as: - `def callback_fn(d: dict) -> None:` - Defaults to None, which means no callback is applied. - """ - - def __init__( - self, - instance: QUBOInstance, - config: SolverConfig, - backend: concepts.Backend, - ): - """Instantiate an `OptimizedDriveShaper`. - - Args: - instance (QUBOInstance): Qubo instance. - config (SolverConfig): Configuration for solving. - backend (Backend): Backend to use during optimization. - - """ - super().__init__(instance, config, backend) - - self.drive = None - self.best_cost = None - self.best_bitstring = None - self.best_params = None - self.bitstrings = None - self.counts = None - self.probabilities = None - self.costs = None - self.optimized_custom_qubo_cost = self.config.drive_shaping.optimized_custom_qubo_cost - self.optimized_custom_objective_fn = self.config.drive_shaping.optimized_custom_objective - self.optimized_callback_objective = self.config.drive_shaping.optimized_callback_objective - - def generate( - self, - register: Register, - ) -> tuple[Drive, QUBOSolution]: - """ - Generate a drive via optimization. - - Args: - register (Register): The physical register layout. - - Returns: - Drive: A generated Drive. - QUBOSolution: An instance of the qubo solution - """ - # TODO: Harmonize the output of the pulse_shaper generate - QUBO = self.qubo_coefficients - self.register = register - - self.norm_weights_list = self._compute_norm_weights() - - n_amp = 3 - n_det = 3 - - eps = 0.0001 - zero = eps - one = 1.0 - eps - - bounds = ( - [(zero, one)] * n_amp + [(-one, -zero)] + [(-one, one)] * (n_det - 2) + [(zero, one)] - ) - x0 = ( - self.config.drive_shaping.optimized_initial_omega_parameters - + self.config.drive_shaping.optimized_initial_detuning_parameters - ) - - def objective(x: list[float]) -> float: - drive = self.build_drive(x) - - try: - bitstrings, counts, probabilities, costs, cost_eval, best_bitstring = ( - self.run_simulation( - self.register, - drive, - QUBO, - convert_to_tensor=False, - ) - ) - if self.optimized_custom_objective_fn is not None: - cost_eval = self.optimized_custom_objective_fn( - bitstrings, - counts, - probabilities, - costs, - cost_eval, - best_bitstring, - ) - if not np.isfinite(cost_eval): - print(f"[Warning] Non-finite cost encountered: {cost_eval} at x={x}") - cost_eval = 1e4 - - except Exception as e: - print(f"[Exception] Error during simulation at x={x}: {e}") - cost_eval = 1e4 - - if self.optimized_callback_objective is not None: - self.optimized_callback_objective({"x": x, "cost_eval": cost_eval}) - return float(cost_eval) - - opt_result = gp_minimize( - objective, - bounds, - x0=x0, - n_calls=self.config.drive_shaping.optimized_n_calls, - random_state=self.config.drive_shaping.optimized_seed, - ) - - if opt_result and opt_result.x: - self.best_params = opt_result.x - self.drive = self.build_drive(self.best_params) # type: ignore[arg-type] - - ( - self.bitstrings, - self.counts, - self.probabilities, - self.costs, - self.best_cost, - self.best_bitstring, - ) = self.run_simulation(self.register, self.drive, QUBO, convert_to_tensor=True) - - if self.bitstrings is None or self.counts is None: - # TODO: what needs to be returned here? - # the generate function should always return a drive - even if it is not good. - # we need to return a drive (self.drive) - which is none here. - # return self.drive, QUBOSolution(None, None) - raise RuntimeError("No solution found") - - assert self.costs is not None - solution = QUBOSolution( - bitstrings=self.bitstrings, - counts=self.counts, - probabilities=self.probabilities, - costs=self.costs, - ) - assert self.drive is not None - return self.drive, solution - - def build_drive(self, params: list) -> Drive: - """Build the drive from a list of parameters for the objective. - - Args: - params (list): List of parameters. - - Returns: - Drive: Drive sequence. - """ - specs = self.device.specs - max_seq_duration: float = specs["max_duration"] or 1e3 - max_amplitude: float = specs["max_amplitude"] or 1e4 - max_detuning: float = specs["max_abs_detuning"] or 1e4 - - amp_params = [1e-9] + list(params[:3]) + [1e-9] - det_params = list(params[3:]) - amp_params = [p * max_amplitude for p in amp_params] - det_params = [p * max_detuning for p in det_params] - - amp_wave = InterpolatedWaveform(max_seq_duration, amp_params) - det_wave = InterpolatedWaveform(max_seq_duration, det_params) - - wdetunings = None - final_detuning = det_params[-1] - if self.dmm and final_detuning > 0: - wdetunings = weighted_detunings( - self.register, - max_seq_duration, - self.norm_weights_list, - final_detuning=-final_detuning, - ) - - shaped_drive = Drive(amplitude=amp_wave, detuning=det_wave, dmm=wdetunings) - - return shaped_drive - - def compute_qubo_cost(self, bitstring: str, QUBO: torch.Tensor) -> float: - """The qubo cost for a single bitstring to apply during optimization. - - Args: - bitstring (str): candidate bitstring. - QUBO (torch.Tensor): qubo coefficients. - - Returns: - float: respective cost of bitstring. - """ - if self.optimized_custom_qubo_cost is None: - return calculate_qubo_cost(bitstring, QUBO) - - return cast(float, self.optimized_custom_qubo_cost(bitstring, QUBO)) - - def run_simulation( - self, - register: Register, - drive: Drive, - QUBO: torch.Tensor, - convert_to_tensor: bool = True, - ) -> tuple: - """Run a quantum program using backend and returns - a tuple of (bitstrings, counts, probabilities, costs, best cost, best bitstring). - - Args: - register (Register): register of quantum program. - drive (Drive): drive to run on backend. - QUBO (torch.Tensor): Qubo coefficients. - convert_to_tensor (bool, optional): Convert tuple components to tensors. - Defaults to True. - - Returns: - tuple: tuple of (bitstrings, counts, probabilities, costs, best cost, best bitstring) - """ - try: - program = QuantumProgram(register=register, drive=drive) - program.compile_to( - self.device, - profile=compiler_profile(self.config), - device_max_duration_ratio=max_duration_ratio(self.config), - ) - job = self.backend.run(program) - bitstring_counts = job.results().final_bitstrings - - cost_dict = {b: self.compute_qubo_cost(b, QUBO) for b in bitstring_counts.keys()} - - best_bitstring = min(cost_dict, key=cost_dict.get) # type: ignore[arg-type] - best_cost = cost_dict[best_bitstring] - - if convert_to_tensor: - keys = list(bitstring_counts.keys()) - values = list(bitstring_counts.values()) - - bitstrings_tensor = torch.tensor( - [[int(b) for b in bitstr] for bitstr in keys], dtype=torch.int32 - ) - counts_tensor = torch.tensor(values, dtype=torch.int32) - probabilities_tensor = counts_tensor.float() / counts_tensor.sum() - - costs_tensor = torch.tensor( - [self.compute_qubo_cost(b, QUBO) for b in keys], dtype=torch.float32 - ) - - return ( - bitstrings_tensor, - counts_tensor, - probabilities_tensor, - costs_tensor, - best_cost, - best_bitstring, - ) - else: - counts = list(bitstring_counts.values()) - nsamples = float(sum(counts)) - return ( - list(bitstring_counts.keys()), - counts, - [c / nsamples for c in counts], - list(cost_dict.values()), - best_cost, - best_bitstring, - ) - - except Exception as e: - print(f"Simulation failed: {e}") - return ( - torch.tensor([]), - torch.tensor([]), - torch.tensor([]), - torch.tensor([]), - float("inf"), - None, - ) - - -def get_drive_shaper( - instance: QUBOInstance, - config: SolverConfig, - backend: concepts.Backend, -) -> BaseDriveShaper: - """ - Method that returns the correct DriveShaper based on configuration. - The correct drive shaping method can be identified using the config, and an - object of this driveshaper can be returned using this function. - - Args: - instance (QUBOInstance): The QUBO problem to embed. - config (SolverConfig): The solver configuration used. - backend (Backend): Backend to extract device from or to use - during drive shaping. - - Returns: - (BaseDriveShaper): The representative Drive Shaper object. - """ - if config.drive_shaping.drive_shaping_method == DriveType.HEURISTIC: - return HeuristicDriveShaper(instance, config, backend) - elif config.drive_shaping.drive_shaping_method == DriveType.OPTIMIZED: - return OptimizedDriveShaper(instance, config, backend) - elif issubclass(config.drive_shaping.drive_shaping_method, BaseDriveShaper): - return cast( - BaseDriveShaper, - config.drive_shaping.drive_shaping_method(instance, config, backend), - ) - else: - raise NotImplementedError diff --git a/qubosolver/pipeline/embedder.py b/qubosolver/pipeline/embedder.py deleted file mode 100644 index b2732679..00000000 --- a/qubosolver/pipeline/embedder.py +++ /dev/null @@ -1,220 +0,0 @@ -from __future__ import annotations - -import typing -from abc import ABC, abstractmethod -import numpy as np -import torch -import warnings - -from qoolqit import Register -from qoolqit.devices import Device -from qoolqit.embedding.matrix_embedder import Blade, BladeConfig - - -from qubosolver import QUBOInstance, concepts -from qubosolver.algorithms.greedy.greedy import Greedy -from qubosolver.config import EmbedderType, SolverConfig -from qubosolver.utils.density import calculate_density - -warnings.filterwarnings("ignore", module="pulser") - - -class BaseEmbedder(ABC): - """ - Abstract base class for all embedders. - - Prepares the geometry (register) of atoms based on the QUBO instance. - Returns a Register compatible with Pasqal/Pulser devices. - """ - - def __init__(self, instance: QUBOInstance, config: SolverConfig, backend: concepts.Backend): - """ - Args: - instance (QUBOInstance): The QUBO problem to embed. - config (SolverConfig): The Solver Configuration. - """ - self.instance: QUBOInstance = instance - self.config: SolverConfig = config - self.register: Register | None = None - self.backend = backend - - # TODO: remove when bumping to qoolqit v1 - # for converting to qoolqit - self._distance_conversion = self.config.device.converter.factors[2] - - @abstractmethod - def embed(self) -> Register: - """ - Creates a layout of atoms as the register. - - Returns: - Register: The register. - """ - ... - - -class BLaDEmbedder(BaseEmbedder): - - def embed(self) -> Register: - - embed_config = self.config.embedding - default = BladeConfig() - step_per_round = embed_config.blade_steps_per_round - if step_per_round is None: - step_per_round = default.steps_per_round - if embed_config.blade_starting_positions is not None: - starting_positions = embed_config.blade_starting_positions.numpy() - else: - starting_positions = None - - min_distance = self.config.embedding.min_distance - max_radial_distance = self.config.device.specs["max_radial_distance"] - if min_distance is None or max_radial_distance is None: - device = self.config.device - max_min_dist_ratio = None - else: - device = None - max_min_dist_ratio = max_radial_distance / min_distance - - config = BladeConfig( - steps_per_round=step_per_round, - starting_positions=starting_positions, - dimensions=tuple(embed_config.blade_dimensions), - max_min_dist_ratio=max_min_dist_ratio, - device=device, - ) - - _blade = Blade(config) - graph = _blade.embed(self.instance.coefficients.numpy()) - if min_distance is not None: - graph.rescale_coords(spacing=min_distance) - register = Register.from_graph(graph) - - return register - - -class GreedyEmbedder(BaseEmbedder): - """Create an embedding in a greedy fashion. - - At each step, place one logical node onto one trap to minimize the - incremental mismatch between the logical QUBO matrix Q and the physical - interaction matrix U (approx. C / ||r_i - r_j||^6). - """ - - @staticmethod - def _number_of_traps_from_device(device: Device) -> int: - """Determine the number of traps to use based on the device constraints. - - Inspects the device's layout and atom number limits to derive an - appropriate trap count. The resolution order is: - - 1. ``max_layout_traps`` – if the device exposes a hard trap limit, use it directly. - 2. ``max_atom_num`` / ``max_layout_filling`` – if only an atom-number limit is - available, derive the minimum number of traps needed to accommodate that - many atoms at the device's maximum filling ratio. - 3. Fallback – return ``200`` when neither property is set. - - Args: - device (Device): The quantum device whose constraints are inspected. - - Returns: - int: The number of traps to allocate for the embedding. - """ - - if device._device.max_layout_traps: - return device._device.max_layout_traps - - if device._device.max_atom_num: - return int(np.ceil(device._device.max_atom_num / device._device.max_layout_filling)) - - return 200 - - @typing.no_type_check - def embed(self) -> Register: - """ - Creates a layout of atoms as the register. - - Returns: - Register: The register. - """ - if self.config.embedding.greedy_traps == -1: - self.config.embedding.greedy_traps = self._number_of_traps_from_device( - self.config.device - ) - - if self.config.embedding.greedy_traps < self.instance.size: - raise ValueError( - "Number of traps must be at least equal to the number of atoms on the register." - ) - - # compute density (unchanged) - self.config.embedding.greedy_density = calculate_density( - self.instance.coefficients, self.instance.size - ) - - # build params for the Greedy algorithm - params = { - "device": self.config.device._device, - "layout": self.config.embedding.greedy_layout, - "traps": int(self.config.embedding.greedy_traps), - "spacing": float(self.config.embedding.greedy_spacing), - # animation controls (all read by Greedy) - "draw_steps": bool(self.config.embedding.draw_steps), # collect per-step data - "animation": bool(self.config.embedding.draw_steps), # render animation after run - "animation_save_path": self.config.embedding.animation_save_path, # optional export - } - - # --- DEBUG / INFO: show where Greedy comes from + the params we’ll pass - dev = params["device"] - dev_str = ( - getattr(dev, "name", None) - or getattr(dev, "device_name", None) - or dev.__class__.__name__ - ) - printable = dict(params) - printable["device"] = dev_str # avoid dumping the whole object - # --- Call Greedy (unchanged public signature) - best, _, coords, _, _ = Greedy().launch_greedy( - Q=self.instance.coefficients, - params=params, - # no extra kwargs; Greedy reads animation/draw/save_path from params - ) - min_distance = self.config.embedding.min_distance - if min_distance is not None: - min_reg_distance = torch.cdist(coords, coords).fill_diagonal_(float("inf")).min() - coords *= min_distance / min_reg_distance - else: - coords /= self._distance_conversion - - # build the register (unchanged) - qubits = {f"q{i}": coord for i, coord in enumerate(coords)} - register = Register(qubits) - return register - - -def get_embedder( - instance: QUBOInstance, config: SolverConfig, backend: concepts.Backend -) -> BaseEmbedder: - """ - Method that returns the correct embedder based on configuration. - The correct embedding method can be identified using the config, and an - object of this embedding can be returned using this function. - - Args: - instance (QUBOInstance): The QUBO problem to embed. - config (Device): The quantum device to target. - - Returns: - (BaseEmbedder): The representative embedder object. - """ - - if config.embedding.embedding_method == EmbedderType.BLADE: - return BLaDEmbedder(instance, config, backend) - elif config.embedding.embedding_method == EmbedderType.GREEDY: - return GreedyEmbedder(instance, config, backend) - elif issubclass(config.embedding.embedding_method, BaseEmbedder): - return typing.cast( - BaseEmbedder, config.embedding.embedding_method(instance, config, backend) - ) - else: - raise NotImplementedError diff --git a/qubosolver/pipeline/fixtures.py b/qubosolver/pipeline/fixtures.py deleted file mode 100644 index cf25da33..00000000 --- a/qubosolver/pipeline/fixtures.py +++ /dev/null @@ -1,307 +0,0 @@ -from __future__ import annotations - -from copy import deepcopy -from typing import Callable, Dict, List, cast - -import numpy as np -import torch - -from qubosolver import QUBOInstance, QUBOSolution -from qubosolver.config import SolverConfig -from qubosolver.qubo_types import SolutionStatusType - - -def bit_flip_local_search( - qubo_func: Callable[[np.ndarray], float], s: np.ndarray, shuffle: bool = True -) -> tuple[np.ndarray, float]: - """ - Performs a local search by flipping bits to improve the objective value. - - Args: - qubo_func: Function that computes the objective value for a solution array. - s (np.ndarray): Binary array representing a candidate solution. - shuffle (bool, optional): Shuffle to diversify - - Returns: - tuple[np.ndarray, float]: The improved solution and its objective value. - """ - s_current = s.copy() - current_objective = qubo_func(s_current) - while True: - best_idx = None - best_obj = current_objective - indices = list(range(len(s_current))) - if shuffle: - np.random.shuffle(indices) # option to diversify - # Evaluate all possible flips, keep best - for i in indices: - s_new = s_current.copy() - s_new[i] = 1 - s_new[i] - new_objective = qubo_func(s_new) - if new_objective < best_obj: - best_obj = new_objective - best_idx = i - # If no improvements, stop - if best_idx is None: - break - # Apply best flip - s_current[best_idx] = 1 - s_current[best_idx] - current_objective = best_obj - return s_current, current_objective - - -def hansen_fixing(qubo: QUBOInstance) -> Dict[int, int]: - """ - Identifies and fixes variables in a QUBO instance based on threshold conditions. - - This method determines whether a variable should be fixed to 0 or 1 by computing - lower and upper bounds from the diagonal and off-diagonal elements of the QUBO matrix. - - Args: - qubo (QUBOInstance): The QUBO instance containing the coefficients matrix. - - Returns: - Dict[int, int]: A dictionary mapping variable indices to fixed values (0 or 1). - """ - if qubo.coefficients is None: - raise ValueError("QUBO coefficients are not initialized.") - - fixed_dict: Dict[int, int] = {} - size_raw = qubo.size - size: int = cast(int, size_raw) - epsilon: float = 1e-8 # Tolerance to avoid floating-point precision issues - - for i in range(size): - ci = qubo.coefficients[i, i].item() # Diagonal element - - q_minus = sum(min(0, qubo.coefficients[i, j].item()) for j in range(size) if j != i) - q_plus = sum(max(0, qubo.coefficients[i, j].item()) for j in range(size) if j != i) - - if ci + q_minus * 2 >= -epsilon: - fixed_dict[i] = 0 - elif ci + q_plus * 2 <= epsilon: - fixed_dict[i] = 1 - - return fixed_dict - - -class Fixtures: - """ - Handles all preprocessing and postprocessing logic for QUBO problems. - - This class centralizes the transformation or validation of the problem instance - before solving and modification or annotation of the solution after solving. - """ - - def __init__(self, instance: QUBOInstance, config: SolverConfig): - """ - Initialize the fixture handler with the QUBO instance and solver configuration. - - Args: - instance (QUBOInstance): The QUBO problem instance to process. - config (SolverConfig): Solver configuration, which may include flags for enabling - postprocessing (e.g., do_postprocessing). - """ - self.instance = instance - self.config = config - - self.reduced_qubo: QUBOInstance = deepcopy(instance) - self.fixation_rule_list: List[Callable[[QUBOInstance], Dict[int, int]]] = [ - hansen_fixing, - ] - self.fixed_var_dict_list: List[Dict[int, int]] = [] - - @property - def n_fixed_variables(self) -> int: - """Returns the number of fixed variables. - - Returns: - int: The number of fixed variables. - """ - return sum([len(fixed) for fixed in self.fixed_var_dict_list]) - - def preprocess(self) -> QUBOInstance: - """ - Apply preprocessing steps to the QUBO instance before solving. - - Returns: - QUBOInstance: The processed or annotated instance. - """ - - # Check if preprocessing is enabled via the configuration. - if not hasattr(self.config, "do_preprocessing") or not self.config.do_preprocessing: - return self.instance - - # Apply every rules until exhaustion - self.apply_full_fixation_exhaust() - - return self.instance - - def postprocess(self, solution: QUBOSolution) -> QUBOSolution: - """ - Apply postprocessing steps to the QUBO solution after solving. - - This method iterates over all solutions in the bitstrings tensor and, for each solution, - performs a local bit-flip search to attempt to improve its objective value. The objective - is computed using the QUBOInstance's evaluate_solution method. - - Args: - solution (QUBOSolution): The raw solution from the solver. - - Returns: - QUBOSolution: The updated solution with improved bitstrings and costs. - """ - - if not hasattr(self.config, "do_postprocessing") or not self.config.do_postprocessing: - return solution - - # If there are no bitstrings, return the solution unchanged. - if solution.bitstrings.numel() == 0: - return solution - - # Define an objective function that uses the existing evaluate_solution method. - def qubo_objective(s_arr: np.ndarray) -> float: - # Convert the solution array to a list of integers - return self.instance.evaluate_solution(s_arr.tolist()) - - improved_bitstrings = [] - improved_costs = [] - num_solutions = solution.bitstrings.shape[0] - - for idx in range(num_solutions): - # Get the current solution (row) as a numpy array of integers. - s_orig = solution.bitstrings[idx].detach().cpu().numpy().astype(int) - # Apply bit-flip local search to improve the solution. - s_improved, new_cost = bit_flip_local_search(qubo_objective, s_orig) - improved_bitstrings.append(s_improved) - improved_costs.append(new_cost) - - # Create new tensors for the improved solutions and their costs. - new_bitstrings_tensor = torch.tensor(np.array(improved_bitstrings), dtype=torch.float32) - new_costs_tensor = torch.tensor(improved_costs, dtype=torch.float32) - - # Update the solution object. - solution.bitstrings = new_bitstrings_tensor - solution.costs = new_costs_tensor - - if self.config.do_preprocessing and self.config.do_postprocessing: - solution.solution_status = SolutionStatusType.PREPOSTPROCESSED - else: - solution.solution_status = SolutionStatusType.POSTPROCESSED - - return solution - - def reduce_qubo(self, fixed_dict: Dict[int, int]) -> None: - """ - Applies variable fixation to reduce the size of the QUBO problem. - - This function modifies the QUBO coefficient matrix by: - - Removing rows and columns corresponding to fixed variables. - - Adjusting diagonal elements to account for fixed variables. - - Args: - fixed_dict (Dict[int, int]): A dictionary of fixed variable assignments. - - Keys are variable indices. - - Values are fixed binary values (0 or 1). - - Returns: - None: Modifies `self.reduced_qubo` in place. - """ - Q = self.reduced_qubo.coefficients.clone() - - fixed_to_0 = {i for i, v in fixed_dict.items() if v == 0} - fixed_to_1 = {i for i, v in fixed_dict.items() if v == 1} - fixed_vars = sorted(fixed_to_0 | fixed_to_1, reverse=True) - - for i in fixed_vars: - if i >= Q.shape[0]: - continue - - if i in fixed_to_1: - for j in range(Q.shape[0]): - if j != i: - Q[j, j] += Q[i, j] * 2 - - Q = torch.cat((Q[:i, :], Q[i + 1 :, :]), dim=0) - Q = torch.cat((Q[:, :i], Q[:, i + 1 :]), dim=1) - - self.reduced_qubo.coefficients = Q - self.reduced_qubo.update_metrics() - - def apply_rule(self, fixation_rule: Callable[[QUBOInstance], Dict[int, int]]) -> int: - """ - Applies a given variable fixation rule to the reduced QUBO instance. - - Args: - fixation_rule (Callable[[QUBOInstance], Dict[int, int]]): A function that - returns a dictionary mapping variable indices to fixed values. - - Returns: - int: The number of variables fixed by this rule. - """ - fixed = fixation_rule(self.reduced_qubo) - self.reduce_qubo(fixed) - if fixed: - self.fixed_var_dict_list.append(fixed) - - return len(fixed) - - def apply_full_fixation_exhaust(self) -> None: - """ - Iteratively applies all fixation rules until no more variables can be fixed. - - This function repeatedly applies all rules in `self.fixation_rule_list` - until no further reduction is possible. - """ - fixed_sum = 1 - while fixed_sum > 0: - fixed_sum = 0 - for fixation_rule in self.fixation_rule_list: - fixed_var_number = self.apply_rule(fixation_rule) - fixed_sum += fixed_var_number - - def post_process_fixation(self, solution: QUBOSolution) -> QUBOSolution: - """ - Restores fixed variables in the solution bitstrings after QUBO reduction. - - This method reconstructs the full-length bitstrings by reinserting the fixed - variables at their original positions. - - Args: - solution (QUBOSolution): The solution object from the reduced QUBO problem. - - Returns: - QUBOSolution: A solution object with bitstrings restored to their original size. - """ - if not getattr(self.config, "do_preprocessing", False): - return solution - - bitstring_list = solution.bitstrings.tolist() or [[]] - - def reinsert_fixed_variables(bitstring: List[int]) -> List[int]: - for fixation_dict in reversed(self.fixed_var_dict_list): - for position, bit_value in sorted(fixation_dict.items()): - bitstring.insert(position, bit_value) - return bitstring - - should_restore = not self.config.use_quantum or (self.instance.size is not None) - assert self.instance.size - should_restore = should_restore and (len(bitstring_list[0]) < self.instance.size) - - if should_restore: - bitstring_list = [reinsert_fixed_variables(bitstring) for bitstring in bitstring_list] - - bitstrings = torch.tensor(bitstring_list, dtype=torch.float32) - - costs = torch.tensor( - [self.instance.evaluate_solution(sample.tolist()) for sample in bitstrings], - dtype=torch.float32, - ) - - return QUBOSolution( - bitstrings=bitstrings, - costs=costs, - counts=solution.counts, - probabilities=solution.probabilities, - solution_status=SolutionStatusType.PREPROCESSED, - ) diff --git a/qubosolver/saveload.py b/qubosolver/saveload.py deleted file mode 100644 index 69badcd2..00000000 --- a/qubosolver/saveload.py +++ /dev/null @@ -1,73 +0,0 @@ -from __future__ import annotations - -from pathlib import Path -from typing import Any -import torch - -from qubosolver.data import QUBODataset, QUBOSolution - -# Modules to be automatically added to the qubosolver.utils namespace -__all__ = [ - "save_qubo_dataset", - "load_qubo_dataset", -] - - -def save_qubo_dataset(dataset: QUBODataset, filepath: str | Path) -> None: - """ - Saves a QUBODataset to a file. - - Args: - dataset (QUBODataset): - The QUBODataset object to save. - filepath (str | Path): - Path to the file where the QUBODataset will be saved. - - Notes: - The saved data includes: - - Coefficients (size x size x num_instances tensor) - - Solutions (optional, includes bitstrings, counts, probabilities, and costs) - """ - data: dict[str, Any] = {"coefficients": dataset.coefficients, "solutions": None} - if dataset.solutions is not None: - data["solutions"] = [ - { - "bitstrings": solution.bitstrings, - "counts": solution.counts, - "probabilities": solution.probabilities, - "costs": solution.costs, - } - for solution in dataset.solutions - ] - torch.save(data, filepath) - - -def load_qubo_dataset(filepath: str | Path) -> QUBODataset: - """ - Loads a QUBODataset from a file. - Notes: - The file should contain data saved in the format used by `save_qubo_dataset`. - - Args: - filepath (str | Path): - Path to the file from which the QUBODataset will be loaded. - - Returns: - QUBODataset: - The loaded QUBODataset object. - - - """ - data = torch.load(filepath) - solutions = None - if data["solutions"] is not None: - solutions = [ - QUBOSolution( - bitstrings=solution["bitstrings"], - counts=solution["counts"], - probabilities=solution["probabilities"], - costs=solution["costs"], - ) - for solution in data["solutions"] - ] - return QUBODataset(coefficients=data["coefficients"], solutions=solutions) diff --git a/qubosolver/types/__init__.py b/qubosolver/types/__init__.py new file mode 100644 index 00000000..014e1923 --- /dev/null +++ b/qubosolver/types/__init__.py @@ -0,0 +1,36 @@ +from __future__ import annotations + +from qubosolver.types import bitstring, bitstrings, matrix, tensor, vector, vectori +from qubosolver.types._linalg import Bitstring, Bitstrings, Matrix, Tensor, Vector, Vectori +from qubosolver.types._solution import QUBOSolution +from qubosolver.types._analyzer import QUBOAnalyzer +from qubosolver.types._instance import QUBOInstance +from qubosolver.types._enums import EmbedderType, LayoutType, DriveType, SolutionStatusType, DensityType + +__all__ = [ + # Submodules + "bitstring", + "bitstrings", + "matrix", + "tensor", + "vector", + "vectori", + # Type Aliases + "Bitstring", + "Bitstrings", + "Matrix", + "Tensor", + "Vector", + "Vectori", + # Classes + "QUBOSolution", + "QUBOAnalyzer", + "QUBOInstance", + # Enums + "EmbedderType", + "LayoutType", + "DriveType", + "SolutionStatusType", + "DensityType", + # Functions +] diff --git a/qubosolver/qubo_analyzer.py b/qubosolver/types/_analyzer.py similarity index 99% rename from qubosolver/qubo_analyzer.py rename to qubosolver/types/_analyzer.py index e1df3e46..c3e904aa 100644 --- a/qubosolver/qubo_analyzer.py +++ b/qubosolver/types/_analyzer.py @@ -4,10 +4,8 @@ import seaborn as sns import torch -from .data import QUBOSolution -from .qubo_instance import QUBOInstance - -__all__ = ["QUBOAnalyzer"] +from ._solution import QUBOSolution +from ._instance import QUBOInstance _BITSTRINGS = "bitstrings" _COSTS = "costs" diff --git a/qubosolver/types/_checks.py b/qubosolver/types/_checks.py new file mode 100644 index 00000000..953bbaed --- /dev/null +++ b/qubosolver/types/_checks.py @@ -0,0 +1,18 @@ +from __future__ import annotations + +import os +from typing import TypeVar + +_T = TypeVar("_T") + +# Tensor Types +_RUNTIME_CHECKS: bool = os.getenv("QUBO_SOLVER_RUNTIME_CHECKS", "0") == "1" + + +def debug_runtime_typecheck(target: _T) -> _T: + """Applies @beartype if QUBO_SOLVER_RUNTIME_CHECKS is enabled, otherwise a no-op.""" + if _RUNTIME_CHECKS: + from beartype import beartype + + return beartype(target) # type: ignore[no-any-return] + return target diff --git a/qubosolver/data.py b/qubosolver/types/_dataset.py similarity index 76% rename from qubosolver/data.py rename to qubosolver/types/_dataset.py index 9f2a1ef1..a2020d77 100644 --- a/qubosolver/data.py +++ b/qubosolver/types/_dataset.py @@ -1,13 +1,12 @@ from __future__ import annotations -from dataclasses import dataclass from typing import TYPE_CHECKING, Any, Iterator import torch from torch.utils.data import Dataset from qubosolver.data_utils import generate_symmetric_mask -from qubosolver.qubo_types import SolutionStatusType + if TYPE_CHECKING: pass @@ -16,93 +15,6 @@ __all__ = ["QUBOSolution", "QUBODataset"] -@dataclass -class QUBOSolution: - """ - Represents a solution to a QUBO problem. - - Attributes: - bitstrings (torch.Tensor): - Tensor of shape (num_solutions, bitstring_length), containing the bitstring solutions. - Each entry is an integer tensor with 0s and 1s. - counts (torch.Tensor | None): - Tensor of shape (num_solutions,), containing the count of occurrences of each bitstring. - Optional, can be None. - probabilities (torch.Tensor | None): - Tensor of shape (num_solutions,), containing the probability of each bitstring solution. - Optional, can be None. - costs (torch.Tensor): - Tensor of shape (num_solutions,), containing the cost associated with each - bitstring solution. - """ - - bitstrings: torch.Tensor - costs: torch.Tensor - counts: torch.Tensor | None = None - probabilities: torch.Tensor | None = None - solution_status: SolutionStatusType = SolutionStatusType.UNPROCESSED - - def compute_costs(self, instance: Any) -> torch.Tensor: - """ - Computes the cost for each bitstring solution based on the provided QUBO instance. - - Args: - instance (QUBOInstance): The QUBO instance containing the QUBO matrix. - - Returns: - torch.Tensor: A tensor of costs for each bitstring. - """ - # Retrieve the QUBO matrix from the QUBOInstance - QUBO = instance.coefficients # Assuming `coefficients` holds the QUBO matrix - - costs = [] - for bitstring in self.bitstrings: - if isinstance(bitstring, str): - z = torch.tensor([int(b) for b in bitstring], dtype=torch.float32) - elif isinstance(bitstring, torch.Tensor): - z = bitstring.detach().clone() - else: - z = torch.tensor(bitstring, dtype=torch.float32) - cost = torch.matmul( - z.permute(*torch.arange(z.ndim - 1, -1, -1)), torch.matmul(QUBO, z) - ).item() # Use the QUBO matrix from the instance - costs.append(cost) - - return torch.tensor(costs) - - def compute_probabilities(self) -> torch.Tensor: - """ - Computes the probabilities of each bitstring solution based on their counts. - - Returns: - torch.Tensor: A tensor of probabilities for each bitstring. - """ - if self.counts is None: - raise ValueError("Counts are required to compute probabilities.") - - total_counts = self.counts.sum().item() - probabilities = ( - self.counts / total_counts if total_counts > 0 else torch.zeros_like(self.counts) - ) - return probabilities - - def sort_by_cost(self) -> None: - """ - Sorts the QUBOSolution in-place by increasing cost. - - Reorders bitstrings, costs, counts, and probabilities (if available) - based on the ascending order of the costs. - """ - - sorted_indices = torch.argsort(self.costs) - self.bitstrings = self.bitstrings[sorted_indices] - self.costs = self.costs[sorted_indices] - if self.counts is not None: - self.counts = self.counts[sorted_indices] - if self.probabilities is not None: - self.probabilities = self.probabilities[sorted_indices] - - class QUBODataset(Dataset): """ Represents a dataset for Quadratic Unconstrained Binary Optimization (QUBO) problems. @@ -337,3 +249,111 @@ def from_random( # Step 4: Return the dataset. return cls(coefficients=coefficients) + + +def save_qubo_dataset(dataset: QUBODataset, filepath: str | Path) -> None: + """ + Saves a QUBODataset to a file. + + Args: + dataset (QUBODataset): + The QUBODataset object to save. + filepath (str | Path): + Path to the file where the QUBODataset will be saved. + + Notes: + The saved data includes: + - Coefficients (size x size x num_instances tensor) + - Solutions (optional, includes bitstrings, counts, probabilities, and costs) + """ + data: dict[str, Any] = {"coefficients": dataset.coefficients, "solutions": None} + if dataset.solutions is not None: + data["solutions"] = [ + { + "bitstrings": solution.bitstrings, + "counts": solution.counts, + "probabilities": solution.probabilities, + "costs": solution.costs, + } + for solution in dataset.solutions + ] + torch.save(data, filepath) + + +def load_qubo_dataset(filepath: str | Path) -> QUBODataset: + """ + Loads a QUBODataset from a file. + Notes: + The file should contain data saved in the format used by `save_qubo_dataset`. + + Args: + filepath (str | Path): + Path to the file from which the QUBODataset will be loaded. + + Returns: + QUBODataset: + The loaded QUBODataset object. + + + """ + data = torch.load(filepath) + solutions = None + if data["solutions"] is not None: + solutions = [ + QUBOSolution( + bitstrings=solution["bitstrings"], + counts=solution["counts"], + probabilities=solution["probabilities"], + costs=solution["costs"], + ) + for solution in data["solutions"] + ] + return QUBODataset(coefficients=data["coefficients"], solutions=solutions) + + +def generate_symmetric_mask( + size: int, target: int, device: str, generator: torch.Generator +) -> torch.Tensor: + """Generate a symmetric boolean mask with an exact number of True elements + to match a certain density of QUBO. + Used in the `from_random` method of `QUBODataset`. + + Args: + size (int): Size of problem. + target (int): Target number of elements. + device (str): Torch device. + generator (torch.Generator): generator for randomness. + + Returns: + torch.Tensor: Mask. + """ + possible_x = [] + for x in range(1, min(size, target) + 1): + if (target - x) % 2 == 0: + y = (target - x) // 2 + if y <= (size * (size - 1)) // 2: + possible_x.append(x) + if not possible_x: + x, y = 1, 0 + else: + x = possible_x[ + int(torch.randint(0, len(possible_x), (1,), device=device, generator=generator).item()) + ] + y = (target - x) // 2 + + mask = torch.zeros((size, size), dtype=torch.bool, device=device) + diag_indices = torch.randperm(size, device=device, generator=generator)[:x] + for i in diag_indices.tolist(): + mask[i, i] = True + + upper_indices = torch.tensor( + [(i, j) for i in range(size) for j in range(i + 1, size)], + device=device, + ) + if upper_indices.size(0) > 0 and y > 0: + perm = torch.randperm(upper_indices.size(0), device=device, generator=generator)[:y] + chosen_upper = upper_indices[perm] + for i, j in chosen_upper.tolist(): + mask[i, j] = True + mask[j, i] = True + return mask diff --git a/qubosolver/qubo_types.py b/qubosolver/types/_enums.py similarity index 100% rename from qubosolver/qubo_types.py rename to qubosolver/types/_enums.py diff --git a/qubosolver/qubo_instance.py b/qubosolver/types/_instance.py similarity index 62% rename from qubosolver/qubo_instance.py rename to qubosolver/types/_instance.py index 93509910..8686ff41 100644 --- a/qubosolver/qubo_instance.py +++ b/qubosolver/types/_instance.py @@ -2,21 +2,19 @@ import torch import io -import qubosolver.io.utils as io_utils -from numpy.typing import ArrayLike +from typing import TYPE_CHECKING -from .data import QUBOSolution -from .data_utils import convert_to_tensor -from .qubo_types import DensityType -from .utils.density import ( - calculate_density, - classify_density, -) +from ._checks import debug_runtime_typecheck +from . import matrix +from ._solution import QUBOSolution +from ._enums import DensityType +from qubosolver._io import utils as io_utils -# Modules to be automatically added to the qubosolver namespace -__all__: list[str] = ["QUBOInstance"] +if TYPE_CHECKING: + from ._linalg import Matrix +@debug_runtime_typecheck class QUBOInstance: """ Represents a single instance of a Quadratic Unconstrained Binary Optimization (QUBO) problem. @@ -38,9 +36,7 @@ class QUBOInstance: def __init__( self, - coefficients: dict | ArrayLike | None = None, - device: str = "cpu", - dtype: torch.dtype = torch.float32, + coefficients: Matrix = matrix.zeros(0), ): """ Initializes a QUBOInstance. @@ -53,15 +49,11 @@ def __init__( dtype (torch.dtype): Data type of the tensors (default: torch.float32). """ - self._coefficients: torch.Tensor = torch.zeros([0, 0], dtype=dtype, device=device) - self.solution: QUBOSolution | None = None - self.density: float | None = None - self.density_type: DensityType | None = None - - if coefficients is None: - self.coefficients = self._coefficients - else: - self.coefficients = coefficients + self._coefficients: Matrix = coefficients + self.solution: QUBOSolution = QUBOSolution() + self.density: float = 0.0 + self.density_type: DensityType = DensityType.SPARSE + self.update_metrics() @property def size(self) -> int: @@ -97,48 +89,6 @@ def coefficients(self) -> torch.Tensor: ) return self._coefficients - @coefficients.setter - def coefficients(self, coeffs: dict | ArrayLike) -> None: - """ - Setter for the QUBO coefficient matrix, with validation checks. - - Exits the program with an error message if a check fails. - """ - # Convert input to tensor - tensor = convert_to_tensor(coeffs, device=self.device, dtype=self.dtype) # type: ignore[arg-type] - - # All checks passed, assign the tensor and update metrics - self._coefficients = tensor - self.update_metrics() - - def set_coefficients( - self, new_coefficients: dict[tuple[int, int], float] | None = None - ) -> None: - """ - Updates the coefficients of the QUBO problem. - - Args: - new_coefficients (dict[tuple[int, int], float] | None): - Dictionary of new coefficients to set. Keys are (row, column) tuples. - """ - if not new_coefficients: - return - - max_index = max(max(i, j) for i, j in new_coefficients.keys()) - if max_index >= self.size: - self._expand_size(max_index + 1) - - indices = torch.tensor(list(new_coefficients.keys()), dtype=torch.long, device=self.device) - values = torch.tensor(list(new_coefficients.values()), dtype=self.dtype, device=self.device) - self._coefficients[indices[:, 0], indices[:, 1]] = values - off_diagonal_mask = indices[:, 0] != indices[:, 1] - symmetric_indices = indices[off_diagonal_mask] - self._coefficients[symmetric_indices[:, 1], symmetric_indices[:, 0]] = values[ - off_diagonal_mask - ] - - self.update_metrics() - @property def _max_off_diag(self) -> float: mask = ~torch.eye(self.size, dtype=torch.bool, device=self.device) @@ -148,29 +98,12 @@ def _max_off_diag(self) -> float: def normalized_coefficients(self) -> torch.Tensor: return self.coefficients / self._max_off_diag - def _expand_size(self, new_size: int) -> None: - """ - Expands the size of the coefficient matrix to accommodate larger indices. - - Args: - new_size (int): - New size of the coefficient matrix. - """ - expanded_coefficients = torch.zeros( - (new_size, new_size), dtype=self.dtype, device=self.device - ) - expanded_coefficients[: self.size, : self.size] = self._coefficients - self._coefficients = expanded_coefficients - def update_metrics(self) -> None: """ Updates the density metrics of the QUBO problem. """ - if self.size > 0: - self.density = calculate_density(self._coefficients, self.size) - self.density_type = classify_density(self.density) - else: - self.density = self.density_type = None + self.density = calculate_density(self._coefficients) + self.density_type = classify_density(self.density) def evaluate_solution(self, solution: list | tuple | ArrayLike) -> float: """ @@ -247,3 +180,52 @@ def load(file_like: io_utils.FileLike[bytes]) -> QUBOInstance: Q = torch.load(buffer, weights_only=True) return QUBOInstance(Q) + + +# Density thresholds +_SPARSE_THRESHOLD: tuple[float, float] = (0.0, 0.3) +_MEDIUM_THRESHOLD: tuple[float, float] = (0.3, 0.7) +_HIGH_THRESHOLD: tuple[float, float] = (0.7, 1.0) + + +def classify_density(density: float) -> DensityType: + """ + Classifies the density of a QUBO problem based on predefined thresholds. + + Args: + density (float): + The density value to classify. Should be in the range [0.0, 1.0]. + + Returns: + DensityType: + The classification of the density (SPARSE, MEDIUM, or HIGH). + """ + if _SPARSE_THRESHOLD[0] <= density < _SPARSE_THRESHOLD[1]: + return DensityType.SPARSE + elif _MEDIUM_THRESHOLD[0] <= density < _MEDIUM_THRESHOLD[1]: + return DensityType.MEDIUM + elif _HIGH_THRESHOLD[0] <= density <= _HIGH_THRESHOLD[1]: + return DensityType.HIGH + else: + raise ValueError(f"Density {density} is outside the defined thresholds.") + + +def calculate_density(m: Matrix) -> float: + """ + Calculates the density of a QUBO coefficient matrix. + + Density is defined as the fraction of non-zero elements in the matrix. + + Args: + coefficients (torch.Tensor | None): + Tensor representing the QUBO coefficients. If None, density is 0. + size (int | None): + Size of the QUBO matrix (number of rows/columns). + + Returns: + float: + The density value, ranging from 0.0 (completely sparse) to 1.0 (completely dense). + """ + if m.numel() == 0: + return 0.0 + return torch.count_nonzero(m).item() / m.numel() diff --git a/qubosolver/types/_linalg.py b/qubosolver/types/_linalg.py new file mode 100644 index 00000000..afe1a3a2 --- /dev/null +++ b/qubosolver/types/_linalg.py @@ -0,0 +1,168 @@ +from __future__ import annotations + +import os +import torch +from typing import TYPE_CHECKING +from ._checks import _RUNTIME_CHECKS + +_FLOAT_DTYPE_MAP: dict[str, torch.dtype] = { + "float16": torch.float16, + "float32": torch.float32, + "float64": torch.float64, + "bfloat16": torch.bfloat16, +} + + +def _device_from_env() -> torch.device: + """Returns the torch device based on environment variables. + + Checks in order: + 1. QUBO_SOLVER_DEVICE — explicit device string (e.g. "cpu", "cuda", "cuda:1", "mps") + 2. USE_GPU — if set to "1" or "true", use "cuda"; otherwise "cpu" + 3. Defaults to "cpu" if neither is set. + """ + device_str = os.getenv("QUBO_SOLVER_DEVICE") + if device_str is not None: + try: + return torch.device(device_str) + except RuntimeError: + raise ValueError(f"Invalid QUBO_SOLVER_DEVICE={device_str!r}.") + + use_gpu = os.getenv("USE_GPU") + if use_gpu is not None: + return torch.device("cuda" if use_gpu.lower() in ("1", "true") else "cpu") + + return torch.device("cpu") + + +def _use_double_precision_from_env() -> bool: + """Returns whether double precision (float64) is enabled based on environment variables. + + Checks in order: + 1. QUBO_SOLVER_FLOAT_DTYPE — if set to "float64", returns True; "float32" returns False + 2. USE_DOUBLE_PRECISION — if set to "1" or "true", returns True; otherwise False + 3. Defaults to False if neither is set (i.e. float32 is the default, matching PyTorch). + """ + float_dtype_str = os.getenv("QUBO_SOLVER_FLOAT_DTYPE") + if float_dtype_str is not None: + return float_dtype_str == "float64" + + use_double = os.getenv("USE_DOUBLE_PRECISION") + if use_double is not None: + return use_double.lower() in ("1", "true") + + return False + + +def _float_type_from_env() -> torch.dtype: + """Returns the float dtype based on environment variables. + + Checks in order: + 1. QUBO_SOLVER_FLOAT_DTYPE — explicit dtype name (float16, bfloat16, float32, float64) + 2. USE_DOUBLE_PRECISION — if set to "1" or "true", use float64; otherwise float32 + 3. Defaults to float32 if neither is set (matching PyTorch default). + """ + float_dtype_str = os.getenv("QUBO_SOLVER_FLOAT_DTYPE") + if float_dtype_str is not None: + if float_dtype_str not in _FLOAT_DTYPE_MAP: + raise ValueError( + f"Invalid QUBO_SOLVER_FLOAT_DTYPE={float_dtype_str!r}. " + f"Valid options: {list(_FLOAT_DTYPE_MAP)}" + ) + return _FLOAT_DTYPE_MAP[float_dtype_str] + + return torch.float64 if _use_double_precision_from_env() else torch.float32 + + +class GlobalConfig: + _float_dtype = _float_type_from_env() + _device = _device_from_env() + + @classmethod + def use_double_precision(cls, enable: bool = True) -> None: + """Allows the user to easily toggle float64 on or off.""" + cls._float_dtype = torch.float64 if enable else torch.float32 + + @classmethod + def set_float_precision(cls, dtype: torch.dtype) -> None: + """Set the global float dtype. + + Args: + dtype: A torch float dtype. Valid options: float16, bfloat16, float32, float64. + """ + if dtype not in _FLOAT_DTYPE_MAP.values(): + raise ValueError( + f"Invalid dtype {dtype!r}. Valid options: {list(_FLOAT_DTYPE_MAP.values())}" + ) + cls._float_dtype = dtype + + @classmethod + def use_gpu(cls, enable: bool = True) -> None: + """Allows the user to easily toggle GPU (cuda) on or off.""" + cls._device = torch.device("cuda" if enable else "cpu") + + @classmethod + def set_device(cls, device: torch.device) -> None: + """Set the global torch device. + + Args: + device: A torch.device instance (e.g. torch.device("cuda:1")). + """ + cls._device = device + + +def dtype() -> torch.dtype: + """Returns the currently configured global float dtype.""" + return GlobalConfig._float_dtype + + +def device() -> torch.device: + """Returns the currently configured global torch device.""" + return GlobalConfig._device + + +if TYPE_CHECKING or _RUNTIME_CHECKS: + import jaxtyping + from typing_extensions import TypeAlias + + USE_DOUBLE_PRECISION: bool = _use_double_precision_from_env() + + Vectorf: TypeAlias = jaxtyping.Float32[torch.Tensor, "n"] # noqa: F821 + Matrixf: TypeAlias = jaxtyping.Float32[torch.Tensor, "n n"] # noqa: F821, F722 + Tensorf: TypeAlias = jaxtyping.Float32[torch.Tensor, "..."] # noqa: F821 + + Vectord: TypeAlias = jaxtyping.Float64[torch.Tensor, "n"] # noqa: F821 + Matrixd: TypeAlias = jaxtyping.Float64[torch.Tensor, "n n"] # noqa: F821, F722 + Tensord: TypeAlias = jaxtyping.Float64[torch.Tensor, "..."] # noqa: F821 + + Vectori: TypeAlias = jaxtyping.Int64[torch.Tensor, "n"] # noqa: F821 + + if USE_DOUBLE_PRECISION: + Vector = Vectord + Matrix = Matrixd + Tensor = Tensord + else: + Vector = Vectorf + Matrix = Matrixf + Tensor = Tensorf + + Bitstring: TypeAlias = jaxtyping.Int8[torch.Tensor, "n"] # noqa: F821 + Bitstrings: TypeAlias = jaxtyping.Int8[torch.Tensor, "n m"] # noqa: F821, F722 + +else: + Vectorf: TypeAlias = torch.Tensor + Matrixf: TypeAlias = torch.Tensor + Tensorf: TypeAlias = torch.Tensor + + Vectord: TypeAlias = torch.Tensor + Matrixd: TypeAlias = torch.Tensor + Tensord: TypeAlias = torch.Tensor + + Vectori: TypeAlias = torch.Tensor + + Vector: TypeAlias = torch.Tensor + Matrix: TypeAlias = torch.Tensor + Tensor: TypeAlias = torch.Tensor + + Bitstring: TypeAlias = torch.Tensor + Bitstrings: TypeAlias = torch.Tensor diff --git a/qubosolver/concepts/__init__.py b/qubosolver/types/_protocols/__init__.py similarity index 100% rename from qubosolver/concepts/__init__.py rename to qubosolver/types/_protocols/__init__.py diff --git a/qubosolver/types/_protocols/backend.py b/qubosolver/types/_protocols/backend.py new file mode 100644 index 00000000..5af21eba --- /dev/null +++ b/qubosolver/types/_protocols/backend.py @@ -0,0 +1,12 @@ +from __future__ import annotations + +from typing import Self +import qoolqit +from qoolqit.execution import job +from pulser.backend import Results +from .protocol import Protocol + + +class Backend(Protocol): + + def run(self: Self, program: qoolqit.QuantumProgram) -> job.Job[Results]: ... diff --git a/qubosolver/types/_protocols/protocol.py b/qubosolver/types/_protocols/protocol.py new file mode 100644 index 00000000..b7f6fb0e --- /dev/null +++ b/qubosolver/types/_protocols/protocol.py @@ -0,0 +1,11 @@ +from __future__ import annotations + +from typing import TYPE_CHECKING +from qubosolver.types.types import _RUNTIME_CHECKS + +if _RUNTIME_CHECKS and not TYPE_CHECKING: + from beartype.typing import Protocol as Protocol +else: + from typing import Protocol as Protocol + +__all__ = ["Protocol"] diff --git a/qubosolver/types/_solution.py b/qubosolver/types/_solution.py new file mode 100644 index 00000000..cc8cd266 --- /dev/null +++ b/qubosolver/types/_solution.py @@ -0,0 +1,116 @@ +from __future__ import annotations + +from dataclasses import dataclass +import torch +from typing import TYPE_CHECKING + +from ._checks import debug_runtime_typecheck +from . import vector, vectori +from . import bitstrings as _bitstrings +from ._enums import SolutionStatusType + +from qubosolver import _utils + +if TYPE_CHECKING: + from ._linalg import Bitstrings, Vector, Vectori + from ._instance import QUBOInstance + from pulser.backend.results import Results + + +@debug_runtime_typecheck +@dataclass +class QUBOSolution: + """ + Represents a solution to a QUBO problem. + + Attributes: + bitstrings (torch.Tensor): + Tensor of shape (num_solutions, bitstring_length), containing the bitstring solutions. + Each entry is an integer tensor with 0s and 1s. + counts (torch.Tensor | None): + Tensor of shape (num_solutions,), containing the count of occurrences of each bitstring. + Optional, can be None. + probabilities (torch.Tensor | None): + Tensor of shape (num_solutions,), containing the probability of each bitstring solution. + Optional, can be None. + costs (torch.Tensor): + Tensor of shape (num_solutions,), containing the cost associated with each + bitstring solution. + """ + + bitstrings: Bitstrings = _bitstrings.zeros(0, 0) + costs: Vector = vector.zeros(0) + counts: Vectori = vectori.zeros(0) + probabilities: Vector = vector.zeros(0) + solution_status: SolutionStatusType = SolutionStatusType.UNPROCESSED + + def empty(self) -> bool: + r = self.bitstrings.numel() == 0 + if not r: + return False + assert self.costs.numel() == 0 + assert self.counts.numel() == 0 + assert self.probabilities.numel() == 0 + + return True + + def __bool__(self) -> bool: + return not self.empty() + + def compute_costs(self, instance: QUBOInstance) -> Vector: + """ + Computes the cost for each bitstring solution based on the provided QUBO instance. + + Args: + instance (QUBOInstance): The QUBO instance containing the QUBO matrix. + + Returns: + torch.Tensor: A tensor of costs for each bitstring. + """ + dtype = instance.dtype + costs = _utils.costs.batched_qubo_cost(self.bitstrings.to(dtype), instance.coefficients) + return costs + + def compute_probabilities(self) -> Vector: + """ + Computes the probabilities of each bitstring solution based on their counts. + + Returns: + torch.Tensor: A tensor of probabilities for each bitstring. + """ + total_counts = self.counts.sum().item() + probabilities = ( + self.counts / total_counts if total_counts > 0 else torch.zeros_like(self.counts, dtype=vector.dtype()) + ) + return probabilities + + def sort_by_cost(self) -> None: + """ + Sorts the QUBOSolution in-place by increasing cost. + + Reorders bitstrings, costs, counts, and probabilities (if available) + based on the ascending order of the costs. + """ + + sorted_indices = torch.argsort(self.costs) + self.bitstrings = self.bitstrings[sorted_indices] + self.costs = self.costs[sorted_indices] + if self.counts.numel() > 0: + self.counts = self.counts[sorted_indices] + if self.probabilities.numel() > 0: + self.probabilities = self.probabilities[sorted_indices] + + @staticmethod + def from_results(results: Results) -> QUBOSolution: + counter = results.final_bitstrings + bitstrings = torch.tensor( + [list(map(int, list(b))) for b in list(counter.keys())], dtype=torch.int8 + ) + if bitstrings.numel() == 0: + bitstrings = torch.empty((0, 0), dtype=torch.int8) + counts = torch.tensor(list(map(int, list(counter.values()))), dtype=torch.int64) + + return QUBOSolution( + bitstrings=bitstrings, + counts=counts, + ) diff --git a/qubosolver/types/bitstring.py b/qubosolver/types/bitstring.py new file mode 100644 index 00000000..410f232c --- /dev/null +++ b/qubosolver/types/bitstring.py @@ -0,0 +1,23 @@ +from __future__ import annotations + +from typing import Any +import torch +from . import _linalg, vector + + +def dtype() -> torch.dtype: + return torch.int8 + + +def device() -> torch.device: + return _linalg.device() + + +def zeros(n: int, *, device: torch.device = device()) -> _linalg.Bitstring: + """Creates a vector of the given size with the specified dtype and device.""" + return vector.zeros(n, dtype=dtype(), device=device) + + +def tensor(data: Any, *, device: torch.device = device(), **kwargs) -> _linalg.Bitstring: + """Creates a vector from the given data.""" + return vector.tensor(data, dtype=dtype(), device=device, **kwargs) diff --git a/qubosolver/types/bitstrings.py b/qubosolver/types/bitstrings.py new file mode 100644 index 00000000..7881418a --- /dev/null +++ b/qubosolver/types/bitstrings.py @@ -0,0 +1,23 @@ +from __future__ import annotations + +from typing import Any +import torch +from . import _linalg, bitstring + + +def dtype() -> torch.dtype: + return bitstring.dtype() + + +def device() -> torch.device: + return bitstring.device() + + +def zeros(m: int, n: int, *, device: torch.device = device()) -> _linalg.Bitstrings: + """Creates a vector of the given size with the specified dtype and device.""" + return torch.zeros((m, n), dtype=dtype(), device=device) + + +def tensor(data: Any, *, device: torch.device = device(), **kwargs) -> _linalg.Bitstrings: + """Creates a vector from the given data.""" + return torch.tensor(data, dtype=dtype(), device=device, **kwargs) diff --git a/qubosolver/types/matrix.py b/qubosolver/types/matrix.py new file mode 100644 index 00000000..436c23ca --- /dev/null +++ b/qubosolver/types/matrix.py @@ -0,0 +1,27 @@ +from __future__ import annotations + +from typing import Any +import torch +from . import _linalg + + +def dtype() -> torch.dtype: + return _linalg.dtype() + + +def device() -> torch.device: + return _linalg.device() + + +def zeros( + n: int, *, dtype: torch.dtype = dtype(), device: torch.device = device() +) -> _linalg.Matrix: + """Creates a vector of the given size with the specified dtype and device.""" + return torch.zeros((n, n), dtype=dtype, device=device) + + +def tensor( + data: Any, *, dtype: torch.dtype = dtype(), device: torch.device = device(), **kwargs +) -> _linalg.Matrix: + """Creates a vector from the given data.""" + return torch.tensor(data, dtype=dtype, device=device, **kwargs) diff --git a/qubosolver/types/tensor.py b/qubosolver/types/tensor.py new file mode 100644 index 00000000..ee4ca687 --- /dev/null +++ b/qubosolver/types/tensor.py @@ -0,0 +1,27 @@ +from __future__ import annotations + +from typing import Any +import torch +from . import _linalg + + +def dtype() -> torch.dtype: + return _linalg.dtype() + + +def device() -> torch.device: + return _linalg.device() + + +def zeros( + *args, dtype: torch.dtype = dtype(), device: torch.device = device(), **kwargs +) -> _linalg.Tensor: + """Creates a vector of the given size with the specified dtype and device.""" + return torch.zeros(*args, dtype=dtype, device=device, **kwargs) + + +def tensor( + data: Any, *, dtype: torch.dtype = dtype(), device: torch.device = device(), **kwargs +) -> _linalg.Tensor: + """Creates a vector from the given data.""" + return torch.tensor(data, dtype=dtype, device=device, **kwargs) diff --git a/qubosolver/types/vector.py b/qubosolver/types/vector.py new file mode 100644 index 00000000..5a2fdeb8 --- /dev/null +++ b/qubosolver/types/vector.py @@ -0,0 +1,27 @@ +from __future__ import annotations + +from typing import Any +import torch +from . import _linalg + + +def dtype() -> torch.dtype: + return _linalg.dtype() + + +def device() -> torch.device: + return _linalg.device() + + +def zeros( + n: int, *, dtype: torch.dtype = dtype(), device: torch.device = device() +) -> _linalg.Vector: + """Creates a vector of the given size with the specified dtype and device.""" + return torch.zeros(n, dtype=dtype, device=device) + + +def tensor( + data: Any, *, dtype: torch.dtype = dtype(), device: torch.device = device(), **kwargs +) -> _linalg.Vector: + """Creates a vector from the given data.""" + return torch.tensor(data, dtype=dtype, device=device, **kwargs) diff --git a/qubosolver/types/vectori.py b/qubosolver/types/vectori.py new file mode 100644 index 00000000..fe27282e --- /dev/null +++ b/qubosolver/types/vectori.py @@ -0,0 +1,23 @@ +from __future__ import annotations + +from typing import Any +import torch +from . import _linalg, vector + + +def dtype() -> torch.dtype: + return torch.int64 + + +def device() -> torch.device: + return _linalg.device() + + +def zeros(n: int, *, device: torch.device = device()) -> _linalg.Vectori: + """Creates a vector of the given size with the specified dtype and device.""" + return vector.zeros(n, dtype=dtype(), device=device) + + +def tensor(data: Any, *, device: torch.device = device(), **kwargs) -> _linalg.Vectori: + """Creates a vector from the given data.""" + return vector.tensor(data, dtype=dtype(), device=device, **kwargs) diff --git a/qubosolver/utils/__init__.py b/qubosolver/utils/__init__.py deleted file mode 100644 index 247e209a..00000000 --- a/qubosolver/utils/__init__.py +++ /dev/null @@ -1,14 +0,0 @@ -from __future__ import annotations - -from .density import ( - calculate_density, - classify_density, -) -from .qubo_eval import calculate_qubo_cost - -# Modules to be automatically added to the qubosolver.utils namespace -__all__ = [ - "classify_density", - "calculate_density", - "calculate_qubo_cost", -] diff --git a/qubosolver/utils/density.py b/qubosolver/utils/density.py deleted file mode 100644 index 8b2a1114..00000000 --- a/qubosolver/utils/density.py +++ /dev/null @@ -1,56 +0,0 @@ -from __future__ import annotations - -import torch - -from qubosolver.qubo_types import DensityType - -# Density thresholds -SPARSE_THRESHOLD: tuple[float, float] = (0.0, 0.3) -MEDIUM_THRESHOLD: tuple[float, float] = (0.3, 0.7) -HIGH_THRESHOLD: tuple[float, float] = (0.7, 1.0) - - -def classify_density(density: float) -> DensityType: - """ - Classifies the density of a QUBO problem based on predefined thresholds. - - Args: - density (float): - The density value to classify. Should be in the range [0.0, 1.0]. - - Returns: - DensityType: - The classification of the density (SPARSE, MEDIUM, or HIGH). - """ - if SPARSE_THRESHOLD[0] <= density < SPARSE_THRESHOLD[1]: - return DensityType.SPARSE - elif MEDIUM_THRESHOLD[0] <= density < MEDIUM_THRESHOLD[1]: - return DensityType.MEDIUM - elif HIGH_THRESHOLD[0] <= density <= HIGH_THRESHOLD[1]: - return DensityType.HIGH - else: - raise ValueError(f"Density {density} is outside the defined thresholds.") - - -def calculate_density(coefficients: torch.Tensor | None, size: int | None) -> float: - """ - Calculates the density of a QUBO coefficient matrix. - - Density is defined as the fraction of non-zero elements in the matrix. - - Args: - coefficients (torch.Tensor | None): - Tensor representing the QUBO coefficients. If None, density is 0. - size (int | None): - Size of the QUBO matrix (number of rows/columns). - - Returns: - float: - The density value, ranging from 0.0 (completely sparse) to 1.0 (completely dense). - """ - if coefficients is None or coefficients.numel() == 0: - return 0.0 - - total_elements = size**2 if size else 0 - non_zero_elements = torch.count_nonzero(coefficients).item() - return float(non_zero_elements / total_elements) diff --git a/tests/conftest.py b/tests/conftest.py index c45c9304..eafa17f4 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -11,27 +11,22 @@ from pulser_simulation import QutipBackendV2 from emu_sv import SVBackend from emu_mps import MPSBackend +import qoolqit -from qoolqit.devices import DigitalAnalogDevice, AnalogDevice, Device, MockDevice -from pulser_pasqal import PasqalCloud -from qubosolver import QUBOInstance, QUBOSolution -from qubosolver.qubo_analyzer import QUBOAnalyzer -from qubosolver.config import ( +from mock.connection import MockConnection + +from qubosolver import ( + QUBOInstance, + QUBOSolution, + QUBOAnalyzer, + SolverConfig, + EmbedderType, + LayoutType, + DriveType, EmbeddingConfig, DriveShapingConfig, - SolverConfig, - LocalEmulator, ) -from qubosolver.qubo_types import EmbedderType, LayoutType, DriveType -from mock.connection import MockConnection - -from typing import TYPE_CHECKING - -if TYPE_CHECKING: - from typing import Generator -connection = PasqalCloud() -connection.fetch_available_devices() def pytest_collection_modifyitems(items: list[pytest.Item]) -> None: @@ -96,8 +91,8 @@ def classical_solver_config() -> SolverConfig: return SolverConfig(use_quantum=False) -locals_bkds: list[LocalEmulator] = [ - LocalEmulator(backend_type=btype, num_shots=500) +locals_bkds: list[qoolqit.execution.LocalEmulator] = [ + qoolqit.execution.LocalEmulator(backend_type=btype, num_shots=500) for btype in [ QutipBackendV2, SVBackend, @@ -109,25 +104,23 @@ def classical_solver_config() -> SolverConfig: @pytest.fixture( params=locals_bkds, ) -def local_backend(request: pytest.FixtureRequest) -> LocalEmulator: +def local_backend(request: pytest.FixtureRequest) -> qoolqit.execution.LocalEmulator: return request.param # type: ignore[no-any-return] @pytest.fixture( params=[ - DigitalAnalogDevice(), - AnalogDevice(), - MockDevice(), - Device(pulser_device=connection.fetch_available_devices()["FRESNEL"]), + qoolqit.DigitalAnalogDevice(), + qoolqit.AnalogDevice(), + qoolqit.MockDevice(), ], ids=[ "DigitalAnalogDevice", "AnalogDevice", "MockDevice", - "FRESNEL", ], ) -def local_device(request: pytest.FixtureRequest) -> Device: +def local_device(request: pytest.FixtureRequest) -> qoolqit.Device: return request.param # type: ignore[no-any-return] @@ -145,7 +138,7 @@ def embedding_method(request: pytest.FixtureRequest) -> EmbedderType: def qutip_solver_config() -> SolverConfig: return SolverConfig( use_quantum=True, - backend=LocalEmulator(backend_type=QutipBackendV2, num_shots=500), + backend=qoolqit.execution.LocalEmulator(backend_type=QutipBackendV2, num_shots=500), ) diff --git a/tests/test_driveshaping.py b/tests/drive_shaping/test_drive_shaper.py similarity index 92% rename from tests/test_driveshaping.py rename to tests/drive_shaping/test_drive_shaper.py index 51cc532a..60d80dca 100644 --- a/tests/test_driveshaping.py +++ b/tests/drive_shaping/test_drive_shaper.py @@ -13,7 +13,8 @@ get_drive_shaper, ) from qubosolver.qubo_instance import QUBOInstance -from qubosolver.qubo_types import DriveType, SolutionStatusType +from qubosolver.types.enums import DriveType, SolutionStatusType +from qubosolver.types import Bitstring @pytest.fixture @@ -80,16 +81,9 @@ def test_generate_optimized_drive_shaper( if isinstance(solution.counts, torch.Tensor): assert solution.counts.numel() > 0 - # try with custom objective_fn + # try with custom objective - def custom_ojective( - bitstrings: list, - counts: list, - probabilities: list, - costs: list, - best_cost: float, - best_bitstring: str, - ) -> float: + def custom_ojective(solution: QUBOSolution) -> float: return float(1e4) opt_res = [] @@ -97,7 +91,7 @@ def custom_ojective( def callback_fn(d: dict) -> None: opt_res.append(d) - def custom_qubo(bitstring: str, QUBO: torch.Tensor) -> float: + def custom_qubo(bitstring: Bitstring, QUBO: torch.Tensor) -> float: return 1.0 custom_fn_ps = DriveShapingConfig( @@ -113,9 +107,6 @@ def custom_qubo(bitstring: str, QUBO: torch.Tensor) -> float: backend, ) assert isinstance(shaper, OptimizedDriveShaper) - assert shaper.optimized_custom_objective_fn is not None - assert shaper.optimized_callback_objective is not None - assert shaper.optimized_custom_qubo_cost is not None drive, solution = shaper.generate(dummy_register) assert len(opt_res) > 0 assert opt_res[-1]["cost_eval"] == float(1e4) diff --git a/tests/drive_shaping/test_heuristic.py b/tests/drive_shaping/test_heuristic.py index 19a8187b..a5354d2f 100644 --- a/tests/drive_shaping/test_heuristic.py +++ b/tests/drive_shaping/test_heuristic.py @@ -1,7 +1,7 @@ from __future__ import annotations from dataclasses import dataclass -from typing import List, Iterable +from typing import Iterable import torch import itertools import numpy as np @@ -27,7 +27,7 @@ def to_solutions( bitstrings: Iterable[str | torch.Tensor], costs: Iterable[float] = itertools.repeat(float("inf")), probabilities: Iterable[float] = itertools.repeat(0.0), -) -> List[Solution]: +) -> list[Solution]: def to_string(b: str | torch.Tensor) -> str: if isinstance(b, torch.Tensor): return "".join(str(int(i)) for i in b) @@ -40,7 +40,7 @@ def to_string(b: str | torch.Tensor) -> str: def gather_optimal_solutions( data: Iterable[Solution], min_cost: float | None = None -) -> List[Solution]: +) -> list[Solution]: if min_cost is None: min_cost = min(d.cost for d in data) return [d for d in data if np.allclose(d.cost, min_cost)] diff --git a/tests/drive_shaping/test_optimized.py b/tests/drive_shaping/test_optimized.py index 6fb55360..1c42213a 100644 --- a/tests/drive_shaping/test_optimized.py +++ b/tests/drive_shaping/test_optimized.py @@ -8,15 +8,14 @@ import pytest import pytest_check as check from unittest.mock import MagicMock, patch -from typing import List, Iterable, Dict, Any +from typing import Iterable, Any from qoolqit.devices.device import DigitalAnalogDevice, AnalogDevice from qoolqit.register import Register from qoolqit.execution import LocalEmulator from qubosolver.pipeline.drive import OptimizedDriveShaper -from qubosolver.qubo_instance import QUBOInstance -from qubosolver.config import SolverConfig, DriveShapingConfig +from qubosolver import QUBOInstance, QUBOSolution, SolverConfig, DriveShapingConfig from qubosolver.qubo_analyzer import QUBOAnalyzer @@ -34,37 +33,29 @@ class Solution: def to_solutions( - bitstrings: Iterable[str | torch.Tensor], - costs: Iterable[float] = itertools.repeat(float("inf")), - probabilities: Iterable[float] = itertools.repeat(0.0), -) -> List[Solution]: - def to_string(b: str | torch.Tensor) -> str: - if isinstance(b, torch.Tensor): - return "".join(str(int(i)) for i in b) - if isinstance(b, str): - return b - raise ValueError() + solution: QUBOSolution, +) -> list[Solution]: + def to_string(b: torch.Tensor) -> str: + return "".join(str(int(i)) for i in b) - return [Solution(to_string(b), c, p) for b, c, p in zip(bitstrings, costs, probabilities)] + return [ + Solution(to_string(b), c.item(), p.item()) + for b, c, p in zip(solution.bitstrings, solution.costs, solution.probabilities) + ] def gather_optimal_solutions( data: Iterable[Solution], min_cost: float | None = None -) -> List[Solution]: +) -> list[Solution]: if min_cost is None: min_cost = min(d.cost for d in data) return [d for d in data if np.allclose(d.cost, min_cost)] def probability_based_ojective( - bitstrings: list, - counts: list, - probabilities: list, - costs: list, - best_cost: float, - best_bitstring: str, + solution: QUBOSolution, ) -> float: - optimal_solutions = gather_optimal_solutions(to_solutions(bitstrings, costs, probabilities)) + optimal_solutions = gather_optimal_solutions(to_solutions(solution)) check.is_not(optimal_solutions, []) min_cost = optimal_solutions[0].cost total_prob = sum(s.probability for s in optimal_solutions) @@ -135,9 +126,7 @@ def test_equilateral_triangular_qubo(seed: int, use_probability_based_objective: print(f"{analyzer.df}") assert isinstance(qubo_solution.probabilities, torch.Tensor) - optimal_solutions = gather_optimal_solutions( - to_solutions(qubo_solution.bitstrings, qubo_solution.costs, qubo_solution.probabilities) - ) + optimal_solutions = gather_optimal_solutions(to_solutions(qubo_solution)) check.is_not(optimal_solutions, []) min_cost = optimal_solutions[0].cost @@ -212,9 +201,7 @@ def test_triangular_qubo(seed: int, use_probability_based_objective: bool) -> No print(f"{analyzer.df}") assert isinstance(qubo_solution.probabilities, torch.Tensor) - optimal_solutions = gather_optimal_solutions( - to_solutions(qubo_solution.bitstrings, qubo_solution.costs, qubo_solution.probabilities) - ) + optimal_solutions = gather_optimal_solutions(to_solutions(qubo_solution)) check.is_not(optimal_solutions, []) min_cost = optimal_solutions[0].cost @@ -249,20 +236,12 @@ def test_errors(raise_exception: bool) -> None: register = Register.from_coordinates(vertices.tolist()) - def error( - bitstrings: list, - counts: list, - probabilities: list, - costs: list, - best_cost: float, - best_bitstring: str, - ) -> float: - + def error(solution: QUBOSolution) -> float: if raise_exception: raise RuntimeError("Error occurred") return float("inf") - def optimized_callback_objective(d: Dict[Any, Any]) -> None: + def optimized_callback_objective(d: dict[Any, Any]) -> None: check.almost_equal(d["cost_eval"], 1e4) mock_error = MagicMock(wraps=error) @@ -335,11 +314,11 @@ def test_failed_simulation_2() -> None: drive_shaper = OptimizedDriveShaper(QUBOInstance(Q), config, config.backend) with patch( - "qubosolver.pipeline.drive.OptimizedDriveShaper.run_simulation", - return_value=(None, None, None, None, None, None), + "qubosolver.drive_shaping.optimized._run_simulation", + return_value=QUBOSolution(), ): - with pytest.raises(RuntimeError, match="No solution found"): - drive, qubo_solution = drive_shaper.generate(register) + drive, qubo_solution = drive_shaper.generate(register) + check.is_true(qubo_solution.empty()) def test_failed_skopt() -> None: @@ -367,5 +346,5 @@ def test_failed_skopt() -> None: drive_shaper = OptimizedDriveShaper(QUBOInstance(Q), config, config.backend) with patch("qubosolver.pipeline.drive.gp_minimize", return_value=None): - with pytest.raises(RuntimeError, match="No solution found"): - drive, qubo_solution = drive_shaper.generate(register) + drive, qubo_solution = drive_shaper.generate(register) + check.is_true(qubo_solution.empty()) diff --git a/tests/embeddings/test_embedding.py b/tests/embedding/test_embedder.py similarity index 100% rename from tests/embeddings/test_embedding.py rename to tests/embedding/test_embedder.py diff --git a/tests/embeddings/test_greedy.py b/tests/embedding/test_greedy.py similarity index 98% rename from tests/embeddings/test_greedy.py rename to tests/embedding/test_greedy.py index e1730861..ea35e19b 100644 --- a/tests/embeddings/test_greedy.py +++ b/tests/embedding/test_greedy.py @@ -5,8 +5,8 @@ import numpy as np from qoolqit.devices.device import DigitalAnalogDevice -from qubosolver.algorithms.greedy.greedy import Greedy -from qubosolver.qubo_types import LayoutType +from qubosolver.embedding.algorithms.greedy.greedy import Greedy +from qubosolver.types.enums import LayoutType from qubosolver.data import QUBODataset import pytest import pytest_check as check diff --git a/tests/embeddings/test_greedy_animation.py b/tests/embedding/test_greedy_animation.py similarity index 86% rename from tests/embeddings/test_greedy_animation.py rename to tests/embedding/test_greedy_animation.py index c3eecf67..ffe05e42 100644 --- a/tests/embeddings/test_greedy_animation.py +++ b/tests/embedding/test_greedy_animation.py @@ -2,13 +2,13 @@ from __future__ import annotations -from typing import Any, Dict, List, Optional +from typing import Any, Optional import torch from qoolqit.devices.device import DigitalAnalogDevice -from qubosolver.algorithms.greedy.greedy import Greedy -from qubosolver.qubo_types import LayoutType +from qubosolver.embedding.algorithms.greedy.greedy import Greedy +from qubosolver.types.enums import LayoutType def _toy_qubo() -> torch.Tensor: @@ -25,7 +25,7 @@ def _toy_qubo() -> torch.Tensor: ) -def _base_params(n: int) -> Dict[str, Any]: +def _base_params(n: int) -> dict[str, Any]: # Paramètres de base, sans animation par défaut return { "layout": LayoutType.TRIANGULAR, @@ -63,11 +63,11 @@ def test_greedy_animation_calls_renderer(monkeypatch: Any) -> None: params["animation_save_path"] = None # ne rien écrire # Compteur d’appels - called: Dict[str, int] = {"count": 0} + called: dict[str, int] = {"count": 0} def _fake_render( self: Greedy, - frames: List[Dict[str, Any]], + frames: list[dict[str, Any]], all_coords_np: Any, spacing: float, layout_name: str, @@ -79,7 +79,7 @@ def _fake_render( return None # 1) S'assurer que le code rentre bien dans le bloc d'animation - import qubosolver.algorithms.greedy.greedy as greedy_mod + import qubosolver.embedding.algorithms.greedy.greedy as greedy_mod monkeypatch.setattr(greedy_mod, "_VIZ_OK", True, raising=False) diff --git a/tests/integration/functional.py b/tests/integration/functional.py new file mode 100644 index 00000000..15a8e4b7 --- /dev/null +++ b/tests/integration/functional.py @@ -0,0 +1,248 @@ +from __future__ import annotations + +from dataclasses import dataclass +import itertools +import numpy as np +import pytest +import pytest_check as check +import random +import torch +from typing import Iterable + +import qoolqit + +from qubosolver import QUBOInstance, QUBOSolution +from qubosolver.qubo_analyzer import QUBOAnalyzer +from qubosolver import solvers, transforms, embedding, drive_shaping +from qubosolver.types import Bitstring + + +@dataclass +class Solution: + bitstring: str + cost: float = float("inf") + probability: float = 0.0 + + +def to_solutions( + bitstrings: Iterable[str | torch.Tensor], + costs: Iterable[float] = itertools.repeat(float("inf")), + probabilities: Iterable[float] = itertools.repeat(0.0), +) -> list[Solution]: + def to_string(b: str | torch.Tensor) -> str: + if isinstance(b, torch.Tensor): + return "".join(str(int(i)) for i in b) + if isinstance(b, str): + return b + raise ValueError() + + return [Solution(to_string(b), c, p) for b, c, p in zip(bitstrings, costs, probabilities)] + + +def gather_optimal_solutions( + data: Iterable[Solution], min_cost: float | None = None +) -> list[Solution]: + if min_cost is None: + min_cost = min(d.cost for d in data) + return [d for d in data if np.allclose(d.cost, min_cost)] + + +def interaction_matrix_from_vertices(vertices: torch.Tensor) -> torch.Tensor: + U = 1.0 / torch.cdist(vertices, vertices) ** 6 + U.fill_diagonal_(0.0) + return U + + +def simple_qubo() -> tuple[QUBOInstance, list[Solution]]: + + sqrt3 = np.sqrt(3.0) + vertices = torch.tensor( + [ + [0.0, 0.0], + [-1.0, 0.0], + [-1.5, -0.5 * sqrt3], + [-0.5, -0.5 * sqrt3], + [4.0, 0.0], + ], + dtype=torch.float32, + ) + n_qubits = vertices.shape[0] + diagonal_scale = -2.0 + diagonal = torch.ones(n_qubits, dtype=torch.float32) + Q = interaction_matrix_from_vertices(vertices) + diagonal_scale * torch.diag(diagonal) + Q /= Q.max() + + results = [] + for bits in itertools.product([0, 1], repeat=n_qubits): + z = torch.tensor(bits, dtype=torch.float32) + cost = (z @ Q @ z).item() + results.append(Solution("".join(str(int(b)) for b in z.flatten()), cost)) + + # Get all bitstrings with minimum cost + expected_optimal_solutions = gather_optimal_solutions(results) + check.is_not(expected_optimal_solutions, []) + + print(f"\nExpected Minimum cost: {expected_optimal_solutions[0].cost}") + print(f"All expected optimal bitstrings: {[s.bitstring for s in expected_optimal_solutions]}") + print(f"Number of expected optimal solutions: {len(expected_optimal_solutions)}\n") + + return QUBOInstance(coefficients=Q), expected_optimal_solutions + + +def manual_seed(seed: int) -> None: + np.random.seed(seed) + torch.manual_seed(seed) + random.seed(seed) + + +def check_solution(qubo: QUBOInstance, expected_optimal_solutions: list[Solution]) -> None: + + solution = qubo.solution + solution.costs = solution.compute_costs(qubo) + solution.probabilities = solution.compute_probabilities() + solution.sort_by_cost() + + # Solutions are not duplicated + check.equal(solution.bitstrings.unique(dim=0).shape[0], solution.bitstrings.shape[0]) + + analyzer = QUBOAnalyzer([solution]) + print(f"{analyzer.df}") + + assert isinstance(solution.probabilities, torch.Tensor) + optimal_solutions = gather_optimal_solutions( + to_solutions(solution.bitstrings, solution.costs, solution.probabilities) + ) + check.is_not(optimal_solutions, []) + + min_cost = optimal_solutions[0].cost + + print(f"\nMinimum cost: {min_cost}") + print(f"All optimal bitstrings: {[s.bitstring for s in optimal_solutions]}") + print(f"Number of optimal solutions: {len(optimal_solutions)}\n") + + check.almost_equal(min_cost, expected_optimal_solutions[0].cost) + expected_optimal_bistrings = [s.bitstring for s in expected_optimal_solutions] + for s in optimal_solutions: + check.is_in(s.bitstring, expected_optimal_bistrings) + + cumulated_probability = sum(s.probability for s in optimal_solutions) + check.greater(cumulated_probability, 0.75) + + +@pytest.mark.usefixtures("restore_rng_state") +@pytest.mark.parametrize("drive_shaping_method", ["heuristic", "optimized"]) +@pytest.mark.parametrize("embedding_method", ["greedy", "blade"]) +@pytest.mark.parametrize("postprocessing", [True, False], ids=["post", "no-post"]) +@pytest.mark.parametrize("preprocessing", [True, False], ids=["pre", "no-pre"]) +def test_quantum_solve( + preprocessing: bool, + postprocessing: bool, + embedding_method: str, + drive_shaping_method: str, +) -> None: + + manual_seed(16844214) + qubo, expected_optimal_solutions = simple_qubo() + + device = qoolqit.AnalogDevice() + + # Try trivial solution (none here) + trivial_solution = solvers.trivial_solution_search(qubo) + check.is_false(trivial_solution) + + effective_qubo = qubo + if preprocessing: + effective_qubo = transforms.variable_fixing.apply_recursively(qubo) + + if embedding_method == "blade": + blade_config = embedding.blade.Config(device=device) + register = embedding.blade.embed(effective_qubo, blade_config, normalize=True) + elif embedding_method == "greedy": + greedy_config = embedding.greedy.Config(traps=100, spacing=0.1) + register = embedding.greedy.embed(effective_qubo, device, greedy_config, normalize=True) + else: + raise ValueError(f"Invalid embedding method: {embedding_method}") + print(f"Register: {register.qubits}") + print(f"Distances: {register.distances()}") + + emulator = qoolqit.execution.LocalEmulator() + + if drive_shaping_method == "heuristic": + drive = drive_shaping.heuristic.generate( + register, effective_qubo, device, dmm=False, kappa=0.5 + ) + elif drive_shaping_method == "optimized": + drive, _ = drive_shaping.optimized.generate( + register, + effective_qubo, + device, + dmm=False, + backend=emulator, + config=drive_shaping.optimized.Config(n_calls=11), + ) + else: + raise ValueError(f"Invalid drive shaping method: {drive_shaping_method}") + + job = solvers.analog_quantum_sample(emulator, device, drive, register) + effective_qubo.solution = QUBOSolution.from_results(job.results()) + + # Post-process fixations of the preprocessing and restore the original QUBO + if preprocessing: + assert isinstance(effective_qubo, transforms.variable_fixing.QUBOInstance) + qubo.solution = transforms.variable_fixing.unapply(effective_qubo) + else: + qubo.solution = effective_qubo.solution + + if postprocessing: + qubo.solution = solvers.iterative_bitflip_local_search(qubo.solution, qubo) + + check_solution(qubo, expected_optimal_solutions) + + +@pytest.mark.usefixtures("restore_rng_state") +@pytest.mark.parametrize("solving_method", ["cplex", "tabu", "sa", "sa+tabu"]) +@pytest.mark.parametrize("postprocessing", [True, False], ids=["post", "no-post"]) +@pytest.mark.parametrize("preprocessing", [True, False], ids=["pre", "no-pre"]) +def test_classical_solve( + preprocessing: bool, + postprocessing: bool, + solving_method: str, +) -> None: + + seed = 16844214 + manual_seed(seed) + qubo, expected_optimal_solutions = simple_qubo() + + # Try trivial solution (none here) + trivial_solution = solvers.trivial_solution_search(qubo) + check.is_false(trivial_solution) + + effective_qubo = qubo + if preprocessing: + effective_qubo = transforms.variable_fixing.apply_recursively(qubo) + + if solving_method == "cplex": + effective_qubo.solution = solvers.cplex(effective_qubo) + elif solving_method == "tabu": + x0: Bitstring = torch.randint(0, 2, size=(effective_qubo.size,)).to(torch.int8) + effective_qubo.solution = solvers.tabu_search(effective_qubo, x0) + elif solving_method == "sa": + effective_qubo.solution = solvers.simulated_annealing(effective_qubo, seed=seed, top_k=1) + elif solving_method == "sa+tabu": + effective_qubo.solution = solvers.simulated_annealing(effective_qubo, seed=seed, top_k=1) + effective_qubo.solution = solvers.tabu_search( + effective_qubo, effective_qubo.solution.bitstrings[0] + ) + else: + raise ValueError(f"Invalid solving method: {solving_method}") + + if preprocessing: + assert isinstance(effective_qubo, transforms.variable_fixing.QUBOInstance) + qubo.solution = transforms.variable_fixing.unapply(effective_qubo) + else: + qubo.solution = effective_qubo.solution + + if postprocessing: + qubo.solution = solvers.iterative_bitflip_local_search(qubo.solution, qubo) + + check_solution(qubo, expected_optimal_solutions) diff --git a/tests/mock/connection.py b/tests/mock/connection.py index 280fdcb9..83a1773a 100644 --- a/tests/mock/connection.py +++ b/tests/mock/connection.py @@ -13,13 +13,13 @@ if TYPE_CHECKING: from pulser import Sequence as PulserSequence from pulser.backend.results import Results - from typing import Any, Dict, Mapping + from typing import Any, Mapping class MockConnection(RemoteConnection): def __init__(self, result: Results, running_iterations: int = 0) -> None: self.result: Results = result - self.results: Dict[str, Results] = dict() + self.results: dict[str, Results] = dict() self._status_calls = 0 self._support_open_batch = True diff --git a/tests/test_config.py b/tests/pipeline/test_config.py similarity index 99% rename from tests/test_config.py rename to tests/pipeline/test_config.py index d4f4f10e..ddaae1ec 100644 --- a/tests/test_config.py +++ b/tests/pipeline/test_config.py @@ -14,7 +14,7 @@ DecompositionConfig, SolverConfig, ) -from qubosolver.qubo_types import ( +from qubosolver.types.enums import ( EmbedderType, LayoutType, DriveType, diff --git a/tests/test_postprocessing.py b/tests/pipeline/test_postprocessing.py similarity index 71% rename from tests/test_postprocessing.py rename to tests/pipeline/test_postprocessing.py index 23f480a2..780c275f 100644 --- a/tests/test_postprocessing.py +++ b/tests/pipeline/test_postprocessing.py @@ -8,8 +8,9 @@ from qubosolver.config import SolverConfig from qubosolver.qubo_instance import QUBOInstance from qubosolver.qubo_analyzer import QUBOAnalyzer -from qubosolver.pipeline.fixtures import Fixtures, bit_flip_local_search from qubosolver.data import QUBODataset +from qubosolver.solvers.bitflip import _bit_flip_local_search +from qubosolver import solvers @pytest.mark.parametrize("postprocessing", [True, False]) @@ -46,8 +47,7 @@ def test_basic_qubo_2d_integration(postprocessing: bool) -> None: print(f"\n{df}") -@pytest.mark.parametrize("postprocessing", [True, False]) -def test_basic_qubo_2d_fixture(postprocessing: bool) -> None: +def test_basic_qubo_2d() -> None: # fmt: off Q = torch.tensor([ @@ -57,23 +57,16 @@ def test_basic_qubo_2d_fixture(postprocessing: bool) -> None: # fmt: on instance = QUBOInstance(coefficients=Q) - fixture = Fixtures(instance, SolverConfig(do_postprocessing=postprocessing)) solution = QUBOSolution(bitstrings=torch.tensor([[0, 0]]), costs=torch.tensor([0.0])) - pp_solution = fixture.postprocess(solution) + pp_solution = solvers.iterative_bitflip_local_search(solution, instance) pp_solution.sort_by_cost() pp_solution.bitstrings = pp_solution.bitstrings.int() - if postprocessing: - torch.testing.assert_close( - pp_solution.bitstrings[0, :], torch.tensor([1, 1], dtype=torch.int32) - ) - torch.testing.assert_close(pp_solution.costs, torch.tensor([-18.0])) - else: - torch.testing.assert_close( - pp_solution.bitstrings[0, :], torch.tensor([0, 0], dtype=torch.int32) - ) - torch.testing.assert_close(pp_solution.costs, torch.tensor([0.0])) + torch.testing.assert_close( + pp_solution.bitstrings[0, :], torch.tensor([1, 1], dtype=torch.int32) + ) + torch.testing.assert_close(pp_solution.costs, torch.tensor([-18.0])) analyzer = QUBOAnalyzer(pp_solution) df = analyzer.df @@ -81,9 +74,8 @@ def test_basic_qubo_2d_fixture(postprocessing: bool) -> None: @pytest.mark.usefixtures("restore_rng_state") -@pytest.mark.parametrize("postprocessing", [True, False]) @pytest.mark.parametrize("density", [0.2, 0.5, 0.8]) -def test_random_qubos(postprocessing: bool, density: float) -> None: +def test_random_qubos(density: float) -> None: size = 5 @@ -94,9 +86,8 @@ def test_random_qubos(postprocessing: bool, density: float) -> None: instance = QUBOInstance(coefficients=Q) bitstring = (torch.rand(size) > 0.5).int() cost = instance.evaluate_solution(bitstring) - fixture = Fixtures(instance, SolverConfig(do_postprocessing=postprocessing)) solution = QUBOSolution(bitstrings=bitstring.unsqueeze(0), costs=torch.tensor([cost])) - pp_solution = fixture.postprocess(solution) + pp_solution = solvers.iterative_bitflip_local_search(solution, instance) pp_solution.sort_by_cost() pp_solution.bitstrings = pp_solution.bitstrings.int() @@ -104,11 +95,7 @@ def test_random_qubos(postprocessing: bool, density: float) -> None: df = analyzer.df print(f"\n{df}") - if postprocessing: - check.less_equal(pp_solution.costs[0], cost) - else: - check.almost_equal(pp_solution.costs[0], cost) - torch.testing.assert_close(pp_solution.bitstrings[0, :], bitstring) + check.less_equal(pp_solution.costs[0], cost) def test_no_solution() -> None: @@ -120,12 +107,11 @@ def test_no_solution() -> None: # fmt: on instance = QUBOInstance(coefficients=Q) - fixture = Fixtures(instance, SolverConfig(do_postprocessing=True)) solution = QUBOSolution(bitstrings=torch.tensor([]), costs=torch.tensor([])) check.equal(solution.bitstrings.numel(), 0) - # Post-processing doesn't find new solutions if there were none to begin with - pp_solution = fixture.postprocess(solution) + # Bitflip doesn't find new solutions if there were none to begin wit + pp_solution = solvers.iterative_bitflip_local_search(solution, instance) check.equal(pp_solution.bitstrings.numel(), 0) @@ -145,7 +131,7 @@ def cost_function(bitstring: np.ndarray) -> float: initial_cost = cost_function(s) check.almost_equal(initial_cost, 0.0) - best_bitstring, best_cost = bit_flip_local_search(cost_function, s, shuffle=shuffle) + best_bitstring, best_cost = _bit_flip_local_search(cost_function, s, shuffle=shuffle) np.testing.assert_allclose(best_bitstring, np.array([1, 1])) check.almost_equal(best_cost, -18.0) @@ -168,5 +154,5 @@ def cost_function(bitstring: np.ndarray) -> float: return float(bitstring.T @ Q.numpy() @ bitstring) initial_cost = cost_function(s) - _, best_cost = bit_flip_local_search(cost_function, s, shuffle=shuffle) + _, best_cost = _bit_flip_local_search(cost_function, s, shuffle=shuffle) check.less_equal(best_cost, initial_cost) diff --git a/tests/test_preprocessing.py b/tests/pipeline/test_preprocessing.py similarity index 82% rename from tests/test_preprocessing.py rename to tests/pipeline/test_preprocessing.py index eee4d958..0e1c5450 100644 --- a/tests/test_preprocessing.py +++ b/tests/pipeline/test_preprocessing.py @@ -4,6 +4,7 @@ import pytest import pytest_check as check import torch +from copy import deepcopy from qoolqit import Drive, Ramp, Register, Constant from qubosolver import QUBOInstance, QUBOSolution @@ -14,13 +15,10 @@ SolverConfig, LocalEmulator, ) -from qubosolver.pipeline.fixtures import ( - Fixtures, - hansen_fixing, -) -from qubosolver.qubo_types import SolutionStatusType, EmbedderType +from qubosolver.types.enums import SolutionStatusType, EmbedderType from qubosolver.solver import QuboSolver from qubosolver.pipeline.drive import BaseDriveShaper +from qubosolver import transforms def test_apply_full_and_post_process_fixation() -> None: @@ -30,20 +28,15 @@ def test_apply_full_and_post_process_fixation() -> None: dtype=torch.int32, ) - qubo = QUBOInstance(matrix) - - config = SolverConfig(do_preprocessing=True) - fix_class = Fixtures(qubo, config) - fix_class.apply_full_fixation_exhaust() + full_qubo = QUBOInstance(matrix) + reduced_qubo = transforms.variable_fixing.apply_recursively(full_qubo) - assert fix_class.fixed_var_dict_list == [{0: 1, 3: 1}, {0: 0, 1: 0}] - assert fix_class.n_fixed_variables > 0 + assert reduced_qubo._fixed_indices == [{0: 1, 3: 1}, {0: 0, 1: 0}] + assert reduced_qubo.n_fixed_indices == 4 - assert isinstance(fix_class.reduced_qubo, QUBOInstance) + reduced_qubo.solution = QUBOSolution(torch.empty(0, 0), torch.empty(0, 0)) - sol = QUBOSolution(torch.empty(0, 0), torch.empty(0, 0)) - - sol_reconstructed = fix_class.post_process_fixation(sol) + sol_reconstructed = transforms.variable_fixing.unapply(reduced_qubo) assert isinstance(sol_reconstructed, QUBOSolution) @@ -57,7 +50,7 @@ def test_hansen_fixing() -> None: qubo_reducible = QUBOInstance(matrix_reducible) - fixed_var = hansen_fixing(qubo_reducible) + fixed_var = transforms.variable_fixing.hansen_fixing(qubo_reducible) assert isinstance(fixed_var, dict) @@ -67,7 +60,7 @@ def test_hansen_fixing() -> None: qubo_reducible = QUBOInstance(matrix_not_reducible) - empty_fixed_var = hansen_fixing(qubo_reducible) + empty_fixed_var = transforms.variable_fixing.hansen_fixing(qubo_reducible) assert empty_fixed_var == {} @@ -76,50 +69,34 @@ def test_reduce_qubo() -> None: matrix = torch.tensor( [[-98, 2, 13, 1], [2, -12, 20, 15], [13, 20, -34, 7], [1, 15, 7, -57]], - dtype=torch.int32, + dtype=torch.float32, ) - qubo = QUBOInstance(matrix) - config = SolverConfig(do_preprocessing=True) - fix_class = Fixtures(qubo, config) - - fix_class.reduce_qubo({0: 1, 3: 1}) - - assert torch.equal( - fix_class.reduced_qubo.coefficients, - torch.tensor([[22, 20], [20, 6]], dtype=torch.int32), - ) + reduced_qubo = deepcopy(qubo) - fix_class_not_reduced = Fixtures(qubo, config) + reduced_qubo = transforms.variable_fixing.reduce_qubo(qubo, {0: 1, 3: 1}) - fix_class_not_reduced.reduce_qubo({}) + expected_coefficients = torch.tensor([[22, 20], [20, 6]], dtype=torch.float32) + torch.testing.assert_close(reduced_qubo.coefficients, expected_coefficients) - assert torch.equal( - fix_class_not_reduced.reduced_qubo.coefficients, - torch.tensor( - [[-98, 2, 13, 1], [2, -12, 20, 15], [13, 20, -34, 7], [1, 15, 7, -57]], - dtype=torch.int32, - ), - ) + non_reduced_qubo = transforms.variable_fixing.reduce_qubo(qubo, {}) + torch.testing.assert_close(non_reduced_qubo.coefficients, qubo.coefficients) def test_apply_rule() -> None: + matrix = torch.tensor( [[-98, 2, 13, 1], [2, -12, 20, 15], [13, 20, -34, 7], [1, 15, 7, -57]], - dtype=torch.int32, + dtype=torch.float32, ) - qubo = QUBOInstance(matrix) - config = SolverConfig(do_preprocessing=True) - fix_class = Fixtures(qubo, config) - - fix_class.apply_rule(hansen_fixing) - - assert torch.equal( - fix_class.reduced_qubo.coefficients, - torch.tensor([[22, 20], [20, 6]], dtype=torch.int32), + reduced_qubo = transforms.variable_fixing.apply( + qubo, [transforms.variable_fixing.hansen_fixing] ) + expected_coefficients = torch.tensor([[22, 20], [20, 6]], dtype=torch.float32) + torch.testing.assert_close(reduced_qubo.coefficients, expected_coefficients) + def test_classical_unprocessed(qubo_instance_for_preprocessing: QUBOInstance) -> None: """ @@ -244,18 +221,18 @@ def test_reduce_qubo_2() -> None: ) qubo = QUBOInstance(matrix) - config = SolverConfig(do_preprocessing=True) - fix_class = Fixtures(qubo, config) - fix_class.preprocess() - check.equal(fix_class.fixed_var_dict_list, [{0: 0}, {2: 1, 3: 1}]) + + reduced_qubo = transforms.variable_fixing.apply_recursively(qubo) + + check.equal(reduced_qubo._fixed_indices, [{0: 0}, {2: 1, 3: 1}]) # Hardcode solution - reduced_solution = QUBOSolution( + reduced_qubo.solution = QUBOSolution( bitstrings=torch.tensor([[0, 1]]), costs=torch.tensor([0.0]), ) - solution = fix_class.post_process_fixation(reduced_solution) + solution = transforms.variable_fixing.unapply(reduced_qubo) check.equal(solution.bitstrings.shape, (1, 5)) bitstring = QUBOAnalyzer.tensor_to_bitstrings(solution.bitstrings.to(torch.int64))[0] @@ -305,7 +282,7 @@ def test_quantum_prepostprocessing_2( ] ) - instance = QUBOInstance(Q) + instance = QUBOInstance(torch.from_numpy(Q)) config = SolverConfig(use_quantum=True, do_preprocessing=preprocessing) config.embedding = EmbeddingConfig( diff --git a/tests/test_batch_id.py b/tests/pipeline/test_remote_job.py similarity index 93% rename from tests/test_batch_id.py rename to tests/pipeline/test_remote_job.py index 5a334663..dcc03799 100644 --- a/tests/test_batch_id.py +++ b/tests/pipeline/test_remote_job.py @@ -1,10 +1,11 @@ from __future__ import annotations +import io import pytest import pytest_check as check import numpy as np import torch -import io + from pulser.backend.remote import RemoteConnection from pulser.backend.results import Results @@ -17,7 +18,7 @@ LocalEmulator, RemoteEmulator, ) -from qubosolver.qubo_types import EmbedderType, DriveType +from qubosolver.types.enums import EmbedderType, DriveType from qubosolver.solver import QUBOInstance, QuboSolverQuantum, QUBOSolution import qubosolver.io.utils as io_utils @@ -30,7 +31,7 @@ @pytest.mark.parametrize("embedding_method", list(EmbedderType)) @pytest.mark.parametrize("preprocessing", [True, False], ids=["pre", "no_pre"]) @pytest.mark.parametrize("dmm", [True, False], ids=["dmm", "no_dmm"]) -def test_quantum_batch_id( +def test_quantum_remote_job( make_mock_connection: type[MockConnection], drive_method: str, embedding_method: str, @@ -55,7 +56,7 @@ def test_quantum_batch_id( def pre( connection: RemoteConnection | None = None, ) -> tuple[job.Job[Results], QuboSolverQuantum]: - instance = QUBOInstance(Q) + instance = QUBOInstance(torch.from_numpy(Q)) min_distance = 1.001 if drive_method == DriveType.HEURISTIC else None @@ -88,14 +89,7 @@ def pre( return job, solver def post(job: job.Job[Results], solver: QuboSolverQuantum) -> QUBOSolution: - bitstrings, counts = QuboSolverQuantum.parse_results(job.results()) - - solution = QUBOSolution( - bitstrings=bitstrings.float(), - counts=counts, - costs=torch.Tensor(), - probabilities=None, - ) + solution = QUBOSolution.from_results(job.results()) # Post-process fixations of the preprocessing and restore the original QUBO solution = solver.post_process_fixation(solution) diff --git a/tests/solvers/test_heuristics.py b/tests/solvers/test_classical_solver.py similarity index 100% rename from tests/solvers/test_heuristics.py rename to tests/solvers/test_classical_solver.py diff --git a/tests/solvers/test_decomposition.py b/tests/solvers/test_decomposition.py index bd150d57..adf5446b 100644 --- a/tests/solvers/test_decomposition.py +++ b/tests/solvers/test_decomposition.py @@ -6,7 +6,6 @@ import itertools import numpy as np import random -from typing import Tuple, List from qoolqit.devices import DigitalAnalogDevice from qoolqit import Register @@ -15,7 +14,7 @@ from qubosolver.config import SolverConfig, DecompositionConfig, EmbeddingConfig from qubosolver.solver import DecomposeQuboSolver, QuboSolver from qubosolver.data import QUBODataset -from qubosolver.algorithms.decompose import compute_distance_interaction_matrix +from qubosolver.transforms.algorithms.decompose import compute_distance_interaction_matrix @pytest.mark.priority(120) @@ -32,7 +31,7 @@ def test_initial_steps_solver(decomposable_qubo: QUBOInstance, use_quantum: bool seed = 79450 random.seed(seed) - from qubosolver.algorithms.decompose import ( + from qubosolver.transforms.algorithms.decompose import ( compute_distance_interaction_matrix, geometric_search, interaction_matrix_from_placed, @@ -233,7 +232,7 @@ def test_compute_distance_interaction_diagonal() -> None: @pytest.mark.usefixtures("restore_rng_state") @pytest.mark.parametrize("dims", [(4,), (3,), (3, 3), (2, 3, 2), (4, 3, 2, 3)]) @pytest.mark.parametrize("seed", [1935225697, 1547, 66987, 55571, 998618750]) -def test_decompose_and_solve_block_qubo(seed: int, dims: Tuple[int]) -> None: +def test_decompose_and_solve_block_qubo(seed: int, dims: tuple[int]) -> None: """Test that the decomposition solver correctly identifies and solves block-diagonal QUBO matrices. The test constructs a block-diagonal QUBO matrix from smaller sub-problems, runs the @@ -255,7 +254,7 @@ def test_decompose_and_solve_block_qubo(seed: int, dims: Tuple[int]) -> None: Args: seed (int): Random seed for reproducibility (controls ``random``, ``torch``, and ``numpy``). - dims (Tuple[int, ...]): Dimensions of the individual QUBO blocks that form the + dims (tuple[int, ...]): Dimensions of the individual QUBO blocks that form the block-diagonal matrix. """ @@ -342,7 +341,7 @@ def test_decompose_and_solve_block_qubo(seed: int, dims: Tuple[int]) -> None: check.equal(block_indices, list(range(N))) # Assume that A and B are partitions of range(N) - def is_refinement_of(A: List[List[int]], B: List[List[int]]) -> bool: + def is_refinement_of(A: list[list[int]], B: list[list[int]]) -> bool: for a in A: if not any(set(a).issubset(b) for b in B): return False diff --git a/tests/solvers/test_solver.py b/tests/solvers/test_solver.py index 747b7104..05cdfac6 100644 --- a/tests/solvers/test_solver.py +++ b/tests/solvers/test_solver.py @@ -11,6 +11,7 @@ from qoolqit.devices import Device from qoolqit.execution import JobStatus +from qubosolver import transforms from qubosolver.config import ( EmbeddingConfig, DriveShapingConfig, @@ -18,7 +19,7 @@ LocalEmulator, RemoteEmulator, ) -from qubosolver.qubo_types import EmbedderType +from qubosolver.types.enums import EmbedderType from qubosolver.solver import ( QUBOInstance, QuboSolver, @@ -27,7 +28,6 @@ QUBOSolution, ) from qubosolver.qubo_analyzer import QUBOAnalyzer -from qubosolver.pipeline.basesolver import BaseSolver from mock.connection import MockConnection from pulser.backend.remote import ( @@ -66,7 +66,7 @@ def test_different_shots(simple_qubo_instance: QUBOInstance) -> None: ), ) solutions = default_solver.solve() - assert solutions.counts.sum() == 500 # type: ignore[union-attr] + assert solutions.counts.sum() == 500 lessshots_solver = QuboSolver( simple_qubo_instance, @@ -75,7 +75,7 @@ def test_different_shots(simple_qubo_instance: QUBOInstance) -> None: ), ) solutions = lessshots_solver.solve() - assert solutions.counts.sum() == 100 # type: ignore[union-attr] + assert solutions.counts.sum() == 100 @pytest.mark.priority(40) @@ -127,13 +127,13 @@ def test_parse_results() -> None: mock_result = Mock() mock_result.final_bitstrings = {"001": 10, "110": 5, "010": 3} - bitstrings, counts = BaseSolver.parse_results(mock_result) + solution = QUBOSolution.from_results(mock_result) expected_bitstrings = torch.tensor([[0, 0, 1], [1, 1, 0], [0, 1, 0]]) expected_counts = torch.tensor([10, 5, 3]) - torch.testing.assert_close(bitstrings, expected_bitstrings) - torch.testing.assert_close(counts, expected_counts) + torch.testing.assert_close(solution.bitstrings, expected_bitstrings) + torch.testing.assert_close(solution.counts, expected_counts) def test_parse_results_empty_final_bitstrings() -> None: @@ -141,12 +141,12 @@ def test_parse_results_empty_final_bitstrings() -> None: mock_result = Mock() mock_result.final_bitstrings = {} - bitstrings, counts = BaseSolver.parse_results(mock_result) + solution = QUBOSolution.from_results(mock_result) - check.equal(bitstrings.shape, (0, 0)) - check.equal(bitstrings.dtype, torch.int64) - check.equal(counts.shape, (0,)) - check.equal(counts.dtype, torch.int64) + check.equal(solution.bitstrings.shape, (0, 0)) + check.equal(solution.bitstrings.dtype, torch.int64) + check.equal(solution.counts.shape, (0,)) + check.equal(solution.counts.dtype, torch.int64) def test_parse_results_binary_string_conversion() -> None: @@ -154,13 +154,13 @@ def test_parse_results_binary_string_conversion() -> None: mock_result = Mock() mock_result.final_bitstrings = {"0101": 8, "1010": 12, "1111": 4} - bitstrings, counts = BaseSolver.parse_results(mock_result) + solution = QUBOSolution.from_results(mock_result) expected_bitstrings = torch.tensor([[0, 1, 0, 1], [1, 0, 1, 0], [1, 1, 1, 1]]) expected_counts = torch.tensor([8, 12, 4]) - torch.testing.assert_close(bitstrings, expected_bitstrings) - torch.testing.assert_close(counts, expected_counts) + torch.testing.assert_close(solution.bitstrings, expected_bitstrings) + torch.testing.assert_close(solution.counts, expected_counts) def test_parse_results_single_bitstring() -> None: @@ -168,13 +168,13 @@ def test_parse_results_single_bitstring() -> None: mock_result = Mock() mock_result.final_bitstrings = {"101": 25} - bitstrings, counts = BaseSolver.parse_results(mock_result) + solution = QUBOSolution.from_results(mock_result) expected_bitstrings = torch.tensor([[1, 0, 1]]) expected_counts = torch.tensor([25]) - torch.testing.assert_close(bitstrings, expected_bitstrings) - torch.testing.assert_close(counts, expected_counts) + torch.testing.assert_close(solution.bitstrings, expected_bitstrings) + torch.testing.assert_close(solution.counts, expected_counts) def test_parse_results_string_counts_to_integer_tensor() -> None: @@ -182,13 +182,13 @@ def test_parse_results_string_counts_to_integer_tensor() -> None: mock_result = Mock() mock_result.final_bitstrings = {"101": "15", "010": "8", "111": "12"} - bitstrings, counts = BaseSolver.parse_results(mock_result) + solution = QUBOSolution.from_results(mock_result) expected_bitstrings = torch.tensor([[1, 0, 1], [0, 1, 0], [1, 1, 1]]) expected_counts = torch.tensor([15, 8, 12]) - torch.testing.assert_close(bitstrings, expected_bitstrings) - torch.testing.assert_close(counts, expected_counts) + torch.testing.assert_close(solution.bitstrings, expected_bitstrings) + torch.testing.assert_close(solution.counts, expected_counts) def trivial_triangular_qubo(connection: Optional[RemoteConnection] = None) -> QuboSolverQuantum: @@ -199,7 +199,7 @@ def trivial_triangular_qubo(connection: Optional[RemoteConnection] = None) -> Qu [6.0, 6.0, -10.0], ] ) - qubo = QUBOInstance(Q) + qubo = QUBOInstance(torch.from_numpy(Q)) config = SolverConfig(use_quantum=True, do_preprocessing=False) @@ -232,15 +232,8 @@ def test_submit_integration(make_mock_connection: type[MockConnection], wait: bo job = solver.submit(drive, embedding) results = job.results() - bitstrings_local, counts_local = QuboSolverQuantum.parse_results(results) - - solution = QUBOSolution( - bitstrings=bitstrings_local.float(), - counts=counts_local, - costs=torch.Tensor(), - probabilities=None, - ) + solution = QUBOSolution.from_results(results) solution.costs = solution.compute_costs(solver.instance) solution.probabilities = solution.compute_probabilities() @@ -253,7 +246,6 @@ def test_submit_integration(make_mock_connection: type[MockConnection], wait: bo torch.testing.assert_close(bitstrings, torch.tensor([[0, 0, 1], [0, 1, 0], [1, 0, 0]])) - solution.bitstrings = solution.bitstrings.int() analyzer = QUBOAnalyzer([solution]) print(f"\n{analyzer.df}") @@ -277,16 +269,16 @@ def test_submit_integration(make_mock_connection: type[MockConnection], wait: bo assert isinstance(results_remote, Results) check.equal(remote_job.get_status(), JobStatus.DONE) - bitstrings_remote, counts_remote = QuboSolverQuantum.parse_results(results_remote) - torch.testing.assert_close(bitstrings_remote, bitstrings_local) - torch.testing.assert_close(counts_remote, counts_local) + solution_remote = QUBOSolution.from_results(results) + torch.testing.assert_close(solution_remote.bitstrings, solution.bitstrings) + torch.testing.assert_close(solution_remote.counts, solution.counts) -@pytest.mark.parametrize("preprocessing", [True, False]) -@pytest.mark.parametrize("postprocessing", [True, False]) +@pytest.mark.parametrize("postprocessing", [True, False], ids=["post", "no-post"]) +@pytest.mark.parametrize("preprocessing", [True, False], ids=["pre", "no-pre"]) def test_save_load_qubo_solver_quantum( - preprocessing: bool, postprocessing: bool, + preprocessing: bool, ) -> None: Q = torch.tensor( @@ -309,10 +301,12 @@ def test_save_load_qubo_solver_quantum( solver.preprocess() if preprocessing: + assert isinstance(solver.instance, transforms.variable_fixing.QUBOInstance) torch.testing.assert_close(solver.instance.coefficients, expected_preprocessed_Q) + torch.testing.assert_close(solver.instance._parent_instance.coefficients, Q) else: + check.is_not_instance(solver.instance, transforms.variable_fixing.QUBOInstance) torch.testing.assert_close(solver.instance.coefficients, Q) - torch.testing.assert_close(solver.fixtures.instance.coefficients, Q) # Save the solver file = io.BytesIO() @@ -323,12 +317,19 @@ def test_save_load_qubo_solver_quantum( loaded_solver = QuboSolverQuantum.load(file) # Verify the loaded solver has the same properties - # No need to have saved the preprocessed Q - torch.testing.assert_close(loaded_solver.instance.coefficients, Q) - torch.testing.assert_close(loaded_solver.fixtures.instance.coefficients, Q) + torch.testing.assert_close(loaded_solver.instance.coefficients, solver.instance.coefficients) + if preprocessing: + assert isinstance(solver.instance, transforms.variable_fixing.QUBOInstance) + assert isinstance(loaded_solver.instance, transforms.variable_fixing.QUBOInstance) + torch.testing.assert_close( + loaded_solver.instance._parent_instance.coefficients, + solver.instance._parent_instance.coefficients, + ) + check.equal(loaded_solver.instance._fixed_indices, solver.instance._fixed_indices) + else: + check.is_not_instance(loaded_solver.instance, transforms.variable_fixing.QUBOInstance) check.equal(loaded_solver.config.do_preprocessing, solver.config.do_preprocessing) check.equal(loaded_solver.config.do_postprocessing, solver.config.do_postprocessing) - check.equal(loaded_solver.fixtures.fixed_var_dict_list, solver.fixtures.fixed_var_dict_list) for method in [ "solve", diff --git a/tests/solvers/test_trival_qubo_solutions.py b/tests/solvers/test_trival_solutions.py similarity index 92% rename from tests/solvers/test_trival_qubo_solutions.py rename to tests/solvers/test_trival_solutions.py index 27488ee6..84e6999b 100644 --- a/tests/solvers/test_trival_qubo_solutions.py +++ b/tests/solvers/test_trival_solutions.py @@ -14,7 +14,7 @@ def test_classical_all_positive_trivial() -> None: should return a batch of one all-zero bitstring with solution_status 'trivial-zero'. """ - coeffs = [[1.0, 0.5], [0.5, 2.0]] + coeffs = torch.tensor([[1.0, 0.5], [0.5, 2.0]]) instance = QUBOInstance(coefficients=coeffs) config = SolverConfig(use_quantum=False) @@ -39,7 +39,7 @@ def test_quantum_all_negative_trivial(local_backend: LocalEmulator) -> None: config = SolverConfig(use_quantum=True, backend=local_backend) - coeffs = [[-1.0, -0.2], [-0.2, -3.0]] + coeffs = torch.tensor([[-1.0, -0.2], [-0.2, -3.0]]) instance = QUBOInstance(coefficients=coeffs) # check if value error is raised. with pytest.raises( @@ -47,7 +47,7 @@ def test_quantum_all_negative_trivial(local_backend: LocalEmulator) -> None: ): solver = QuboSolverQuantum(instance, config) - coeffs = [[-1.0, 0.0], [0.0, -3.0]] + coeffs = torch.tensor([[-1.0, 0.0], [0.0, -3.0]]) instance = QUBOInstance(coefficients=coeffs) solver = QuboSolverQuantum(instance, config) @@ -63,7 +63,7 @@ def test_quantum_all_negative_trivial(local_backend: LocalEmulator) -> None: def test_diagonal_trivial(local_backend: LocalEmulator) -> None: - coeffs = [[-1.0, 0.0], [0.0, 3.0]] + coeffs = torch.tensor([[-1.0, 0.0], [0.0, 3.0]]) instance = QUBOInstance(coefficients=coeffs) config = SolverConfig(use_quantum=True, backend=local_backend) diff --git a/tests/test_notebooks.py b/tests/test_notebooks.py index baa58048..1be61703 100644 --- a/tests/test_notebooks.py +++ b/tests/test_notebooks.py @@ -6,7 +6,6 @@ import subprocess import sys from pathlib import Path -from typing import List import pytest @@ -20,7 +19,7 @@ } -def get_ipynb_files(dir: Path) -> List[Path]: +def get_ipynb_files(dir: Path) -> list[Path]: files = [] for it in dir.iterdir(): diff --git a/tests/concepts/test_concepts.py b/tests/types/protocols/test_concepts.py similarity index 77% rename from tests/concepts/test_concepts.py rename to tests/types/protocols/test_concepts.py index eaf3d34d..1fc58986 100644 --- a/tests/concepts/test_concepts.py +++ b/tests/types/protocols/test_concepts.py @@ -1,6 +1,6 @@ from __future__ import annotations -from qubosolver import concepts +from qubosolver.types import protocols from typing import TYPE_CHECKING @@ -16,7 +16,7 @@ # to detect whether a class matches the concept when it's expected not to. -def _test_backend(backend: concepts.Backend) -> None: ... +def _test_backend(backend: protocols.Backend) -> None: ... def test_not_backend() -> None: @@ -24,28 +24,28 @@ def test_not_backend() -> None: class NotBackend0: ... _test_backend(NotBackend0()) # type: ignore[arg-type] - _b0: concepts.Backend = NotBackend0() # type: ignore[assignment] + _b0: protocols.Backend = NotBackend0() # type: ignore[assignment] class NotBackend1: def run(self: Self, program: QuantumProgram) -> None: return _test_backend(NotBackend1()) # type: ignore[arg-type] - _b1: concepts.Backend = NotBackend1() # type: ignore[assignment] + _b1: protocols.Backend = NotBackend1() # type: ignore[assignment] class NotBackend2: def run(self: Self, program: int) -> job.Job[Results]: return job._LocalJob(Results(atom_order=(), total_duration=0)) _test_backend(NotBackend2()) # type: ignore[arg-type] - _b2: concepts.Backend = NotBackend2() # type: ignore[assignment] + _b2: protocols.Backend = NotBackend2() # type: ignore[assignment] class NotBackend3: def run(self: Self, program: QuantumProgram) -> job.Job[Results] | None: return job._LocalJob(Results(atom_order=(), total_duration=0)) _test_backend(NotBackend3()) # type: ignore[arg-type] - _b3: concepts.Backend = NotBackend3() # type: ignore[assignment] + _b3: protocols.Backend = NotBackend3() # type: ignore[assignment] def test_backend() -> None: @@ -55,11 +55,11 @@ def run(self: Self, program: QuantumProgram) -> job.Job[Results]: return job._LocalJob(Results(atom_order=(), total_duration=0)) _test_backend(Backend0()) - _b0: concepts.Backend = Backend0() + _b0: protocols.Backend = Backend0() class Backend1: def run(self: Self, program: QuantumProgram | None) -> job.Job[Results]: return job._LocalJob(Results(atom_order=(), total_duration=0)) _test_backend(Backend1()) - _b1: concepts.Backend = Backend1() + _b1: protocols.Backend = Backend1() diff --git a/tests/data/test_qubo_analyzer.py b/tests/types/test_analyzer.py similarity index 93% rename from tests/data/test_qubo_analyzer.py rename to tests/types/test_analyzer.py index d3f8f13b..075c71fb 100644 --- a/tests/data/test_qubo_analyzer.py +++ b/tests/types/test_analyzer.py @@ -8,7 +8,7 @@ from qubosolver.qubo_instance import QUBOInstance from qubosolver.solver import QuboSolver from qubosolver.config import SolverConfig, ClassicalConfig -from qubosolver.qubo_types import ClassicalSolverType +from qubosolver.types.enums import ClassicalSolverType def test_init_single_solution(basic_solution: QUBOSolution) -> None: @@ -113,12 +113,8 @@ def test_analyzer_classical(simple_qubo_instance: QUBOInstance, classical_method analyzer = QUBOAnalyzer([solution], labels=["sol1"]) assert len(analyzer.df) == len(solution.bitstrings) - if classical_method == ClassicalSolverType.CPLEX: - assert "counts" not in analyzer.df.columns - assert "probs" not in analyzer.df.columns - else: - assert "counts" in analyzer.df.columns - assert "probs" in analyzer.df.columns + assert "counts" in analyzer.df.columns + assert "probs" in analyzer.df.columns def test_analyzer_quantum(simple_qubo_instance: QUBOInstance) -> None: diff --git a/tests/data/test_qubo_dataset.py b/tests/types/test_dataset.py similarity index 95% rename from tests/data/test_qubo_dataset.py rename to tests/types/test_dataset.py index 5ab00752..d7738870 100644 --- a/tests/data/test_qubo_dataset.py +++ b/tests/types/test_dataset.py @@ -44,7 +44,7 @@ def test_dataset_generation(negative_offdiag_rate: float | None) -> None: for qubo, _ in dataset: assert qubo.shape[0] == size - assert np.isclose(calculate_density(qubo, size), density, atol=1e-1) + assert np.isclose(calculate_density(qubo), density, atol=1e-1) assert torch.all(qubo >= coefficient_bounds[0]) assert torch.all(qubo <= coefficient_bounds[1]) if negative_offdiag_rate: diff --git a/tests/data/test_qubo_instance.py b/tests/types/test_instance.py similarity index 100% rename from tests/data/test_qubo_instance.py rename to tests/types/test_instance.py diff --git a/tests/types/test_linalg.py b/tests/types/test_linalg.py new file mode 100644 index 00000000..d5d3705d --- /dev/null +++ b/tests/types/test_linalg.py @@ -0,0 +1,79 @@ +from __future__ import annotations + +import os +import pytest +import torch +from unittest.mock import patch + +from qubosolver.types._linalg import ( + _device_from_env, + _float_type_from_env, + GlobalConfig, + device, +) + + +def test_device_from_env_defaults_to_cpu() -> None: + with patch.dict(os.environ, {}, clear=True): + result = _device_from_env() + assert result == torch.device("cpu") + + +def test_device_from_env_use_gpu_1() -> None: + with patch.dict(os.environ, {"USE_GPU": "1"}, clear=False): + result = _device_from_env() + assert result == torch.device("cuda") + + +def test_device_from_env_use_gpu_false() -> None: + with patch.dict(os.environ, {"USE_GPU": "false"}, clear=False): + result = _device_from_env() + assert result == torch.device("cpu") + + +def test_device_from_env_invalid_qubo_solver_device_raises() -> None: + with patch.dict(os.environ, {"QUBO_SOLVER_DEVICE": "invalid_device_string"}, clear=False): + with pytest.raises(ValueError, match="Invalid QUBO_SOLVER_DEVICE"): + _device_from_env() + + +def test_float_type_from_env_float64_dtype() -> None: + with patch.dict(os.environ, {"QUBO_SOLVER_FLOAT_DTYPE": "float64"}, clear=False): + result = _float_type_from_env() + assert result == torch.float64 + + +def test_float_type_from_env_defaults_to_float32() -> None: + with patch.dict(os.environ, {}, clear=True): + result = _float_type_from_env() + assert result == torch.float32 + + +def test_float_type_from_env_invalid_dtype_raises() -> None: + with patch.dict(os.environ, {"QUBO_SOLVER_FLOAT_DTYPE": "float128"}, clear=False): + with pytest.raises(ValueError, match="Invalid QUBO_SOLVER_FLOAT_DTYPE"): + _float_type_from_env() + + +def test_global_config_use_double_precision_sets_float64() -> None: + original = GlobalConfig._float_dtype + try: + GlobalConfig.use_double_precision(enable=True) + assert GlobalConfig._float_dtype == torch.float64 + finally: + GlobalConfig._float_dtype = original + + +def test_global_config_set_float_precision_invalid_dtype_raises() -> None: + with pytest.raises(ValueError): + GlobalConfig.set_float_precision(torch.int32) + + +def test_global_config_use_gpu_sets_cuda_device() -> None: + original = GlobalConfig._device + try: + GlobalConfig.use_gpu(enable=True) + assert GlobalConfig._device == torch.device("cuda") + assert device() == torch.device("cuda") + finally: + GlobalConfig._device = original