# Photonic Circuits (Perceval)

Experimental

The photonics package is experimental. Its public interface may change, and it currently evaluates one sample and one scalar observable at a time. The [roadmap](#roadmap) lists what comes next.

eigenkit models **linear-optics (photonic) circuits** in a separate package, with its own model of computation. They are not a backend under the qubit `Provider`, because the two models of computation are different:

| Qubit circuits (`core`)                                 | Photonic circuits (`photonics`)                                              |
| ------------------------------------------------------- | ---------------------------------------------------------------------------- |
| A statevector of size `2ⁿ`.                             | A Fock space, made of the configurations of photon numbers over the *modes*. |
| Gates act on qubits.                                    | Beam splitters and phase shifters act on *modes*.                            |
| The readout is an expectation value of Pauli operators. | The readout is a distribution of photon numbers, or detection statistics.    |
| The parameter-shift rule shifts by `±π/2`.              | A photonic parameter-shift rule uses more terms.                             |

Placing photonics under the qubit `Provider` would work only through dual-rail encoding. That encoding is probabilistic, because its two-qubit gates rely on post-selection, and it gives up everything that makes photonics interesting. So photonics has its own intermediate representation (IR) and its own executor. It reuses the parts of the core that do not depend on the model of computation: the `Parameter` expression tree, the separation between compiling and running, and the bridge to PyTorch through `autograd.Function`.

## Install

```bash
pip install eigenkit[perceval]
```

This installs [Perceval](https://github.com/Quandela/Perceval), the photonic software development kit of Quandela, under the MIT licence. It is not part of the `all` extra, because the Perceval simulator pins `protobuf` to an exact version. Install it explicitly, for example with `eigenkit[all,perceval]`. Importing `eigenkit`, or even `eigenkit.photonics`, never imports Perceval. Perceval is loaded only when you create or run a Perceval executor.

## Quickstart: the Hong-Ou-Mandel dip

A circuit is a `PhotonicCircuitSpec`: a number of modes `n_modes`, an ordered tuple of components, and an input Fock state, which gives the number of photons in each mode. The angles of the components are the **same** `Parameter` expressions as in the qubit intermediate representation, or plain numbers.

```python
import math
import numpy as np

from eigenkit.photonics import BeamSplitter, PercevalBackend, PhotonicCircuitSpec

# Two indistinguishable photons, one in each mode, into a balanced beam-splitter.
spec = PhotonicCircuitSpec(
    n_modes=2,
    n_weights=0,
    n_inputs=0,
    components=(BeamSplitter(modes=(0, 1), theta=math.pi / 2),),  # theta=pi/2 => 50:50
    input_state=(1, 1),
)

compiled = PercevalBackend().compile(spec)  # compile once
dist = compiled.run(np.array([]), np.array([]))  # {fock_tuple: probability}
print(dist)  # {(0, 2): 0.5, (2, 0): 0.5}: the photons bunch, and P(1, 1) = 0
```

The photons *bunch*: `(1, 1)` has probability `0`. This result is not classical, because it comes from the permanent of the transfer matrix of the circuit. It is a good check that the executor computes real boson statistics.

The components available today are `BeamSplitter(modes, theta, phi=0)`, `PhaseShifter(mode, phi)` and `Permutation(perm)`. An angle can be a plain number (stored as a `Constant`), a `WeightRef` (trainable), an `InputRef` (data), or any combination with `Scale`, `Sum` and `Multiply`. This is exactly the `Parameter` tree of the qubit circuits, and it accepts the same values as the qubit `CircuitBuilder`.

## Observables

The readouts are observables that are **diagonal in the Fock basis**, which are the ones that photonic quantum machine learning uses in practice. A `FockProbability` target must hold as many photons as the input state, because linear optics conserves the number of photons.

```python
from eigenkit.photonics import FockProbability, PhotonNumber

w, x = np.array([]), np.array([])  # the circuit above has no weights and no inputs
compiled.expectation(w, x, PhotonNumber(mode=0))  # <n̂_0>, expected photons in mode 0
compiled.expectation(w, x, FockProbability(target=(2, 0)))  # P(detect exactly |2,0>)
```

## Gradients: the photonic parameter-shift rule

The qubit rule with shifts of `±π/2` does **not** apply to linear optics. The generator of a phase shifter is the number operator, whose spectrum is made of integers. So an expectation that is diagonal in the Fock basis is a trigonometric polynomial in the phase, of degree at most `N`, where `N` is the total number of photons. Its exact derivative needs the general equidistant shift rule with `2N` terms. The rule with `±π/2` is the case `N = 1`.

The angle `theta` of a beam splitter is generated by an operator of the special unitary group SU(2) on its two modes. The eigenvalues of this operator are integers or half-integers, depending on how many photons the two modes hold. In a circuit with two modes, this count is always `N`, and the same rule with `2N` terms applies. With more modes, the expectation also has half-integer frequencies in `theta`, so the rule is applied in `theta / 2` with order `2N`, which takes `4N` evaluations per beam-splitter angle. Both rules are exact. They reuse `chain_factor` from the qubit intermediate representation to map the derivatives with respect to the angles onto the weights and the inputs.

```python
from eigenkit.core.ir import WeightRef
from eigenkit.photonics import PhaseShifter

# One photon through a Mach-Zehnder interferometer with a trainable phase.
spec = PhotonicCircuitSpec(
    n_modes=2,
    n_weights=1,
    n_inputs=0,
    components=(
        BeamSplitter(modes=(0, 1), theta=math.pi / 2),
        PhaseShifter(mode=0, phi=WeightRef(0)),
        BeamSplitter(modes=(0, 1), theta=math.pi / 2),
    ),
    input_state=(1, 0),
)
compiled = PercevalBackend().compile(spec)
w = np.array([0.3])
weight_jac, input_jac = compiled.jacobian(w, x, PhotonNumber(0), method="parameter_shift")
```

`method="finite_difference"` gives a gradient by central differences, as a check. With two or more photons, the naive rule with `±π/2` is *wrong*, and only the rule with more terms matches the finite differences (see `tests/photonics/test_gradients.py`). Pass `need_input_grad=False`, or `need_weight_grad=False`, to skip the shift evaluations for a side that you do not need. That side is then returned as `None`. `PhotonicCircuit` does this automatically for the inputs that do not require a gradient.

## Training with PyTorch

`PhotonicCircuit` is an `nn.Module` that holds the weights as an `nn.Parameter` and passes the photonic gradient to PyTorch autograd. So `loss.backward()` and any `torch.optim` optimiser work as usual:

```python
import torch
from eigenkit.photonics import FockProbability
from eigenkit.photonics import PhotonicCircuit

circuit = PhotonicCircuit(spec, FockProbability(target=(0, 1)), init_weights=[1.5])
opt = torch.optim.Adam(circuit.parameters(), lr=0.1)
for _ in range(120):
    opt.zero_grad()
    loss = (1.0 - circuit()) ** 2  # for example, send the photon to mode 1
    loss.backward()
    opt.step()
```

A complete worked classifier is in [`examples/photonics_classification.py`](https://github.com/eigenorbit/eigenkit/blob/main/examples/photonics_classification.py). It encodes 2-D features as phases, trains a shallow interferometer with Adam, and reaches an accuracy of about `0.97`.

Place the encodings of a photonic feature map so that they do not commute

When two features are encoded as phases on *adjacent* single modes of a one-photon state, they only change their *relative* phase. The model then sees `x0 − x1` and cannot see `x0 + x1`. Separate the two encodings with a beam splitter, which does not commute with them, so that both features enter independently. This is the photonic version of ordering the rotations in a qubit feature map.

## Reservoir computing

Photonics also supports **reservoir computing**, in the form of a quantum extreme learning machine. A *fixed* random interferometer is the reservoir. Each input is encoded in its phases and read out as photon-number statistics, which gives a high-dimensional nonlinear map of the features. Only a classical linear head is trained. No gradient through the optics is needed, only the forward path of the executor. It reuses the `Stateless` windowing policy of the qubit reservoir package without changes, and the readout can be the closed-form `LinearReadout` of the same package.

```python
import numpy as np
from eigenkit.photonics.reservoir import random_reservoir_spec, run_photonic_reservoir
from eigenkit.reservoir import LinearReadout, Stateless

rng = np.random.default_rng(0)
sequences = rng.uniform(0, np.pi, size=(8, 10, 1))  # (B, T, F), one feature per step
targets = sequences[..., 0] * np.roll(sequences[..., 0], 1, axis=1)  # x_t * x_{t-1}

spec, weights = random_reservoir_spec(n_modes=4, n_inputs=2, input_state=(1, 1, 0, 0))
feats = run_photonic_reservoir(
    spec, weights, sequences, policy=Stateless(window=2)
)  # (B, T, n_fock)
head = LinearReadout(ridge=1e-6).fit(
    feats.reshape(-1, feats.shape[-1]), targets.reshape(-1)
)
```

A worked example is in [`examples/photonics_reservoir.py`](https://github.com/eigenorbit/eigenkit/blob/main/examples/photonics_reservoir.py). It learns the bilinear memory task `y_t = x_t·x_{t-1}`, which a linear readout on the raw inputs cannot learn.

No memory yet

Only reservoirs without a state, or with a classical window, are supported. A photonic carry-over of the state, which would be a quantum memory across time steps, needs Fock states to be injected into the executor, and it is not supported yet.

## Running on Quandela Cloud

Besides local simulation, a circuit can run on **Quandela Cloud**, on a real photonic quantum processing unit (QPU) or on a cloud simulator, through `QuandelaCloudBackend`. The results are *sampled* (counts of shots), so `run` returns an estimated Fock distribution, and gradients are not available, as with the qubit QPU backends. The backend follows the QPU conventions of the library. It issues a `CostWarning` on the first submission, `submit()` returns a handle with `.result()` and `.status()` without waiting, and `dry_run_info()` works offline.

```python
import numpy as np
from eigenkit.photonics import QuandelaCloudBackend

backend = QuandelaCloudBackend("qpu:ascella", shots=20000)  # token from $QUANDELA_TOKEN
compiled = backend.compile(spec)  # `spec` and `weights` from the reservoir example above
print(compiled.dry_run_info())  # inspect offline, before submitting
handle = compiled.submit(weights, np.zeros(spec.n_inputs))
print(handle.status())
dist = handle.result()  # blocks until the cloud job completes
```

The access token is read from the `token=` argument, then from `$QUANDELA_TOKEN` or `$PERCEVAL_TOKEN`, and then from the token that Perceval itself would use (saved with `pcvl.RemoteConfig`). It is never written in the code. The processor connects once per compiled handle, and the connection is reused for later runs. `run` and `result()` wait at most `poll_timeout_seconds` (`3600` by default). They raise an error if the job fails, or if the post-selection on the number of detected photons (`min_detected_photons`) discards every shot. See [`examples/quandela_cloud_hardware.py`](https://github.com/eigenorbit/eigenkit/blob/main/examples/quandela_cloud_hardware.py), and the live test marked `qpu` in `tests/photonics/test_quandela_cloud.py`, which you run with `pytest -m qpu` and a token set.

## How it maps to the qubit core

The photonics package has the same layers as the core, for a different model of computation:

| Qubit core                                                                                                     | Photonics package                                     |
| -------------------------------------------------------------------------------------------------------------- | ----------------------------------------------------- |
| `core/ir.py`                                                                                                   | `photonics/ir.py`                                     |
| `core/provider.py` (the protocol, and the general parameter-shift rule used when a backend has nothing faster) | the `photonics` executor and `photonics/gradients.py` |
| `core/backends/*_backend.py`                                                                                   | the Perceval executor, which is a photonic backend    |
| `torch/` (the PyTorch interface)                                                                               | `PhotonicCircuit` (the photonic PyTorch interface)    |

The Perceval executor is a *photonic* backend, not a qubit `Provider`. The PyTorch adapter belongs to the photonics package, not to the qubit `torch/` package. The `Parameter` tree, the separation between compiling and running, and the `autograd.Function` bridge are reused from the core without changes.

## Status and limitations

- **Local simulation** uses the strong linear optical simulation (SLOS) algorithm of Perceval. It is exact within the sizes it can handle. Its cost grows with the numbers of photons and modes, and the sizes it can handle are smaller than for a qubit statevector of size `2ⁿ`.
- **Execution on Quandela Cloud**, on a real QPU or a cloud simulator, goes through `QuandelaCloudBackend`. The results are sampled, so there are no gradients.
- **One scalar observable is evaluated at a time, for one sample.** There is no batching yet.
- **Only observables that are diagonal in the Fock basis are supported.**

## Roadmap

These items are planned, in no fixed order:

- A `PhotonicMeasurement` output type will allow several observables per evaluation and readouts that are not diagonal. It will also report the success probability of the post-selection with a sampled distribution.
- `run`, `expectation` and `PhotonicCircuit` will accept a leading batch axis.
- The Perceval circuit will be built once, with symbolic `pcvl.P` parameters that are bound at each call, in place of being rebuilt at every call.
- Photonic reservoirs will get a quantum memory (an analogue of `Carry`), which needs Fock states to be injected between time steps.
- `QuandelaCloudBackend` will move off `RemoteProcessor` and `RemoteJob`, which Perceval 1.3 keeps under `perceval.runtime.legacy` (they are still re-exported as `pcvl.*`).
- **MerLin**, the PyTorch framework of Quandela for quantum machine learning, will be evaluated as an alternative to the parameter-shift bridge, with native autograd.