Metadata-Version: 2.2
Name: bqs
Version: 0.2.1
Summary: Better Quantum Software
Home-page: https://github.com/MichiCatte/bqs.git
Author: Michele Cattelan
Author-email: michelecat97@gmail.com
License: License :: OSI Approved :: MIT License
Classifier: Development Status :: 2 - Pre-Alpha
Classifier: Intended Audience :: Developers
Classifier: Intended Audience :: Science/Research
Classifier: Programming Language :: Python :: 3.9
Classifier: Programming Language :: Python :: 3.10
Description-Content-Type: text/markdown
License-File: LICENSE
Requires-Dist: cirq>=1.4.1
Requires-Dist: networkx>=3.0
Requires-Dist: numpy>=1.24
Requires-Dist: scipy>=1.10
Dynamic: author
Dynamic: author-email
Dynamic: classifier
Dynamic: description
Dynamic: description-content-type
Dynamic: home-page
Dynamic: license
Dynamic: requires-dist
Dynamic: summary

# bqs — Better Quantum Software

`bqs` is a Python framework for building **scalable, parallelized Variational Quantum Algorithms (VQAs)**, with a particular focus on **hardware-efficient implementations in the NISQ era**.

The design and methodology implemented in this package are described in
**[Parallel Circuit Implementation of Variational Quantum Algorithms](https://www.nature.com/articles/s41534-025-00982-6)**
(preprint: [arXiv:2304.03037](https://arxiv.org/pdf/2304.03037)).

---

## Installation

```bash
pip install bqs
```

---

## Quick start

```python
import numpy as np
from bqs import QAOA, COBYLA

# Define a symmetric QUBO matrix
Q = np.array([[0.0, 0.5, 0.5],
              [0.5, 0.0, 0.5],
              [0.5, 0.5, 0.0]])

alg = QAOA(Q, p=2)
alg.set_optimizer(COBYLA)
alg.optimize(num_samples_training=100, num_iterations=50, opt_arguments={"display": False})

samples = alg.circuit.sample(num_samples=200)   # (200, 3) array of ±1 spins
energy  = alg.objective_function.evaluate_samples(samples)
```

---

## Background

### Problem formulation

`bqs` works with optimization problems expressed as **QUBOs** (Quadratic Unconstrained Binary Optimization):

$$E(\mathbf{x}) = \mathbf{x}^T Q \mathbf{x} + \text{offset}, \qquad x_i \in \{0, 1\}$$

Internally, all algorithms convert the QUBO to its equivalent **Ising Hamiltonian** via the substitution $x_i = (1 - s_i)/2$, $s_i \in \{-1, +1\}$:

$$H(\mathbf{s}) = \mathbf{s}^T J \mathbf{s} + \mathbf{h}^T \mathbf{s} + \text{offset}$$

where $J = \frac{1}{2}\text{triu}(Q, k=1)$ and $h_i = \frac{1}{2}\sum_j Q_{ij}$.

Users can supply either a QUBO matrix $Q$ to the constructor, or the Ising triple $(J, h, \text{offset})$ directly via the `from_ising` classmethod.

### QAOA

The Quantum Approximate Optimization Algorithm (QAOA) is a variational hybrid algorithm. It applies $p$ alternating layers of a cost unitary $e^{-i\gamma_l H}$ and a mixer unitary $e^{-i\beta_l \sum_i X_i}$ to an initial uniform superposition:

$$U(\boldsymbol{\gamma}, \boldsymbol{\beta}) = \prod_{l=1}^{p} e^{-i\beta_l H_\text{mix}} e^{-i\gamma_l H} \cdot H^{\otimes n}|0\rangle^{\otimes n}$$

The $2p$ real parameters $(\boldsymbol{\gamma}, \boldsymbol{\beta})$ are tuned by a classical optimizer to minimize the expected energy $\langle H \rangle$ estimated from circuit samples.

### Problem decomposition (pQAOA)

For large problems, the full QUBO may require more qubits than available hardware supports. `bqs` implements a **parallel circuit decomposition** strategy: the global problem $P$ is split into $k$ smaller sub-problems $\{SP_i\}$, each implemented on a smaller QAOA circuit. Samples from the sub-circuits are **glued back together** and evaluated against the global cost function.

This approach preserves solution quality while drastically reducing the per-circuit qubit count.

---

## Algorithm classes

All algorithm classes live in `bqs` and share a common interface.

### `QAOA`

Vanilla QAOA on the full problem.

```python
QAOA(
    Q,                        # (n, n) symmetric QUBO matrix
    offset=0.0,               # constant energy offset
    p=1,                      # number of QAOA layers
    hyperparameters=None,     # (2*p,) array; random in [0, 2π) if None
    qubit_order=None,         # (n,) int array mapping variable i → physical qubit index
    quantum_objects=None,     # AbstractQuantumObjects backend; CirqQuantumObjects if None
)
```

**Named constructor** — bypass the internal QUBO→Ising conversion when the Ising form is already available:

```python
QAOA.from_ising(
    J,                        # (n, n) upper-triangular Ising coupling matrix
    h,                        # (n,) Ising local-field vector
    offset=0.0,
    p=1,
    hyperparameters=None,
    qubit_order=None,
    quantum_objects=None,
)
```

---

### `pQAOA`

Parallel QAOA: each sub-problem slice runs its own QAOA circuit with **independent** $2p$ hyperparameters. The full problem must be decomposable as a concatenation of the slice variables.

```python
pQAOA(
    Q,                        # (n_full, n_full) full-model QUBO
    slice_Qs,                 # list of per-slice QUBO matrices
    offset=0.0,
    slice_offsets=None,       # per-slice offsets; all 0 if None
    p=1,
    hyperparameters=None,     # (n_slices, 2*p) array; random if None
    qubit_orders=None,        # list of per-slice qubit order arrays
    cartesian=True,           # True → Cartesian-product gluing; False → zip gluing
    quantum_objects=None,
)
```

**Named constructor:**

```python
pQAOA.from_ising(
    J, h, offset=0.0,
    slice_isings,             # list of (J_k, h_k) or (J_k, h_k, offset_k) tuples
    p=1,
    hyperparameters=None,
    qubit_orders=None,
    cartesian=True,
    quantum_objects=None,
)
```

**Sample gluing:**
With `cartesian=True`, if each slice circuit produces $m$ samples, the glued output has $m^k$ rows — every combination of slice outcomes. With `cartesian=False`, slices are zipped row-by-row, producing $m$ rows (each slice must produce the same number of samples).

---

### `pQAOASingleParameters`

Same as `pQAOA` but all slices **share a single set of $2p$ hyperparameters**.

```python
pQAOASingleParameters(
    Q, slice_Qs,
    offset=0.0, slice_offsets=None,
    p=1,
    hyperparameters=None,     # (2*p,) array; random if None
    qubit_orders=None,
    cartesian=True,
    quantum_objects=None,
)
```

**Named constructor:**

```python
pQAOASingleParameters.from_ising(
    J, h, offset=0.0,
    slice_isings,
    p=1,
    hyperparameters=None,
    qubit_orders=None,
    cartesian=True,
    quantum_objects=None,
)
```

---

### `SingleSliceQAOA`

A single shared circuit is run repeatedly for each copy of the slice. A **variable map** then remaps each copy's outputs back onto the full-model columns.

This is most useful when the global problem has a repetitive structure (e.g., identical sub-blocks), so one circuit suffices for all slices.

```python
SingleSliceQAOA(
    Q,                        # (n_full, n_full) full-model QUBO
    slice_Q,                  # (n_slice, n_slice) QUBO for the repeated circuit
    var_map,                  # int array (n_slice, num_slices):
                              #   var_map[i, k] = full-model column for slice var i in copy k
    offset=0.0,
    slice_offset=0.0,
    p=1,
    hyperparameters=None,
    qubit_order=None,
    quantum_objects=None,
)
```

**Named constructor:**

```python
SingleSliceQAOA.from_ising(
    J, h, offset=0.0,
    slice_J, slice_h, slice_offset=0.0,
    var_map,
    p=1,
    hyperparameters=None,
    qubit_order=None,
    quantum_objects=None,
)
```

**Example `var_map`** — 6-variable problem with 2-variable slice repeated 3 times:

```python
var_map = np.array([[0, 2, 4],   # slice var 0 → full-model columns 0, 2, 4
                    [1, 3, 5]])  # slice var 1 → full-model columns 1, 3, 5
```

---

## Running the optimization

All algorithm classes share the same three-step interface:

```python
# 1. Instantiate
alg = QAOA(Q, p=2)

# 2. Attach an optimizer
alg.set_optimizer(COBYLA)

# 3. Optimize
alg.optimize(
    num_samples_training=100,     # samples per cost-function evaluation
    num_iterations=200,           # max optimizer iterations
    opt_arguments={"display": False},
    # circuit_options={"parallel": True},  # run slice circuits in parallel (pQAOA)
)
```

After optimization the best hyperparameters are stored in `alg.hyperparameters` and `alg.circuit.hyperparameters`.

### Sampling results

```python
samples = alg.circuit.sample(num_samples=500)
# → np.ndarray (500, n) of ±1 spin values

samples_binary = alg.circuit.sample(num_samples=500, original_basis=True)
# → np.ndarray (500, n) of {0, 1} binary values

energy = alg.objective_function.evaluate_samples(samples)
# → float: mean Ising energy over the sample batch
```

---

## Complete example

```python
import numpy as np
from bqs import QAOA, pQAOA, pQAOASingleParameters, SingleSliceQAOA, COBYLA
from bqs.utils.bqm import qubo_to_ising

# 6-variable chain: (0,0)-(0,1), (0,0)-(1,0), (1,0)-(1,1), (1,1)-(2,0), (2,0)-(2,1)
N = 6
Q_full = np.zeros((N, N))
for i, j in [(0,1),(0,2),(2,3),(3,4),(4,5)]:
    Q_full[i,j] = Q_full[j,i] = 0.5

# Sub-problem slices
Q_slice1 = np.zeros((2, 2)); Q_slice1[0,1] = Q_slice1[1,0] = 0.5
Q_slice2 = np.zeros((4, 4))
for i, j in [(0,1),(1,2),(2,3)]:
    Q_slice2[i,j] = Q_slice2[j,i] = 0.5

def run(alg, num_samples=50, num_iter=20):
    alg.set_optimizer(COBYLA)
    alg.optimize(num_samples_training=num_samples, num_iterations=num_iter,
                 opt_arguments={"display": False})
    return alg.objective_function.evaluate_samples(
        alg.circuit.sample(num_samples=num_samples))

# Vanilla QAOA
print("QAOA energy:", run(QAOA(Q_full, p=2)))

# pQAOA: per-slice hyperparameters, Cartesian-product gluing
print("pQAOA energy:", run(pQAOA(Q_full, [Q_slice1, Q_slice2])))

# pQAOA: zip gluing
print("pQAOA (zip) energy:", run(pQAOA(Q_full, [Q_slice1, Q_slice2], cartesian=False)))

# Shared hyperparameters across slices
print("pQAOASingleParameters energy:", run(pQAOASingleParameters(Q_full, [Q_slice1, Q_slice2])))

# Single shared circuit for all slices
Q_ss = np.zeros((2,2)); Q_ss[0,1] = Q_ss[1,0] = 0.5
var_map = np.array([[0,2,4],[1,3,5]])
print("SingleSliceQAOA energy:", run(SingleSliceQAOA(Q_full, Q_ss, var_map)))

# Using Ising inputs directly (skips the internal QUBO→Ising conversion)
J, h, off = qubo_to_ising(Q_full)
print("QAOA (from_ising):", run(QAOA.from_ising(J, h, off, p=2)))
```

---

## Quantum backend

The `quantum_objects` parameter accepts any subclass of `AbstractQuantumObjects`. This allows you to plug in a different quantum simulator or hardware SDK without changing any algorithm code.

The default backend is `CirqQuantumObjects`, which uses Google Cirq's statevector simulator.

### Using the default backend explicitly

```python
from bqs import QAOA, CirqQuantumObjects

alg = QAOA(Q, quantum_objects=CirqQuantumObjects())
```

### Writing a custom backend

Subclass `AbstractQuantumObjects` and implement the required interface:

```python
from bqs import AbstractQuantumObjects
import numpy as np

class MyBackend(AbstractQuantumObjects):

    def qubit(self, index, **kwargs):
        """Return the backend's qubit object for physical qubit `index`."""
        ...

    def circuit(self, n_qubits):
        """Create and return an empty circuit for `n_qubits` qubits."""
        ...

    def append(self, circuit, operation):
        """Append `operation` to `circuit` in-place."""
        ...

    def sample(self, circuit, qubit_name_to_object, num_samples, **kwargs) -> np.ndarray:
        """
        Run `circuit` for `num_samples` shots.

        Returns
        -------
        np.ndarray of shape (num_samples, n) with values in {-1, +1}.
        qubit_name_to_object maps integer variable index i → qubit object.
        """
        ...

    def H(self, qubit):    ...   # Hadamard gate
    def X(self, qubit):    ...   # Pauli-X gate
    def Y(self, qubit):    ...   # Pauli-Y gate
    def Z(self, qubit):    ...   # Pauli-Z gate
    def CNOT(self, q1, q2):  ... # CNOT gate
    def rx(self, qubit, angle): ...   # Rx(angle) rotation
    def ry(self, qubit, angle): ...   # Ry(angle) rotation
    def rz(self, qubit, angle): ...   # Rz(angle) rotation

alg = QAOA(Q, quantum_objects=MyBackend())
```

---

## Optimizers

Three optimizers are included out of the box:

| Class | Description |
|---|---|
| `COBYLA` | Gradient-free; wraps `scipy.optimize.minimize` with method `'COBYLA'`. Good default choice. |
| `SPSA` | Simultaneous Perturbation Stochastic Approximation. Useful for noisy landscapes. |
| `MonteCarlo` | Random hyperparameter search; keeps the best set found. Useful as a baseline. |

All are passed as classes (not instances) to `set_optimizer`:

```python
alg.set_optimizer(COBYLA)
alg.set_optimizer(SPSA)
alg.set_optimizer(MonteCarlo)
```

### Writing a custom optimizer

Subclass `Optimizer` and implement `_run`:

```python
from bqs.utils.vqa_utils import CostFunction, VariationalCircuit
from bqs.utils.optimizers import Optimizer
from scipy.optimize import minimize, OptimizeResult

class NelderMead(Optimizer):

    def __init__(self, cost_function: CostFunction, circuit: VariationalCircuit):
        super().__init__(cost_function, circuit)

    def _run(self, objective_function, num_iterations: int, display: bool = False, **kwargs):
        result = minimize(
            objective_function,
            self.hyperparameters,
            method="Nelder-Mead",
            options={"maxiter": num_iterations, "disp": display},
            **kwargs,
        )
        return result   # must be a scipy OptimizeResult (has .x and .fun)

alg.set_optimizer(NelderMead)
```

---

## Building new algorithm classes

Beyond swapping the backend or the optimizer, you can add an entirely new *algorithm variant* to the framework — e.g. your own flavor of QAOA — by subclassing `BaseQAOA` (or `BaseVQA` directly, if you're not doing QAOA at all). This is the most powerful extension point, but also the least "plug-and-play" one, so read this section fully before starting.

### The three interfaces involved

Every algorithm in `bqs` is assembled from three pieces:

| Abstract class | Role | Required methods |
|---|---|---|
| `BaseVQA` / `BaseQAOA` | The object you instantiate and call `.optimize()` on | `_set_circuit`, `_set_objective_function`, `get_circuit` |
| `VariationalCircuit` | Builds the circuit from hyperparameters and samples it | `get_circuit_from_hyperparameters`, `sample` |
| `CostFunction` | Scores a batch of samples | `evaluate_samples` |

A new algorithm variant usually means writing a `BaseQAOA` subclass paired with a matching `VariationalCircuit` subclass (and, only if your scoring differs from plain Ising energy, a `CostFunction` subclass). `optimize()` itself is already implemented on `BaseQAOA` and rarely needs to be overridden.

### How construction works

`BaseVQA.__init__` takes the constructor arguments as keyword attributes and sets them on `self` *before* calling the hooks that need them:

```python
class BaseVQA(ABC):
    def __init__(self, **attributes):
        for name, value in attributes.items():
            setattr(self, name, value)
        self.optimizer = None
        self._validate_hyperparameters()
        self._set_circuit()
        self._set_objective_function()
```

So a subclass never assigns anything to `self` before calling `super().__init__()` — it just forwards everything as keywords, and by the time `_validate_hyperparameters`/`_set_circuit`/`_set_objective_function` run, every attribute is already there:

```python
class MyQAOAVariant(BaseQAOA):
    def __init__(self, Q, my_extra_param, offset=0.0, p=1,
                 hyperparameters=None, qubit_order=None, quantum_objects=None):
        super().__init__(Q=Q, my_extra_param=my_extra_param, offset=offset, p=p,
                          hyperparameters=hyperparameters, qubit_order=qubit_order,
                          quantum_objects=quantum_objects)
```

### Adding a `from_ising` constructor

If your variant should also accept `(J, h, offset)` directly, don't branch on it inside `__init__`. Build the attribute dict for the Ising case and hand it to `cls._from_attributes(...)`, which allocates the instance and runs the same `BaseVQA.__init__` path — no `cls.__new__` boilerplate to write yourself:

```python
@classmethod
def from_ising(cls, J, h, offset=0.0, p=1, my_extra_param=None, **kwargs):
    return cls._from_attributes(
        Q=None,
        _J_in=np.asarray(J, dtype=float),   # stash Ising inputs under a private name
        _h_in=np.asarray(h, dtype=float),
        offset=offset, p=p, my_extra_param=my_extra_param, **kwargs,
    )
```

Your `_set_circuit()` then checks `if self.Q is not None` and falls back to `self._J_in` / `self._h_in` otherwise — that's the branch point, not the constructor.

### Non-flat hyperparameters

`BaseQAOA` assumes a flat `(2*p,)` hyperparameter array by default, via the `_validate_hyperparameters` hook:

```python
class BaseQAOA(BaseVQA):
    def _validate_hyperparameters(self):
        self.hyperparameters = self._validate_hp_1d(self.hyperparameters, self.p)
```

If your variant's hyperparameters don't fit that shape — e.g. `pQAOA`'s per-slice `(n_slices, 2*p)` array — override `_validate_hyperparameters` to normalize/validate whatever shape you need; everything else about construction stays the same:

```python
class pQAOA(BaseQAOA):
    def _validate_hyperparameters(self):
        n_slices = len(self.slice_Qs) if self.Q is not None else len(self._slice_isings)
        self.hyperparameters = self._validate_hp_2d(self.hyperparameters, n_slices, self.p)
```

### Worked example

A small variant of `QAOA` that adds a caller-supplied bias to the local fields before sampling:

```python
import numpy as np
from bqs.algorithms.qaoa import BaseQAOA
from bqs.utils import QAOACircuit, QAOACostFunction

class WarmStartQAOA(BaseQAOA):
    """QAOA that nudges the Ising local fields by a fixed 'warm_start' bias
    before the circuit is built — a simple way to bias sampling toward a
    known good solution."""

    def __init__(self, Q, warm_start, offset=0.0, p=1,
                 hyperparameters=None, qubit_order=None, quantum_objects=None):
        super().__init__(Q=Q, warm_start=np.asarray(warm_start, dtype=float),
                          offset=offset, p=p, hyperparameters=hyperparameters,
                          qubit_order=qubit_order, quantum_objects=quantum_objects)

    def _set_circuit(self):
        self.circuit = QAOACircuit(
            Q=self.Q, offset=self.offset, p=self.p,
            hyperparameters=self.hyperparameters, qubit_order=self.qubit_order,
            quantum_objects=self.quantum_objects,
        )
        self.circuit._h = self.circuit._h + self.warm_start   # apply the bias

    def _set_objective_function(self):
        self.objective_function = QAOACostFunction(
            self.circuit._J, self.circuit._h, self.circuit._ising_offset
        )

    def get_circuit(self):
        return self.circuit.qasm


alg = WarmStartQAOA(Q, warm_start=np.array([0.1, -0.1, 0.0]), p=2)
```

Reusing `QAOACircuit` and `QAOACostFunction` as-is (rather than writing new circuit/cost classes) is the cheapest way to prototype a variant — write your own `VariationalCircuit` subclass only once you actually need different gates or a different sampling scheme.

### Writing a custom cost function

If your variant's scoring isn't plain Ising energy, subclass `CostFunction`:

```python
from bqs.utils.vqa_utils import CostFunction

class MyCostFunction(CostFunction):
    def evaluate_samples(self, solutions):
        """
        solutions: list of dicts, {variable_name: measured_value}, or an
        equivalent batch of samples.
        Returns a single float — the mean energy/cost over the batch.
        """
        ...
```

### Summary checklist for a new `BaseQAOA` subclass

1. Pass every subclass-specific value as a keyword to `super().__init__(...)` — there's no attribute-ordering rule to follow.
2. Decide whether `BaseQAOA`'s 1-D hyperparameter validation fits your data. If not (e.g. 2-D per-slice hyperparameters), override `_validate_hyperparameters`.
3. Implement `_set_circuit`, `_set_objective_function`, and `get_circuit`.
4. Reuse existing `VariationalCircuit`/`CostFunction` classes where possible; only write new ones when gates or scoring genuinely differ.
5. If adding `from_ising`, build the Ising-case attribute dict and pass it to `cls._from_attributes(...)`.

---

## BQM utilities

```python
from bqs.utils.bqm import qubo_to_ising, ising_to_qubo, qubo_energy, ising_energy

# QUBO → Ising
J, h, offset = qubo_to_ising(Q)
# J: (n, n) upper-triangular; h: (n,) local fields; offset: float

# Ising → QUBO
Q2, new_offset = ising_to_qubo(J, h, offset)

# Evaluate a single sample
x = np.array([0, 1, 1, 0, 1, 0])          # binary {0,1}
s = np.array([-1, 1, 1, -1, 1, -1])       # spin {-1,+1}
e_qubo  = qubo_energy(Q, 0.0, x)
e_ising = ising_energy(J, h, offset, s)   # equals e_qubo
```

---

## Architecture overview

```
BaseVQA (abstract)
├── BaseQAOA (abstract)
│   ├── QAOA
│   ├── pQAOA
│   │   └── pQAOASingleParameters
│   └── SingleSliceQAOA

VariationalCircuit (abstract)
├── QAOACircuit
├── pQAOACircuit
│   └── pQAOASingleParametersCircuit
└── SingleSliceQAOACircuit  (inherits QAOACircuit)

CostFunction (abstract)
└── QAOACostFunction

Optimizer (abstract)
├── COBYLA
├── SPSA
└── MonteCarlo

AbstractQuantumObjects (abstract)
└── CirqQuantumObjects
```

Each `BaseQAOA` subclass wires together one `VariationalCircuit`, one `QAOACostFunction`, and an `Optimizer`. Users interact only with the `BaseQAOA` layer; the circuit and cost-function objects are set up automatically on construction — see [Building new algorithm classes](#building-new-algorithm-classes) if you're adding a branch to this tree yourself.

---

## Reference

> Michele Cattelan and Sheir Yarkoni,
> *Parallel Circuit Implementation of Variational Quantum Algorithms*,
> npj Quantum Information, 2025.
> [https://www.nature.com/articles/s41534-025-00982-6](https://www.nature.com/articles/s41534-025-00982-6)
