# This code is part of qtealeaves.
#
# This code is licensed under the Apache License, Version 2.0. You may
# obtain a copy of this license in the LICENSE.txt file in the root directory
# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
#
# Any modifications or derivative works of this code must retain this
# copyright notice, and modified files need to carry a notice indicating
# that they have been altered from the originals.
"""
This submodule contains both exact and tensor-network based solvers
for a generic QUBO problem in a standard format.
Exact solvers are limited to small problem sizes due to exponential
computational complexity, whereas tensor network–based methods can
handle larger and more complex QUBO matrices, making them suitable
for more demanding optimization tasks.
"""
# pylint: disable=invalid-name
# pylint: disable=too-many-lines
# pylint: disable=too-many-locals
# pylint: disable=too-many-arguments
# pylint: disable=too-many-statements
# pylint: disable=too-many-instance-attributes
# Modules
# ----------------------------------------------------------------------------
# for sphinx resolving | in typing
from __future__ import annotations
import codecs
import heapq
import json
import logging
import re
from collections import defaultdict
from dataclasses import dataclass, field
from pathlib import Path
from time import time
from typing import Self
import numpy as np
import numpy.typing as npt
from qtealeaves import __version__
from qtealeaves.convergence_parameters import TNConvergenceParameters
from qtealeaves.modeling import QuantumModel, RandomizedLocalTerm, TwoBodyAllToAllTerm1D
from qtealeaves.observables import Local, State2File, TNObservables
from qtealeaves.operators import TNSpin12Operators
from qtealeaves.simulation import QuantumGreenTeaSimulation
from qtealeaves.tensors import TensorBackend
from qtealeaves.tooling import QTeaLeavesError
# isort complains, but a clean solution would do two things
# 1) move torch up to third-party (after numpy) or down to optional
# 2) qredtea cannot be imported here (cyclic dependency as qtealeaves
# is imported in qredtea.
# isort: off
try:
import torch as to
from qredtea.torchapi import default_pytorch_backend
except ImportError as exc:
to = exc
default_pytorch_backend = None
try:
import cupy as cp
except ImportError as exc:
cp = exc
# isort: on
from .opt_tooling import (
construct_exact_observables,
construct_exact_spinglass_hamiltonian,
construct_spinglass_driver_state,
generate_perturbation_transverse_fields,
get_exact_spectrum,
measure_exact_observables,
read_tn_state,
spins_to_binaries,
)
__all__ = ["QUBOSolver", "QUBOConvergenceParameters"]
# ----------------------------------------------------------------------------
# ----------------------------------------------------------------------------
logger = logging.getLogger(__name__)
class QUBOProblemError(QTeaLeavesError):
"""
Custom derived exception class for errors in the QUBO solver
instantiation process.
"""
class QUBOSetUpError(QTeaLeavesError):
"""
Custom derived exception class for errors during the set-up of
the QUBO solver simulation.
"""
class QUBOSolverError(QTeaLeavesError):
"""
Custom derived exception class for errors when or after
solving an input QUBO problem.
"""
[docs]
@dataclass
class QUBOConvergenceParameters(TNConvergenceParameters):
# pylint: disable=anomalous-backslash-in-string
r"""
Convergence parameters for the tensor-network based QUBO solver.
Attributes
----------
Same as the base class
:class:`qtealeaves.convergence_parameters.TNConvergenceParameters`.
Attributes (ground-state search)
--------------------------------
max_bond_dimension : int, optional
Maximum bond dimension of the tensor network during the
ground-state search.
The maximum value of the bond dimension is upper bounded by
:math:`2^{n \\mathbin{//} 2}`, where :math:`n` is the number
of spins in the spin glass system.
Default to 8.
min_bond_dimension : int, optional
Minimum bond dimension of the tensor network during the
ground-state search.
(Suggested for simulation which use single tensor updates)
Default to :obj:`None`, which means, it is set to the same
value of the maximum bond dimension.
This default behavior is specific to the solving of classical
optimization problems.
cut_ratio : float, optional
Control which singular values are truncated during splitting tensors.
Either relative values or truncated norm can be considered.
This fine tuning can be set via ``trunc_method``.
Default to 1e-9
trunc_method : str, optional
Method use to truncate the singular values.
Available:
- ``R``: use cut_ratio
Cut ratio :math:`\\epsilon` after which the singular values are
neglected, i.e. if :math:`\\lambda_1` is the bigger singular values
then after an SVD we neglect all the singular values such that
:math:`\\lambda_i/\\lambda_1\\leq\\epsilon`.
- ``N``: use maximum norm
Maximum value of the norm neglected for the singular values during
the trunctation.
Default to 'R'.
max_iter : int, optional
The maximum number of sweeps for the ground-state search.
After ``max_iter`` the ground-state search ends.
Default to 20 sweeps.
abs_deviation : float, optional
Exit criterion for the ground-state search.
If the energy of the current sweep has an absolute per-sweep
deviation from the previous sweeps below this threshold, the
tensor network optimization ends.
Default to 4e-12.
rel_deviation : float, optional
Exit criterion for ground-state search.
If the energy of the current sweep has a relative per-sweep
deviation from the previous sweeps below this threshold, the
tensor network optimization ends.
Default to 1e-12.
n_points_conv_check : int, optional
The number of sweeps used when checking convergence for the
ground state search.
Exit criteria are not checked until this many sweeps have been
performed.
Default to 6.
random_sweeps : bool, optional
If :obj:`True`, perform random sweeps to optimize the tensors.
It could help in escaping local minima for classical combinatorial
optimization problems.
Default to :obj:`False`.
enable_spacelink_expansion : bool, optional
Mode used for sweeping along the tensors during the optimization.
If :obj:`False`, standard sweeps without space link expansion are
performed to obtain the ground state; if :obj:`True`, the first
half of the sweeps use space link expansion, while in the second
half the space link expansion is disabled.
Note that this second option could automatically change some other
convergence parameters, for example the
:attr:`~TensorNetworkConvergenceParameters.n_points_conv_check` for
consistency. In that case a warning to the user will be thrown.
Default to :obj:`False`.
arnoldi_maxiter : int, optional
Maximum number of Lanczos vectors to be used in the eigensolver.
Default to 32.
psi0 : str, optional
The initial tensor network state from which to start the ground-state
search for the spin glass Hamiltonian representing the QUBO problem.
The available options are:
- ``random``: the initial state of the tensors is initialized
randomly;
- ``driver_product_state``: the tensors are initialized in the
product state
.. math::
\\left\\vert{+}\\right\\rangle^{n}
with
.. math::
\\left\\vert{+}\\right\\rangle
=
\\dfrac{1}{\\sqrt{2}}
\\left(\\left\\vert{0}\\right\\rangle
+
\\left\\vert{1}\\right\\rangle\\right)
i.e., in the ground state of the driver Hamiltonian
.. math::
\\mathcal{H}_{\\scriptscriptstyle driver}
=
\\sum\\limits_{j=1}^n \\sigma_j^x .
Within this option, the product state is exactly constructed
with the given *initial bond dimension*
:math:`\\chi =` ``ini_bond_dimension``;
Default to 'random'.
ini_bond_dimension : int, optional
The bond dimension used to construct the initial tensor-network
state. Using subspace expansion the bond dimension can grow.
Default to :obj:`None`, i.e., the initial bond dimension is
initialized to ``max_bond_dimension``.
psi0_noise : float, optional
Numerical noise to construct the exact initial product state with a
given bond dimension.
More details can be found in the :mod:`opt_tooling` submodule.
Default to 1e-7.
transverse_field_ratio : float, optional
The ratio of the magnitude of the strength of the classical
longitudinal terms (Z-Pauli) to the quantum off-diagonal
transverse field terms (X-Pauli) to be added.
This ratio controls the magnitude of the transverse field
perturbation relative to the classical Hamiltonian, whose ground
state must be determined.
See the function :func:`generate_perturbation_transverse_fields`
in the :mod:`opt_tooling` submodule.
Default to :obj:`None`, i.e., no transverse fields will be
added to the classical Hamiltonian.
randomseed : list[int], optional
Random seeds for reprodubibility and solution diversity.
As for :mod:`qtealeaves` compatibility, it has to be a
list of 4 integers.
Default to [1907, 79, 17, 299].
Attributes (ansatz)
-------------------
data_type : str | list[str], optional
Data type of the tensors in the chosen ansatz.
Available options are
- "Z": double precision complex (.complex128)
- "C": single precision complex (.complex64)
- "D": double precision real (.float64)
- "S": single precision real (.float32)
- ["DT", "DT", ..., "DT"]: allows selecting a
specific data type for each sweep, where "DT"
is a character in {"Z", "C", "D", "S"}.
Default to ``Z``.
tn_ansatz : str, optional
Ansatz to be used for approximating the wave function of
the mapped MP4EO Hamiltonian.
Available topologies are ``TTN`` for Tree Tensor Networks and
``MPS`` for Matrix Product States.
Default to ``TTN``.
mpo_mode : int, optional
MPO representation inside the tensor network simulation.
Other options are:
- 0 --> TPO via iTPO
- 1 --> iTPO
- 3 --> iTPO with compression
- 10 --> sparse MPO
- 11 --> indexed sparse MPO
Default to -1, which enables an auto-selection..
use_spinglass_sparse_mpo : bool, optional
Use a sparse MPO for spin-glass systems.
It can be used only with mpo_mode = 10, 11.
Default to :obj:`False`.
linalg_backend_library : str, optional
Library to handle linear algebra operations at the level of
the tensors.
Available libraries are:
- ``numpy``, both on CPU and GPU;
- ``pytorch`` (pass **torch**), both on CPU and GPU;
Note that :mod:`qredtea` is required for non-trivial backend.
Default to 'numpy'.
device : str, optional
Device where to run the optimization. Available options are
``cpu`` to run the ground-state search on the CPUs, and ``gpu``
to run the ground-state search on the GPUs. Also the mix device
``cgpu`` is available.
Default to ``cpu``.
Attributes (logging)
--------------------
base_log_dir : str | Path, optional
Base directory for saving log files associated with
the :class:`QUBOSolver`.
Default to ``QUBO_RESULTS``.
log_dir : str | Path, optional
Sub-directory in :attr:`~QUBOConvergenceParameters.base_log_dir`
for saving log files associated with :mod:`qtealeaves` I/O processes.
Default to ``qtealeaves``.
tn_io_tag : str, optional
Extra common tag to be added at the end of simulation log
folders in order to distinguish different simulations.
Default to :obj:`None`.
_tn_io_tag_history : List[str]
A list to store tensor-network log-file I/O tag history.
tn_input_folder : str, optional
Name of the folder where to save the input data and parameters
for the tensor network ground-state search.
If the folder does not exist, it will be automatically created.
Default to ``tn_inputs``.
tn_input_filepath : str | Path
Path to the folder for tensor-network input files and logs.
tn_output_folder : str, optional
Name of the folder where to save the output data and log files
for the tensor network ground-state search.
If the folder does not exist, it will be automatically created.
Default to ``tn_outputs``.
tn_output_filepath : str | Path
Path to the folder for tensor-network output files and logs.
tn_states_folder : str, optional
Name of the folder where to save the tensor network states
produced.
If the folder does not exist, it will be automatically created.
Default to ``tn_states``.
tn_state_file : str | Path
Full path to the file where the optimized tensor
network state will be stored after optimization.
tn_states_filepath : str | Path
Path to the folder for tensor-network state files
and logs.
tn_state_filename : str
The name of the file where the optimized tensor
network state will be stored after optimization.
tn_initial_state_filename : str
The name of the file where the initial tensor
network state will be stored before optimization.
tn_initial_state_file : str | Path
Full path to the file where the initial tensor
network state will be stored before optimization.
Notes
-----
An instance of this class must be created and passed to the QUBO
solver whenever the QUBO problem is to be solved using tensor
networks.
"""
# Instance attribute initialization
# ------------------------------------------------------------------------
# TN convergence attributes: ground-state search
max_bond_dimension: int = 8
min_bond_dimension: int | None = None
cut_ratio: float | None = 1e-9
trunc_method: str = "R"
max_iter: int = 20
abs_deviation: float = 4e-12
rel_deviation: float = 1e-12
n_points_conv_check: int = 6
arnoldi_maxiter: int = 32
ini_bond_dimension: int | None = None
psi0: str = "random"
psi0_noise: float = 1e-7
random_sweeps: bool = False
enable_spacelink_expansion: bool = False
transverse_field_ratio: float | None = None
randomseed: list[int] = field(default_factory=lambda: [1907, 79, 17, 299])
# TN convergence attributes: ansatz
data_type: str | list[str] = "Z"
tn_ansatz: str = "TTN"
mpo_mode: int = -1
use_spinglass_sparse_mpo: bool = False
linalg_backend_library: str = "numpy"
device: str = "cpu"
# TN convergence attributes: logging and I/O
tn_io_tag: str | None = None
base_log_dir: str | Path = Path("QUBO_RESULTS")
log_dir: str | Path = Path("qtealeaves")
tn_input_folder: str = "inputs"
tn_output_folder: str = "outputs"
tn_states_folder: str = "states"
# Internals (not user-supplied)
tn_type: int = field(init=False)
_tn_io_tag_history: list[str] = field(init=False, default_factory=list)
tn_input_filepath: Path = field(init=False)
tn_output_filepath: Path = field(init=False)
tn_states_filepath: Path = field(init=False)
tn_state_filename: str = field(init=False, default="psi")
tn_state_file: str = field(init=False)
tn_initial_state_file: str | None = field(init=False)
tn_initial_state_filename: str | None = field(init=False)
# ------------------------------------------------------------------------
# ------------------------------------------------------------------------
# Dunder methods
# ------------------------------------------------------------------------
# pylint: disable=too-many-branches
def __post_init__(self) -> None:
"""
Initializes the convergence parameters for tensor-network
simulations.
"""
# Validate user parameters
if not (
isinstance(self.randomseed, list)
and all(isinstance(seed, int) for seed in self.randomseed)
):
raise QUBOSetUpError(
"Random seeds must be a list of (at least) 4 positive integers."
)
# Set default minimum bond dimension
if self.min_bond_dimension is None:
self.min_bond_dimension = self.max_bond_dimension
if self.min_bond_dimension > self.max_bond_dimension:
raise QUBOSetUpError(
f"Minimum bond dimension ({self.min_bond_dimension}) cannot "
f"exceed the maximum one {self.max_bond_dimension}."
)
self.ini_bond_dimension = self.ini_bond_dimension or self.max_bond_dimension
# Choose method for sweeping
if self.enable_spacelink_expansion:
half_nsweeps = self.max_iter // 2
sweep_mode = [0] * (half_nsweeps)
sweep_mode += [1] * (self.max_iter - half_nsweeps)
## Re-define convergence criteria for consistency
self.n_points_conv_check += half_nsweeps
logger.warning(
"Enabling spacelink expansion for the tensor-network "
"ground-state search...\nRedefining the number of sweeps "
"after which to check convergence criteria accordingly to %d...",
self.n_points_conv_check,
)
else:
sweep_mode = 1
# Check tensors data type
valid_data_type = {"Z", "C", "S", "D"}
if isinstance(self.data_type, list):
if len(self.data_type) != self.max_iter:
raise QUBOSetUpError(
"Mismatch between list of data types and"
" number of sweeps. Reset one of the two."
)
if any(dt not in valid_data_type for dt in self.data_type):
raise QUBOSetUpError("Invalid data type for tensors.")
elif self.data_type not in valid_data_type:
raise QUBOSetUpError(f"Invalid data type '{self.data_type}' for tensors.")
# Inherited attributes from qtealeaves
super().__init__(
max_iter=self.max_iter,
abs_deviation=self.abs_deviation,
rel_deviation=self.rel_deviation,
n_points_conv_check=self.n_points_conv_check,
max_bond_dimension=self.max_bond_dimension,
cut_ratio=self.cut_ratio,
min_bond_dimension=self.min_bond_dimension,
ini_bond_dimension=self.ini_bond_dimension,
svd_ctrl="A",
random_sweep=self.random_sweeps,
arnoldi_maxiter=self.arnoldi_maxiter,
statics_method=sweep_mode,
data_type=self.data_type,
device=self.device,
)
# Validate ansatz type
if self.tn_ansatz not in {"TTN", "MPS"}:
raise QUBOSetUpError(f"Invalid '{self.tn_ansatz}' ansatz-type.")
self.tn_type = 5 if self.tn_ansatz == "TTN" else 6
# Validate the linear algebra library for the tensors
if self.linalg_backend_library not in {"numpy", "torch"}:
raise QUBOSetUpError(
f"Invalid linear algebra library {self.linalg_backend_library}"
" for tensors."
)
self.tensor_backend = self._get_tensor_backend(
linalg_library=self.linalg_backend_library,
tensor_dtype=self.data_type,
device=self.device,
)
# Validate initial tensor-network state
if self.psi0 not in {"random", "driver_product_state"}:
raise QUBOSetUpError(
"Tensor network state initialization option"
f" '{self.psi0}' not available."
)
# Choose the MPO representation
mpo_mode = self._get_mpo_mode(self.mpo_mode, self.tn_ansatz)
if self.use_spinglass_sparse_mpo and mpo_mode not in {10, 11}:
raise QUBOSetUpError(
"Cannot use SPO for spin glasses when MPO mode is "
f"{mpo_mode}. Use instead MPO mode '10' or '11'."
)
self.mpo_mode = mpo_mode
# Prepare log files for qtealeaves I/O
self.base_log_dir = Path(self.base_log_dir)
self.log_dir = Path(self.log_dir)
self.base_dir = self.base_log_dir / self.log_dir
self.base_dir.mkdir(parents=True, exist_ok=True)
self.tn_io_tag = f"_{self.tn_io_tag}" if self.tn_io_tag else ""
self._tn_io_tag_history.append(self.tn_io_tag)
self.tn_initial_state_filename = (
None
if self.psi0 == "random"
else f"psi0_{self.psi0}.pkl{self.tn_ansatz.lower()}"
)
self._update_tn_folder_paths(create=True)
# ------------------------------------------------------------------------
# ------------------------------------------------------------------------
# Properties
# ------------------------------------------------------------------------
@property
def bond_dim(self) -> int:
"""Read-only access to the bond dimension."""
return self.max_bond_dimension
@bond_dim.setter
def bond_dim(self, new_bond_dimension: int) -> None:
"""Sets a new value for the bond dimension."""
# Validate new bond dimension
if new_bond_dimension < 1:
raise ValueError("Maximum bond dimension must be positive.")
if (
self.min_bond_dimension is not None
and new_bond_dimension < self.min_bond_dimension
):
raise ValueError(
f"Maximum bond dimension ({new_bond_dimension}) cannot be "
f"less than minimum bond dimension ({self.min_bond_dimension})."
)
self.max_bond_dimension = new_bond_dimension
self.sim_params["max_bond_dimension"] = new_bond_dimension
@property
def min_bond_dim(self) -> int:
"""Read-only access to the minimum bond dimension."""
return self.min_bond_dimension
@min_bond_dim.setter
def min_bond_dim(self, new_bond_dimension: int) -> None:
"""Sets a new value for the minimum bond dimension."""
# Validate new bond dimension
if new_bond_dimension < 1:
raise ValueError("Minimum bond dimension must be positive.")
if new_bond_dimension > self.max_bond_dimension:
raise ValueError(
f"Minimum bond dimension ({new_bond_dimension}) cannot be "
"greater than maximum bond dimension ({self.max_bond_dimension})."
)
self.min_bond_dimension = new_bond_dimension
self.sim_params["min_bond_dimension"] = new_bond_dimension
# ------------------------------------------------------------------------
# ------------------------------------------------------------------------
# Static methods
# ------------------------------------------------------------------------
@staticmethod
def _get_mpo_mode(input_mode: int, tn_ansatz: str) -> int:
"""
Selects the correct MPO mode in :mod:`qtealeaves` based on the
selected level from the user, e.g., coming from CLI options.
Arguments
---------
input_mode : int
The selected representation for the MPOs.
tn_ansatz : str
The chosen tensor-network ansatz (MPS or TTN).
Returns
-------
mpo_mode : int
The correct MPO mode to be used in :mod:`qtealeaves`.
"""
# Validate user's choice
if input_mode not in {-1, 0, 3, 10, 11}:
raise QUBOSetUpError(f"Invalid MPO mode '{input_mode}' selected.")
# Enable auto-selection from the user
mpo_mode = {("TTN", -1): 3, ("MPS", -1): 0}.get(
(tn_ansatz, input_mode), input_mode
)
# Output
return mpo_mode
@staticmethod
def _get_tensor_backend(
linalg_library: str, tensor_dtype: str, device: str
) -> TensorBackend:
"""
Selects the correct tensor backend from :mod:`qtealeaves` and
:mod:`qredtea` based on user simulation parameter choices.
Arguments
---------
linalg_library : str
The chosen linear algebra library to use as a backend
for tensor operations.
tensor_dtype : str
Character that identifies the numerical data-type of
the tensors.
device : str
The device to use for executing the tensor operations
and the computations.
Returns
-------
qtea_backend : qtealeaves.tensors.TensorBackend
The selected tensor-network backend.
"""
# Validate and select the correct linear algebra library
if linalg_library == "numpy":
if device == "gpu":
if isinstance(cp, Exception):
raise ImportError(
"CuPy is required for 'numpy' with GPU device."
) from cp
xp = cp
else:
xp = np
qtea_backend = TensorBackend
elif linalg_library == "torch":
if isinstance(to, Exception):
raise ImportError("PyTorch is required for 'torch' backend.") from to
xp = to
qtea_backend = default_pytorch_backend
else:
raise QUBOSetUpError(
f"Unknown linear algebra library {linalg_library}"
" for the Quantum TEA backend."
)
# Validate and select the correct data type for tensors
dtype_opts = {
"Z": xp.complex128,
"C": xp.complex64,
"D": xp.float64,
"S": xp.float32,
}
backend_dtype = dtype_opts.get(tensor_dtype)
# Output
return qtea_backend(dtype=backend_dtype, device=device)
# ------------------------------------------------------------------------
# ------------------------------------------------------------------------
# Instance methods
# ------------------------------------------------------------------------
[docs]
def print_details(self) -> None:
"""
Prints details about convergence parameters for tensor-network
based simulations.
Arguments
---------
:obj:`None`
Returns
-------
:obj:`None`
"""
# Header
head = "Convergence parameters for the tensor-network QUBO solver"
sep_line = "-" * len(head)
print(f"{sep_line}\n{head}\n{sep_line}\n")
# Details
self._print_topology_details()
self._print_initial_state_details()
self._print_ground_state_search_details()
def _print_topology_details(self) -> None:
"""
Prints details about the tensor-network topology.
Arguments
---------
:obj:`None`
Returns
-------
:obj:`None`
"""
print("Chosen tensor-network topology details")
print(f"\tAnsatz <-- {self.tn_ansatz}")
data_type_details = {
"Z": "double-precision complex (xp.complex128)",
"C": "single-precision complex (xp.complex64)",
"D": "double-precision real (xp.float64)",
"S": "single-precision real (xp.float32)",
}
data_type_repr = (
self.data_type
if isinstance(self.data_type, list)
else data_type_details.get(self.data_type, "Unknown data-type")
)
print(f"\tData type of the tensors <-- using {data_type_repr}")
mpo_mode_details = {
0: "SPO via iTPO",
1: "iTPO",
3: "iTPO with compression",
10: "sparse MPO",
11: "indexed sparse MPO",
}
print(
"\tMPO representation <-- using "
f"{mpo_mode_details.get(self.mpo_mode, 'Unknown representation')}"
)
if self.use_spinglass_sparse_mpo:
print(
"\t"
+ " " * len("MPO representation ")
+ "<-- with efficient spin-glass SPO"
)
print(
"\tChosen random seeds for diversity and reproducibility"
f" <-- {self.randomseed}"
)
def _print_ground_state_search_details(self) -> None:
"""
Prints details about the ground-state search.
Arguments
---------
:obj:`None`
Returns
-------
:obj:`None`
"""
print("Parameters of the ground-state search")
print(f"\tMaximum bond dimension <-- {self.max_bond_dimension}")
print(f"\tMinimum bond dimension <-- {self.min_bond_dimension}")
print(f"\tCut ratio <-- {self.cut_ratio}")
trunc_method_details = {
"R": "by relative cut ratio w.r.t. the largest singular value",
"N": "by maximum norm",
}
print(
"\tRelated truncation method <-- "
f"{trunc_method_details.get(self.trunc_method, 'Unknown strategy')}"
)
print(
"\tMaximum number of iteration for the local solver <-- "
f"{self.arnoldi_maxiter}"
)
print(f"\tNumber of sweeps <-- {self.max_iter}")
print(f"\tAbsolute deviation <-- {self.abs_deviation}")
print(f"\tRelative deviation <-- {self.rel_deviation}")
print(f"\tNumber of convergence points <-- {self.n_points_conv_check}")
print(f"\tUse random sweep strategy <-- {self.random_sweeps}")
print(
f"\tUse spacelink expansion strategy <-- {self.enable_spacelink_expansion}"
)
if self.transverse_field_ratio:
print(
"\tApply random transverse X-Pauli like fields "
f"with perturbation ratio <-- {self.transverse_field_ratio:.2e}"
)
print(
"\tLinear algebra library for the backend "
f"<-- {self.linalg_backend_library}"
)
print(f"\tSelected device for the computations <-- {self.device}\n")
def _print_initial_state_details(self) -> None:
"""
Prints details about the initial tensor-network state
before starting the ground-state search.
"""
print("Initial tensor-network state details")
print(f"\tInitial state mode <-- {self.psi0}")
print(f"\tInitial bond dimension <-- {self.ini_bond_dimension}")
if self.psi0 != "random":
print(f"\tNoise for padding <-- {self.psi0_noise}\n")
[docs]
def update_tn_io_tag(self, new_info: str) -> None:
"""
Updates :attr:`~QUBOConvergenceParameters.tn_io_tag` with
``new_info``, keeps track of the history and updates I/O
folder paths dynamically.
Arguments
---------
new_info : str
The new information to append at the end of
:attr:`~QUBOConvergenceParameters.tn_io_tag`.
Returns
-------
:obj:`None`
"""
# Ensure string data type
if not isinstance(new_info, str):
new_info = str(new_info)
# Add underscore if not already present
if not new_info.startswith("_"):
new_info = f"_{new_info}"
# Save the current tag in history
self._tn_io_tag_history.append(self.tn_io_tag)
# Update the tag
self.tn_io_tag += new_info
# Update folder paths based on the new tag
self._update_tn_folder_paths(True)
[docs]
def revert_to_previous_tn_io_tag(self) -> None:
"""
Reverts :attr:`~QUBOConvergenceParameters.tn_io_tag` to its
previous value stored in the history stack, and updates folder
paths accordingly.
Arguments
---------
:obj:`None`
Returns
-------
:obj:`None`
Raises
------
:exc:`ValueError`
If there is no previous tag to revert to.
"""
if len(self._tn_io_tag_history) <= 1:
raise ValueError("No previous tensor-network log-file tag to revert to.")
# Pop the last tag from the history and revert
self.tn_io_tag = self._tn_io_tag_history.pop()
# Update folder paths based on the new tag
self._update_tn_folder_paths()
[docs]
def reset_tn_io_tag(self) -> None:
"""
Resets :attr:`~QUBOConvergenceParameters.tn_io_tag` to its
original value (first value in history), and updates folder
paths accordingly.
Arguments
---------
:obj:`None`
Returns
-------
:obj:`None`
"""
if self._tn_io_tag_history:
## Reset to the first tag in history
self.tn_io_tag = self._tn_io_tag_history[0]
## Clear history except the original tag
self._tn_io_tag_history = [self.tn_io_tag]
# Update folder paths based on the reset tag
self._update_tn_folder_paths()
def _update_tn_folder_paths(self, create: bool = False) -> None:
"""
Updates tensor-network I/O folder paths based on the current
:attr:`~QUBOConvergenceParameters.tn_io_tag`.
Also ensures directories exist.
Arguments
---------
create : bool, optional
Whether to create again the folder after updating
or not.
Default to :obj:`False`.
Returns
-------
:obj:`None`
"""
self.tn_input_filepath = self.base_dir / Path(
f"{self.tn_input_folder}{self.tn_io_tag}/"
)
self.tn_output_filepath = self.base_dir / Path(
f"{self.tn_output_folder}{self.tn_io_tag}/"
)
self.tn_states_filepath = self.base_dir / Path(
f"{self.tn_states_folder}{self.tn_io_tag}/"
)
# Ensure the directories exist
if create:
self.tn_input_filepath.mkdir(parents=True, exist_ok=True)
self.tn_output_filepath.mkdir(parents=True, exist_ok=True)
self.tn_states_filepath.mkdir(parents=True, exist_ok=True)
self.tn_state_file = str(self.tn_states_filepath / self.tn_state_filename)
self.tn_initial_state_file = str(
None
if not self.tn_initial_state_filename
else self.tn_input_filepath / self.tn_initial_state_filename
)
# ------------------------------------------------------------------------
# ------------------------------------------------------------------------
[docs]
class QUBOSolver:
r"""
Solver for a generic QUBO problem
.. math::
\min_{x \in \{0, 1\}^n} f(\boldsymbol{x})
.. math::
f(\boldsymbol{x}) = \boldsymbol{x}^T \mathcal{Q} \boldsymbol{x}
where :math:`\boldsymbol{x}` is the :math:`n`-dimensional binary
vector describing the binary configuration of the QUBO decision
variables. The QUBO problem is uniquely described by the QUBO matrix
:math:`\mathcal{Q}`, which must be a symmetric or upper triangular
real matrix.
This QUBO solver implements different methods to obtain a solution:
- ``brute-force`` obtains the exact solution (global optimum)
but it's limited by the number of binaries (:math:`n \leq 30`);
- ``exact ground-state search`` maps the QUBO problem into an
equivalent spin glass Hamiltonian and encodes the QUBO solution
into the ground state of this Hamiltonian. The ground state is
exactly reconstructed by allocating the full Hamiltonian matrix
and sorting its diagonal. This method obtains the exact solution
(global optimum), but it's limited by the number of binaries
(:math:`n \leq 30`, depending on the available RAM);
- ``tensor-network ground-state search`` also maps the QUBO problem
into a corresponding ground-state search of a spin glass Hamiltonian,
but the solution is found by applying tensor-network optimization
using tools from qtealeaves.
This solver is not limited by the number of binaries but by the
amount of entanglement needed for the computation.
Parameters
----------
qubo_matrix : numpy.ndarray[float]
The input QUBO problem represented by its QUBO matrix.
The input matrix must be either a symmetric or an
upper triangular 2D :mod:`numpy` array of real numbers.
Only the upper triangular part of the input matrix will
be considered, i.e., if :math:`Q_{ij}` are the matrix
elements of the input QUBO problem, only those for
:math:`i \leq j` will be used by the solver.
qubo_offset : float, optional
The QUBO problem energy cost offset.
Ignored during the optimization.
Default to 0.0.
qubo_variable_labeling : dict[str, int], optional
Labeling scheme from the original problem variable names
to the index assign in the QUBO problem reformulation.
Keys are variable names, values are the corresponding indices
within the one-dimensional embedding.
Default to :obj:`None`, i.e., a trivial labeling scheme
where the :math:`j`-th decision variable takes QUBO index
:math:`j`.
Attributes
----------
qubo_matrix : numpy.ndarray[float]
The QUBO matrix that defines the input QUBO problem.
n_binaries int
Number of decision variables in the QUBO problem.
n_interactions : int
Number of two-body couplings between pairs of decision
variables.
spinglass_couplings : dict[str, float | numpy.ndarray[float]]
Dictionary that contains the Hamiltonian couplings of
mapped 1-dimensional spin glass model corresponding
to the input QUBO problem.
cost : list[float]
QUBO cost function values related to the obtained solution
space after optimization.
solution : list[numpy.ndarray[int]]
Solution space obtained through the QUBO solver.
_variables_to_qubo_index_map : dict[str, int]
Map from the original problem variable names to the index
assign in the QUBO problem reformulation.
Keys are variable names, values are the corresponding
indices within the one-dimensional embedding.
If empty, the trivial mapping has been considered, i.e.,
variable :math:`j` has index :math:`j` in the embedded
1D chain.
status : str
Status of the solver.
runtimes : dict[str, float]
Dictionary that contains runtime information for the profiling
of the solver algorithm.
Notes
-----
**A very important note for users**: if the QUBO matrix describing
the problem passed to this class is **symmetric**, the solver
automatically considers only its upper triangular part, automatically
**doubling** the remaining off-diagonal elements, while leaving the
diagonal elements unchanged.
If instead the matrix is already **upper triangular**, the solver assumes
that the off-diagonal elements **have already been multiplied by 2**.
The definition of the QUBO problem cost function is thus as follows:
.. math::
f(\boldsymbol{x}) = \sum_{j \leq k} \mathcal{Q}_{jk} x_{j} x_{k}
"""
# Define and compile generic QUBO matrix structure
_COEFF_PATTERN = r"((?:\d+(?:\.\d*)?|\.\d+)(?:[eE][+-]?\d+)?)"
_LIN_REGEX = re.compile(
rf"""
([+\-])\s* # sign: + or -
{_COEFF_PATTERN}\s+ # local bias strength
([^\s]+) # variable name (no spaces)
""",
re.VERBOSE,
)
_SELF_LIN_REGEX = re.compile(
rf"""
([+\-])\s* # sign: + or -
{_COEFF_PATTERN}\s+ # auto-interaction strength
([^\s\*]+?) # first variable name (non-greedy, no '*')
\^2\b # literal exponent '^2' followed by word boundary
""",
re.VERBOSE,
)
_QUAD_REGEX = re.compile(
rf"""
([+\-])\s* # sign
{_COEFF_PATTERN}\s+ # interaction strength
([^\s\*^]+)\s* # first variable (anything except space, *, ^)
\*\s* # * with optional spaces
([^\s\*^]+) # second variable
""",
re.VERBOSE,
)
_OFFSET_REGEX = re.compile(
rf"""
([+\-])\s* # sign: + or -
{_COEFF_PATTERN} # constant value
(?=\s*[+\-]|$) # lookahead for next term or end of string
""",
re.VERBOSE,
)
# Constructor for a QUBO solver instance
# ------------------------------------------------------------------------
def __init__(
self,
qubo_matrix: npt.NDArray[np.floating] | None = None,
qubo_offset: float | None = 0.0,
qubo_variable_labeling: dict[str, int] | None = None,
) -> None:
"""
Initializes the QUBO solver.
Notes
-----
Users can optionally provide a QUBO matrix at initialization if
they are solving a single problem instance.
For workflows involving multiple QUBO matrices (e.g., in loops
or batch evaluations), defer setting the matrix and use
:meth:`~QUBOSolver.set_problem()` at a later stage in the solver
pipeline.
"""
self.qubo_matrix = None
self.offset = 0.0
self.n_binaries = 0
self.n_interactions = 0
self.spinglass_couplings = {}
self._variables_to_qubo_index_map = {}
self.cost = []
self.solution = []
self.status = "\n".join(
[
"QUBO solver object from qtealeaves has been initialized.",
"Still missing the QUBO problem instance to be solved.",
]
)
self.runtimes = {}
if qubo_matrix is not None:
self.set_qubo_problem(
qubo_matrix=qubo_matrix,
qubo_offset=qubo_offset,
qubo_variable_labeling=qubo_variable_labeling,
)
# ------------------------------------------------------------------------
# ------------------------------------------------------------------------
# Properties
# ------------------------------------------------------------------------
@property
def solution_space_dim(self) -> int:
"""
Returns the dimension of the solution space of all possible
QUBO decision-variable configurations.
"""
return 2**self.n_binaries
@property
def best_config(self) -> npt.NDArray[np.integer]:
"""
Returns the optimal QUBO problem solution found by the solver.
"""
if not self.solution:
raise QUBOSolverError("Empty solution space!")
return (self.solution)[0]
@property
def best_cost(self) -> float:
"""Returns the optimal QUBO cost value found by the solver."""
if not self.cost:
raise QUBOSolverError("Empty solution space!")
return (self.cost)[0]
@property
def variables_to_qubo_index_map(self) -> dict[str, int]:
"""
Returns the mapping between orginal problem variables and QUBO
indices.
"""
if not self._variables_to_qubo_index_map:
logger.warning(
"Trivial mapping to relate decision variables to QUBO indices."
"\nNo knowledge about variable name in the original optimization problem."
)
return self._variables_to_qubo_index_map
# ------------------------------------------------------------------------
# ------------------------------------------------------------------------
# Instance methods
# ------------------------------------------------------------------------
[docs]
def set_qubo_problem(
self,
qubo_matrix: npt.NDArray[np.floating],
qubo_offset: float | None = 0.0,
qubo_variable_labeling: dict[str, int] | None = None,
) -> None:
"""
Sets or updates the QUBO matrix, with validation and internal
preprocessing.
Arguments
---------
qubo_matrix : numpy.ndarray[float]
The QUBO matrix to be solved.
qubo_offset : float, optional
The QUBO problem energy cost offset.
Ignored during the optimization.
Default to 0.0.
qubo_variable_labeling : dict[str, int], optional
Labeling scheme from the original problem variable names
to the index assign in the QUBO problem reformulation.
Keys are variable names, values are the corresponding indices
within the one-dimensional embedding.
Default to :obj:`None`, i.e., a trivial labeling scheme
where the :math:`j`-th decision variable takes QUBO index
:math:`j`.
Returns
-------
:obj:`None`
"""
qubo_matrix = np.asarray(qubo_matrix, dtype=float)
# Validate input problem
if np.iscomplexobj(qubo_matrix):
raise QUBOProblemError("QUBO matrix must not be a complex matrix!")
# Adjust the QUBO problem for qtealeaves
# Already upper triangular – use as-is
if np.allclose(qubo_matrix, np.triu(qubo_matrix), rtol=1e-9, atol=1e-12):
self.qubo_matrix = qubo_matrix
# Symmetric – convert to upper triangular
elif np.allclose(qubo_matrix, qubo_matrix.T, rtol=1e-9, atol=1e-12):
upper = np.triu(qubo_matrix)
diag_indices = np.diag_indices_from(upper)
off_diag = upper.copy()
off_diag[diag_indices] = 0
self.qubo_matrix = upper.copy()
self.qubo_matrix[~np.eye(upper.shape[0], dtype=bool)] *= 2
logger.warning(
"Input QUBO matrix is symmetric. It has been converted to upper"
" triangular form with off-diagonal elements multiplied by 2.",
)
# Non-valid input from the user
else:
raise QUBOProblemError("QUBO matrix must be symmetric or upper triangular!")
# Other specific attributes for the solver
self.offset = qubo_offset
self._variables_to_qubo_index_map = qubo_variable_labeling or {}
self.n_binaries = qubo_matrix.shape[0]
self.n_interactions = np.count_nonzero(self.qubo_matrix)
self.n_interactions -= np.count_nonzero(np.diagonal(qubo_matrix).copy())
self.status = " ".join(
[
"QUBO solver object from qtealeaves has been initialized.",
"QUBO problem instance to be solved has also been initialized.",
]
)
[docs]
def solve(
self,
solver: str | None = "tensor-network",
n_solutions: int | None = 1,
**kwargs,
) -> None:
r"""
Solves the QUBO problem using a specific ``solver``.
Arguments
---------
solver : str, optional
Method to solve the QUBO problem.
Options are 'tensor-network', 'exact' and 'brute-force'.
Default to 'tensor-network'.
n_solutions : int, optional
Number of optimal/near-optimal solutions to compute with
the given ``solver``.
Default to 1, i.e., the solver searches only for the
optimal QUBO solution.
**kwargs
Additional optional arguments, depending on the chosen solver.
List of optional arguments can be found in the *Details*
section.
Returns
-------
:obj:`None`
Raises
------
:exc:`ValueError`
If ``n_solutions`` is not a positive integer.
:exc:`NotImplementedError`
If ``solver`` is not recognized.
:exc:`QUBOSetUpError`
If tensor-network solver is selected but missing
required parameters.
Details
-------
Additional arguments can be passed to this method, depending
on the chosen ``solver``.
List of available optional arguments:
- ``renormalize_couplings`` : bool (quantum-inspired only)
Whether to renormalize the spin glass Hamiltonian couplings
within :math:`\left[-1, 1\right]`.
Default to :obj:`True`.
- ``tn_conv_params`` : QUBOConvergenceParameter (tensor-network only)
The convergence parameters for the QUBO solver simulation
based on tensor-network methods.
It contains the main parameters for the ansatz, the ground-state
search and the backend at the tensor level.
It is required to run the tensor-network solver.
"""
# Validate arguments
if not isinstance(n_solutions, int) or n_solutions < 1:
raise ValueError(
f"'nsolutions' must be an integer >= 1, got {n_solutions}."
)
if solver not in {"brute-force", "exact", "tensor-network"}:
raise NotImplementedError(f"QUBO solver {solver} not available.")
# Optimization starts here
self.status = "FAILED"
# Brute force
if solver == "brute-force":
self._solve_via_brute_force(n_solutions)
self._check_spectrum_degeneracy()
return
# Quantum-inspired method starts here
# Mapping into spin glass
self._to_spinglass_couplings()
if kwargs.get("renormalize_couplings", True):
self._renormalize_spinglass_couplings()
# Ground-state search: tensor networks
if solver == "tensor-network":
tn_conv_params = kwargs.get("tn_conv_params")
if tn_conv_params is None:
raise QUBOSetUpError(
"Cannot run tensor-network based optimization. "
"Pass hyper-parameters!"
)
self._solve_via_tn_ground_state_search(tn_conv_params, n_solutions)
# Ground-state search: exact
elif solver == "exact":
self._solve_via_exact_sorting(n_solutions)
# Check ground-state degeneracy
self._check_spectrum_degeneracy()
[docs]
def evaluate(self, binary_configuration: npt.NDArray[np.integer]) -> float:
r"""
Evaluates the QUBO cost function for a given binaries configuration.
The QUBO cost function is defined as
.. math::
f_{\mathrm{\scriptscriptstyle QUBO}}
=
\sum_{j < j' = 1}^n
Q_{j j'} x_j x_{j'}
+
\sum_{j=1}^n Q_{jj} x_j
where :math:`n` is the number of binary decision variables in the
QUBO problem and :math:`\mathcal{Q}` is the QUBO matrix.
Arguments
---------
binary_configuration : numpy.ndarray[int]
The :math:`n`-dimensional bitstring representing the values of
the QUBO problem decision variables.
Returns
-------
cost : float
The value of the QUBO cost function evaluated on the given
binaries configuration ``binary_configuration``.
Raises
------
:exc:`QUBOSolverError`
If the input configuration does not match the problem dimension
or contains invalid values.
"""
# Validate instance initialization
if self.qubo_matrix is None:
raise QUBOSetUpError("QUBO problem matrix has not been initialized.")
# Validate the input bit-string
if binary_configuration.shape != (self.n_binaries,):
raise QUBOSolverError(
f"Expected bitstring of shape ({self.n_binaries},), "
f"got {binary_configuration.shape} instead."
)
if not np.all((binary_configuration == 0) | (binary_configuration == 1)):
raise QUBOSolverError("Input must be a binary array (only 0s and 1s).")
# Evaluate the cost function
# NOTE: we don't need to transpose the first array
# NOTE: because bit-strings are one-dimensional objects
qubo_cost = float(
binary_configuration @ self.qubo_matrix @ binary_configuration
)
qubo_cost += self.offset
return qubo_cost
def _to_spinglass_couplings(self) -> None:
r"""
Derives the corresponding spin glass model of the QUBO
problem.
Arguments
---------
:obj:`None`
Returns
-------
:obj:`None`
Details
-------
The couplings of the spin glass Hamiltonian are obtained from
the QUBO matrix elements by transforming binary variables to
spin-:math:`1 \mathbin{/} 2` variables via the following linear
transformation:
.. math::
x = \dfrac{1 + \sigma}{2}
where :math:`x \\in \\left\{0, 1\\right\}` is the binary variable
and :math:`\\sigma \\in \\left\{-1, 1\\right\}` is the classical
spin variable.
The transformation implementing the mapping between the QUBO
cost function and the spin glass Hamiltonian is as follows:
.. math::
\\begin{align}
& \\mathcal{E}_0 = Cte
+ \\dfrac{1}{2} \\sum_{j=1}^{n_Q} \\mathcal{Q}_{jj}
+ \\dfrac{1}{4} \\sum{k > j}^{n_Q} \\mathcal{Q}_{jk} \\
& \\mathrm{h}_j = \\dfrac{1}{2} \\mathcal{Q}_{jj}
+ \\dfrac{1}{4} \\sum_{k > j}^{n_Q} \\mathcal{Q}_{jk}
+ \\dfrac{1}{4} \\sum_{k = 1}^{j - 1} \\mathcal{Q}_{kj} \\
& \\mathrm{J}_{jk} = \\dfrac{1}{4} \\mathcal{Q}_{jk}
\\end{align}
and the spin glass Hamiltonian reads as
.. math::
\\mathcal{H}
=
\\mathcal{E}_0
+ \\sum_{j=1}^{n} \\mathrm{h}_j \\sigma_{j}
+ \\sum_{j < k}^{n} \\mathrm{J}_{jk} \\sigma_j \\sigma_k
Specifically, the spin glass coupling defining the model are
arranged as the following python dictionary:
- ``offset``: the constant term :math:`\\mathcal{E}_0` proportional
to the identity, which is just an energy shift and will be ignored
during the optimization;
- ``one-body``: the set of single-body couplings :math:`\\mathrm{h}_j`,
i.e., the set of local longitudinal magnetic fields (biases), one
for each spin-:math:`1 \\mathbin{/} 2` in the spin glass system;
- ``two-body``: the set of two-body (in general all-to-all) couplings
:math:`\\mathrm{J}_{jk}` describing the interactions between pairs
of spin-:math:`1 \\mathbin{/} 2`.
"""
# Validate instance initialization
if self.qubo_matrix is None:
raise QUBOSetUpError("QUBO problem matrix has not been initialized.")
qjj = np.diagonal(self.qubo_matrix).copy()
qij_upper = np.triu(self.qubo_matrix, k=1)
# Energy offset
cte = self.offset + 0.5 * np.sum(qjj) + 0.25 * np.sum(qij_upper)
# One-body couplings, aka, spin glass local magnetic fields
magnetic_fields = 0.5 * qjj
for jj in range(self.n_binaries):
magnetic_fields[jj] += 0.25 * (
np.sum(qij_upper[jj, jj + 1 :]) + np.sum(qij_upper[:jj, jj])
)
# Two-body J matrix: upper triangular only
Jmat_upper = 0.25 * qij_upper
# Update the corresponding class instance attribute
self.spinglass_couplings = {
"offset": cte,
"one-body": magnetic_fields,
"two-body": Jmat_upper,
}
def _solve_via_brute_force(self, n_solutions: int | None = 1) -> None:
"""
Solves the QUBO problem via brute-force by generating all the
possible decision-variable configurations.
Arguments
---------
nsolutions : int, optional
Number of brute-force solutions to compute.
If greater than 1, the method returns the first ``n_solutions``
bit-strings with the lowest QUBO cost values.
Default to 1.
Returns
-------
:obj:`None`
Notes
-------
This solver is limited to problem sizes less than 30 binaries.
No mapping to a spin glass is needed for this solver.
"""
if self.n_binaries > 30:
raise RuntimeError(
"\n\n**********************************************************\n"
"**********************************************************\n"
"With great power comes great responsibility!\n"
"The space of all the possible binaries configuration\n"
"is huge, but sadly, your poor old CLASSICAL computer's\n"
"memory just ain't up to the task!\n"
f"Crazy number of binaries --> {self.n_binaries} :O."
"\n**********************************************************\n"
"**********************************************************\n"
)
self.status = "BRUTE-FORCE METHOD STARTED"
st_to_sol = time()
# Generate all the possible configurations
solution_pool = []
for bitstring in range(self.solution_space_dim):
tmp_config = np.array(
[(bitstring >> i) & 1 for i in range(self.n_binaries)], dtype=np.uint8
)
tmp_cost = self.evaluate(tmp_config)
## Push to the pool and mantain size
## NOTE: heapq always pops the smallest item
## NOTE: So, put a - sign to track the largest
## NOTE: cost currently in the heap
if len(solution_pool) < n_solutions:
heapq.heappush(solution_pool, (-tmp_cost, tmp_config.copy()))
else:
heapq.heappushpop(solution_pool, (-tmp_cost, tmp_config.copy()))
# Sort back the results
solution_pool = sorted(
[(-cost, config) for cost, config in solution_pool], key=lambda x: x[0]
)
self.runtimes["time-to-solution"] = time() - st_to_sol
self.status = "WHOLE SOLUTION SPACE EXPLORED"
# Collect results
self.cost = [float(cost) for cost, _ in solution_pool]
self.solution = [config for _, config in solution_pool]
self.status = "SUCCESS"
def _solve_via_exact_sorting(self, n_solutions: int | None = 1) -> None:
"""
Finds the optimal solution(s) to the QUBO problem by mapping it
onto a spin glass Hamiltonian.
Then, finds the ground-state (or low-energy eigenstates) of this
Hamiltonian by performing an exact sorting of its diagonal.
Arguments
---------
n_solutions : int, optional
Number of QUBO problem optimal solutions to computed, i.e.,
the number of eigenstates to be extracted from the ordered
spin glass spectrum.
Returns
-------
:obj:`None`
"""
if self.n_binaries > 30:
raise RuntimeError(
"\n\n**********************************************************\n"
"**********************************************************\n"
"With great power comes great responsibility!\n"
"The Hilbert space is huge, but sadly, your\n"
"poor old CLASSICAL computer's memory just\n"
"ain't up to the task!\n"
f"Crazy number of qubits --> {self.n_binaries} :O."
"\n**********************************************************\n"
"**********************************************************\n"
)
st_to_sol = time()
# Construct the spin glass Hamiltonian
st_sg_mapping = time()
sg_ham = construct_exact_spinglass_hamiltonian(self.spinglass_couplings)
self.runtimes["hamiltonian mapping"] = time() - st_sg_mapping
# Solve the spectrum
st_diagonalization = time()
sg_spectrum = get_exact_spectrum(sg_ham, n_solutions)
self.status = "EXACT DIAGONALIZATION EXECUTED"
self.runtimes["hamiltonian diagonalization"] = time() - st_diagonalization
st_post_processing = time()
# Measure observables on the exact spectrum
obs = construct_exact_observables(self.n_binaries)
sg_config = measure_exact_observables(obs, sg_spectrum)["sz"]
self.status = "PROJECTIVE MEASUREMENTS PERFORMED"
# Map back to QUBO solutions
self.solution = [
np.array(spins_to_binaries(config), dtype=np.uint8) for config in sg_config
]
self.cost = [self.evaluate(config) for config in self.solution]
self.runtimes["post-processing"] = time() - st_post_processing
self.status = "SUCCESS"
self.runtimes["time-to-solution"] = time() - st_to_sol
def _renormalize_spinglass_couplings(
self, max_diff: float | None = 1e2, oom_threshold: float | None = 5e2
) -> None:
r"""
Renormalizes the couplings of the quantum Hamiltonian model
that encodes the QUBO problem instance to be solved to
ensure numerical stability.
The couplings are renormalized in the range
:math:`\\left[-1, 1\\right]` whenever the difference between
the largest and the smallest value of the couplings exceeds
``max_diff``, or their order of magnitude exceeds ``oom_threshold``.
Arguments
---------
max_diff : float, optional
If the difference between the largest and the smallest
values of the coupling strengths (zero-couplings excluded)
exceeds this number, the model couplings are enormalized
in :math:`\\left[-1, 1\\right]`.
Default to 100.
oom_threshold : float, optional
Threshold for the overall order of magnitude of couplings.
If any coupling exceeds this threshold, all couplings are
renormalized.
Default is 500.
Returns
-------
:obj:`None`
Notes
-----
The original spin glass model is stored in the class instance
attribute :attr:`~QUBOSolver.spinglass_couplings`.
The renormalization process also works in the presence of other
terms beyond pairwise interaction, which may arise in more
sophisticated quantum-inspired mappings.
The set of renormalized couplings overwrites the class instance
attribute :attr:`~QUBOSolver.spinglass_couplings`, with the addition
of the following keys:
- **``max``**: the rescaling factor to get back the original
Hamiltonian couplings.
"""
# Look-up for the original couplings
couplings = self.spinglass_couplings
# Ensure that not all couplings are zero
assert not all(
np.all(coupling == 0)
for coupling in couplings.values()
if isinstance(coupling, np.ndarray)
and np.issubdtype(coupling.dtype, np.number)
), QUBOProblemError(
"The spin glass Hamiltonian resulting from the "
"input QUBO is zero. Check-out your input QUBO problem."
)
# Compute the absolute values of the original couplings
abs_values = [
np.abs(coupling[np.nonzero(coupling)])
for coupling in couplings.values()
if isinstance(coupling, np.ndarray)
and np.issubdtype(coupling.dtype, np.number)
]
# Compute max and min of the absolute couplings
min_abs_coupling = min(value.min() for value in abs_values if value.size > 0)
max_abs_coupling = max(value.max() for value in abs_values if value.size > 0)
# Rescale the couplings if needed
rescaled_couplings = {"max": max_abs_coupling}
if (
max_abs_coupling >= max_diff * min_abs_coupling
or max_abs_coupling >= oom_threshold
):
rescaled_couplings.update(
{
key: (
coupling / max_abs_coupling
if isinstance(coupling, np.ndarray)
else coupling
)
for key, coupling in couplings.items()
}
)
else:
rescaled_couplings.update(couplings)
rescaled_couplings["max"] = 1
# Update class instance attibute
self.spinglass_couplings = rescaled_couplings
def _solve_via_tn_ground_state_search(
self,
convergence_params: QUBOConvergenceParameters,
n_eigenstates: int | None = 1,
) -> None:
r"""
Tensor-network based solver for QUBO problems.
Arguments
---------
convergence_params : QUBOConvergenceParameters
The convergence parameters for the QUBO solver simulation, for
example the bond dimension, the number of sweeps, and the
initial tensor network state.
n_eigenstates : int, optional
The number of optimal solutions to the QUBO problem to compute,
i.e., the number of eigenstates to be extracted from the mapped
spin glass spectrum.
Default to 1, i.e., only the ground state.
Returns
-------
:obj:`None`
Details
-------
Finds the optimal solution(s) to the QUBO problem by mapping it onto
a spin glass Hamiltonian.
Then, find the ground-state (or low-energy eigenstates) of this
Hamiltonian by performing a variational optimization via tensor
network methods.
Specifically, the wave function of the spin glass system is represented
as a tensor network (TN) with a given topology and bond dimension, and
the energy is optimized using numerical methods from :mod:`qtealeaves`.
Notes
-----
If ``n_eigenstates`` is greater than 1, excited states are obtained
via sampling of spin-glass configurations from the converged tensor
network state.
In this case, there is no guarantee that the sampled bit-strings are
exactly the first ``n_eigenstates`` eigenvectors in the Hamiltonian
spectrum (in general, they don't).
"""
# Check method arguments type
if not isinstance(convergence_params, QUBOConvergenceParameters):
raise TypeError(
"Expected an instance of QUBOConvergenceParameters, "
f"but got {type(convergence_params).__name__}"
)
# Look-ups from convergence parameters
tn_state_filename = convergence_params.tn_state_file
psi0_filename = convergence_params.tn_initial_state_filename
tn_ansatz = convergence_params.tn_ansatz
backend = convergence_params.tensor_backend
tfields_ratio = convergence_params.transverse_field_ratio
st_to_sol = time()
# Construct the spin glass Hamiltonian as a TN operator
st_sg_mapping = time()
hbiases = self.spinglass_couplings["one-body"]
sg_ints = self.spinglass_couplings["two-body"]
tn_ops = TNSpin12Operators()
tn_obs = TNObservables()
tn_obs += Local(name="sz", operator="sz")
tn_obs += State2File(tn_state_filename, "U")
sg_ham = QuantumModel(1, "L", name="QUBO-Spin Glass Model")
sg_ham += RandomizedLocalTerm(operator="sz", coupling_entries=hbiases)
sg_ham += TwoBodyAllToAllTerm1D(operators=["sz", "sz"], coupling_matrix=sg_ints)
# Add transverse quantum perturbation if requested
transverse_hfields = generate_perturbation_transverse_fields(
couplings=self.spinglass_couplings, ratio=tfields_ratio
)
if transverse_hfields is not None:
sg_ham += RandomizedLocalTerm(
operator="sx", coupling_entries=transverse_hfields
)
self.runtimes["hamiltonian mapping"] = time() - st_sg_mapping
# Setup simulation parameter dictionary
tn_params = {"L": self.n_binaries, "seed": convergence_params.randomseed}
# Initial state preparation
st_init_state = time()
psi0_state = None
if convergence_params.psi0 == "driver_product_state":
psi0_state = construct_spinglass_driver_state(
nsites=self.n_binaries,
bond_dimension=convergence_params.ini_bond_dimension,
numerical_noise=convergence_params.psi0_noise,
tn_ansatz=tn_ansatz,
tn_backend=backend,
)
self.runtimes["initial-state preparation"] = time() - st_init_state
if psi0_state:
psi0_state.save_pickle(psi0_filename)
tn_params["continue_file"] = psi0_filename
# Instantiate the tensor network simulator
simulator = QuantumGreenTeaSimulation(
model=sg_ham,
operators=tn_ops,
convergence=convergence_params,
observables=tn_obs,
folder_name_output=convergence_params.tn_output_filepath,
tn_type=convergence_params.tn_type,
mpo_mode=convergence_params.mpo_mode,
has_log_file=True,
store_checkpoints=False,
py_tensor_backend=backend,
)
# Run the ground-state search
st_ground_state_search = time()
simulator.run(params=tn_params, delete_existing_folder=True)
gs_search_result = simulator.get_static_obs(tn_params)
self.status = "GROUND-STATE SEARCH EXECUTED"
self.runtimes["ground-state search"] = time() - st_ground_state_search
# Get QUBO solution(s)
st_post_processing = time()
if n_eigenstates == 1: # Best QUBO bit-string only
sg_config = [round(ss) for ss in gs_search_result["sz"]]
try:
self.solution = [np.array(spins_to_binaries(sg_config), dtype=np.uint8)]
self.status = "NON-DEGENERATE GROUND STATE"
except ValueError:
## Get the optimized tensor-network state
psi_gs = read_tn_state(
filepath=tn_state_filename,
tn_ansatz=tn_ansatz,
tn_backend=backend,
)
psi_gs.convergence_parameters = convergence_params
## Sample and returns spin glass configurations
sg_configs = [
np.array(list(config), dtype=np.uint8)
for config in psi_gs.sample_n_unique_states(4).keys()
]
# Flipping spin glass configuration to match qtealeaves
sg_configs = [1 - bitstring for bitstring in sg_configs]
self.solution = [
np.array(config, dtype=np.uint8) for config in sg_configs
]
self.status = "EIGENSTATES REQUIRED"
else: # Get more low-energy QUBO bit-strings
## Get the optimized tensor-network state
psi_gs = read_tn_state(
filepath=tn_state_filename,
tn_ansatz=tn_ansatz,
tn_backend=backend,
)
## Sample and returns spin glass configurations
sg_configs = [
np.array(list(config), dtype=np.uint8)
for config in psi_gs.sample_n_unique_states(n_eigenstates).keys()
]
# Flipping spin glass configuration to match qtealeaves
sg_configs = [1 - bitstring for bitstring in sg_configs]
self.solution = [np.array(config, dtype=np.uint8) for config in sg_configs]
self.status = "EIGENSTATES REQUIRED"
self.cost = [self.evaluate(config) for config in self.solution]
self.runtimes["post-processing"] = time() - st_post_processing
# Finalize
self.status = "SUCCESS"
self.runtimes["time-to-solution"] = time() - st_to_sol
def _check_spectrum_degeneracy(self) -> None:
"""
Checks degeneracy in the ground state of the
solution space.
Arguments
---------
:obj:`None`
Returns
-------
:obj:`None`
"""
if len(self.solution) == 1:
return
if self.cost.count(self.best_cost) > 1:
self.status += " | DEGENERATE GROUND STATE"
else:
self.status = " | NON-DEGENERATE GROUND STATE"
[docs]
def prettyprint(self) -> None:
"""
Pretty prints method to visualize the QUBO matrix to standard
output.
Arguments
---------
:obj:`None`
Returns
-------
:obj:`None`
"""
# Validate instance initialization
if self.qubo_matrix is None:
raise QUBOSetUpError("QUBO problem matrix has not been initialized.")
# Determine the maximum width of each column
col_widths = [
max(len(str(item)) for item in col) for col in zip(*self.qubo_matrix)
]
# Construct the format string for each row
row_format = " | ".join(["{{:{}}}".format(width) for width in col_widths])
# Print each row using the format string
print("QUBO matrix\n-------------")
for row in self.qubo_matrix:
print(row_format.format(*row))
# pylint: disable=too-many-branches
def _save_to_lp_file(
self,
filename: str | Path,
problem_name: str | None = "Quantum TEA QUBO Problem",
) -> None:
"""
Saves the QUBO problem to a standard LP file.
Arguments
---------
filename : str | pathlib.Path
Full path to the output LP file.
If the correct extension ``.lp``is not provided,
it will be added automatically.
problem_name : str, optional
Name of the exported QUBO problem in the LP file.
Default to 'Quantum TEA QUBO Problem'.
Returns
-------
:obj:`None`
"""
# Validate instance initialization
if self.qubo_matrix is None:
raise QUBOSetUpError("QUBO problem has not been initialized.")
if self._variables_to_qubo_index_map:
variable_names = [
vname
for vname, _ in sorted(
self._variables_to_qubo_index_map.items(), key=lambda item: item[1]
)
]
else:
variable_names = [f"x_{jj}" for jj in range(self.n_binaries)]
assert (
len(variable_names) == self.n_binaries
), "QUBO variable name list length must match QUBO matrix size."
filename = Path(filename)
if filename.suffix.lower() != ".lp":
filename = filename.with_suffix(".lp")
lp_lines = [
r"\ This file has been generated by QUBOSolver from qtealeaves",
r"\ ENCODING=UTF-8",
rf"\ Problem name: {problem_name}",
rf"\ qtealeaves version (from the Quantum TEA suite): {__version__}",
"",
"Minimize",
]
obj_linear_terms = []
obj_quadratic_terms = []
# Linear terms (QUBO matrix diagonal)
for jj in range(self.n_binaries):
coeff = self.qubo_matrix[jj, jj]
if abs(coeff) > 1e-15:
sign = "+" if coeff >= 0 else "-"
qubo_term = f"{sign} {abs(coeff):.15f} {variable_names[jj]}"
obj_linear_terms.append(qubo_term)
# Quadratic terms (QUBO matrix off-diagonal)
for jj in range(self.n_binaries):
for kk in range(jj + 1, self.n_binaries):
coeff = self.qubo_matrix[jj, kk]
if abs(coeff) > 1e-15:
sign = "+" if coeff >= 0 else "-"
qubo_term = (
f"{abs(coeff):.15f} {variable_names[jj]}*{variable_names[kk]}"
)
obj_quadratic_terms.append((sign, qubo_term))
# Offset
offset_str = ""
if abs(self.offset) > 1e-15:
sign = "+" if self.offset >= 0 else "-"
offset_str = f"{sign} {abs(self.offset):.15f}"
# Compose the objective function
obj_lines = obj_linear_terms[:]
if obj_quadratic_terms:
first_sign, first_qubo_term = obj_quadratic_terms[0]
obj_lines.append(f"{first_sign} [ {first_qubo_term}")
for sign, qubo_term in obj_quadratic_terms[1:]:
obj_lines.append(f"{sign} {qubo_term}")
obj_lines.append(" ] / 2")
obj_lines.append(offset_str)
# Break long lines for LP format readability
wrapped_lines = []
current_line = " obj: "
for line in obj_lines:
if len(current_line) + len(line) + 1 > 80:
wrapped_lines.append(current_line.rstrip())
current_line = " " + line + " "
else:
current_line += line + " "
wrapped_lines.append(current_line.rstrip())
lp_lines.extend(wrapped_lines)
# Put QUBO variables
lp_lines.extend(["", "Subject To", "", "Bounds"]) # No constraints for QUBO!
for var in variable_names:
lp_lines.append(f" 0 <= {var} <= 1")
lp_lines.append("\nBinaries")
line = " "
for var in variable_names:
next_entry = f"{var} "
if len(line) + len(next_entry) > 80:
lp_lines.append(line.rstrip())
line = " "
line += next_entry
if line.strip():
lp_lines.append(line.rstrip())
# End
lp_lines.append("End")
# Write to file
with filename.open("w", encoding="utf-8") as lp_file:
lp_file.write("\n".join(lp_lines) + "\n")
[docs]
def save_to_file(
self, filename: str | Path, qubo_format: str | None = "lp"
) -> None:
"""
Saves the QUBO problem to a log file in a given ``qubo_format``.
Arguments
---------
filename : str | pathlib.Path
Full path to the output log file.
If the correct extension ``qubo_format``is not provided,
it will be added automatically.
qubo_format : str, optional
The (allowed) extension of the output log file.
Currently available extensions are 'lp' (suggested), 'json'
and 'text'.
Default to 'lp'.
Returns
-------
:obj:`None`
"""
# Validate instance initialization
if self.qubo_matrix is None:
raise QUBOSetUpError("QUBO problem has not been initialized.")
qubo_format = (qubo_format or "lp").lower()
if qubo_format not in ["json", "txt", "lp"]:
raise NotImplementedError(f"Format '{qubo_format}' not supported.")
filename = Path(filename)
filename.parent.mkdir(parents=True, exist_ok=True)
expected_ext = f".{qubo_format}"
if filename.suffix.lower() != expected_ext:
filename = filename.with_suffix(expected_ext)
# Save to standard LP file
if qubo_format == "lp":
self._save_to_lp_file(filename)
# JSONify the QUBO matrix and the offset
if qubo_format == "json":
jsonified_qubo = self.qubo_matrix.tolist()
with codecs.open(filename, "w", encoding="utf-8") as out_file:
json.dump(
obj={"qubo_matrix": jsonified_qubo, "offset": self.offset},
fp=out_file,
separators=(",", ":"),
sort_keys=True,
indent=4,
)
# Textify the QUBO matrix and the offset
elif qubo_format == "txt":
with open(filename, "w", encoding="utf-8") as out_file:
out_file.write("# " + "-" * 79)
out_file.write(
f"\n# [n = {self.n_binaries}] \t n: QUBO problem binaries\n"
)
out_file.write("# " + "-" * 79 + "\n#\n")
out_file.writelines(
f"{row_idx}\t{col_idx}\t{value:.15e}\n"
for (row_idx, col_idx), value in np.ndenumerate(self.qubo_matrix)
if col_idx >= row_idx
)
out_file.write(f"# Offset: {self.offset:.15e}\n")
out_file.write("# END\n")
# ------------------------------------------------------------------------
# ------------------------------------------------------------------------
# Static methods
# ------------------------------------------------------------------------
# pylint: disable=too-many-branches
@staticmethod
def _read_from_lp_file(
filename: str | Path,
) -> tuple[npt.NDArray[np.floating], float, dict[str, int]]:
"""
Reads a QUBO problem from a standard LP file and construct
the corresponding QUBO matrix.
Arguments
---------
filename : str | pathlib.Path
Full path (with the ".lp" extension) to the LP log file where
the QUBO matrix is stored.
Returns
-------
qubo_matrix : numpy.ndarray[float]
The QUBO matrix constructed from the QUBO problem.
offset : float
The QUBO energy cost offset.
variable_qubo_idx_map : dict[str, int]
Mapping dictionary from variable names to QUBO matrix
indices.
Keys are the variable names read from the LP ``filename``,
values are the corresponding QUBO indices within the one
dimensional embedding of the problem.
Notes
-----
At the moment, the LP-reader output dictionary ``variable_qubo_idx_map``
is not used to perform any calculations within the solver pipeline.
Anyway, users can request it to the solver in order to perform
post-processing and map back the QUBO solution to their original
combinatorial optimization problem.
See :meth:`~QUBOSolver.solve` method of this class.
When reading from an LP file, variables in the QUBO problem are
assigned indices in the order they are encountered in the linear
terms of the QUBO objective function.
This ordering remains consistent within the solver, but it may
differ from the original problem formulation the user intended to
solve.
For this reason, the ``variable_qubo_idx_map`` is stored as a solver
instance attribute and can be post-processed to recover the position
of each binary variable within the 1D chain where the problem is
embedded.
Finally, the construction of the QUBO matrix assumes that interaction
terms appear only once in the LP file, meaning that symmetric terms
have already been accounted for by doubling the interaction
coefficient.
"""
filename = Path(filename)
variable_qubo_idx_map = {}
qubo_sign = 1.0
offset = 0.0
linear_terms = defaultdict(float)
quadratic_terms = defaultdict(float)
next_var_index = 0
objective_lines = []
in_objective = False
# Read LP file lines
with filename.open("r", encoding="utf-8") as file:
for line in file:
stripped_line = line.strip().lower()
if stripped_line.startswith(("minimize", "maximize")):
in_objective = True
qubo_sign = -1.0 if stripped_line.startswith("maximize") else 1.0
continue
if in_objective and stripped_line.startswith(
("subject to", "bounds", "binary", "end")
):
break
if in_objective:
objective_lines.append(line.strip())
# Concatenate and clean up multiline objective
full_read_obj = QUBOSolver._flatten_lp_qubo_obj(" ".join(objective_lines))
# Match and extract linear terms h_j * x_j
# NOTE: self-interacting terms counts as linear h_j * (x_j)^2
for match in list(QUBOSolver._LIN_REGEX.finditer(full_read_obj)):
sign, coeff, binary = match.groups()
if "*" in binary or "^" in binary:
continue
if binary not in variable_qubo_idx_map:
variable_qubo_idx_map[binary] = next_var_index
next_var_index += 1
linear_terms[binary] += float(coeff) * (-1 if sign == "-" else 1)
full_read_obj = full_read_obj.replace(match.group(0), "", 1)
full_read_obj = re.sub(r"\s{2,}", " ", full_read_obj).strip()
for match in list(QUBOSolver._SELF_LIN_REGEX.finditer(full_read_obj)):
sign, coeff, binary = match.groups()
if binary not in variable_qubo_idx_map:
variable_qubo_idx_map[binary] = next_var_index
next_var_index += 1
linear_terms[binary] += float(coeff) * (-1 if sign == "-" else 1)
full_read_obj = full_read_obj.replace(match.group(0), "", 1)
full_read_obj = re.sub(r"\s{2,}", " ", full_read_obj).strip()
# Match quadratic terms J_jk * x_j * x_k
for match in list(QUBOSolver._QUAD_REGEX.finditer(full_read_obj)):
sign, coeff, binary1, binary2 = match.groups()
for bin_var in (binary1, binary2):
if bin_var not in variable_qubo_idx_map:
variable_qubo_idx_map[bin_var] = next_var_index
next_var_index += 1
# Order does not matter
quadratic_terms[tuple(sorted([binary1, binary2]))] += float(coeff) * (
-1 if sign == "-" else 1
)
full_read_obj = full_read_obj.replace(match.group(0), "", 1)
full_read_obj = re.sub(r"\s{2,}", " ", full_read_obj).strip()
# Match offset constant term
for sign, cte in QUBOSolver._OFFSET_REGEX.findall(full_read_obj):
offset += float(cte) * (-1 if sign == "-" else 1)
# Final QUBO matrix build
n_vars = len(variable_qubo_idx_map)
qubo_matrix = np.zeros((n_vars, n_vars), dtype=np.float64)
qubo_idx_map = dict(variable_qubo_idx_map)
diag_idxs = [qubo_idx_map[vname] for vname in linear_terms]
qubo_matrix[diag_idxs, diag_idxs] += np.fromiter(
linear_terms.values(), dtype=np.float64
)
row_idx_list = []
col_idx_list = []
couplings = []
for (var1, var2), value in quadratic_terms.items():
row_idx, col_idx = qubo_idx_map[var1], qubo_idx_map[var2]
if row_idx > col_idx:
row_idx, col_idx = col_idx, row_idx
row_idx_list.append(row_idx)
col_idx_list.append(col_idx)
couplings.append(value)
qubo_matrix[
np.array(row_idx_list, dtype=int), np.array(col_idx_list, dtype=int)
] += np.array(couplings, dtype=np.float64)
return qubo_sign * qubo_matrix, qubo_sign * offset, variable_qubo_idx_map
@staticmethod
def _flatten_lp_qubo_obj(obj_string_from_lp: str) -> str:
"""
Normalizes the objective string by expanding/removing bracketed
quadratic terms and eliminating any trailing '/2' division.
Arguments
---------
obj_string_from_lp : str
The native parsed LP file string containing the
full QUBO objective function.
Returns
-------
processed_obj_string : str
The flattened and processed full QUBO objective function
string coming from the LP file.
Notes
-----
This could be the situation with standard LP file, for example
generated through the **CPLEX** Optimization suite.
"""
# Define a regular expression to capture standard bracket blocks
bracket_block_regex = re.compile(
r"""
# opening bracket
\[
# non-lazy capture of inner content,
# because [...] could be on multiple lines
# This also defines the 'inner' group
(?P<inner>.*?)
# closing bracket
\]
# optional '/2' with spaces
(?:\s*/\s*2)?
""",
re.DOTALL | re.VERBOSE,
)
# Define how to replace the matched bracket blocks
def _bracket_replacement(match: re.Match) -> str:
inner = match.group("inner")
clean = inner.replace("[", "").replace("]", "")
return f" {clean} "
# Substitute all bracketed expressions
processed_obj_string = bracket_block_regex.sub(
_bracket_replacement, obj_string_from_lp
)
# Clean-up and output
processed_obj_string = processed_obj_string.replace("[", "").replace("]", "")
return processed_obj_string
# ------------------------------------------------------------------------
# ------------------------------------------------------------------------
# Class methods
# ------------------------------------------------------------------------
# pylint: disable=too-many-nested-blocks
[docs]
@classmethod
def read(cls, filename: str | Path) -> Self:
"""
Initializes the QUBO problem from an input file.
Arguments
---------
filename : str | pathlib.Path
The full path to the file containing the QUBO problem as
a matrix in a given (allowed) format.
A check on file's extension is performed before reading.
Returns
-------
obj : :class:`qtealeaves.optimization.QUBOSolver`
An instance of QUBOSolver initialized with the read QUBO
matrix.
Notes
-----
**Important note for the users regarding LP file**: this method
allows you to read QUBO problems from standard LP files, such as
those generated by classical solvers like CPLEX.
If you don't trust the input file or the correctness of the QUBO
obtained this way, we recommend providing the QUBO matrix directly
through the object constructor or via the
:meth:`~QUBOSolver.set_qubo_problem` method.
"""
# Get the extension
filename = Path(filename)
ext = filename.suffix.lower()
offset = 0.0
var_index_map = None
# JSON
if ext == ".json":
try:
with filename.open("r", encoding="utf-8") as in_file:
qubo_data = json.load(in_file)
if isinstance(qubo_data, dict) and "qubo_matrix" in qubo_data:
qubo_matrix = np.array(qubo_data["qubo_matrix"])
offset = float(qubo_data.get("offset", 0.0))
else:
# Fallback for raw matrix
qubo_matrix = np.array(qubo_data)
except Exception as err:
raise ValueError(f"Failed to parse JSON QUBO file: {filename}") from err
# Plain TEXT
elif ext == ".txt":
try:
with filename.open("r", encoding="utf-8") as in_file:
lines = in_file.readlines()
# Extract numerical lines and offset
data_lines = []
for line in lines:
line = line.strip()
if not line or line.startswith("#"):
# Try to parse the offset if available
if line.lower().startswith("# offset:"):
try:
offset = float(line.split(":", 1)[1].strip())
except ValueError as exc:
raise ValueError(
f"Invalid offset format in file: {line}"
) from exc
continue
data_lines.append(line)
row_idxs, col_idxs, qij = np.loadtxt(
data_lines, delimiter="\t", usecols=(0, 1, 2), unpack=True
)
row_idxs = row_idxs.astype(int)
col_idxs = col_idxs.astype(int)
qubo_size = int(np.max([row_idxs.max(), col_idxs.max()]) + 1)
qubo_matrix = np.zeros((qubo_size, qubo_size))
qubo_matrix[row_idxs, col_idxs] = qij
except Exception as err:
raise ValueError(
f"Failed to parse .txt QUBO format in {filename}"
) from err
# CPLEX like LP files
elif ext == ".lp":
try:
qubo_matrix, offset, var_index_map = cls._read_from_lp_file(filename)
except Exception as err:
raise ValueError(f"Failed to read .lp file: {filename}") from err
# Fallback
else:
raise NotImplementedError(
f"Unsupported file extension '{ext}' for QUBO input."
)
# Output
return cls(qubo_matrix, offset, var_index_map)
# ------------------------------------------------------------------------
# ------------------------------------------------------------------------