Metadata-Version: 2.2
Name: bqs
Version: 0.2.0
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 solved 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)
```

---

## 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.

---

## 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)
