diff --git a/.github/workflows/lint.yml b/.github/workflows/lint.yml new file mode 100644 index 00000000..7ab03801 --- /dev/null +++ b/.github/workflows/lint.yml @@ -0,0 +1,30 @@ +name: Linting + +on: + push: + branches: + - main + pull_request: {} + workflow_dispatch: {} + +jobs: + lint: + name: Linting + runs-on: ubuntu-latest + strategy: + matrix: + python-version: ["3.10"] + steps: + - name: Checkout main code + uses: actions/checkout@v4 + - name: Set up Python ${{ matrix.python-version }} + uses: actions/setup-python@v5 + with: + python-version: ${{ matrix.python-version }} + - name: Install pre-commit + run: | + pip install pre-commit + pre-commit install + - name: Check files and lint + run: | + pre-commit run --all-files diff --git a/.github/workflows/run-tests-and-mypy.yml b/.github/workflows/test.yml similarity index 62% rename from .github/workflows/run-tests-and-mypy.yml rename to .github/workflows/test.yml index 87d2bd5a..b394b13d 100644 --- a/.github/workflows/run-tests-and-mypy.yml +++ b/.github/workflows/test.yml @@ -1,4 +1,4 @@ -name: Linting / Tests / Documentation +name: Tests on: push: @@ -16,29 +16,6 @@ concurrency: cancel-in-progress: true jobs: - lint: - name: Linting - runs-on: ubuntu-latest - strategy: - matrix: - python-version: ["3.10"] - steps: - - name: Checkout main code and submodules - uses: actions/checkout@v4 - - name: Set up Python ${{ matrix.python-version }} - uses: actions/setup-python@v4 - with: - python-version: ${{ matrix.python-version }} - - name: Install hatch - run: | - python -m pip install hatch - - name: Install main dependencies - run: | - python -m hatch -v -e tests - - name: Lint - run: | - python -m hatch -e tests run pre-commit run --all-files - unit_tests: name: Unit testing runs-on: ubuntu-latest @@ -49,21 +26,19 @@ jobs: - name: Checkout main code and submodules uses: actions/checkout@v4 - name: Set up Python ${{ matrix.python-version }} - uses: actions/setup-python@v4 + uses: actions/setup-python@v5 with: python-version: ${{ matrix.python-version }} - name: Install hatch run: | - python -m pip install hatch - - name: Install main dependencies - run: | - python -m hatch -v -e tests - - name: Lint + pip install --upgrade pip + pip install hatch + - name: Install main dependencies in test env run: | - python -m hatch -e tests run pre-commit run --all-files + hatch -v -e tests - name: Perform unit tests run: | - python -m hatch -e tests run pytest -n auto + hatch -e tests run pytest -n auto test_docs: name: Documentation @@ -75,41 +50,46 @@ jobs: - name: Checkout main code and submodules uses: actions/checkout@v4 - name: Set up Python ${{ matrix.python-version }} - uses: actions/setup-python@v4 + uses: actions/setup-python@v5 with: python-version: ${{ matrix.python-version }} - name: Install hatch run: | - python -m pip install hatch + pip install --upgrade pip + pip install hatch - name: Install main dependencies run: | - python -m hatch -v -e docs + hatch -v -e docs - name: Test docs run: | - python -m hatch run docs:build + hatch run docs:build publish: name: Publish to PyPI if: startsWith(github.ref, 'refs/tags/v') needs: unit_tests runs-on: ubuntu-latest + permissions: + # IMPORTANT: this permission is mandatory for trusted publishing + id-token: write steps: - - name: Checkout main code and submodules + - name: Checkout main code uses: actions/checkout@v4 with: ref: ${{ github.ref }} - name: Set up Python - uses: actions/setup-python@v4 + uses: actions/setup-python@v5 with: python-version: "3.10" - name: Install Python dependencies run: | - python -m pip install --upgrade pip + pip install --upgrade pip pip install hatch - - name: Build and publish package + - name: Build package run: | hatch build - hatch publish -u __token__ -a ${{ secrets.PYPI_API_TOKEN }} + - name: Publish to PyPI + uses: pypa/gh-action-pypi-publish@release/v1 - name: Confirm deployment timeout-minutes: 5 run: | @@ -123,17 +103,28 @@ jobs: deploy_docs: name: Deploy documentation if: startsWith(github.ref, 'refs/tags/v') - needs: unit_tests + needs: [unit_tests, test_docs] runs-on: ubuntu-latest steps: - - name: Checkout main code and submodules + - name: Checkout main code uses: actions/checkout@v4 + - name: Install JetBrains Mono font + run: | + sudo apt install -y wget unzip fontconfig + wget https://download.jetbrains.com/fonts/JetBrainsMono-2.304.zip + unzip JetBrainsMono-2.304.zip -d JetBrainsMono + mkdir -p /usr/share/fonts/truetype/jetbrains + cp JetBrainsMono/fonts/ttf/*.ttf /usr/share/fonts/truetype/jetbrains/ + fc-cache -f -v + - name: Install graphviz + run: sudo apt-get install -y graphviz - name: Set up Python 3.10 - uses: actions/setup-python@v4 + uses: actions/setup-python@v5 with: python-version: '3.10' - name: Install Hatch run: | + pip install --upgrade pip pip install hatch - name: Deploy docs run: | diff --git a/README.md b/README.md index 0d7961d8..9afde80e 100644 --- a/README.md +++ b/README.md @@ -4,7 +4,7 @@ It acts as the main backend for [`Qadence`](https://github.com/pasqal-io/qadence), a digital-analog quantum programming interface. `pyqtorch` allows for writing fully differentiable quantum programs using both digital and analog operations; enabled via a intuitive, torch-based syntax. -[![Linting / Tests/ Documentation](https://github.com/pasqal-io/pyqtorch/actions/workflows/run-tests-and-mypy.yml/badge.svg)](https://github.com/pasqal-io/pyqtorch/actions/workflows/run-tests-and-mypy.yml) +[![Linting / Tests/ Documentation](https://github.com/pasqal-io/pyqtorch/actions/workflows/test.yml/badge.svg)](https://github.com/pasqal-io/pyqtorch/actions/workflows/test.yml) [![License](https://img.shields.io/badge/License-Apache_2.0-blue.svg)](https://opensource.org/licenses/Apache-2.0) [![Pypi](https://badge.fury.io/py/pyqtorch.svg)](https://pypi.org/project/pyqtorch/) diff --git a/docs/analog.md b/docs/analog.md index d9e4cf25..0ab89c72 100644 --- a/docs/analog.md +++ b/docs/analog.md @@ -2,6 +2,9 @@ An analog operation is one whose unitary is best described by the evolution of some hermitian generator, or Hamiltonian, acting on an arbitrary number of qubits. For a time-independent generator $\mathcal{H}$ and some time variable $t$, the evolution operator is $\exp(-i\mathcal{H}t)$. `pyqtorch` provides the HamiltonianEvolution class to initialize analog operations. There exists several ways to pass a generator, and we present them next. +!!! warning "Dimensionless units" + The quantity $\mathcal{H}t$ has to be **dimensionless** for exponentiation in PyQTorch. + ### Tensor generator The first case of generator we can provide is simply an arbitrary hermitian tensor. diff --git a/docs/differentiation.md b/docs/differentiation.md index bb0f23e3..359e8416 100644 --- a/docs/differentiation.md +++ b/docs/differentiation.md @@ -13,12 +13,11 @@ The [adjoint differentiation mode](https://arxiv.org/abs/2009.02823) computes fi The Generalized parameter shift rule (GPSR mode) is an extension of the well known [parameter shift rule (PSR)](https://arxiv.org/abs/1811.11184) algorithm [to arbitrary quantum operations](https://arxiv.org/abs/2108.01218). Indeed, PSR only works for quantum operations whose generator has a single gap in its eigenvalue spectrum, GPSR extending to multi-gap. !!! warning "Usage restrictions" - At the moment, only operations with two distinct eigenvalues - from their generator (single gap shift rule) are supported. The multi gap case - will be supported in a later release. - Circuits with one or more Scale or HamiltonianEvolution operations are not supported. - Finally, circuits with operations sharing a same parameter name - are also not supported. + At the moment, circuits with one or more Scale or HamiltonianEvolution operations are not supported. + They should be handled differently as GPSR requires operations to be of the form presented below. + Also, circuits with operations sharing a same parameter name are also not supported + as such cases are handled by our other Python package for differentiable digital-analog quantum programs Qadence + which uses pyqtorch as a backend. Qadence convert circuits to use different parameter names when applying GPSR. For this, we define the differentiable function as quantum expectation value @@ -56,11 +55,11 @@ batch_size = 1 ry = pyq.RY(0, param_name="x") cnot = pyq.CNOT(1, 2) -ops = [ry] +ops = [ry, cnot] n_qubits = 3 circ = pyq.QuantumCircuit(n_qubits, ops) -obs = pyq.QuantumCircuit(n_qubits, [pyq.Z(0)]) +obs = pyq.Observable(pyq.Z(0)) state = pyq.zero_state(n_qubits) values_ad = {"x": torch.tensor([torch.pi / 2], requires_grad=True)} diff --git a/docs/dropout.md b/docs/dropout.md index 3ba5b13c..1b2f4210 100644 --- a/docs/dropout.md +++ b/docs/dropout.md @@ -1 +1,401 @@ -Coming soon. +## Fitting a noisy sinusoid with quantum dropout + +Here we will demonstrate an implemention [quantum dropout](https://arxiv.org/abs/2310.04120), for the case of fitting a noisy sine function. To show the usefulness of dropout for quantum neural networks (QNNs), we shall compare the performance of a QNN with dropout and one without. + +Firstly, we define the dataset that we will perform regression on, this function is $sin(\pi x)+\epsilon$, where $x\in\reals$ and $\epsilon$ is noise sampled from a normal distribution which is then added to each point. + +```python exec="on" source="material-block" session="dropout" +from __future__ import annotations + +import matplotlib.pyplot as plt +from matplotlib import ticker +import numpy as np +import torch +from torch import manual_seed, optim, tensor + +import pyqtorch as pyq +from pyqtorch.circuit import DropoutQuantumCircuit +from pyqtorch.primitives import Parametric +from pyqtorch.utils import DropoutMode + +seed = 70 +manual_seed(seed) +np.random.seed(seed) + +# choose device and hyperparameters + +device = torch.device("cpu") +n_qubits = 2 # a greater performance difference is observed with 5 or more qubits +depth = 5 # a greater performance difference is observed at depth 10 +lr = 0.01 +n_points = 20 +epochs = 100 # a greater performance difference is observed at 200-250 epochs of training +dropout_prob = 0.06 +noise = 0.4 + +def sin_dataset(dataset_size: int = 100, test_size: float = 0.4, noise: float = 0.4): + """Generates points (x,y) which follow sin(πx)+ϵ, + where epsilon is noise randomly sampled from the normal + distribution for each datapoint. + Args: + dataset_size (int): total number of points. Defaults to 100. + test_size (float): fraction of points for testing. Defaults to 0.4. + noise (float): standard deviation of added noise. Defaults to 0.4. + Returns: + data (tuple): data divided into training and testing + """ + x_ax = np.linspace(-1, 1, dataset_size) + y = np.sin(x_ax * np.pi) + noise = np.random.normal(0, 0.5, y.shape) * noise + y += noise + + # permuting the points around before dividing into train and test sets. + rng = np.random.default_rng(seed) + indices = rng.permutation(dataset_size) + n_test = int(dataset_size * test_size) + n_train = int(dataset_size * (1 - test_size)) + test_indices = indices[:n_test] + train_indices = indices[n_test : (n_test + n_train)] + x_train, x_test, y_train, y_test = ( + x_ax[train_indices], + x_ax[test_indices], + y[train_indices], + y[test_indices], + ) + + x_train = x_train.reshape(-1, 1) + x_test = x_test.reshape(-1, 1) + y_train = y_train.reshape(-1, 1) + y_test = y_test.reshape(-1, 1) + return x_train, x_test, y_train, y_test + + +# generates points following sin(πx)+ϵ, split into training and testing sets. +x_train, x_test, y_train, y_test = sin_dataset(dataset_size=n_points, test_size=0.25, noise=0.4) +``` + +We can now visualise the function we will train QNNs to learn. + +```python exec="on" source="material-block" html="1" session="dropout" +fig, ax = plt.subplots() +plt.plot(x_train, y_train, "o", label="training") +plt.plot(x_test, y_test, "o", label="testing") +plt.plot( + np.linspace(-1, 1, 100), + [np.sin(x * np.pi) for x in np.linspace(-1, 1, 100)], + linestyle="dotted", + label=r"$\sin(x)$", +) +plt.ylabel(r"$y = \sin(\pi\cdot x) + \epsilon$") +plt.xlabel(r"$x$") +ax.xaxis.set_major_locator(ticker.MultipleLocator(0.5)) +ax.yaxis.set_major_locator(ticker.MultipleLocator(0.5)) +plt.legend() +from io import StringIO # markdown-exec: hide +from matplotlib.figure import Figure # markdown-exec: hide +def fig_to_html(fig: Figure) -> str: # markdown-exec: hide + buffer = StringIO() # markdown-exec: hide + fig.savefig(buffer, format="svg") # markdown-exec: hide + return buffer.getvalue() # markdown-exec: hide +print(fig_to_html(plt.gcf())) # markdown-exec: hide +``` + +Since our QNN can only output values between -1 and 1 we need to scale the data to be between these values. + +```python exec="on" source="material-block" session="dropout" +class MinMaxScaler: + """A class which scales data to be within a chosen range""" + + def __init__(self, feature_range: tuple =(0, 1)): + self.feature_range = feature_range + self.min = None + self.scale = None + + def fit(self, X): + self.min = X.min(axis=0) + self.scale = X.max(axis=0) - self.min + self.scale[self.scale == 0] = 1 # Avoid division by zero for constant features + + def transform(self, X): + X_scaled = (X - self.min) / self.scale + X_scaled = ( + X_scaled * (self.feature_range[1] - self.feature_range[0]) + self.feature_range[0] + ) + return X_scaled + + def fit_transform(self, X): + self.fit(X) + return self.transform(X) + + def inverse_transform(self, X_scaled): + X = (X_scaled - self.feature_range[0]) / (self.feature_range[1] - self.feature_range[0]) + X = X * self.scale + self.min + return X + +scaler = MinMaxScaler(feature_range=(-1, 1)) +y_train = scaler.fit_transform(y_train) +y_test = scaler.transform(y_test) +``` + +Now we need to construct our QNNs, they are comprised of feature maps and an ansatz. The ansatz contains CNOT gates with nearest neighbour entanglement after every rotation gate. The feature map contains two rotation gates RY and RZ, taking as rotation angles $\arcsin(x)$ and $\arccos(x^2)$ respectively. + +```python exec="on" source="material-block" session="dropout" + +def hea_ansatz(n_qubits, layer): + """creates an ansatz which performs RX, RZ, RX rotations on each qubit, + which nearest neighbour CNOT gates interpersed between each rotational gate.""" + ops = [] + for i in range(n_qubits): + ops.append(pyq.RX(i, param_name=f"theta_{i}{layer}{0}")) + + for j in range(n_qubits - 1): + ops.append(pyq.CNOT(control=j, target=j + 1)) + + for i in range(n_qubits): + ops.append(pyq.RZ(i, param_name=f"theta_{i}{layer}{1}")) + + for j in range(n_qubits - 1): + ops.append(pyq.CNOT(control=j, target=j + 1)) + + for i in range(n_qubits): + ops.append(pyq.RX(i, param_name=f"theta_{i}{layer}{2}")) + + for j in range(n_qubits - 1): + ops.append(pyq.CNOT(control=j, target=j + 1)) + return ops + +# the two feature maps we will be using +def fm1(n_qubits): + return [pyq.RY(i, "x1") for i in range(n_qubits)] + + +def fm2(n_qubits): + return [pyq.RZ(i, "x2") for i in range(n_qubits)] + +# The class which constructs QNNs +class QuantumModelBase(torch.nn.Module): + def __init__(self, n_qubits, n_layers, device): + super().__init__() + self.n_qubits = n_qubits + self.n_layers = n_layers + self.device = device + + self.embedding1 = fm1(n_qubits=n_qubits) + self.embedding2 = fm2(n_qubits=n_qubits) + + self.params = torch.nn.ParameterDict() + operations = self.build_operations() + self.circuit = self.build_circuit(operations) + self.observable = pyq.Observable([pyq.Z(i) for i in range(n_qubits)]).to( + device, dtype=torch.complex64 + ) + self.params = self.params.to(device=device, dtype=torch.float32) + + def build_operations(self): + """defines operations for the quantum circuit and trainable parameters.""" + operations = [] + for i in range(self.n_layers): + operations += self.embedding1 + self.embedding2 + layer_i_ansatz = hea_ansatz(n_qubits=n_qubits, layer=i) + operations += layer_i_ansatz + for op in layer_i_ansatz: + if isinstance(op, Parametric): + self.params[f"{op.param_name}"] = torch.randn(1, requires_grad=True) + + return operations + + def build_circuit(self, operations): + """constructs a QuantumCircuit object and puts it on device.""" + return pyq.QuantumCircuit( + n_qubits=self.n_qubits, + operations=operations, + ).to(device=self.device, dtype=torch.complex64) + + def forward(self, x): + """the forward pass for the QNN""" + x = x.flatten() + x_1 = {"x1": torch.asin(x)} + x_2 = {"x2": torch.acos(x**2)} + state = self.circuit.init_state(batch_size=int(x.shape[0])) + + out = pyq.expectation( + circuit=self.circuit, + state=state, + values={**self.params, **x_1, **x_2}, + observable=self.observable, + ) + + return out + +# first we will define a QNN which is overparameterised with no dropout. +model = QuantumModelBase(n_qubits=n_qubits, n_layers=depth, device=device) +# define the corresponding optimizer for the problem +opt = optim.Adam(model.parameters(), lr=lr) +``` + +We now wish to train the QNN to learn the noisy sinusoidal function we have defined. This function will return the loss curves for the train and test sets. + +```python exec="on" source="material-block" session="dropout" +def train(model, opt, x_train, y_train, x_test, y_test, epochs, device): + # lists which will store losses for train and tests sets as we train. + train_loss_history = [] + validation_loss_history = [] + + x_test = tensor(x_test).to(device, dtype=torch.float32) + y_test = ( + tensor(y_test) + .to(device, dtype=torch.float32) + .reshape( + -1, + ) + ) + + x_train = tensor(x_train).to(device, dtype=torch.float32) + y_train = ( + tensor(y_train) + .to(device, dtype=torch.float32) + .reshape( + -1, + ) + ) + + # we will be using the mean squared error as our loss function. + cost_fn = torch.nn.MSELoss() + + for epoch in range(epochs): + model.train() + + opt.zero_grad() + y_preds = model(x_train) + train_loss = cost_fn(y_preds, y_train.flatten()) + train_loss.backward() + opt.step() + + # no dropout is performed during evaluation of the model. + model.eval() + train_preds = model(x_train) + train_loss = cost_fn(train_preds, y_train.flatten()).detach().cpu().numpy() + + test_preds = model(x_test) + test_loss = cost_fn(test_preds, y_test.flatten()).detach().cpu().numpy() + + train_loss_history.append(train_loss) + validation_loss_history.append(test_loss) + + # log performance every 100 epochs. + if epoch % 100 == 0: + print(f"epoch: {epoch}, train loss {train_loss}, val loss: {test_loss}") + + return train_loss_history, validation_loss_history + +# train the vanilla QNN, extracting the training and testing loss curves +no_dropout_train_loss_hist, no_dropout_test_loss_hist = train( + model=model, + opt=opt, + x_train=x_train, + y_train=y_train, + x_test=x_test, + y_test=y_test, + epochs=epochs, + device=device, +) +``` + +Next we define a QNN with the same archecture as before but now includes rotational dropout. Rotional dropout will randomly drop single qubit trainable parameterised gates with some probability dropout_prob. + +```python exec="on" source="material-block" session="dropout" + +class DropoutModel(QuantumModelBase): + """Inherits from QuantumModelBase but the build_circuit function now creates a circuit which drops certain gates with some probability during training.""" + def __init__( + self, n_qubits, n_layers, device, dropout_mode="rotational_dropout", dropout_prob=0.03 + ): + self.dropout_mode = dropout_mode + self.dropout_prob = dropout_prob + super().__init__(n_qubits, n_layers, device) + + def build_circuit(self, operations): + return DropoutQuantumCircuit( + n_qubits=self.n_qubits, + operations=operations, + dropout_mode=self.dropout_mode, + dropout_prob=self.dropout_prob, + ).to(device=self.device, dtype=torch.complex64) + +# Define the QNN with rotational quantum dropout. +model = DropoutModel( + n_qubits=n_qubits, + n_layers=depth, + device=device, + dropout_mode=DropoutMode.ROTATIONAL, + dropout_prob=dropout_prob, +) + +# Define the corresponding optimiser +opt = optim.Adam(model.parameters(), lr=lr) + +# train the QNN which contains rotational dropout. +dropout_train_loss_hist, dropout_test_loss_hist = train( + model=model, + opt=opt, + x_train=x_train, + y_train=y_train, + x_test=x_test, + y_test=y_test, + epochs=epochs, + device=device, +) +``` + +Now that we have trained both a regular QNN and one with rotational quantum dropout, we can visualise the training and testing loss curves. What we observe is the regular QNN performing better on the train set but poorer on the test set. We can attribute this discrepency to the regular QNN overfitting to the noise rather than the true underlying function. The QNN with rotational dropout does no exhibit overfitting and adheres to the true function better, thus performs better on the test set. + +```python exec="on" source="material-block" html="1" session="dropout" +def plot_results( + no_dropout_train_history, + no_dropout_test_history, + dropout_train_history, + dropout_test_history, + epochs, +): + fig, ax = plt.subplots(1, 2) + plt.subplots_adjust(wspace=0.05) + ax[0].set_title("MSE train") + + ax[0].plot(range(epochs), no_dropout_train_history, label="no dropout") + ax[0].plot(range(epochs), dropout_train_history, label="dropout") + + ax[1].set_title("MSE test") + ax[1].plot(range(epochs), no_dropout_test_history, label="no dropout") + ax[1].plot(range(epochs), dropout_test_history, label="dropout") + + ax[0].legend( + loc="upper center", bbox_to_anchor=(1.01, 1.25), ncol=4, fancybox=True, shadow=True + ) + + for ax in ax.flat: + ax.set_xlabel("Epochs") + ax.set_ylabel("MSE") + ax.set_yscale("log") + ax.set_ylim([1e-3, 0.6]) + ax.label_outer() + + plt.subplots_adjust(bottom=0.3) + + from io import StringIO # markdown-exec: hide + from matplotlib.figure import Figure # markdown-exec: hide + def fig_to_html(fig: Figure) -> str: # markdown-exec: hide + buffer = StringIO() # markdown-exec: hide + fig.savefig(buffer, format="svg") # markdown-exec: hide + return buffer.getvalue() # markdown-exec: hide + print(fig_to_html(plt.gcf())) # markdown-exec: hide + + +# finally we compare the +plot_results( + no_dropout_train_history=no_dropout_train_loss_hist, + dropout_train_history=dropout_train_loss_hist, + no_dropout_test_history=no_dropout_test_loss_hist, + dropout_test_history=dropout_test_loss_hist, + epochs=epochs, +) +``` diff --git a/docs/embed.md b/docs/embed.md index bcbb4694..4479f6b2 100644 --- a/docs/embed.md +++ b/docs/embed.md @@ -30,7 +30,8 @@ circ = pyq.QuantumCircuit(1, [pyq.RX(0, mul_sinx_y)]) state= pyq.zero_state(1) y = torch.rand(1, requires_grad=True) values = {'x': x, 'y': y} -expval = pyq.expectation(circuit=circ, state=state, values=values, observable= pyq.Observable(1, [pyq.Z(0)]),diff_mode=pyq.DiffMode.AD,embedding=embedding) +obs = pyq.Observable([pyq.Z(0)]) +expval = pyq.expectation(circuit=circ, state=state, values=values, observable=obs, diff_mode=pyq.DiffMode.AD, embedding=embedding) print(torch.autograd.grad(expval, (x, y), torch.ones_like(expval))) ``` diff --git a/docs/fitting_a_function.md b/docs/fitting_a_function.md index 4cf708e2..42704782 100644 --- a/docs/fitting_a_function.md +++ b/docs/fitting_a_function.md @@ -7,9 +7,9 @@ from operator import add from functools import reduce import torch import pyqtorch as pyq -from pyqtorch.circuit import hea +from pyqtorch.composite import hea from pyqtorch.utils import DiffMode -from pyqtorch.parametric import Parametric +from pyqtorch.primitives import Parametric import matplotlib.pyplot as plt from torch.nn.functional import mse_loss @@ -35,7 +35,7 @@ feature_map = [pyq.RX(i, f'x') for i in range(N_QUBITS)] ansatz, params = hea(N_QUBITS, DEPTH, 'theta') # Lets move all necessary components to the DEVICE circ = pyq.QuantumCircuit(N_QUBITS, feature_map + ansatz).to(device=DEVICE, dtype=COMPLEX_DTYPE) -observable = pyq.DiagonalObservable(N_QUBITS, pyq.Z(0)).to(device=DEVICE, dtype=COMPLEX_DTYPE) +observable = pyq.Observable(pyq.Z(0)).to(device=DEVICE, dtype=COMPLEX_DTYPE) params = params.to(device=DEVICE, dtype=REAL_DTYPE) x, y = x.to(device=DEVICE, dtype=REAL_DTYPE), y.to(device=DEVICE, dtype=REAL_DTYPE) state = circ.init_state() diff --git a/docs/index.md b/docs/index.md index 93b40513..37994dd8 100644 --- a/docs/index.md +++ b/docs/index.md @@ -93,6 +93,9 @@ psi_end = hamiltonian_evolution( assert is_normalized(psi_end, atol=1e-05) ``` +!!! warning "Dimensionless units" + The quantity $\mathcal{H}t$ has to be considered **dimensionless** for exponentiation in `pyqtorch`. + ## Circuits Using digital and analog operations, you can can build fully differentiable quantum circuits using the `QuantumCircuit` class; note that the default differentiation mode in pyqtorch is using torch.autograd. @@ -109,7 +112,7 @@ n_qubits = 2 circ = pyq.QuantumCircuit(n_qubits, ops) state = pyq.random_state(n_qubits) theta = torch.rand(1, requires_grad=True) -obs = pyq.QuantumCircuit(n_qubits, [pyq.Z(0)]) +obs = pyq.Observable(pyq.Z(0)) expval = pyq.expectation(circ, state, {"theta": theta}, obs) dfdtheta = torch.autograd.grad(expval, theta, torch.ones_like(expval)) ``` diff --git a/docs/noise.md b/docs/noise.md index 631c5f37..1f4aa84a 100644 --- a/docs/noise.md +++ b/docs/noise.md @@ -116,7 +116,7 @@ import torch from pyqtorch.circuit import QuantumCircuit from pyqtorch.noise import BitFlip -from pyqtorch.primitive import X +from pyqtorch.primitives import X from pyqtorch.utils import product_state diff --git a/docs/pde.md b/docs/pde.md index bced8693..e87ba627 100644 --- a/docs/pde.md +++ b/docs/pde.md @@ -13,9 +13,9 @@ import numpy as np import torch from torch import Tensor, exp, linspace, ones_like, optim, rand, sin, tensor from torch.autograd import grad -from pyqtorch.circuit import hea -from pyqtorch import CNOT, RX, RY, QuantumCircuit, Z, expectation, DiagonalObservable, Sequence, Merge, Add -from pyqtorch.parametric import Parametric +from pyqtorch.composite import hea +from pyqtorch import CNOT, RX, RY, QuantumCircuit, Z, expectation, Sequence, Merge, Add, Observable +from pyqtorch.primitives import Parametric from pyqtorch.utils import DiffMode DIFF_MODE = DiffMode.AD @@ -90,7 +90,7 @@ feature_map = [RX(i, VARIABLES[X_POS]) for i in range(N_QUBITS // 2)] + [ ] ansatz, params = hea(N_QUBITS, DEPTH, "theta") circ = QuantumCircuit(N_QUBITS, feature_map + ansatz).to(device=DEVICE, dtype=COMPLEX_DTYPE) -sumZ_obs = DiagonalObservable(N_QUBITS, Add(Sequence([Z(i) for i in range(N_QUBITS)]))).to(device=DEVICE, dtype=COMPLEX_DTYPE) +sumZ_obs = Observable([Z(i) for i in range(N_QUBITS)]).to(device=DEVICE, dtype=COMPLEX_DTYPE) params = params.to(device=DEVICE, dtype=REAL_DTYPE) state = circ.init_state() diff --git a/pyproject.toml b/pyproject.toml index 43bbc31e..e7c42857 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -5,6 +5,11 @@ build-backend = "hatchling.build" [project] name = "pyqtorch" description = "An efficient, large-scale emulator designed for quantum machine learning, seamlessly integrated with a PyTorch backend. Please refer to https://pyqtorch.readthedocs.io/en/latest/ for setup and usage info, along with the full documentation." +readme = "README.md" +version = "1.4.1" +requires-python = ">=3.8,<3.13" +license = { text = "Apache 2.0" } +keywords = ["quantum"] authors = [ { name = "Slimane Thabet", email = "slimane.thabet@pasqal.com" }, { name = "Aleksander Wennersteen", email = "aleksander.wennersteen@pasqal.com" }, @@ -15,11 +20,9 @@ authors = [ { name = "Joao P. Moutinho", email = "joao.moutinho@pasqal.com"}, { name = "Vytautas Abramavicius", email = "vytautas.abramavicius@pasqal.com" }, { name = "Anton Quelle", email = "anton.quelle@pasqal.com" }, - { name = "Charles Moussa", email = "charles.moussa@pasqal.com" } + { name = "Charles Moussa", email = "charles.moussa@pasqal.com" }, + { name = "Callum Duffy", email = "callum.duffy@pasqal.com" }, ] -requires-python = ">=3.8,<3.13" -license = {text = "Apache 2.0"} -version = "1.3.1" classifiers=[ "License :: OSI Approved :: Apache Software License", "Programming Language :: Python", @@ -32,19 +35,68 @@ classifiers=[ "Programming Language :: Python :: Implementation :: CPython", "Programming Language :: Python :: Implementation :: PyPy", ] -dependencies = ["torch", "numpy"] +# always specify a version for each package +# to maintain consistency +dependencies = [ + "torch", + "numpy" +] + +[tool.hatch.metadata] +allow-direct-references = true +allow-ambiguous-features = true + +# add one or more extras in the dependencies [project.optional-dependencies] -dev = ["flaky","black", "pytest", "pytest-xdist", "pytest-cov", "flake8", "mypy", "pre-commit", "ruff", "matplotlib"] -testlibs = [ "jax", "qutip~=4.7.5", "jaxlib"] +dev = [ + "flaky", + "black", + "pytest", + "pytest-xdist", + "pytest-cov", + "flake8", + "mypy", + "pre-commit", + "ruff", + "matplotlib" +] +testlibs = [ + "jax", + "qutip~=4.7.5", + "jaxlib" +] + +[project.urls] +Documentation = "https://pasqal-io.github.io/pyqtorch/latest/" +Issues = "https://github.com/pasqal-io/pyqtorch/-/issues" +Source = "https://github.com/pasqal-io/pyqtorch" + [tool.hatch.envs.tests] -features=["testlibs", "dev"] +features = [ + "testlibs", + "dev" +] [tool.hatch.envs.tests.scripts] test = "pytest -n auto {args}" test-docs = "hatch -e docs run mkdocs build --clean --strict" test-cov = "pytest -n auto --cov=pyqtorch tests/" +[tool.pytest.ini_options] +testpaths = ["tests"] +addopts = """-vvv --cov-report=term-missing --cov-config=pyproject.toml --cov=template_python --cov=tests""" +xfail_strict = true +filterwarnings = [ + "ignore:Call to deprecated create function FieldDescriptor", + "ignore:Call to deprecated create function Descriptor", + "ignore:Call to deprecated create function EnumDescriptor", + "ignore:Call to deprecated create function EnumValueDescriptor", + "ignore:Call to deprecated create function FileDescriptor", + "ignore:Call to deprecated create function OneofDescriptor", + "ignore:distutils Version classes are deprecated." +] + [tool.hatch.envs.docs] dependencies = [ "mkdocs", diff --git a/pyqtorch/__init__.py b/pyqtorch/__init__.py index 146cf700..5c47de2f 100644 --- a/pyqtorch/__init__.py +++ b/pyqtorch/__init__.py @@ -46,17 +46,17 @@ logger.info(f"PyQTorch logger successfully setup with log level {LOG_LEVEL}") -from .analog import ( +from .api import expectation, run, sample +from .apply import apply_operator +from .circuit import DropoutQuantumCircuit, QuantumCircuit +from .composite import ( Add, - DiagonalObservable, - HamiltonianEvolution, - Observable, + Merge, Scale, + Sequence, ) -from .api import expectation, run, sample -from .apply import apply_operator -from .circuit import Merge, QuantumCircuit, Sequence from .embed import ConcretizedCallable, Embedding +from .hamiltonians import HamiltonianEvolution, Observable from .noise import ( AmplitudeDamping, BitFlip, @@ -66,12 +66,27 @@ PhaseDamping, PhaseFlip, ) -from .parametric import CPHASE, CRX, CRY, CRZ, PHASE, RX, RY, RZ, U -from .primitive import ( +from .primitives import ( CNOT, + CPHASE, + CRX, + CRY, + CRZ, CSWAP, CY, CZ, + OPS_1Q, + OPS_2Q, + OPS_3Q, + OPS_DIGITAL, + OPS_PARAM, + OPS_PARAM_1Q, + OPS_PARAM_2Q, + OPS_PAULI, + PHASE, + RX, + RY, + RZ, SWAP, H, I, @@ -81,6 +96,7 @@ SDagger, T, Toffoli, + U, X, Y, Z, @@ -111,6 +127,7 @@ "Scale", "Merge", "QuantumCircuit", + "DropoutQuantumCircuit", "Sequence", "CPHASE", "CRX", @@ -148,7 +165,6 @@ "Z", "apply_operator", "Observable", - "DiagonalObservable", "sesolve", "mesolve", ] diff --git a/pyqtorch/analog.py b/pyqtorch/analog.py deleted file mode 100644 index 10e8613c..00000000 --- a/pyqtorch/analog.py +++ /dev/null @@ -1,768 +0,0 @@ -from __future__ import annotations - -import logging -from collections import OrderedDict -from functools import reduce -from logging import getLogger -from operator import add -from typing import Any, Callable, Tuple, Union - -import torch -from torch import Tensor -from torch.nn import Module, ModuleList, ParameterDict - -from pyqtorch.apply import apply_operator -from pyqtorch.circuit import Sequence -from pyqtorch.embed import Embedding -from pyqtorch.matrices import _dagger -from pyqtorch.primitive import Primitive -from pyqtorch.time_dependent.sesolve import sesolve -from pyqtorch.utils import ( - ATOL, - Operator, - SolverType, - State, - StrEnum, - inner_prod, - is_diag, - operator_to_sparse_diagonal, -) - -BATCH_DIM = 2 -TGenerator = Union[Tensor, str, Primitive, Sequence] - -logger = getLogger(__name__) - - -def forward_hook(*args, **kwargs) -> None: # type: ignore[no-untyped-def] - logger.debug("Forward complete") - torch.cuda.nvtx.range_pop() - - -def pre_forward_hook(*args, **kwargs) -> None: # type: ignore[no-untyped-def] - logger.debug("Executing forward") - torch.cuda.nvtx.range_push("HamiltonianEvolution.forward") - - -def backward_hook(*args, **kwargs) -> None: # type: ignore[no-untyped-def] - logger.debug("Backward complete") - torch.cuda.nvtx.range_pop() - - -def pre_backward_hook(*args, **kwargs) -> None: # type: ignore[no-untyped-def] - logger.debug("Executed backward") - torch.cuda.nvtx.range_push("Hamiltonian Evolution.backward") - - -class GeneratorType(StrEnum): - """ - Options for types of generators allowed in HamiltonianEvolution. - """ - - PARAMETRIC_OPERATION = "parametric_operation" - """Generators of type Primitive or Sequence which contain - possibly trainable or non-trainable parameters.""" - OPERATION = "operation" - """Generators of type Primitive or Sequence which do not contain parameters or contain - constants as Parameters for example pyq.Scale(Z(0), torch.tensor([1.])).""" - TENSOR = "tensor" - """Generators of type torch.Tensor in which case a qubit_support needs to be passed.""" - SYMBOL = "symbol" - """Generators which are symbolic, i.e. will be passed via the 'values' dict by the user.""" - TIME_DEPENDENT = "time_dependent" - """Generators which are time dependent.""" - - -class Scale(Sequence): - """ - Generic container for multiplying a 'Primitive' or 'Sequence' instance by a parameter. - - Attributes: - operations: Operations making the Sequence. - param_name: Name of the parameter to multiply operations with. - """ - - def __init__(self, operations: Union[Sequence, Module], param_name: str | Tensor): - """ - Initializes a Scale object. - - Arguments: - operations: Operations making the Sequence. - param_name: Name of the parameter to multiply operations with. - """ - super().__init__( - operations.operations if isinstance(operations, Sequence) else [operations] - ) - self.param_name = param_name - assert len(self.operations) == 1 - - def forward( - self, - state: Tensor, - values: dict[str, Tensor] | ParameterDict = dict(), - embedding: Embedding | None = None, - ) -> State: - """ - Apply the operation(s) multiplying by the parameter value. - - Arguments: - state: Input state. - values: Parameter value. - - Returns: - The transformed state. - """ - return ( - values[self.param_name] * super().forward(state, values) - if isinstance(self.operations, Sequence) - else self._forward(state, values) - ) - - def _forward( - self, - state: Tensor, - values: dict[str, Tensor], - embedding: Embedding | None = None, - ) -> State: - """ - Apply the single operation of Scale multiplied by the parameter value. - - Arguments: - state: Input state. - values: Parameter value. - - Returns: - The transformed state. - """ - return apply_operator( - state, self.unitary(values), self.operations[0].qubit_support - ) - - def unitary( - self, values: dict[str, Tensor], embedding: Embedding | None = None - ) -> Operator: - """ - Get the corresponding unitary. - - Arguments: - values: Parameter value. - - Returns: - The unitary representation. - """ - thetas = ( - values[self.param_name] - if isinstance(self.param_name, str) - else self.param_name - ) - return thetas * self.operations[0].unitary(values, embedding) - - def dagger( - self, values: dict[str, Tensor], embedding: Embedding | None = None - ) -> Operator: - """ - Get the corresponding unitary of the dagger. - - Arguments: - values: Parameter value. - - Returns: - The unitary representation of the dagger. - """ - return _dagger(self.unitary(values, embedding)) - - def jacobian( - self, values: dict[str, Tensor], embedding: Embedding | None = None - ) -> Operator: - """ - Get the corresponding unitary of the jacobian. - - Arguments: - values: Parameter value. - - Returns: - The unitary representation of the jacobian. - - Raises: - NotImplementedError - """ - raise NotImplementedError( - "The Jacobian of `Scale` is done via decomposing it into the gradient w.r.t\ - the scale parameter and the gradient w.r.t to the scaled block." - ) - # TODO make scale a primitive block with an additional parameter - # So you can do the following: - # thetas = values[self.param] if isinstance(self.param, str) else self.param_name - # return thetas * ones_like(self.unitary(values)) - - def tensor( - self, - values: dict[str, Tensor] = dict(), - n_qubits: int | None = None, - diagonal: bool = False, - ) -> Operator: - """ - Get the corresponding unitary over n_qubits. - - Arguments: - values: Parameter value. - n_qubits: The number of qubits the unitary is represented over. - Can be higher than the number of qubit support. - diagonal: Whether the operation is diagonal. - - - Returns: - The unitary representation. - Raises: - NotImplementedError for the diagonal case. - """ - thetas = ( - values[self.param_name] - if isinstance(self.param_name, str) - else self.param_name - ) - return thetas * self.operations[0].tensor(values, n_qubits, diagonal) - - def flatten(self) -> list[Scale]: - """This method should only be called in the AdjointExpectation, - where the `Scale` is only supported for Primitive (and not Sequences) - so we don't want to flatten this to preserve the scale parameter. - - Returns: - The Scale within a list. - """ - return [self] - - def to(self, *args: Any, **kwargs: Any) -> Scale: - """Perform conversions for dtype or device. - - Returns: - Converted Scale. - """ - super().to(*args, **kwargs) - if not isinstance(self.param_name, str): - self.param_name = self.param_name.to(*args, **kwargs) - - return self - - -class Add(Sequence): - """ - The 'add' operation applies all 'operations' to 'state' and returns the sum of states. - - Attributes: - operations: List of operations to add up. - """ - - def __init__(self, operations: list[Module]): - - super().__init__(operations=operations) - - def forward( - self, - state: State, - values: dict[str, Tensor] | ParameterDict = dict(), - embedding: Embedding | None = None, - ) -> State: - """ - Apply the operations multiplying by the parameter values. - - Arguments: - state: Input state. - values: Parameter value. - - Returns: - The transformed state. - """ - return reduce(add, (op(state, values) for op in self.operations)) - - def tensor( - self, values: dict = {}, n_qubits: int | None = None, diagonal: bool = False - ) -> Tensor: - """ - Get the corresponding sum of unitaries over n_qubits. - - Arguments: - values: Parameter value. - n_qubits: The number of qubits the unitary is represented over. - Can be higher than the number of qubit support. - diagonal: Whether the operation is diagonal. - - - Returns: - The unitary representation. - Raises: - NotImplementedError for the diagonal case. - """ - if n_qubits is None: - n_qubits = max(self.qubit_support) + 1 - mat = torch.zeros((2, 2, 1), device=self.device) - for _ in range(n_qubits - 1): - mat = torch.kron(mat, torch.zeros((2, 2, 1), device=self.device)) - return reduce( - add, (op.tensor(values, n_qubits, diagonal) for op in self.operations), mat - ) - - -class Observable(Sequence): - """ - The Observable :math:`O` represents an operator from which - we can extract expectation values from quantum states. - - Given an input state :math:`\\ket\\rangle`, the expectation value with :math:`O` is defined as - :math:`\\langle\\bra|O\\ket\\rangle` - - Attributes: - operations: List of operations. - n_qubits: Number of qubits it is defined on. - """ - - def __init__( - self, - n_qubits: int | None, - operations: list[Module] | Primitive | Sequence, - ): - super().__init__(operations) - if n_qubits is None: - n_qubits = max(self.qubit_support) + 1 - self.n_qubits = n_qubits - - def run( - self, - state: Tensor, - values: dict[str, Tensor], - embedding: Embedding | None = None, - ) -> State: - """ - Apply the observable onto a state to obtain :math:`\\|O\\ket\\rangle`. - - Arguments: - state: Input state. - values: Values of parameters. - - Returns: - The transformed state. - """ - for op in self.operations: - state = op(state, values) - return state - - def forward( - self, - state: Tensor, - values: dict[str, Tensor] | ParameterDict = dict(), - embedding: Embedding | None = None, - ) -> Tensor: - """Calculate the inner product :math:`\\langle\\bra|O\\ket\\rangle` - - Arguments: - state: Input state. - values: Values of parameters. - - Returns: - The expectation value. - """ - return inner_prod(state, self.run(state, values)).real - - -class DiagonalObservable(Primitive): - """ - Special case of diagonal observables where computation is simpler. - We simply do a element-wise vector-product instead of a tensordot. - - Attributes: - pauli: The tensor representation from Primitive. - qubit_support: Qubits the operator acts on. - n_qubits: Number of qubits the operator is defined on. - """ - - def __init__( - self, - n_qubits: int | None, - operations: list[Module] | Primitive | Sequence, - to_sparse: bool = False, - ): - """Initializes the DiagonalObservable. - - Arguments: - n_qubits: Number of qubits the operator is defined on. - operations: Operations defining the observable. - to_sparse: Whether to convert the operator to its sparse representation or not. - """ - if isinstance(operations, list): - operations = Sequence(operations) - if n_qubits is None: - n_qubits = max(operations.qubit_support) + 1 - hamiltonian = operations.tensor({}, n_qubits).squeeze(2) - if to_sparse: - operator = operator_to_sparse_diagonal(hamiltonian) - else: - operator = torch.diag(hamiltonian).reshape(-1, 1) - super().__init__(operator, operations.qubit_support[0]) - self.qubit_support = operations.qubit_support - self.n_qubits = n_qubits - - def run( - self, - state: Tensor, - values: dict[str, Tensor], - embedding: Embedding | None = None, - ) -> Tensor: - """ - Apply the observable onto a state to obtain :math:`\\|O\\ket\\rangle`. - - We flatten the state, do a element-wise multiplication with the diagonal hamiltonian - and reshape it back to pyq-shape. - - - Arguments: - state: Input state. - values: Values of parameters. Unused here. - - Returns: - The transformed state. - """ - return torch.einsum( - "ij,ib->ib", self.pauli, state.flatten(start_dim=0, end_dim=-2) - ).reshape([2] * self.n_qubits + [state.shape[-1]]) - - def forward( - self, - state: Tensor, - values: dict[str, Tensor] | ParameterDict = dict(), - embedding: Embedding | None = None, - ) -> Tensor: - """Calculate the inner product :math:`\\langle\\bra|O\\ket\\rangle` - - Arguments: - state: Input state. - values: Values of parameters. - - Returns: - The expectation value. - """ - return inner_prod(state, self.run(state, values)).real - - -def is_diag_hamiltonian(hamiltonian: Operator, atol: Tensor = ATOL) -> bool: - """ - Returns True if the batched tensors H are diagonal. - - Arguments: - H: Input tensors. - atol: Tolerance for near-zero values. - - Returns: - True if diagonal, else False. - """ - diag_check = torch.tensor( - [ - is_diag(hamiltonian[..., i], atol) - for i in range(hamiltonian.shape[BATCH_DIM]) - ], - device=hamiltonian.device, - ) - return bool(torch.prod(diag_check)) - - -def evolve(hamiltonian: Operator, time_evolution: Tensor) -> Operator: - """Get the evolved operator. - - For a hamiltonian :math:`H` and a time evolution :math:`t`, returns :math:`exp(-i H, t)` - - Arguments: - hamiltonian: The operator :math:`H` for evolution. - time_evolution: The evolution time :math:`t`. - - Returns: - The evolution operator. - """ - if is_diag_hamiltonian(hamiltonian): - evol_operator = torch.diagonal(hamiltonian) * (-1j * time_evolution).view( - (-1, 1) - ) - evol_operator = torch.diag_embed(torch.exp(evol_operator)) - else: - evol_operator = torch.transpose(hamiltonian, 0, -1) * ( - -1j * time_evolution - ).view((-1, 1, 1)) - evol_operator = torch.linalg.matrix_exp(evol_operator) - return torch.transpose(evol_operator, 0, -1) - - -class HamiltonianEvolution(Sequence): - """ - The HamiltonianEvolution corresponds to :math:`t`, returns :math:`exp(-i H, t)` where - a hamiltonian/generator :math:`H` and a time evolution :math:`t` are given. - - We can create such operation by passing different generator types: - - A tensor representation of the generator, - - A string when we consider the generator as a symbol. - - Operations as a single primitive or a sequence, parameterized or not. - - Attributes: - generator_type: This attribute informs on how the generator is inputed - and sets the logic for applying hamiltonian evolution. - time: The evolution time :math:`t`. - operations: List of operations. - cache_length: LRU cache cache_length evolution operators for given set - of parameter values. - """ - - def __init__( - self, - generator: TGenerator, - time: Tensor | str, - qubit_support: Tuple[int, ...] | None = None, - generator_parametric: bool = False, - cache_length: int = 1, - is_time_dependent: bool = False, - duration: int = 100, - steps: int = 10, - solver=SolverType.KRYLOV_SE, - ): - """Initializes the HamiltonianEvolution. - Depending on the generator argument, set the type and set the right generator getter. - - Arguments: - generator: The generator :math:`H`. - time: The evolution time :math:`t`. - qubit_support: The qubits the operator acts on. - generator_parametric: Whether the generator is parametric or not. - """ - self.is_time_dependent = is_time_dependent - self.solver_name = solver - self.duration = duration - self.steps = steps - if isinstance(generator, Tensor): - assert ( - qubit_support is not None - ), "When using a Tensor generator, please pass a qubit_support." - if len(generator.shape) < 3: - generator = generator.unsqueeze(2) - generator = [Primitive(generator, target=-1)] - self.generator_type = GeneratorType.TENSOR - - elif isinstance(generator, str): - assert ( - qubit_support is not None - ), "When using a symbolic generator, please pass a qubit_support." - self.generator_type = GeneratorType.SYMBOL - self.generator_symbol = generator - generator = [] - elif isinstance(generator, (Primitive, Sequence)): - qubit_support = ( - generator.qubit_support - if ( - not qubit_support - or len(qubit_support) <= len(generator.qubit_support) - ) - else qubit_support - ) - if generator_parametric: - generator = [generator] - self.generator_type = GeneratorType.PARAMETRIC_OPERATION - else: - generator = [ - Primitive( - generator.tensor({}, len(qubit_support)), - target=generator.qubit_support[0], - ) - ] - - self.generator_type = GeneratorType.OPERATION - else: - raise TypeError( - f"Received generator of type {type(generator)},\ - allowed types are: [Tensor, str, Primitive, Sequence]" - ) - super().__init__(generator) - self._qubit_support = qubit_support # type: ignore - self.time = time - logger.debug("Hamiltonian Evolution initialized") - if logger.isEnabledFor(logging.DEBUG): - # When Debugging let's add logging and NVTX markers - # WARNING: incurs performance penalty - self.register_forward_hook(forward_hook, always_call=True) - self.register_full_backward_hook(backward_hook) - self.register_forward_pre_hook(pre_forward_hook) - self.register_full_backward_pre_hook(pre_backward_hook) - - self._generator_map: dict[GeneratorType, Callable[[dict], Tensor]] = { - GeneratorType.SYMBOL: self._symbolic_generator, - GeneratorType.TENSOR: self._tensor_generator, - GeneratorType.OPERATION: self._tensor_generator, - GeneratorType.PARAMETRIC_OPERATION: self._parametric_generator, - } - - # to avoid recomputing hamiltonians and evolution - self._cache_hamiltonian_evo: dict[str, Tensor] = dict() - self.cache_length = cache_length - - @property - def generator(self) -> ModuleList: - """Returns the operations making the generator. - - Returns: - The generator as a ModuleList. - """ - return self.operations - - def _symbolic_generator(self, values: dict) -> Operator: - """Returns the generator for the SYMBOL case. - - Returns: - The generator as a tensor. - """ - hamiltonian = values[self.generator_symbol] - # add batch dim - if len(hamiltonian.shape) == 2: - return hamiltonian.unsqueeze(2) - # cases when the batchdim is at index 0 instead of 2 - if len(hamiltonian.shape) == 3 and ( - hamiltonian.shape[0] != hamiltonian.shape[1] - ): - return torch.transpose(hamiltonian, 0, 2) - if len(hamiltonian.shape) == 4 and ( - hamiltonian.shape[0] != hamiltonian.shape[1] - ): - return torch.permute(hamiltonian.squeeze(3), (1, 2, 0)) - return hamiltonian - - def _tensor_generator(self, values: dict = {}) -> Operator: - """Returns the generator for the TENSOR and OPERATION cases. - - Arguments: - values: Non-used argument for consistency with other generator getters. - - Returns: - The generator as a tensor. - """ - return self.generator[0].pauli - - def _parametric_generator(self, values: dict) -> Operator: - """Returns the generator for the PARAMETRIC_OPERATION case. - - Arguments: - values: Values of parameters. - - Returns: - The generator as a tensor. - """ - return self.generator[0].tensor(values, len(self.qubit_support)) - - @property - def create_hamiltonian(self) -> Callable[[dict], Operator]: - """A utility method for setting the right generator getter depending on the init case. - - Returns: - The right generator getter. - """ - return self._generator_map[self.generator_type] - - def forward( - self, - state: Tensor, - values: dict[str, Tensor] | ParameterDict = dict(), - embedding: Embedding | None = None, - ) -> State: - """ - Apply the hamiltonian evolution with input parameter values. - - Arguments: - state: Input state. - values: Values of parameters. - - Returns: - The transformed state. - """ - if self.is_time_dependent and embedding is not None: - n_qubits = len(state.shape) - 1 - batch_size = state.shape[-1] - tsave = torch.linspace(0, float(self.duration), self.steps) - - def Ht(t: torch.Tensor) -> torch.Tensor: - """Accepts a value 't' for time - and returns a (2**n_qubits, 2**n_qubits) Hamiltonian. - """ - # We use the origial embedded params and return a new dict - # where we reembedded all parameters depending on time with value 't' - reembedded_time_values = embedding.reembed_tparam(values, t) - return ( - self.generator[0] - .tensor(reembedded_time_values, n_qubits) - .squeeze(2) - ) # squeeze the batch dim and return shape (2**n_qubits, 2**n_qubits) - - sol = sesolve( - Ht, # returns a tensor of shape (2**n_qubits, 2**n_qubits) - torch.flatten( - state, start_dim=0, end_dim=-2 - ), # reshape to (2**n_qubits, batch_size) - tsave, - self.solver_name, - ) - - state = sol.states[ - -1 - ] # Retrieve the last state of shape (2**n_qubits, batch_size) - - return state.reshape( - [2] * n_qubits + [batch_size] - ) # reshape to [2] * n_qubits + [batch_size] - - else: - evolved_op = self.tensor(values, None) - return apply_operator( - state=state, - operator=evolved_op, - qubits=self.qubit_support, - n_qubits=len(state.size()) - 1, - batch_size=evolved_op.shape[BATCH_DIM], - ) - - def tensor( - self, - values: dict = {}, - n_qubits: int | None = None, - diagonal: bool = False, - ) -> Operator: - """Get the corresponding unitary over n_qubits. - - To avoid computing the evolution operator, we store it in cache wrt values. - - Arguments: - values: Parameter value. - n_qubits: The number of qubits the unitary is represented over. - Can be higher than the number of qubit support. - diagonal: Whether the operation is diagonal. - - - Returns: - The unitary representation. - Raises: - NotImplementedError for the diagonal case. - """ - - values_cache_key = str(OrderedDict(values)) - if self.cache_length > 0 and values_cache_key in self._cache_hamiltonian_evo: - return self._cache_hamiltonian_evo[values_cache_key] - - if diagonal: - raise NotImplementedError - if n_qubits is None: - n_qubits = max(self.qubit_support) + 1 - hamiltonian: torch.Tensor = self.create_hamiltonian(values) - time_evolution: torch.Tensor = ( - values[self.time] if isinstance(self.time, str) else self.time - ) # If `self.time` is a string / hence, a Parameter, - # we expect the user to pass it in the `values` dict - evolved_op = evolve(hamiltonian, time_evolution) - nb_cached = len(self._cache_hamiltonian_evo) - - # LRU caching - if (nb_cached > 0) and (nb_cached == self.cache_length): - self._cache_hamiltonian_evo.pop(next(iter(self._cache_hamiltonian_evo))) - if nb_cached < self.cache_length: - self._cache_hamiltonian_evo[values_cache_key] = evolved_op - return evolved_op diff --git a/pyqtorch/api.py b/pyqtorch/api.py index d58f6dbc..124e7158 100644 --- a/pyqtorch/api.py +++ b/pyqtorch/api.py @@ -1,16 +1,22 @@ from __future__ import annotations from collections import Counter +from functools import partial from logging import getLogger +import torch from torch import Tensor -from pyqtorch.adjoint import AdjointExpectation -from pyqtorch.analog import Observable +from pyqtorch.apply import apply_operator from pyqtorch.circuit import QuantumCircuit +from pyqtorch.differentiation import ( + AdjointExpectation, + PSRExpectation, + check_support_psr, +) from pyqtorch.embed import Embedding -from pyqtorch.gpsr import PSRExpectation, check_support_psr -from pyqtorch.utils import DiffMode, inner_prod +from pyqtorch.hamiltonians import Observable +from pyqtorch.utils import DiffMode, sample_multinomial logger = getLogger(__name__) @@ -96,12 +102,100 @@ def sample( ) +def analytical_expectation( + circuit: QuantumCircuit, + state: Tensor, + observable: Observable, + values: dict[str, Tensor] = dict(), + embedding: Embedding | None = None, +) -> Tensor: + """Compute the analytical expectation value. + + Given an initial state :math:`\\ket\\rangle`, + a quantum circuit :math:`U(\\theta)`, + the analytical expectation value with :math:`O` is defined as + :math:`\\langle\\bra U_{\\dag}(\\theta) |O| U(\\theta) \\ket\\rangle` + + Args: + circuit (QuantumCircuit): Quantum circuit :math:`U(\\theta)`. + state (Tensor): Input state :math:`\\ket\\rangle`. + observable (Observable): Observable O. + values (dict[str, Tensor], optional): Parameter values for the observable if any. + embedding (Embedding | None, optional): An optional instance of `Embedding`. + + Returns: + Tensor: Expectation value. + """ + state = run(circuit, state, values, embedding=embedding) + return observable.expectation(state, values, embedding=embedding) + + +def sampled_expectation( + circuit: QuantumCircuit, + state: Tensor, + observable: Observable, + values: dict[str, Tensor] = dict(), + embedding: Embedding | None = None, + n_shots: int = 1, +) -> Tensor: + """Expectation value approximated via sampling. + + Given an input state :math:`\\ket\\rangle`, + the expectation analytical value with :math:`O` is defined as + :math:`\\langle\\bra|O\\ket\\rangle` + + Args: + state (Tensor): Input state :math:`\\ket\\rangle`. + observable (Observable): Observable O. + values (dict[str, Tensor], optional): Parameter values for the observable if any. + embedding (Embedding | None, optional): An optional instance of `Embedding`. + n_shots: (int, optional): Number of samples to compute expectation on. + + Returns: + Tensor: Expectation value. + """ + state = run(circuit, state, values, embedding=embedding) + n_qubits = circuit.n_qubits + + # batchsize needs to be first dim for eigh + eigvals, eigvecs = torch.linalg.eigh( + observable.tensor( + values=values, embedding=embedding, full_support=tuple(range(n_qubits)) + ).permute((2, 0, 1)) + ) + eigvals = eigvals.squeeze() + eigvec_state_prod = apply_operator( + state, + eigvecs.T.conj(), + tuple(range(n_qubits)), + n_qubits=circuit.n_qubits, + ) + eigvec_state_prod = torch.flatten(eigvec_state_prod, start_dim=0, end_dim=-2).t() + probs = torch.pow(torch.abs(eigvec_state_prod), 2) + batch_sample_multinomial = torch.func.vmap( + lambda p: sample_multinomial( + p, n_qubits, n_shots, return_counter=False, minlength=probs.shape[-1] + ), + randomness="different", + ) + batch_samples = batch_sample_multinomial(probs) + normalized_samples = torch.div( + batch_samples, torch.tensor(n_shots, dtype=probs.dtype) + ) + normalized_samples.requires_grad = True + expectations = torch.einsum( + "i,ji ->j", eigvals.to(dtype=normalized_samples.dtype), normalized_samples # type: ignore[union-attr] + ) + return expectations + + def expectation( circuit: QuantumCircuit, state: Tensor = None, values: dict[str, Tensor] = dict(), observable: Observable = None, # type: ignore[assignment] diff_mode: DiffMode = DiffMode.AD, + n_shots: int | None = None, embedding: Embedding | None = None, ) -> Tensor: """Compute the expectation value of `circuit` given a `state`, @@ -115,6 +209,8 @@ def expectation( denoting the current parameter values for each parameter in `circuit`. observable: A pyq.Observable instance. diff_mode: The differentiation mode. + n_shots: Number of shots for estimating expectation values. + Only used with DiffMode.GPSR or DiffMode.AD. embedding: An optional instance of `Embedding`. Returns: @@ -131,11 +227,12 @@ def expectation( circ = QuantumCircuit(n_qubits, [RY(0, 'theta')]) state = random_state(n_qubits) theta = tensor(pi, requires_grad=True) - observable = Observable(n_qubits, Add([Z(i) for i in range(n_qubits)])) + observable = Observable(Add([Z(i) for i in range(n_qubits)])) expval = expectation(circ, state, {'theta': theta}, observable, diff_mode = DiffMode.ADJOINT) dfdtheta= grad(expval, theta, ones_like(expval))[0] ``` """ + if embedding is not None and diff_mode != DiffMode.AD: raise NotImplementedError("Only diff_mode AD supports embedding") logger.debug( @@ -144,13 +241,19 @@ def expectation( ) if observable is None: logger.error("Please provide an observable to compute expectation.") + if state is None: state = circuit.init_state(batch_size=1) + + expectation_fn = analytical_expectation + if n_shots is not None: + if isinstance(n_shots, int) and n_shots > 0: + expectation_fn = partial(sampled_expectation, n_shots=n_shots) + else: + logger.error("Please provide a 'n_shots' in options of type 'int'.") + if diff_mode == DiffMode.AD: - state = run(circuit, state, values, embedding=embedding) - return inner_prod( - state, run(observable, state, values, embedding=embedding) - ).real + return expectation_fn(circuit, state, observable, values, embedding) elif diff_mode == DiffMode.ADJOINT: return AdjointExpectation.apply( circuit, @@ -163,7 +266,13 @@ def expectation( elif diff_mode == DiffMode.GPSR: check_support_psr(circuit) return PSRExpectation.apply( - circuit, state, observable, embedding, values.keys(), *values.values() + circuit, + state, + observable, + embedding, + expectation_fn, + values.keys(), + *values.values(), ) else: logger.error(f"Requested diff_mode '{diff_mode}' not supported.") diff --git a/pyqtorch/apply.py b/pyqtorch/apply.py index 828e286e..b1539c52 100644 --- a/pyqtorch/apply.py +++ b/pyqtorch/apply.py @@ -2,11 +2,12 @@ from string import ascii_letters as ABC -from numpy import array, log2 +from numpy import array from numpy.typing import NDArray from torch import Tensor, einsum -from pyqtorch.utils import promote_operator +from pyqtorch.matrices import _dagger +from pyqtorch.utils import DensityMatrix ABC_ARRAY: NDArray = array(list(ABC)) @@ -58,30 +59,45 @@ def apply_operator( return einsum(f"{operator_dims},{in_state_dims}->{out_state_dims}", operator, state) -def operator_product(op1: Tensor, op2: Tensor, target: int) -> Tensor: +def apply_density_mat(op: Tensor, density_matrix: DensityMatrix) -> DensityMatrix: + """ + Apply an operator to a density matrix, i.e., compute: + op1 * density_matrix * op1_dagger. + + Args: + op (Tensor): The operator to apply. + density_matrix (DensityMatrix): The density matrix. + + Returns: + DensityMatrix: The resulting density matrix after applying the operator and its dagger. + """ + batch_size_op = op.size(-1) + batch_size_dm = density_matrix.size(-1) + if batch_size_dm > batch_size_op: + # The other condition is impossible because + # operators are always initialized with batch_size = 1. + op = op.repeat(1, 1, batch_size_dm) + return einsum("ijb,jkb,klb->ilb", op, density_matrix, _dagger(op)) + + +def operator_product(op1: Tensor, op2: Tensor) -> Tensor: """ Compute the product of two operators. + `torch.bmm` is not suitable for our purposes because, + in our convention, the batch_size is in the last dimension. Args: op1 (Tensor): The first operator. op2 (Tensor): The second operator. - target (int): The target qubit index. - Returns: Tensor: The product of the two operators. """ - - n_qubits_1 = int(log2(op1.size(1))) - n_qubits_2 = int(log2(op2.size(1))) + # ? Should we continue to adjust the batch here? + # ? as now all gates are init with batch_size=1. batch_size_1 = op1.size(-1) batch_size_2 = op2.size(-1) - if n_qubits_1 > n_qubits_2: - op2 = promote_operator(op2, target, n_qubits_1) - elif n_qubits_1 < n_qubits_2: - op1 = promote_operator(op1, target, n_qubits_2) if batch_size_1 > batch_size_2: op2 = op2.repeat(1, 1, batch_size_1)[:, :, :batch_size_1] elif batch_size_2 > batch_size_1: op1 = op1.repeat(1, 1, batch_size_2)[:, :, :batch_size_2] - return einsum("ijb,jkb->ikb", op1, op2) diff --git a/pyqtorch/circuit.py b/pyqtorch/circuit.py index 3036bad8..ce7f7f9d 100644 --- a/pyqtorch/circuit.py +++ b/pyqtorch/circuit.py @@ -1,158 +1,28 @@ from __future__ import annotations -import logging from collections import Counter from functools import reduce from logging import getLogger from operator import add -from typing import Any, Generator, Iterator, NoReturn import torch -from numpy import int64 -from torch import Tensor, complex128, einsum, rand -from torch import device as torch_device -from torch import dtype as torch_dtype -from torch.nn import Module, ModuleList, ParameterDict +from torch import Tensor, bernoulli, tensor +from torch.nn import Module, ParameterDict -from pyqtorch.apply import apply_operator +from pyqtorch.composite import Sequence from pyqtorch.embed import Embedding -from pyqtorch.matrices import IMAT, add_batch_dim -from pyqtorch.parametric import RX, RY, Parametric -from pyqtorch.primitive import CNOT, Primitive -from pyqtorch.utils import State, product_state, zero_state +from pyqtorch.utils import ( + DensityMatrix, + DropoutMode, + State, + product_state, + sample_multinomial, + zero_state, +) logger = getLogger(__name__) -def forward_hook(*args, **kwargs) -> None: # type: ignore[no-untyped-def] - logger.debug("Forward complete") - torch.cuda.nvtx.range_pop() - - -def pre_forward_hook(*args, **kwargs) -> None: # type: ignore[no-untyped-def] - logger.debug("Executing forward") - torch.cuda.nvtx.range_push("QuantumCircuit.forward") - - -def backward_hook(*args, **kwargs) -> None: # type: ignore[no-untyped-def] - logger.debug("Backward complete") - torch.cuda.nvtx.range_pop() - - -def pre_backward_hook(*args, **kwargs) -> None: # type: ignore[no-untyped-def] - logger.debug("Executed backward") - torch.cuda.nvtx.range_push("QuantumCircuit.backward") - - -class Sequence(Module): - """A generic container for pyqtorch operations""" - - def __init__(self, operations: list[Module]): - super().__init__() - self.operations = ModuleList(operations) - self._device = torch_device("cpu") - self._dtype = complex128 - if len(self.operations) > 0: - try: - self._device = next(iter(set((op.device for op in self.operations)))) - except StopIteration: - pass - logger.debug("QuantumCircuit initialized") - if logger.isEnabledFor(logging.DEBUG): - # When Debugging let's add logging and NVTX markers - # WARNING: incurs performance penalty - self.register_forward_hook(forward_hook, always_call=True) - self.register_full_backward_hook(backward_hook) - self.register_forward_pre_hook(pre_forward_hook) - self.register_full_backward_pre_hook(pre_backward_hook) - - self._qubit_support = tuple( - set( - sum( - [op.qubit_support for op in self.operations], - (), - ) - ) - ) - - self._qubit_support = tuple( - map( - lambda x: x if isinstance(x, (int, int64)) else x[0], - self._qubit_support, - ) - ) - assert all( - [isinstance(q, (int, int64)) for q in self._qubit_support] - ) # TODO fix numpy.int issue - - @property - def qubit_support(self) -> tuple: - return self._qubit_support - - def __iter__(self) -> Iterator: - return iter(self.operations) - - def __len__(self) -> int: - return len(self.operations) - - def __hash__(self) -> int: - return hash(reduce(add, (hash(op) for op in self.operations))) - - def forward( - self, - state: State, - values: dict[str, Tensor] | ParameterDict = {}, - embedding: Embedding | None = None, - ) -> State: - for op in self.operations: - state = op(state, values, embedding) - return state - - @property - def device(self) -> torch_device: - return self._device - - @property - def dtype(self) -> torch_dtype: - return self._dtype - - def to(self, *args: Any, **kwargs: Any) -> Sequence: - self.operations = ModuleList([op.to(*args, **kwargs) for op in self.operations]) - if len(self.operations) > 0: - self._device = self.operations[0].device - self._dtype = self.operations[0].dtype - return self - - def flatten(self) -> ModuleList: - ops = [] - for op in self.operations: - if isinstance(op, Sequence): - ops += op.flatten() - else: - ops.append(op) - return ModuleList(ops) - - def tensor( - self, - values: dict[str, Tensor] = dict(), - n_qubits: int | None = None, - diagonal: bool = False, - ) -> Tensor: - if diagonal: - raise NotImplementedError - if n_qubits is None: - n_qubits = max(self.qubit_support) + 1 - mat = IMAT.clone().unsqueeze(2).to(self.device) - for _ in range(n_qubits - 1): - mat = torch.kron(mat, IMAT.clone().unsqueeze(2).to(self.device)) - - return reduce( - lambda t0, t1: einsum("ijb,jkb->ikb", t1, t0), - (add_batch_dim(op.tensor(values, n_qubits)) for op in self.operations), - mat, - ) - - class QuantumCircuit(Sequence): """A QuantumCircuit defining a register / number of qubits of the full system.""" @@ -162,8 +32,8 @@ def __init__(self, n_qubits: int, operations: list[Module]): def run( self, - state: State = None, - values: dict[str, Tensor] | ParameterDict = {}, + state: Tensor = None, + values: dict[str, Tensor] | ParameterDict = dict(), embedding: Embedding | None = None, ) -> State: if state is None: @@ -193,127 +63,131 @@ def sample( embedding: Embedding | None = None, ) -> list[Counter]: if n_shots < 1: - raise ValueError("You can only call sample with n_shots>0.") - - def sample(probs: Tensor) -> Counter: - return Counter( - { - format(k, "0{}b".format(self.n_qubits)): count.item() - for k, count in enumerate( - torch.bincount( - torch.multinomial( - input=probs, num_samples=n_shots, replacement=True - ) - ) - ) - if count > 0 - } + raise ValueError( + f"You can only call sample with a non-negative value for `n_shots`. Got {n_shots}." ) with torch.no_grad(): - state = torch.flatten( - self.run(state=state, values=values, embedding=embedding), - start_dim=0, - end_dim=-2, - ).t() + state = self.run(state=state, values=values, embedding=embedding) + if isinstance(state, DensityMatrix): + probs = torch.diagonal(state, dim1=0, dim2=1).real + else: + state = torch.flatten( + state, + start_dim=0, + end_dim=-2, + ).t() + probs = torch.abs(torch.pow(state, 2)) + + probs = torch.pow(torch.abs(state), 2) + return list( + map(lambda p: sample_multinomial(p, self.n_qubits, n_shots), probs) + ) - probs = torch.abs(torch.pow(state, 2)) - return list(map(lambda p: sample(p), probs)) +class DropoutQuantumCircuit(QuantumCircuit): + """Creates a quantum circuit able to perform quantum dropout, based on the work of https://arxiv.org/abs/2310.04120. + Args: + dropout_mode (DropoutMode): type of dropout to perform. Defaults to DropoutMode.ROTATIONAL + dropout_prob (float): dropout probability. Defaults to 0.06. + """ -class Merge(Sequence): def __init__( self, + n_qubits: int, operations: list[Module], + dropout_mode: DropoutMode = DropoutMode.ROTATIONAL, + dropout_prob: float = 0.06, ): + super().__init__(n_qubits, operations) + self.dropout_mode = dropout_mode + self.dropout_prob = dropout_prob + + self.dropout_fn = getattr(self, dropout_mode) + + def forward( + self, + state: State, + values: dict[str, Tensor] | ParameterDict = dict(), + embedding: Embedding | None = None, + ) -> State: + if self.training: + state = self.dropout_fn(state, values) + else: + for op in self.operations: + state = op(state, values) + return state + + def rotational_dropout( + self, state: State = None, values: dict[str, Tensor] | ParameterDict = dict() + ) -> State: + """Randomly drops entangling rotational gates. + + Args: + state (State, optional): pure state vector . Defaults to None. + values (dict[str, Tensor] | ParameterDict, optional): gate parameters. Defaults to {}. + + Returns: + State: pure state vector """ - Merge a sequence of single qubit operations acting on the same qubit into a single - einsum operation. + for op in self.operations: + if not ( + (hasattr(op, "param_name")) + and (values[op.param_name].requires_grad) + and not (int(1 - bernoulli(tensor(self.dropout_prob)))) + ): + state = op(state, values) + + return state + + def entangling_dropout( + self, state: State = None, values: dict[str, Tensor] | ParameterDict = dict() + ) -> State: + """Randomly drops entangling gates. - Arguments: - operations: A list of single qubit operations. + Args: + state (State, optional): pure state vector. Defaults to None. + values (dict[str, Tensor] | ParameterDict, optional): gate parameters. Defaults to {}. + Returns: + State: pure state vector """ + for op in self.operations: + has_param = hasattr(op, "param_name") + keep = int(1 - bernoulli(tensor(self.dropout_prob))) - if ( - isinstance(operations, (list, ModuleList)) - and all([isinstance(op, (Primitive, Parametric)) for op in operations]) - and all(list([len(op.qubit_support) == 1 for op in operations])) - and len(list(set([op.qubit_support[0] for op in operations]))) == 1 - ): - # We want all operations to act on the same qubit + if has_param or keep: + state = op(state, values) - super().__init__(operations) - self.qubits = operations[0].qubit_support - else: - raise TypeError( - f"Require all operations to act on a single qubit. Got: {operations}." - ) + return state - def forward( - self, - state: Tensor, - values: dict[str, Tensor] | None = None, - embedding: Embedding | None = None, - ) -> Tensor: - batch_size = state.shape[-1] - if values: - batch_size = max( - batch_size, - max( - list( - map( - lambda t: len(t) if isinstance(t, Tensor) else 1, - values.values(), - ) - ) - ), - ) - return apply_operator( - state, - self.unitary(values, embedding, batch_size), - self.qubits, - ) + def canonical_fwd_dropout( + self, state: State = None, values: dict[str, Tensor] | ParameterDict = dict() + ) -> State: + """Randomly drops rotational gates and next immediate entangling + gates whose target bit is located on dropped rotational gates. - def unitary( - self, - values: dict[str, Tensor] | None, - embedding: Embedding | None, - batch_size: int, - ) -> Tensor: - # We reverse the list of tensors here since matmul is not commutative. - return reduce( - lambda u0, u1: einsum("ijb,jkb->ikb", u0, u1), - ( - add_batch_dim(op.unitary(values, embedding), batch_size) - for op in reversed(self.operations) - ), - ) + Args: + state (State, optional): pure state vector. Defaults to None. + values (dict[str, Tensor] | ParameterDict, optional): gate parameters. Defaults to {}. + Returns: + State: pure state vector + """ + entanglers_to_drop = dict.fromkeys(range(state.ndim - 1), 0) # type: ignore + for op in self.operations: + if ( + hasattr(op, "param_name") + and (values[op.param_name].requires_grad) + and not (int(1 - bernoulli(tensor(self.dropout_prob)))) + ): + entanglers_to_drop[op.target] = 1 + else: + if not hasattr(op, "param_name") and ( + entanglers_to_drop[op.control[0]] == 1 + ): + entanglers_to_drop[op.control[0]] = 0 + else: + state = op(state, values) -def hea(n_qubits: int, depth: int, param_name: str) -> tuple[ModuleList, ParameterDict]: - def _idx() -> Generator[int, Any, NoReturn]: - i = 0 - while True: - yield i - i += 1 - - def idxer() -> Generator[int, Any, None]: - yield from _idx() - - idx = idxer() - ops = [] - for _ in range(depth): - layer = [] - for i in range(n_qubits): - layer += [ - Merge([fn(i, f"{param_name}_{next(idx)}") for fn in [RX, RY, RX]]) - ] - ops += layer - ops += [ - Sequence([CNOT(i % n_qubits, (i + 1) % n_qubits) for i in range(n_qubits)]) - ] - params = ParameterDict( - {f"{param_name}_{n}": rand(1, requires_grad=True) for n in range(next(idx))} - ) - return ops, params + return state diff --git a/pyqtorch/composite/__init__.py b/pyqtorch/composite/__init__.py new file mode 100644 index 00000000..6a66a230 --- /dev/null +++ b/pyqtorch/composite/__init__.py @@ -0,0 +1,4 @@ +from __future__ import annotations + +from .compose import Add, Merge, Scale, hea +from .sequence import Sequence diff --git a/pyqtorch/composite/compose.py b/pyqtorch/composite/compose.py new file mode 100644 index 00000000..4e2a5ade --- /dev/null +++ b/pyqtorch/composite/compose.py @@ -0,0 +1,282 @@ +from __future__ import annotations + +from functools import reduce +from logging import getLogger +from operator import add +from typing import Any, Generator, NoReturn, Union + +import torch +from torch import Tensor, einsum, rand +from torch.nn import Module, ModuleList, ParameterDict + +from pyqtorch.apply import apply_operator +from pyqtorch.embed import Embedding +from pyqtorch.matrices import add_batch_dim +from pyqtorch.primitives import CNOT, RX, RY, Parametric, Primitive +from pyqtorch.utils import ( + Operator, + State, +) + +from .sequence import Sequence + +BATCH_DIM = 2 + +logger = getLogger(__name__) + + +class Scale(Sequence): + """ + Generic container for multiplying a 'Primitive', 'Sequence' or 'Add' instance by a parameter. + + Attributes: + operations: Operations as a Sequence, Add, or a single Primitive operation. + param_name: Name of the parameter to multiply operations with. + """ + + def __init__( + self, operations: Union[Primitive, Sequence, Add], param_name: str | Tensor + ): + """ + Initializes a Scale object. + + Arguments: + operations: Operations as a Sequence, Add, or a single Primitive operation. + param_name: Name of the parameter to multiply operations with. + """ + if not isinstance(operations, (Primitive, Sequence, Add)): + raise ValueError("Scale only supports a single operation, Sequence or Add.") + super().__init__([operations]) + self.param_name = param_name + + def forward( + self, + state: Tensor, + values: dict[str, Tensor] | ParameterDict = dict(), + embedding: Embedding | None = None, + ) -> State: + """ + Apply the operation(s) multiplying by the parameter value. + + Arguments: + state: Input state. + values: Parameter value. + + Returns: + The transformed state. + """ + + scale = ( + values[self.param_name] + if isinstance(self.param_name, str) + else self.param_name + ) + return scale * self.operations[0].forward(state, values, embedding) + + def tensor( + self, + values: dict[str, Tensor] = dict(), + embedding: Embedding | None = None, + full_support: tuple[int, ...] | None = None, + ) -> Operator: + """ + Get the corresponding unitary over n_qubits. + + Arguments: + values: Parameter value. + embedding: An optional embedding. + full_support: Can be higher than the number of qubit support. + + Returns: + The unitary representation. + """ + scale = ( + values[self.param_name] + if isinstance(self.param_name, str) + else self.param_name + ) + return scale * self.operations[0].tensor(values, embedding, full_support) + + def flatten(self) -> list[Scale]: + """This method should only be called in the AdjointExpectation, + where the `Scale` is only supported for Primitive (and not Sequences) + so we don't want to flatten this to preserve the scale parameter. + + Returns: + The Scale within a list. + """ + return [self] + + def to(self, *args: Any, **kwargs: Any) -> Scale: + """Perform conversions for dtype or device. + + Returns: + Converted Scale. + """ + super().to(*args, **kwargs) + if not isinstance(self.param_name, str): + self.param_name = self.param_name.to(*args, **kwargs) + + return self + + +class Add(Sequence): + """ + The 'add' operation applies all 'operations' to 'state' and returns the sum of states. + + Attributes: + operations: List of operations to add up. + """ + + def __init__(self, operations: list[Module]): + + super().__init__(operations=operations) + + def forward( + self, + state: State, + values: dict[str, Tensor] | ParameterDict = dict(), + embedding: Embedding | None = None, + ) -> State: + """ + Apply the operations multiplying by the parameter values. + + Arguments: + state: Input state. + values: Parameter value. + + Returns: + The transformed state. + """ + return reduce(add, (op(state, values, embedding) for op in self.operations)) + + def tensor( + self, + values: dict = dict(), + embedding: Embedding | None = None, + full_support: tuple[int, ...] | None = None, + ) -> Tensor: + """ + Get the corresponding sum of unitaries over n_qubits. + + Arguments: + values: Parameter value. + Can be higher than the number of qubit support. + + + Returns: + The unitary representation. + """ + if full_support is None: + full_support = self.qubit_support + elif not set(self.qubit_support).issubset(set(full_support)): + raise ValueError( + "Expanding tensor operation requires a `full_support` argument " + "larger than or equal to the `qubit_support`." + ) + mat = torch.zeros( + (2 ** len(full_support), 2 ** len(full_support), 1), device=self.device + ) + return reduce( + add, + (op.tensor(values, embedding, full_support) for op in self.operations), + mat, + ) + + +class Merge(Sequence): + def __init__( + self, + operations: list[Module], + ): + """ + Merge a sequence of single qubit operations acting on the same qubit into a single + einsum operation. + + Arguments: + operations: A list of single qubit operations. + + """ + + if ( + isinstance(operations, (list, ModuleList)) + and all([isinstance(op, (Primitive, Parametric)) for op in operations]) + and all(list([len(op.qubit_support) == 1 for op in operations])) + and len(list(set([op.qubit_support[0] for op in operations]))) == 1 + ): + # We want all operations to act on the same qubit + + super().__init__(operations) + self.qubits = operations[0].qubit_support + else: + raise TypeError( + f"Require all operations to act on a single qubit. Got: {operations}." + ) + + def forward( + self, + state: Tensor, + values: dict[str, Tensor] = dict(), + embedding: Embedding | None = None, + ) -> Tensor: + batch_size = state.shape[-1] + if values: + batch_size = max( + batch_size, + max( + list( + map( + lambda t: len(t) if isinstance(t, Tensor) else 1, + values.values(), + ) + ) + ), + ) + return apply_operator( + state, + add_batch_dim(self.tensor(values, embedding), batch_size), + self.qubits, + ) + + def tensor( + self, + values: dict[str, Tensor] = dict(), + embedding: Embedding | None = None, + full_support: tuple[int, ...] | None = None, + ) -> Tensor: + # We reverse the list of tensors here since matmul is not commutative. + return reduce( + lambda u0, u1: einsum("ijb,jkb->ikb", u0, u1), + ( + op.tensor(values, embedding, full_support) + for op in reversed(self.operations) + ), + ) + + +def hea(n_qubits: int, depth: int, param_name: str) -> tuple[ModuleList, ParameterDict]: + def _idx() -> Generator[int, Any, NoReturn]: + i = 0 + while True: + yield i + i += 1 + + def idxer() -> Generator[int, Any, None]: + yield from _idx() + + idx = idxer() + ops = [] + for _ in range(depth): + layer = [] + for i in range(n_qubits): + layer += [ + Merge([fn(i, f"{param_name}_{next(idx)}") for fn in [RX, RY, RX]]) + ] + ops += layer + ops += [ + Sequence([CNOT(i % n_qubits, (i + 1) % n_qubits) for i in range(n_qubits)]) + ] + params = ParameterDict( + {f"{param_name}_{n}": rand(1, requires_grad=True) for n in range(next(idx))} + ) + return ops, params diff --git a/pyqtorch/composite/sequence.py b/pyqtorch/composite/sequence.py new file mode 100644 index 00000000..97174300 --- /dev/null +++ b/pyqtorch/composite/sequence.py @@ -0,0 +1,163 @@ +from __future__ import annotations + +import logging +from functools import reduce +from logging import getLogger +from operator import add +from typing import Any, Iterator + +import torch +from numpy import int64 +from torch import Tensor, complex128, einsum +from torch import device as torch_device +from torch import dtype as torch_dtype +from torch.nn import Module, ModuleList, ParameterDict + +from pyqtorch.embed import Embedding +from pyqtorch.matrices import _dagger, add_batch_dim +from pyqtorch.utils import ( + State, +) + +logger = getLogger(__name__) + + +def forward_hook(*args, **kwargs) -> None: # type: ignore[no-untyped-def] + logger.debug("Forward complete") + torch.cuda.nvtx.range_pop() + + +def pre_forward_hook(*args, **kwargs) -> None: # type: ignore[no-untyped-def] + logger.debug("Executing forward") + torch.cuda.nvtx.range_push("QuantumCircuit.forward") + + +def backward_hook(*args, **kwargs) -> None: # type: ignore[no-untyped-def] + logger.debug("Backward complete") + torch.cuda.nvtx.range_pop() + + +def pre_backward_hook(*args, **kwargs) -> None: # type: ignore[no-untyped-def] + logger.debug("Executed backward") + torch.cuda.nvtx.range_push("QuantumCircuit.backward") + + +class Sequence(Module): + """A generic container for pyqtorch operations""" + + def __init__(self, operations: list[Module]): + super().__init__() + self.operations = ModuleList(operations) + self._device = torch_device("cpu") + self._dtype = complex128 + if len(self.operations) > 0: + try: + self._device = next(iter(set((op.device for op in self.operations)))) + except StopIteration: + pass + logger.debug("QuantumCircuit initialized") + if logger.isEnabledFor(logging.DEBUG): + # When Debugging let's add logging and NVTX markers + # WARNING: incurs performance penalty + self.register_forward_hook(forward_hook, always_call=True) + self.register_full_backward_hook(backward_hook) + self.register_forward_pre_hook(pre_forward_hook) + self.register_full_backward_pre_hook(pre_backward_hook) + + self._qubit_support = tuple( + set( + sum( + [op.qubit_support for op in self.operations], + (), + ) + ) + ) + + self._qubit_support = tuple( + map( + lambda x: x if isinstance(x, (int, int64)) else x[0], + self._qubit_support, + ) + ) + assert all( + [isinstance(q, (int, int64)) for q in self._qubit_support] + ) # TODO fix numpy.int issue + + @property + def qubit_support(self) -> tuple: + return self._qubit_support + + def __iter__(self) -> Iterator: + return iter(self.operations) + + def __len__(self) -> int: + return len(self.operations) + + def __hash__(self) -> int: + return hash(reduce(add, (hash(op) for op in self.operations))) + + def forward( + self, + state: Tensor, + values: dict[str, Tensor] | ParameterDict = dict(), + embedding: Embedding | None = None, + ) -> State: + for op in self.operations: + state = op(state, values, embedding) + return state + + @property + def device(self) -> torch_device: + return self._device + + @property + def dtype(self) -> torch_dtype: + return self._dtype + + def to(self, *args: Any, **kwargs: Any) -> Sequence: + self.operations = ModuleList([op.to(*args, **kwargs) for op in self.operations]) + if len(self.operations) > 0: + self._device = self.operations[0].device + self._dtype = self.operations[0].dtype + return self + + def flatten(self) -> ModuleList: + ops = [] + for op in self.operations: + if isinstance(op, Sequence): + ops += op.flatten() + else: + ops.append(op) + return ModuleList(ops) + + def tensor( + self, + values: dict[str, Tensor] = dict(), + embedding: Embedding | None = None, + full_support: tuple[int, ...] | None = None, + ) -> Tensor: + if full_support is None: + full_support = self.qubit_support + elif not set(self.qubit_support).issubset(set(full_support)): + raise ValueError( + "Expanding tensor operation requires a `full_support` argument " + "larger than or equal to the `qubit_support`." + ) + mat = torch.eye( + 2 ** len(full_support), dtype=self.dtype, device=self.device + ).unsqueeze(2) + return reduce( + lambda t0, t1: einsum("ijb,jkb->ikb", t1, t0), + ( + add_batch_dim(op.tensor(values, embedding, full_support)) + for op in self.operations + ), + mat, + ) + + def dagger( + self, + values: dict[str, Tensor] | Tensor = dict(), + embedding: Embedding | None = None, + ) -> Tensor: + return _dagger(self.tensor(values, embedding)) diff --git a/pyqtorch/differentiation/__init__.py b/pyqtorch/differentiation/__init__.py new file mode 100644 index 00000000..12cc1fe1 --- /dev/null +++ b/pyqtorch/differentiation/__init__.py @@ -0,0 +1,4 @@ +from __future__ import annotations + +from .adjoint import AdjointExpectation +from .gpsr import PSRExpectation, check_support_psr diff --git a/pyqtorch/adjoint.py b/pyqtorch/differentiation/adjoint.py similarity index 92% rename from pyqtorch/adjoint.py rename to pyqtorch/differentiation/adjoint.py index c0471846..1da56896 100644 --- a/pyqtorch/adjoint.py +++ b/pyqtorch/differentiation/adjoint.py @@ -6,12 +6,12 @@ from torch import Tensor, no_grad from torch.autograd import Function -from pyqtorch.analog import Observable, Scale from pyqtorch.apply import apply_operator from pyqtorch.circuit import QuantumCircuit +from pyqtorch.composite import Scale from pyqtorch.embed import Embedding -from pyqtorch.parametric import Parametric -from pyqtorch.primitive import Primitive +from pyqtorch.hamiltonians import Observable +from pyqtorch.primitives import Parametric, Primitive from pyqtorch.utils import inner_prod, param_dict logger = getLogger(__name__) @@ -65,7 +65,7 @@ def forward( logger.error(msg) raise NotImplementedError(msg) ctx.out_state = circuit.run(state, values, embedding) - ctx.projected_state = observable.run(ctx.out_state, values) + ctx.projected_state = observable(ctx.out_state, values) ctx.save_for_backward(*param_values) return inner_prod(ctx.out_state, ctx.projected_state).real @@ -76,7 +76,7 @@ def backward(ctx: Any, grad_out: Tensor) -> Tuple[None, ...]: values = param_dict(ctx.param_names, param_values) grads_dict = {k: None for k in values.keys()} for op in ctx.circuit.flatten()[::-1]: - if isinstance(op, Primitive): + if isinstance(op, (Primitive, Parametric)): ctx.out_state = apply_operator( ctx.out_state, op.dagger(values, ctx.embedding), op.qubit_support ) @@ -98,7 +98,12 @@ def backward(ctx: Any, grad_out: Tensor) -> Tuple[None, ...]: op.dagger(values, ctx.embedding), op.qubit_support, ) + elif isinstance(op, Scale): + # TODO: Properly fix or not support Scale in adjoint at all + logger.error( + f"AdjointExpectation does not support operation: {type(op)}." + ) if not len(op.operations) == 1 and isinstance( op.operations[0], Primitive ): diff --git a/pyqtorch/gpsr.py b/pyqtorch/differentiation/gpsr.py similarity index 59% rename from pyqtorch/gpsr.py rename to pyqtorch/differentiation/gpsr.py index b25fb73d..e784c748 100644 --- a/pyqtorch/gpsr.py +++ b/pyqtorch/differentiation/gpsr.py @@ -1,17 +1,19 @@ from __future__ import annotations from logging import getLogger -from typing import Any, Tuple +from typing import Any, Callable, Tuple import torch from torch import Tensor, no_grad from torch.autograd import Function -from pyqtorch.analog import HamiltonianEvolution, Observable, Scale -from pyqtorch.circuit import QuantumCircuit, Sequence +from pyqtorch.circuit import QuantumCircuit +from pyqtorch.composite import Scale, Sequence from pyqtorch.embed import Embedding -from pyqtorch.parametric import Parametric -from pyqtorch.utils import inner_prod, param_dict +from pyqtorch.hamiltonians import HamiltonianEvolution, Observable +from pyqtorch.matrices import DEFAULT_REAL_DTYPE +from pyqtorch.primitives import Parametric +from pyqtorch.utils import param_dict logger = getLogger(__name__) @@ -68,8 +70,10 @@ class PSRExpectation(Function): Arguments: circuit: A QuantumCircuit instance - observable: A hamiltonian. state: A state in the form of [2 * n_qubits + [batch_size]] + observable: An hermitian operator. + embedding: An optional instance of `Embedding`. + expectation_method: Callable for computing expectations. param_names: A list of parameter names. *param_values: A unpacked tensor of values for each parameter. @@ -87,6 +91,7 @@ def forward( state: Tensor, observable: Observable, embedding: Embedding, + expectation_method: Callable, param_names: list[str], *param_values: Tensor, ) -> Tensor: @@ -98,19 +103,14 @@ def forward( ctx.state = state ctx.embedding = embedding values = param_dict(param_names, param_values) - ctx.out_state = circuit.run(state, values) - ctx.projected_state = observable.run(ctx.out_state, values) ctx.save_for_backward(*param_values) - return inner_prod(ctx.out_state, ctx.projected_state).real + ctx.expectation_method = expectation_method + return expectation_method(circuit, state, observable, values, embedding) @staticmethod def backward(ctx: Any, grad_out: Tensor) -> Tuple[None, ...]: """The PSRExpectation backward call. - Note that only operations with two distinct eigenvalues - from their generator (i.e., compatible with single_gap_shift) - are supported at the moment. - Arguments: ctx (Any): Context object for accessing stored information. grad_out (Tensor): Current jacobian tensor. @@ -123,7 +123,15 @@ def backward(ctx: Any, grad_out: Tensor) -> Tuple[None, ...]: """ values = param_dict(ctx.param_names, ctx.saved_tensors) - shift = torch.tensor(torch.pi) / 2.0 + shift_pi2 = torch.tensor(torch.pi) / 2.0 + shift_multi = 0.5 + + dtype_values = DEFAULT_REAL_DTYPE + device = torch.device("cpu") + try: + dtype_values, device = [(v.dtype, v.device) for v in values.values()][0] + except Exception: + pass def expectation_fn(values: dict[str, Tensor]) -> Tensor: """Use the PSRExpectation for nested grad calls. @@ -139,16 +147,17 @@ def expectation_fn(values: dict[str, Tensor]) -> Tensor: ctx.state, ctx.observable, ctx.embedding, + ctx.expectation_method, values.keys(), *values.values(), ) def single_gap_shift( param_name: str, - values: dict[str, torch.Tensor], - spectral_gap: torch.Tensor, - shift: torch.Tensor = torch.tensor(torch.pi) / 2.0, - ) -> torch.Tensor: + values: dict[str, Tensor], + spectral_gap: Tensor, + shift: Tensor = torch.tensor(torch.pi) / 2.0, + ) -> Tensor: """Implements single gap PSR rule. Args: @@ -160,6 +169,12 @@ def single_gap_shift( Returns: Gradient evaluation for param_name. """ + + # device conversions + spectral_gap = spectral_gap.to(device=device) + shift = shift.to(device=device) + + # apply shift rule shifted_values = values.copy() shifted_values[param_name] = shifted_values[param_name] + shift f_plus = expectation_fn(shifted_values) @@ -171,9 +186,65 @@ def single_gap_shift( / (4 * torch.sin(spectral_gap * shift / 2)) ) - def multi_gap_shift(*args, **kwargs) -> Tensor: - """Implements multi gap PSR rule.""" - raise NotImplementedError("Multi-gap is not yet supported.") + def multi_gap_shift( + param_name: str, + values: dict[str, Tensor], + spectral_gaps: Tensor, + shift_prefac: float = 0.5, + ) -> Tensor: + """Implement multi gap PSR rule. + + See Kyriienko1 and Elfving, 2021 for details: + https://arxiv.org/pdf/2108.01218.pdf + + Args: + param_name: Name of the parameter to apply PSR. + values: Dictionary with parameter values. + spectral_gaps: Spectral gaps value for PSR. + shift_prefac: Shift prefactor value for PSR shifts. + Defaults to torch.tensor(0.5). + + Returns: + Gradient evaluation for param_name. + """ + n_eqs = len(spectral_gaps) + dtype = torch.promote_types(dtype_values, spectral_gaps.dtype) + spectral_gaps = spectral_gaps.to(device=device) + PI = torch.tensor(torch.pi, dtype=dtype) + shifts = shift_prefac * torch.linspace( + PI / 2.0 - PI / 5.0, PI / 2.0 + PI / 5.0, n_eqs, dtype=dtype + ) + shifts = shifts.to(device=device) + + # calculate F vector and M matrix + # (see: https://arxiv.org/pdf/2108.01218.pdf on p. 4 for definitions) + F = [] + M = torch.empty((n_eqs, n_eqs), dtype=dtype).to(device=device) + batch_size = 1 + shifted_params = values.copy() + for i in range(n_eqs): + # + shift + shifted_params[param_name] = shifted_params[param_name] + shifts[i] + f_plus = expectation_fn(shifted_params) + + # - shift + shifted_params[param_name] = shifted_params[param_name] - 2 * shifts[i] + f_minus = expectation_fn(shifted_params) + shifted_params[param_name] = shifted_params[param_name] + shifts[i] + F.append((f_plus - f_minus)) + + # calculate M matrix + for j in range(n_eqs): + M[i, j] = 4 * torch.sin(shifts[i] * spectral_gaps[j] / 2) + + # get number of observables from expectation value tensor + if f_plus.numel() > 1: + batch_size = F[0].shape[0] + + F = torch.stack(F).reshape(n_eqs, -1) + R = torch.linalg.solve(M, F) + dfdx = torch.sum(spectral_gaps * R, dim=0).reshape(batch_size) + return dfdx def vjp(operation: Parametric, values: dict[str, Tensor]) -> Tensor: """Vector-jacobian product between `grad_out` and jacobians of parameters. @@ -185,23 +256,36 @@ def vjp(operation: Parametric, values: dict[str, Tensor]) -> Tensor: Returns: Updated jacobian by PSR. """ - psr_fn = ( - multi_gap_shift if len(operation.spectral_gap) > 1 else single_gap_shift + psr_fn, shift = ( + (multi_gap_shift, shift_multi) + if len(operation.spectral_gap) > 1 + else (single_gap_shift, shift_pi2) ) - return grad_out * psr_fn( # type: ignore[operator] - operation.param_name, values, operation.spectral_gap, shift + operation.param_name, # type: ignore + values, + operation.spectral_gap, + shift, ) grads = {p: None for p in ctx.param_names} for op in ctx.circuit.flatten(): - if isinstance(op, Parametric) and values[op.param_name].requires_grad: # type: ignore[index] - if grads[op.param_name] is not None: - grads[op.param_name] += vjp(op, values) - else: - grads[op.param_name] = vjp(op, values) - - return (None, None, None, None, None, *[grads[p] for p in ctx.param_names]) + if isinstance(op, Parametric) and isinstance(op.param_name, str): + if values[op.param_name].requires_grad: + if grads[op.param_name] is not None: + grads[op.param_name] += vjp(op, values) + else: + grads[op.param_name] = vjp(op, values) + + return ( + None, + None, + None, + None, + None, + None, + *[grads[p] for p in ctx.param_names], + ) def check_support_psr(circuit: QuantumCircuit): @@ -212,7 +296,6 @@ def check_support_psr(circuit: QuantumCircuit): Raises: ValueError: When circuit contains Scale, HamiltonianEvolution, - or one operation has more than two eigenvalues (multi-gap), or a param_name is used multiple times in the circuit. """ @@ -224,15 +307,15 @@ def check_support_psr(circuit: QuantumCircuit): ) if isinstance(op, Sequence): for subop in op.flatten(): + if isinstance(subop, Scale) or isinstance(subop, HamiltonianEvolution): + raise ValueError( + f"PSR is not applicable as circuit contains \ + an operation of type: {type(subop)}." + ) if isinstance(subop, Parametric): if isinstance(subop.param_name, str): - if len(subop.spectral_gap) > 1: - raise NotImplementedError("Multi-gap is not yet supported.") param_names.append(subop.param_name) - elif isinstance(op, Parametric): - if len(op.spectral_gap) > 1: - raise NotImplementedError("Multi-gap is not yet supported.") if isinstance(op.param_name, str): param_names.append(op.param_name) else: diff --git a/pyqtorch/embed.py b/pyqtorch/embed.py index 3cff1913..80f7174a 100644 --- a/pyqtorch/embed.py +++ b/pyqtorch/embed.py @@ -2,7 +2,7 @@ from importlib import import_module from logging import getLogger -from typing import Any, Optional, Tuple +from typing import Any, Tuple from numpy.typing import ArrayLike, DTypeLike @@ -21,7 +21,7 @@ "sub": ("jax.numpy", "subtract"), "div": ("jax.numpy", "divide"), } -DEFAULT_TORCH_MAPPING: dict = {} +DEFAULT_TORCH_MAPPING: dict = dict() DEFAULT_NUMPY_MAPPING = { "mul": ("numpy", "multiply"), "sub": ("numpy", "subtract"), @@ -190,7 +190,7 @@ class Embedding: ops = [rx, cry, phase, ry, cnot] n_qubits = 3 circ = pyq.QuantumCircuit(n_qubits, ops) - obs = pyq.Observable(n_qubits, [pyq.Z(0)]) + obs = pyq.Observable([pyq.Z(0)]) state = pyq.zero_state(n_qubits) @@ -208,7 +208,7 @@ def __init__( vparam_names: list[str] = [], fparam_names: list[str] = [], var_to_call: dict[str, ConcretizedCallable] = dict(), - tparam_name: Optional[str] = None, + tparam_name: str | None = None, engine_name: str = "torch", device: str = "cpu", ) -> None: diff --git a/pyqtorch/hamiltonians/__init__.py b/pyqtorch/hamiltonians/__init__.py new file mode 100644 index 00000000..32f77688 --- /dev/null +++ b/pyqtorch/hamiltonians/__init__.py @@ -0,0 +1,4 @@ +from __future__ import annotations + +from .evolution import GeneratorType, HamiltonianEvolution +from .observable import Observable diff --git a/pyqtorch/hamiltonians/evolution.py b/pyqtorch/hamiltonians/evolution.py new file mode 100644 index 00000000..daa431f8 --- /dev/null +++ b/pyqtorch/hamiltonians/evolution.py @@ -0,0 +1,385 @@ +from __future__ import annotations + +import logging +from collections import OrderedDict +from logging import getLogger +from typing import Callable, Tuple, Union + +import torch +from torch import Tensor +from torch.nn import ModuleList, ParameterDict + +from pyqtorch.apply import apply_operator +from pyqtorch.circuit import Sequence +from pyqtorch.embed import Embedding +from pyqtorch.primitives import Primitive +from pyqtorch.time_dependent.sesolve import sesolve +from pyqtorch.utils import ( + ATOL, + Operator, + SolverType, + State, + StrEnum, + expand_operator, + is_diag, +) + +BATCH_DIM = 2 +TGenerator = Union[Tensor, str, Primitive, Sequence] + +logger = getLogger(__name__) + + +def forward_hook(*args, **kwargs) -> None: # type: ignore[no-untyped-def] + logger.debug("Forward complete") + torch.cuda.nvtx.range_pop() + + +def pre_forward_hook(*args, **kwargs) -> None: # type: ignore[no-untyped-def] + logger.debug("Executing forward") + torch.cuda.nvtx.range_push("HamiltonianEvolution.forward") + + +def backward_hook(*args, **kwargs) -> None: # type: ignore[no-untyped-def] + logger.debug("Backward complete") + torch.cuda.nvtx.range_pop() + + +def pre_backward_hook(*args, **kwargs) -> None: # type: ignore[no-untyped-def] + logger.debug("Executed backward") + torch.cuda.nvtx.range_push("Hamiltonian Evolution.backward") + + +class GeneratorType(StrEnum): + """ + Options for types of generators allowed in HamiltonianEvolution. + """ + + PARAMETRIC_OPERATION = "parametric_operation" + """Generators of type Primitive or Sequence which contain + possibly trainable or non-trainable parameters.""" + OPERATION = "operation" + """Generators of type Primitive or Sequence which do not contain parameters or contain + constants as Parameters for example pyq.Scale(Z(0), torch.tensor([1.])).""" + TENSOR = "tensor" + """Generators of type torch.Tensor in which case a qubit_support needs to be passed.""" + SYMBOL = "symbol" + """Generators which are symbolic, i.e. will be passed via the 'values' dict by the user.""" + + +def is_diag_hamiltonian(hamiltonian: Operator, atol: Tensor = ATOL) -> bool: + """ + Returns True if the batched tensors H are diagonal. + + Arguments: + H: Input tensors. + atol: Tolerance for near-zero values. + + Returns: + True if diagonal, else False. + """ + diag_check = torch.tensor( + [ + is_diag(hamiltonian[..., i], atol) + for i in range(hamiltonian.shape[BATCH_DIM]) + ], + device=hamiltonian.device, + ) + return bool(torch.prod(diag_check)) + + +def evolve(hamiltonian: Operator, time_evolution: Tensor) -> Operator: + """Get the evolved operator. + + For a hamiltonian :math:`H` and a time evolution :math:`t`, returns :math:`exp(-i H, t)` + + Arguments: + hamiltonian: The operator :math:`H` for evolution. + time_evolution: The evolution time :math:`t`. + + Returns: + The evolution operator. + """ + if is_diag_hamiltonian(hamiltonian): + evol_operator = torch.diagonal(hamiltonian) * (-1j * time_evolution).view( + (-1, 1) + ) + evol_operator = torch.diag_embed(torch.exp(evol_operator)) + else: + evol_operator = torch.transpose(hamiltonian, 0, -1) * ( + -1j * time_evolution + ).view((-1, 1, 1)) + evol_operator = torch.linalg.matrix_exp(evol_operator) + return torch.transpose(evol_operator, 0, -1) + + +class HamiltonianEvolution(Sequence): + """ + The HamiltonianEvolution corresponds to :math:`t`, returns :math:`exp(-i H, t)` where + a hamiltonian/generator :math:`H` and a time evolution :math:`t` are given. + + Note that the quantity :math:`H.t` is considered dimensionless. + + We can create such operation by passing different generator types: + - A tensor representation of the generator, + - A string when we consider the generator as a symbol. + - Operations as a single primitive or a sequence, parameterized or not. + + Attributes: + generator_type: This attribute informs on how the generator is inputed + and sets the logic for applying hamiltonian evolution. + time: The evolution time :math:`t`. + operations: List of operations. + cache_length: LRU cache cache_length evolution operators for given set + of parameter values. + """ + + def __init__( + self, + generator: TGenerator, + time: Tensor | str, + qubit_support: Tuple[int, ...] | None = None, + generator_parametric: bool = False, + cache_length: int = 1, + is_time_dependent: bool = False, + duration: int = 100, + steps: int = 10, + solver=SolverType.KRYLOV_SE, + ): + """Initializes the HamiltonianEvolution. + Depending on the generator argument, set the type and set the right generator getter. + + Arguments: + generator: The generator :math:`H`. + time: The evolution time :math:`t`. + qubit_support: The qubits the operator acts on. + generator_parametric: Whether the generator is parametric or not. + """ + + self.is_time_dependent = is_time_dependent + self.solver_type = solver + self.duration = duration + self.steps = steps + + if isinstance(generator, Tensor): + if qubit_support is None: + raise ValueError( + "When using a Tensor generator, please pass a qubit_support." + ) + if len(generator.shape) < 3: + generator = generator.unsqueeze(2) + generator = [Primitive(generator, qubit_support)] + self.generator_type = GeneratorType.TENSOR + + elif isinstance(generator, str): + if qubit_support is None: + raise ValueError( + "When using a symbolic generator, please pass a qubit_support." + ) + self.generator_type = GeneratorType.SYMBOL + self.generator_symbol = generator + generator = [] + elif isinstance(generator, (Primitive, Sequence)): + if qubit_support is not None: + logger.warning( + "Taking support from generator and ignoring qubit_support input." + ) + qubit_support = generator.qubit_support + if generator_parametric: + generator = [generator] + self.generator_type = GeneratorType.PARAMETRIC_OPERATION + else: + generator = [ + Primitive( + generator.tensor(), + generator.qubit_support, + ) + ] + + self.generator_type = GeneratorType.OPERATION + else: + raise TypeError( + f"Received generator of type {type(generator)},\ + allowed types are: [Tensor, str, Primitive, Sequence]" + ) + super().__init__(generator) + self._qubit_support = qubit_support # type: ignore + self.time = time + logger.debug("Hamiltonian Evolution initialized") + if logger.isEnabledFor(logging.DEBUG): + # When Debugging let's add logging and NVTX markers + # WARNING: incurs performance penalty + self.register_forward_hook(forward_hook, always_call=True) + self.register_full_backward_hook(backward_hook) + self.register_forward_pre_hook(pre_forward_hook) + self.register_full_backward_pre_hook(pre_backward_hook) + + self._generator_map: dict[GeneratorType, Callable] = { + GeneratorType.SYMBOL: self._symbolic_generator, + GeneratorType.TENSOR: self._tensor_generator, + GeneratorType.OPERATION: self._tensor_generator, + GeneratorType.PARAMETRIC_OPERATION: self._tensor_generator, + } + + # to avoid recomputing hamiltonians and evolution + self._cache_hamiltonian_evo: dict[str, Tensor] = dict() + self.cache_length = cache_length + + @property + def generator(self) -> ModuleList: + """Returns the operations making the generator. + + Returns: + The generator as a ModuleList. + """ + return self.operations + + def _symbolic_generator( + self, + values: dict, + embedding: Embedding | None = None, + ) -> Operator: + """Returns the generator for the SYMBOL case. + + Arguments: + values: + + Returns: + The generator as a tensor. + """ + hamiltonian = values[self.generator_symbol] + # add batch dim + if len(hamiltonian.shape) == 2: + return hamiltonian.unsqueeze(2) + # cases when the batchdim is at index 0 instead of 2 + if len(hamiltonian.shape) == 3 and ( + hamiltonian.shape[0] != hamiltonian.shape[1] + ): + return torch.transpose(hamiltonian, 0, 2) + if len(hamiltonian.shape) == 4 and ( + hamiltonian.shape[0] != hamiltonian.shape[1] + ): + return torch.permute(hamiltonian.squeeze(3), (1, 2, 0)) + return hamiltonian + + def _tensor_generator( + self, values: dict = dict(), embedding: Embedding | None = None + ) -> Operator: + """Returns the generator for the TENSOR, OPERATION and PARAMETRIC_OPERATION cases. + + Arguments: + values: Values dict with any needed parameters. + + Returns: + The generator as a tensor. + """ + return self.generator[0].tensor(values, embedding) + + @property + def create_hamiltonian(self) -> Callable[[dict], Operator]: + """A utility method for setting the right generator getter depending on the init case. + + Returns: + The right generator getter. + """ + return self._generator_map[self.generator_type] + + def forward( + self, + state: Tensor, + values: dict[str, Tensor] | ParameterDict = dict(), + embedding: Embedding | None = None, + ) -> State: + """ + Apply the hamiltonian evolution with input parameter values. + + Arguments: + state: Input state. + values: Values of parameters. + + Returns: + The transformed state. + """ + if self.is_time_dependent and embedding is not None: + n_qubits = len(state.shape) - 1 + batch_size = state.shape[-1] + t_grid = torch.linspace(0, float(self.duration), self.steps) + + def Ht(t: torch.Tensor) -> torch.Tensor: + """Accepts a value 't' for time and returns + a (2**n_qubits, 2**n_qubits) Hamiltonian evaluated at time 't'. + """ + # We use the origial embedded params and return a new dict + # where we reembedded all parameters depending on time with value 't' + reembedded_time_values = embedding.reembed_tparam(values, t) + return ( + self.generator[0] + .tensor(reembedded_time_values, n_qubits) + .squeeze(2) + ) # squeeze the batch dim and return shape (2**n_qubits, 2**n_qubits) + + sol = sesolve( + Ht, # returns a tensor of shape (2**n_qubits, 2**n_qubits) + torch.flatten( + state, start_dim=0, end_dim=-2 + ), # reshape to (2**n_qubits, batch_size) + t_grid, + self.solver_name, + ) + + # Retrieve the last state of shape (2**n_qubits, batch_size) + state = sol.states[-1] + + return state.reshape( + [2] * n_qubits + [batch_size] + ) # reshape to [2] * n_qubits + [batch_size] + + else: + evolved_op = self.tensor(values, embedding) + return apply_operator( + state=state, + operator=evolved_op, + qubits=self.qubit_support, + n_qubits=len(state.size()) - 1, + batch_size=evolved_op.shape[BATCH_DIM], + ) + + def tensor( + self, + values: dict = dict(), + embedding: Embedding | None = None, + full_support: tuple[int, ...] | None = None, + ) -> Operator: + """Get the corresponding unitary over n_qubits. + + To avoid computing the evolution operator, we store it in cache wrt values. + + Arguments: + values: Parameter value. + Can be higher than the number of qubit support. + + Returns: + The unitary representation. + """ + values_cache_key = str(OrderedDict(values)) + if self.cache_length > 0 and values_cache_key in self._cache_hamiltonian_evo: + evolved_op = self._cache_hamiltonian_evo[values_cache_key] + else: + hamiltonian: torch.Tensor = self.create_hamiltonian(values, embedding) # type: ignore [call-arg] + time_evolution: torch.Tensor = ( + values[self.time] if isinstance(self.time, str) else self.time + ) # If `self.time` is a string / hence, a Parameter, + # we expect the user to pass it in the `values` dict + evolved_op = evolve(hamiltonian, time_evolution) + nb_cached = len(self._cache_hamiltonian_evo) + + # LRU caching + if (nb_cached > 0) and (nb_cached == self.cache_length): + self._cache_hamiltonian_evo.pop(next(iter(self._cache_hamiltonian_evo))) + if nb_cached < self.cache_length: + self._cache_hamiltonian_evo[values_cache_key] = evolved_op + + if full_support is None: + return evolved_op + else: + return expand_operator(evolved_op, self.qubit_support, full_support) diff --git a/pyqtorch/hamiltonians/observable.py b/pyqtorch/hamiltonians/observable.py new file mode 100644 index 00000000..61edd454 --- /dev/null +++ b/pyqtorch/hamiltonians/observable.py @@ -0,0 +1,49 @@ +from __future__ import annotations + +from torch import Tensor +from torch.nn import Module, ParameterDict + +from pyqtorch.circuit import Sequence +from pyqtorch.composite import Add +from pyqtorch.embed import Embedding +from pyqtorch.primitives import Primitive +from pyqtorch.utils import ( + inner_prod, +) + + +class Observable(Add): + """ + The Observable :math:`O` represents an operator from which + we can extract expectation values from quantum states. + + Given an input state :math:`\\ket\\rangle`, the expectation value with :math:`O` is defined as + :math:`\\langle\\bra|O\\ket\\rangle` + + Attributes: + operations: List of operations. + n_qubits: Number of qubits it is defined on. + """ + + def __init__( + self, + operations: list[Module] | Primitive | Sequence, + ): + super().__init__(operations if isinstance(operations, list) else [operations]) + + def expectation( + self, + state: Tensor, + values: dict[str, Tensor] | ParameterDict = dict(), + embedding: Embedding | None = None, + ) -> Tensor: + """Calculate the inner product :math:`\\langle\\bra|O\\ket\\rangle` + + Arguments: + state: Input state. + values: Values of parameters. + + Returns: + The expectation value. + """ + return inner_prod(state, self.forward(state, values, embedding)).real diff --git a/pyqtorch/matrices.py b/pyqtorch/matrices.py index 4fed91c5..d5e94c97 100644 --- a/pyqtorch/matrices.py +++ b/pyqtorch/matrices.py @@ -66,7 +66,7 @@ def PROJMAT(ket: Tensor, bra: Tensor) -> Tensor: } -def _unitary( +def parametric_unitary( theta: torch.Tensor, P: torch.Tensor, I: torch.Tensor, batch_size: int # noqa: E741 ) -> torch.Tensor: cos_t = torch.cos(theta / 2).unsqueeze(0).unsqueeze(1) @@ -98,17 +98,19 @@ def _jacobian( return -1 / 2 * (sin_t * batch_imat + 1j * cos_t * batch_operation_mat) -def _controlled( - unitary: torch.Tensor, batch_size: int, n_control_qubits: int = 1 +def controlled( + operation: torch.Tensor, batch_size: int, n_control_qubits: int = 1 ) -> torch.Tensor: _controlled: torch.Tensor = ( - torch.eye(2 ** (n_control_qubits + 1), dtype=unitary.dtype) + torch.eye( + 2 ** (n_control_qubits + 1), dtype=operation.dtype, device=operation.device + ) .unsqueeze(2) .repeat(1, 1, batch_size) ) _controlled[ 2 ** (n_control_qubits + 1) - 2 :, 2 ** (n_control_qubits + 1) - 2 :, : - ] = unitary + ] = operation return _controlled diff --git a/pyqtorch/noise/__init__.py b/pyqtorch/noise/__init__.py new file mode 100644 index 00000000..aa41334c --- /dev/null +++ b/pyqtorch/noise/__init__.py @@ -0,0 +1,13 @@ +from __future__ import annotations + +from .gates import ( + AmplitudeDamping, + BitFlip, + Depolarizing, + GeneralizedAmplitudeDamping, + Noise, + PauliChannel, + PhaseDamping, + PhaseFlip, +) +from .protocol import NoiseProtocol, _repr_noise diff --git a/pyqtorch/noise.py b/pyqtorch/noise/gates.py similarity index 63% rename from pyqtorch/noise.py rename to pyqtorch/noise/gates.py index 9b048a58..4cb87a57 100644 --- a/pyqtorch/noise.py +++ b/pyqtorch/noise/gates.py @@ -1,44 +1,41 @@ from __future__ import annotations -from math import sqrt +from math import log2, sqrt +from typing import Any import torch from torch import Tensor -from pyqtorch.apply import operator_product +from pyqtorch.apply import apply_density_mat from pyqtorch.embed import Embedding -from pyqtorch.matrices import DEFAULT_MATRIX_DTYPE, IMAT, XMAT, YMAT, ZMAT, _dagger -from pyqtorch.utils import DensityMatrix, density_mat +from pyqtorch.matrices import DEFAULT_MATRIX_DTYPE, IMAT, XMAT, YMAT, ZMAT +from pyqtorch.utils import DensityMatrix, density_mat, qubit_support_as_tuple class Noise(torch.nn.Module): def __init__( - self, kraus: list[Tensor], target: int, probabilities: tuple[float, ...] | float + self, + kraus: list[Tensor], + target: int, + error_probabilities: tuple[float, ...] | float, ) -> None: super().__init__() self.target: int = target - self.qubit_support: tuple[int, ...] = (self.target,) + self.qubit_support: tuple[int, ...] = qubit_support_as_tuple(self.target) for index, tensor in enumerate(kraus): self.register_buffer(f"kraus_{index}", tensor) self._device: torch.device = kraus[0].device - self.probabilities: tuple[float, ...] | float = probabilities - - def __eq__(self, other: object) -> bool: - if isinstance(other, type(self)): - return ( - self.__class__.__name__ == other.__class__.__name__ - and self.probabilities == other.probabilities - ) - return False + self._dtype: torch.dtype = kraus[0].dtype + self.error_probabilities: tuple[float, ...] | float = error_probabilities def extra_repr(self) -> str: - return f"'qubit_support'={self.qubit_support}, 'probabilities'={self.probabilities}" + return f"target: {self.qubit_support}, prob: {self.probabilities}" @property def kraus_operators(self) -> list[Tensor]: return [getattr(self, f"kraus_{i}") for i in range(len(self._buffers))] - def unitary( + def _tensor( self, values: dict[str, Tensor] | Tensor = dict(), embedding: Embedding | None = None, @@ -46,9 +43,6 @@ def unitary( # Since PyQ expects tensor.Size = [2**n_qubits, 2**n_qubits,batch_size]. return [kraus_op.unsqueeze(2) for kraus_op in self.kraus_operators] - def dagger(self, values: dict[str, Tensor] | Tensor = dict()) -> list[Tensor]: - return [_dagger(kraus_op) for kraus_op in self.unitary(values)] - def forward( self, state: Tensor, @@ -75,15 +69,10 @@ def forward( """ if not isinstance(state, DensityMatrix): state = density_mat(state) + n_qubits = int(log2(state.size(1))) rho_evols: list[Tensor] = [] - for kraus_unitary, kraus_dagger in zip( - self.unitary(values), self.dagger(values) - ): - rho_evol: Tensor = operator_product( - kraus_unitary, - operator_product(state, kraus_dagger, self.target), - self.target, - ) + for kraus in self.tensor(values, n_qubits): + rho_evol: Tensor = apply_density_mat(kraus, state) rho_evols.append(rho_evol) rho_final: Tensor = torch.stack(rho_evols, dim=0) rho_final = torch.sum(rho_final, dim=0) @@ -93,11 +82,42 @@ def forward( def device(self) -> torch.device: return self._device - def to(self, device: torch.device) -> Noise: - super().to(device) - self._device = device + @property + def dtype(self) -> torch.dtype: + return self._dtype + + def to(self, *args: Any, **kwargs: Any) -> Noise: + super().to(*args, **kwargs) + self._device = self.kraus_0.device + self._dtype = self.kraus_0.dtype return self + def tensor( + self, values: dict[str, Tensor] = dict(), n_qubits: int = 1 + ) -> list[Tensor]: + block_mats = self._tensor(values) + mats = [] + for blockmat in block_mats: + if n_qubits == 1: + mats.append(blockmat) + else: + full_sup = tuple(i for i in range(n_qubits)) + support = tuple(sorted(self.qubit_support)) + mat = ( + IMAT.clone().to(self.device, self.dtype).unsqueeze(2) + if support[0] != full_sup[0] + else blockmat + ) + for i in full_sup[1:]: + if i == support[0]: + other = blockmat + mat = torch.kron(mat.contiguous(), other.contiguous()) + elif i not in support: + other = IMAT.clone().to(self.device, self.dtype).unsqueeze(2) + mat = torch.kron(mat.contiguous(), other.contiguous()) + mats.append(mat) + return mats + class BitFlip(Noise): """ @@ -110,23 +130,23 @@ class BitFlip(Noise): Args: target (int): The index of the qubit being affected by the noise. - probability (float): The probability of a bit flip error. + error_probability (float): The probability of a bit flip error. Raises: - TypeError: If the probability value is not a float. + TypeError: If the error_probability value is not a float. """ def __init__( self, target: int, - probability: float, + error_probability: float, ): - if probability > 1.0 or probability < 0.0: - raise ValueError("The probability value is not a correct probability") - K0: Tensor = sqrt(1.0 - probability) * IMAT - K1: Tensor = sqrt(probability) * XMAT + if error_probability > 1.0 or error_probability < 0.0: + raise ValueError("The error_probability value is not a correct probability") + K0: Tensor = sqrt(1.0 - error_probability) * IMAT + K1: Tensor = sqrt(error_probability) * XMAT kraus_bitflip: list[Tensor] = [K0, K1] - super().__init__(kraus_bitflip, target, probability) + super().__init__(kraus_bitflip, target, error_probability) class PhaseFlip(Noise): @@ -140,23 +160,23 @@ class PhaseFlip(Noise): Args: target (int): The index of the qubit being affected by the noise. - probability (float): The probability of phase flip error. + error_probability (float): The probability of phase flip error. Raises: - TypeError: If the probability value is not a float. + TypeError: If the error_probability value is not a float. """ def __init__( self, target: int, - probability: float, + error_probability: float, ): - if probability > 1.0 or probability < 0.0: - raise ValueError("The probability value is not a correct probability") - K0: Tensor = sqrt(1.0 - probability) * IMAT - K1: Tensor = sqrt(probability) * ZMAT + if error_probability > 1.0 or error_probability < 0.0: + raise ValueError("The error_probability value is not a correct probability") + K0: Tensor = sqrt(1.0 - error_probability) * IMAT + K1: Tensor = sqrt(error_probability) * ZMAT kraus_phaseflip: list[Tensor] = [K0, K1] - super().__init__(kraus_phaseflip, target, probability) + super().__init__(kraus_phaseflip, target, error_probability) class Depolarizing(Noise): @@ -173,25 +193,25 @@ class Depolarizing(Noise): Args: target (int): The index of the qubit being affected by the noise. - probability (float): The probability of depolarizing error. + error_probability (float): The probability of depolarizing error. Raises: - TypeError: If the probability value is not a float. + TypeError: If the error_probability value is not a float. """ def __init__( self, target: int, - probability: float, + error_probability: float, ): - if probability > 1.0 or probability < 0.0: - raise ValueError("The probability value is not a correct probability") - K0: Tensor = sqrt(1.0 - probability) * IMAT - K1: Tensor = sqrt(probability / 3) * XMAT - K2: Tensor = sqrt(probability / 3) * YMAT - K3: Tensor = sqrt(probability / 3) * ZMAT + if error_probability > 1.0 or error_probability < 0.0: + raise ValueError("The error_probability value is not a correct probability") + K0: Tensor = sqrt(1.0 - error_probability) * IMAT + K1: Tensor = sqrt(error_probability / 3) * XMAT + K2: Tensor = sqrt(error_probability / 3) * YMAT + K3: Tensor = sqrt(error_probability / 3) * ZMAT kraus_depolarizing: list[Tensor] = [K0, K1, K2, K3] - super().__init__(kraus_depolarizing, target, probability) + super().__init__(kraus_depolarizing, target, error_probability) class PauliChannel(Noise): @@ -208,7 +228,8 @@ class PauliChannel(Noise): Args: target (int): The index of the qubit being affected by the noise. - probabilities (Tuple[float, ...]): Tuple containing probabilities of X, Y, and Z errors. + error_probability (Tuple[float, ...]): Tuple containing probabilities + of X, Y, and Z errors. Raises: TypeError: If the probabilities values are not provided in a tuple. @@ -217,21 +238,25 @@ class PauliChannel(Noise): def __init__( self, target: int, - probabilities: tuple[float, ...], + error_probability: tuple[float, ...], ): - sum_prob = sum(probabilities) + sum_prob = sum(error_probability) if sum_prob > 1.0: raise ValueError("The sum of probabilities can't be greater than 1.0") - for probability in probabilities: + for probability in error_probability: if probability > 1.0 or probability < 0.0: raise ValueError("The probability values are not correct probabilities") - px, py, pz = probabilities[0], probabilities[1], probabilities[2] + px, py, pz = ( + error_probability[0], + error_probability[1], + error_probability[2], + ) K0: Tensor = sqrt(1.0 - (px + py + pz)) * IMAT K1: Tensor = sqrt(px) * XMAT K2: Tensor = sqrt(py) * YMAT K3: Tensor = sqrt(pz) * ZMAT kraus_pauli_channel: list[Tensor] = [K0, K1, K2, K3] - super().__init__(kraus_pauli_channel, target, probabilities) + super().__init__(kraus_pauli_channel, target, error_probability) class AmplitudeDamping(Noise): @@ -252,7 +277,7 @@ class AmplitudeDamping(Noise): Args: target (int): The index of the qubit being affected by the noise. - rate (float): The damping rate, indicating the probability of amplitude loss. + error_probability (float): The damping rate, indicating the probability of amplitude loss. Raises: TypeError: If the rate value is not a float. @@ -262,8 +287,9 @@ class AmplitudeDamping(Noise): def __init__( self, target: int, - rate: float, + error_probability: float, ): + rate = error_probability if rate > 1.0 or rate < 0.0: raise ValueError("The damping rate is not a correct probability") K0: Tensor = torch.tensor( @@ -292,7 +318,7 @@ class PhaseDamping(Noise): Args: target (int): The index of the qubit being affected by the noise. - rate (float): The damping rate, indicating the probability of phase damping. + error_probability (float): The damping rate, indicating the probability of phase damping. Raises: TypeError: If the rate value is not a float. @@ -302,8 +328,9 @@ class PhaseDamping(Noise): def __init__( self, target: int, - rate: float, + error_probability: float, ): + rate = error_probability if rate > 1.0 or rate < 0.0: raise ValueError("The damping rate is not a correct probability") K0: Tensor = torch.tensor( @@ -335,8 +362,9 @@ class GeneralizedAmplitudeDamping(Noise): Args: target (int): The index of the qubit being affected by the noise. - probability (float): The probability of amplitude damping error. - rate (float): The damping rate, indicating the probability of generalized amplitude damping. + error_probabilities (tuple[float,...]): The first float must be the probability + of amplitude damping error, and the second float is the damping rate, indicating + the probability of generalized amplitude damping. Raises: TypeError: If the probability or rate values is not float. @@ -346,9 +374,10 @@ class GeneralizedAmplitudeDamping(Noise): def __init__( self, target: int, - probability: float, - rate: float, + error_probability: tuple[float, ...], ): + probability = error_probability[0] + rate = error_probability[1] if probability > 1.0 or probability < 0.0: raise ValueError("The probability value is not a correct probability") if rate > 1.0 or rate < 0.0: @@ -366,4 +395,4 @@ def __init__( [[0, 0], [sqrt(rate), 0]], dtype=DEFAULT_MATRIX_DTYPE ) kraus_generalized_amplitude_damping: list[Tensor] = [K0, K1, K2, K3] - super().__init__(kraus_generalized_amplitude_damping, target, probability) + super().__init__(kraus_generalized_amplitude_damping, target, error_probability) diff --git a/pyqtorch/noise/protocol.py b/pyqtorch/noise/protocol.py new file mode 100644 index 00000000..a3c3782c --- /dev/null +++ b/pyqtorch/noise/protocol.py @@ -0,0 +1,54 @@ +from __future__ import annotations + +import sys + + +class NoiseProtocol: + BITFLIP = "BitFlip" + PHASEFLIP = "PhaseFlip" + DEPOLARIZING = "Depolarizing" + PAULI_CHANNEL = "PauliChannel" + AMPLITUDE_DAMPING = "AmplitudeDamping" + PHASE_DAMPING = "PhaseDamping" + GENERALIZED_AMPLITUDE_DAMPING = "GeneralizedAmplitudeDamping" + + def __init__(self, protocol: str, options: dict = dict()) -> None: + self.protocol: str = protocol + self.options: dict = options + + def __repr__(self) -> str: + if self.target: + return ( + f"{self.protocol}(prob: {self.error_probability}, " + f"target: {self.target})" + ) + return f"{self.protocol}(prob: {self.error_probability})" + + @property + def error_probability(self): + return self.options.get("error_probability") + + @property + def target(self): #! init_state not good size + return self.options.get("target") + + def protocol_to_gate(self): + try: + gate_class = getattr(sys.modules["pyqtorch.noise.gates"], self.protocol) + return gate_class + except AttributeError: + raise ValueError( + f"The protocol {self.protocol} has not been implemented in pyq yet." + ) + + +def _repr_noise(noise: NoiseProtocol | dict[str, NoiseProtocol] | None = None) -> str: + """Returns the string for noise representation in gates.""" + noise_info = "" + if noise is None: + return noise_info + elif isinstance(noise, NoiseProtocol): + noise_info = str(noise) + elif isinstance(noise, dict): + noise_info = ", ".join(str(noise_instance) for noise_instance in noise.values()) + return f", noise: {noise_info}" diff --git a/pyqtorch/primitive.py b/pyqtorch/primitive.py deleted file mode 100644 index b9758083..00000000 --- a/pyqtorch/primitive.py +++ /dev/null @@ -1,303 +0,0 @@ -from __future__ import annotations - -import logging -from functools import cached_property -from logging import getLogger -from typing import Any - -import numpy as np -import torch -from torch import Tensor - -from pyqtorch.apply import apply_operator, operator_product -from pyqtorch.embed import Embedding -from pyqtorch.matrices import ( - IMAT, - OPERATIONS_DICT, - _controlled, - _dagger, -) -from pyqtorch.utils import DensityMatrix, product_state - -logger = getLogger(__name__) - - -def forward_hook(*args, **kwargs) -> None: # type: ignore[no-untyped-def] - torch.cuda.nvtx.range_pop() - - -def pre_forward_hook(*args, **kwargs) -> None: # type: ignore[no-untyped-def] - torch.cuda.nvtx.range_push("Primitive.forward") - - -def backward_hook(*args, **kwargs) -> None: # type: ignore[no-untyped-def] - torch.cuda.nvtx.range_pop() - - -def pre_backward_hook(*args, **kwargs) -> None: # type: ignore[no-untyped-def] - torch.cuda.nvtx.range_push("Primitive.backward") - - -class Primitive(torch.nn.Module): - def __init__( - self, - pauli: Tensor, - target: int | tuple[int, ...], - pauli_generator: Tensor | None = None, - ) -> None: - super().__init__() - self.target: int | tuple[int, ...] = target - - self.qubit_support: tuple[int, ...] = ( - (target,) if isinstance(target, int) else target - ) - if isinstance(target, np.integer): - self.qubit_support = (target.item(),) - self.register_buffer("pauli", pauli) - self.pauli_generator = pauli_generator - self._device = self.pauli.device - self._dtype = self.pauli.dtype - - if logger.isEnabledFor(logging.DEBUG): - # When Debugging let's add logging and NVTX markers - # WARNING: incurs performance penalty - self.register_forward_hook(forward_hook, always_call=True) - self.register_full_backward_hook(backward_hook) - self.register_forward_pre_hook(pre_forward_hook) - self.register_full_backward_pre_hook(pre_backward_hook) - - def __hash__(self) -> int: - return hash(self.qubit_support) - - def extra_repr(self) -> str: - return f"{self.qubit_support}" - - def unitary( - self, - values: dict[str, Tensor] | Tensor = dict(), - embedding: Embedding | None = None, - ) -> Tensor: - return self.pauli.unsqueeze(2) if len(self.pauli.shape) == 2 else self.pauli - - def forward( - self, - state: Tensor, - values: dict[str, Tensor] | Tensor = dict(), - embedding: Embedding | None = None, - ) -> Tensor: - if isinstance(state, DensityMatrix): - # TODO: fix error type int | tuple[int, ...] expected "int" - # Only supports single-qubit gates - return DensityMatrix( - operator_product( - self.unitary(values, embedding), - operator_product(state, self.dagger(values), self.target), # type: ignore [arg-type] - self.target, # type: ignore [arg-type] - ) - ) - else: - return apply_operator( - state, - self.unitary(values, embedding), - self.qubit_support, - len(state.size()) - 1, - ) - - def dagger( - self, - values: dict[str, Tensor] | Tensor = dict(), - embedding: Embedding | None = None, - ) -> Tensor: - return _dagger(self.unitary(values, embedding)) - - @property - def device(self) -> torch.device: - return self._device - - @property - def dtype(self) -> torch.dtype: - return self._dtype - - def to(self, *args: Any, **kwargs: Any) -> Primitive: - super().to(*args, **kwargs) - self._device = self.pauli.device - self._dtype = self.pauli.dtype - return self - - @cached_property - def eigenvals_generator(self) -> Tensor: - """Get eigenvalues of the underlying generator. - - Note that for a primitive, the generator is unclear - so we execute pass. - - Arguments: - values: Parameter values. - - Returns: - Eigenvalues of the generator operator. - """ - if self.pauli_generator is not None: - return torch.linalg.eigvalsh(self.pauli_generator).reshape(-1, 1) - pass - - @cached_property - def spectral_gap(self) -> Tensor: - """Difference between the moduli of the two largest eigenvalues of the generator. - - Returns: - Tensor: Spectral gap value. - """ - spectrum = self.eigenvals_generator - spectral_gap = torch.unique(torch.abs(torch.tril(spectrum - spectrum.T))) - return spectral_gap[spectral_gap.nonzero()] - - def tensor( - self, values: dict[str, Tensor] = {}, n_qubits: int = 1, diagonal: bool = False - ) -> Tensor: - if diagonal: - raise NotImplementedError - blockmat = self.unitary(values) - if n_qubits == 1: - return blockmat - full_sup = tuple(i for i in range(n_qubits)) - support = tuple(sorted(self.qubit_support)) - mat = ( - IMAT.clone().to(self.device).unsqueeze(2) - if support[0] != full_sup[0] - else blockmat - ) - for i in full_sup[1:]: - if i == support[0]: - other = blockmat - mat = torch.kron(mat.contiguous(), other.contiguous()) - elif i not in support: - other = IMAT.clone().to(self.device).unsqueeze(2) - mat = torch.kron(mat.contiguous(), other.contiguous()) - return mat - - -class X(Primitive): - def __init__(self, target: int): - super().__init__(OPERATIONS_DICT["X"], target) - - -class Y(Primitive): - def __init__(self, target: int): - super().__init__(OPERATIONS_DICT["Y"], target) - - -class Z(Primitive): - def __init__(self, target: int): - super().__init__(OPERATIONS_DICT["Z"], target) - - -class I(Primitive): # noqa: E742 - def __init__(self, target: int): - super().__init__(OPERATIONS_DICT["I"], target) - - def forward( - self, - state: Tensor, - values: dict[str, Tensor] = dict(), - embedding: Embedding | None = None, - ) -> Tensor: - return state - - -class H(Primitive): - def __init__(self, target: int): - super().__init__(OPERATIONS_DICT["H"], target) - - -class T(Primitive): - def __init__(self, target: int): - super().__init__(OPERATIONS_DICT["T"], target) - - -class S(Primitive): - def __init__(self, target: int): - super().__init__(OPERATIONS_DICT["S"], target, 0.5 * OPERATIONS_DICT["Z"]) - - -class SDagger(Primitive): - def __init__(self, target: int): - super().__init__( - OPERATIONS_DICT["SDAGGER"], target, -0.5 * OPERATIONS_DICT["Z"] - ) - - -class Projector(Primitive): - def __init__(self, qubit_support: int | tuple[int, ...], ket: str, bra: str): - support = (qubit_support,) if isinstance(qubit_support, int) else qubit_support - if len(ket) != len(bra): - raise ValueError("Input ket and bra bitstrings must be of same length.") - ket_state = product_state(ket).flatten() - bra_state = product_state(bra).flatten() - super().__init__(OPERATIONS_DICT["PROJ"](ket_state, bra_state), support[-1]) - # Override the attribute in AbstractOperator. - self.qubit_support = support - - -class N(Primitive): - def __init__(self, target: int): - super().__init__(OPERATIONS_DICT["N"], target) - - -class SWAP(Primitive): - def __init__(self, control: int, target: int): - super().__init__(OPERATIONS_DICT["SWAP"], target) - self.control = (control,) if isinstance(control, int) else control - self.qubit_support = self.control + (target,) - - -class CSWAP(Primitive): - def __init__(self, control: int | tuple[int, ...], target: tuple[int, ...]): - if not isinstance(target, tuple) or len(target) != 2: - raise ValueError("Target qubits must be a tuple with two qubits") - super().__init__(OPERATIONS_DICT["CSWAP"], target) - self.control = (control,) if isinstance(control, int) else control - self.target = target - self.qubit_support = self.control + self.target - - def extra_repr(self) -> str: - return f"control:{self.control}, target:{self.target}" - - -class ControlledOperationGate(Primitive): - def __init__(self, gate: str, control: int | tuple[int, ...], target: int): - self.control = (control,) if isinstance(control, int) else control - mat = OPERATIONS_DICT[gate] - mat = _controlled( - unitary=mat.unsqueeze(2), - batch_size=1, - n_control_qubits=len(self.control), - ).squeeze(2) - super().__init__(mat, target) - self.qubit_support = self.control + (self.target,) # type: ignore[operator] - - def extra_repr(self) -> str: - return f"control:{self.control}, target:{(self.target,)}" - - -class CNOT(ControlledOperationGate): - def __init__(self, control: int | tuple[int, ...], target: int): - super().__init__("X", control, target) - - -CX = CNOT - - -class CY(ControlledOperationGate): - def __init__(self, control: int | tuple[int, ...], target: int): - super().__init__("Y", control, target) - - -class CZ(ControlledOperationGate): - def __init__(self, control: int | tuple[int, ...], target: int): - super().__init__("Z", control, target) - - -class Toffoli(ControlledOperationGate): - def __init__(self, control: int | tuple[int, ...], target: int): - super().__init__("X", control, target) diff --git a/pyqtorch/primitives/__init__.py b/pyqtorch/primitives/__init__.py new file mode 100644 index 00000000..97ef5841 --- /dev/null +++ b/pyqtorch/primitives/__init__.py @@ -0,0 +1,41 @@ +from __future__ import annotations + +from .parametric import ControlledParametric, ControlledRotationGate, Parametric +from .parametric_gates import ( + CPHASE, + CRX, + CRY, + CRZ, + OPS_PARAM, + OPS_PARAM_1Q, + OPS_PARAM_2Q, + PHASE, + RX, + RY, + RZ, + U, +) +from .primitive import ControlledPrimitive, Primitive +from .primitive_gates import ( + CNOT, + CSWAP, + CY, + CZ, + OPS_1Q, + OPS_2Q, + OPS_3Q, + OPS_DIGITAL, + OPS_PAULI, + SWAP, + H, + I, + N, + Projector, + S, + SDagger, + T, + Toffoli, + X, + Y, + Z, +) diff --git a/pyqtorch/primitives/parametric.py b/pyqtorch/primitives/parametric.py new file mode 100644 index 00000000..157f518f --- /dev/null +++ b/pyqtorch/primitives/parametric.py @@ -0,0 +1,333 @@ +from __future__ import annotations + +from functools import cached_property +from typing import Any, Tuple + +import torch +from torch import Tensor + +from pyqtorch.embed import Embedding +from pyqtorch.matrices import ( + OPERATIONS_DICT, + _jacobian, + controlled, + parametric_unitary, +) +from pyqtorch.noise import NoiseProtocol, _repr_noise +from pyqtorch.quantum_operation import QuantumOperation, Support + +pauli_singleq_eigenvalues = torch.tensor([[-1.0], [1.0]], dtype=torch.cdouble) + + +class Parametric(QuantumOperation): + """ + QuantumOperation taking parameters as input. + + Attributes: + param_name: Name of parameters. + parse_values: Method defining how to handle the values dictionary input. + """ + + n_params = 1 + + def __init__( + self, + generator: str | Tensor, + qubit_support: int | tuple[int, ...] | Support, + param_name: str | int | float | torch.Tensor = "", + noise: NoiseProtocol | dict[str, NoiseProtocol] | None = None, + ): + """Initializes Parametric. + + Arguments: + generator: Generator to use. + qubit_support: Qubits to act on. + param_name: Name of parameters. + noise: Optional noise protocols to apply. + """ + + generator_operation = ( + OPERATIONS_DICT[generator] if isinstance(generator, str) else generator + ) + self.param_name = param_name + + def parse_values( + values: dict[str, Tensor] | Tensor = dict(), + embedding: Embedding | None = None, + ) -> Tensor: + """The legacy way of using parametric gates: + The Parametric gate received a string as a 'param_name' and performs a + a lookup in the passed `values` dict for to retrieve the torch.Tensor passed + under the key `param_name`. + + Arguments: + values: A dict containing param_name:torch.Tensor pairs + embedding: An optional embedding. + Returns: + A Torch Tensor denoting values for the `param_name`. + """ + return Parametric._expand_values(values[self.param_name]) # type: ignore[index] + + def parse_tensor( + values: dict[str, Tensor] | Tensor = dict(), + embedding: Embedding | None = None, + ) -> Tensor: + """Functional version of the Parametric gate: + In case the user did not pass a `param_name`, + pyqtorch assumes `values` will be a torch.Tensor instead of a dict. + + Arguments: + values: A dict containing param_name:torch.Tensor pairs + Returns: + A Torch Tensor with which to evaluate the Parametric Gate. + """ + # self.param_name will be "" + return Parametric._expand_values(values) + + def parse_constant( + values: dict[str, Tensor] | Tensor = dict(), + embedding: Embedding | None = None, + ) -> Tensor: + """Fix a the parameter of a Parametric Gate to a numeric constant + if the user passed a numeric input for the `param_name`. + + Arguments: + values: A dict containing param_name:torch.Tensor pairs + Returns: + A Torch Tensor with which to evaluate the Parametric Gate. + """ + # self.param_name will be a torch.Tensor + return Parametric._expand_values( + torch.tensor(self.param_name, device=self.device, dtype=self.dtype) + ) + + if param_name == "": + self.parse_values = parse_tensor + self.param_name = self.param_name + elif isinstance(param_name, str): + self.parse_values = parse_values + elif isinstance(param_name, (float, int, torch.Tensor)): + self.parse_values = parse_constant + + # Parametric is defined by generator operation and a function + # The function will use parsed parameter values to compute the unitary + super().__init__( + generator_operation, + qubit_support, + operator_function=self._construct_parametric_base_op, + noise=noise, + ) + self.register_buffer("identity", OPERATIONS_DICT["I"]) + + def extra_repr(self) -> str: + """String representation of the operation. + + Returns: + String with information on operation. + """ + return f"target: {self.qubit_support}, param: {self.param_name}" + _repr_noise( + self.noise + ) + + def __hash__(self) -> int: + """Hash qubit support and param_name + + Returns: + Hash value + """ + return hash(self.qubit_support) + hash(self.param_name) + + @staticmethod + def _expand_values(values: Tensor) -> Tensor: + """Expand values if necessary. + + Arguments: + values: Values of parameters + Returns: + Values of parameters expanded. + + """ + return values.unsqueeze(0) if len(values.size()) == 0 else values + + def _construct_parametric_base_op( + self, + values: dict[str, Tensor] | Tensor = dict(), + embedding: Embedding | None = None, + ) -> Tensor: + """ + Get the corresponding unitary with parsed values. + + Arguments: + values: A dict containing a Parameter name and value. + embedding: An optional embedding for parameters. + + Returns: + The unitary representation. + """ + thetas = self.parse_values(values, embedding) + batch_size = len(thetas) + mat = parametric_unitary(thetas, self.operation, self.identity, batch_size) + return mat + + def jacobian( + self, + values: dict[str, Tensor] | Tensor = dict(), + embedding: Embedding | None = None, + ) -> Tensor: + """ + Get the corresponding unitary of the jacobian. + + Arguments: + values: Parameter value. + + Returns: + The unitary representation of the jacobian. + """ + thetas = self.parse_values(values, embedding) + batch_size = len(thetas) + return _jacobian(thetas, self.operation, self.identity, batch_size) + + def to(self, *args: Any, **kwargs: Any) -> Parametric: + super().to(*args, **kwargs) + self._device = self.operation.device + self.param_name = ( + self.param_name.to(*args, **kwargs) + if isinstance(self.param_name, torch.Tensor) + else self.param_name + ) + self._dtype = self.operation.dtype + return self + + +class ControlledParametric(Parametric): + """ + Primitives for controlled parametric operations. + + Attributes: + control: Control qubit(s). + """ + + def __init__( + self, + operation: str | Tensor, + control: int | Tuple[int, ...], + target: int | Tuple[int, ...], + param_name: str | int | float | torch.Tensor = "", + noise: NoiseProtocol | dict[str, NoiseProtocol] | None = None, + ): + """Initializes a ControlledParametric. + + Arguments: + operation: Rotation gate. + control: Control qubit(s). + qubit_targets: Target qubit(s). + param_name: Name of parameters. + """ + support = Support(target, control) + super().__init__(operation, support, param_name, noise=noise) + + def extra_repr(self) -> str: + """String representation of the operation. + + Returns: + String with information on operation. + """ + return ( + f"control: {self.control}, target: {(self.target,)}, param: {self.param_name}" + + _repr_noise(self.noise) + ) + + def _construct_parametric_base_op( + self, values: dict[str, Tensor] = dict(), embedding: Embedding | None = None + ) -> Tensor: + """ + Get the corresponding unitary with parsed values and kronned identities + for control. + + Arguments: + values: A dict containing a Parameter name and value. + embedding: An optional embedding for parameters. + + Returns: + The unitary representation. + """ + thetas = self.parse_values(values, embedding) + batch_size = len(thetas) + mat = parametric_unitary(thetas, self.operation, self.identity, batch_size) + mat = controlled(mat, batch_size, len(self.control)) + return mat + + def jacobian( + self, values: dict[str, Tensor] = dict(), embedding: Embedding | None = None + ) -> Tensor: + """ + Get the corresponding unitary of the jacobian. + + Arguments: + values: Parameter value. + + Returns: + The unitary representation of the jacobian. + """ + thetas = self.parse_values(values, embedding) + batch_size = len(thetas) + n_control = len(self.control) + jU = _jacobian(thetas, self.operation, self.identity, batch_size) + n_dim = 2 ** (n_control + 1) + jC = ( + torch.zeros((n_dim, n_dim), dtype=self.identity.dtype) + .unsqueeze(2) + .repeat(1, 1, batch_size) + ) + unitary_idx = 2 ** (n_control + 1) - 2 + jC[unitary_idx:, unitary_idx:, :] = jU + return jC + + +class ControlledRotationGate(ControlledParametric): + """ + Primitives for controlled rotation operations. + """ + + n_params = 1 + + def __init__( + self, + operation: str | Tensor, + control: int | Tuple[int, ...], + target: int, + param_name: str | int | float | torch.Tensor = "", + noise: NoiseProtocol | dict[str, NoiseProtocol] | None = None, + ): + """Initializes a ControlledRotationGate. + + Arguments: + gate: Rotation gate. + control: Control qubit(s). + qubit_support: Target qubit. + param_name: Name of parameters. + """ + super().__init__(operation, control, target, param_name, noise=noise) + + @cached_property + def eigenvals_generator(self) -> Tensor: + """Get eigenvalues of the underlying generator. + + Arguments: + values: Parameter values. + + Returns: + Eigenvalues of the generator operator. + """ + return torch.cat( + ( + torch.zeros( + 2 ** len(self.qubit_support) - 2, + device=self.device, + dtype=self.dtype, + ), + pauli_singleq_eigenvalues.flatten().to( + device=self.device, dtype=self.dtype + ), + ) + ).reshape(-1, 1) diff --git a/pyqtorch/parametric.py b/pyqtorch/primitives/parametric_gates.py similarity index 58% rename from pyqtorch/parametric.py rename to pyqtorch/primitives/parametric_gates.py index 4db1685a..6d9b0ee3 100644 --- a/pyqtorch/parametric.py +++ b/pyqtorch/primitives/parametric_gates.py @@ -1,7 +1,7 @@ from __future__ import annotations from functools import cached_property -from typing import Any, Tuple +from typing import Any import torch from torch import Tensor @@ -9,190 +9,13 @@ from pyqtorch.embed import Embedding from pyqtorch.matrices import ( DEFAULT_MATRIX_DTYPE, - OPERATIONS_DICT, - _controlled, - _jacobian, - _unitary, + controlled, ) -from pyqtorch.primitive import Primitive -from pyqtorch.utils import Operator +from pyqtorch.noise import NoiseProtocol -pauli_singleq_eigenvalues = torch.tensor([[-1.0], [1.0]], dtype=torch.cdouble) - - -class Parametric(Primitive): - """ - Primitives taking parameters as input. - - Attributes: - param_name: Name of parameters. - parse_values: Method defining how to handle the values dictionary input. - """ - - n_params = 1 - - def __init__( - self, - generator_name: str, - target: int, - param_name: str | int | float | torch.Tensor = "", - ): - """Initializes Parametric. - - Arguments: - generator_name: Name of the operation. - target: Target qubit. - param_name: Name of parameters. - """ - super().__init__(OPERATIONS_DICT[generator_name], target) - self.register_buffer("identity", OPERATIONS_DICT["I"]) - self.param_name = param_name - - def parse_values( - values: dict[str, Tensor] | Tensor = dict(), - embedding: Embedding | None = None, - ) -> Tensor: - """The legacy way of using parametric gates: - The Parametric gate received a string as a 'param_name' and performs a - a lookup in the passed `values` dict for to retrieve the torch.Tensor passed - under the key `param_name`. - - Arguments: - values: A dict containing param_name:torch.Tensor pairs - embedding: An optional embedding. - Returns: - A Torch Tensor denoting values for the `param_name`. - """ - return Parametric._expand_values(values[self.param_name]) # type: ignore[index] - - def parse_tensor( - values: dict[str, Tensor] | Tensor = dict(), - embedding: Embedding | None = None, - ) -> Tensor: - """Functional version of the Parametric gate: - In case the user did not pass a `param_name`, - pyqtorch assumes `values` will be a torch.Tensor instead of a dict. - - Arguments: - values: A dict containing param_name:torch.Tensor pairs - Returns: - A Torch Tensor with which to evaluate the Parametric Gate. - """ - # self.param_name will be "" - return Parametric._expand_values(values) - - def parse_constant( - values: dict[str, Tensor] | Tensor = dict(), - embedding: Embedding | None = None, - ) -> Tensor: - """Fix a the parameter of a Parametric Gate to a numeric constant - if the user passed a numeric input for the `param_name`. - - Arguments: - values: A dict containing param_name:torch.Tensor pairs - Returns: - A Torch Tensor with which to evaluate the Parametric Gate. - """ - # self.param_name will be a torch.Tensor - return Parametric._expand_values( - torch.tensor(self.param_name, device=self.device, dtype=self.dtype) - ) - - if param_name == "": - self.parse_values = parse_tensor - self.param_name = self.param_name - elif isinstance(param_name, str): - self.parse_values = parse_values - elif isinstance(param_name, (float, int, torch.Tensor)): - self.parse_values = parse_constant - - def extra_repr(self) -> str: - """String representation of the operation. - - Returns: - String with information on operation. - """ - return f"target:{self.qubit_support}, param:{self.param_name}" - - @cached_property - def eigenvals_generator(self) -> Tensor: - """Get eigenvalues of the underlying generator. - - Arguments: - values: Parameter values. - - Returns: - Eigenvalues of the generator operator. - """ - return torch.linalg.eigvalsh(self.pauli).reshape(-1, 1) - - def __hash__(self) -> int: - """Hash qubit support and param_name - - Returns: - Hash value - """ - return hash(self.qubit_support) + hash(self.param_name) - - @staticmethod - def _expand_values(values: Tensor) -> Tensor: - """Expand values if necessary. - - Arguments: - values: Values of parameters - Returns: - Values of parameters expanded. - - """ - return values.unsqueeze(0) if len(values.size()) == 0 else values - - def unitary( - self, - values: dict[str, Tensor] | Tensor = dict(), - embedding: Embedding | None = None, - ) -> Operator: - """ - Get the corresponding unitary. - - Arguments: - values: A dict containing a Parameter name and value. - embedding: An optional embedding for parameters. - - Returns: - The unitary representation. - """ - thetas = self.parse_values(values, embedding) - batch_size = len(thetas) - return _unitary(thetas, self.pauli, self.identity, batch_size) - - def jacobian( - self, - values: dict[str, Tensor] | Tensor = dict(), - embedding: Embedding | None = None, - ) -> Operator: - """ - Get the corresponding unitary of the jacobian. +from .parametric import ControlledRotationGate, Parametric - Arguments: - values: Parameter value. - - Returns: - The unitary representation of the jacobian. - """ - thetas = self.parse_values(values, embedding) - batch_size = len(thetas) - return _jacobian(thetas, self.pauli, self.identity, batch_size) - - def to(self, *args: Any, **kwargs: Any) -> Primitive: - super().to(*args, **kwargs) - self._device = self.pauli.device - self.param_name = ( - self.param_name.to(*args, **kwargs) - if isinstance(self.param_name, torch.Tensor) - else self.param_name - ) - self._dtype = self.pauli.dtype - return self +pauli_singleq_eigenvalues = torch.tensor([[-1.0], [1.0]], dtype=torch.cdouble) class RX(Parametric): @@ -211,15 +34,17 @@ class RX(Parametric): def __init__( self, target: int, - param_name: str = "", + param_name: str | int | float | torch.Tensor = "", + noise: NoiseProtocol | dict[str, NoiseProtocol] | None = None, ): """Initializes RX. Arguments: target: Target qubit. param_name: Name of parameters. + noise: Optional noise protocols to apply. """ - super().__init__("X", target, param_name) + super().__init__("X", target, param_name, noise=noise) @cached_property def eigenvals_generator(self) -> Tensor: @@ -250,15 +75,17 @@ class RY(Parametric): def __init__( self, target: int, - param_name: str = "", + param_name: str | int | float | torch.Tensor = "", + noise: NoiseProtocol | dict[str, NoiseProtocol] | None = None, ): """Initializes RY. Arguments: target: Target qubit. param_name: Name of parameters. + noise: Optional noise protocols to apply. """ - super().__init__("Y", target, param_name) + super().__init__("Y", target, param_name, noise=noise) @cached_property def eigenvals_generator(self) -> Tensor: @@ -289,15 +116,17 @@ class RZ(Parametric): def __init__( self, target: int, - param_name: str = "", + param_name: str | int | float | torch.Tensor = "", + noise: NoiseProtocol | dict[str, NoiseProtocol] | None = None, ): """Initializes RZ. Arguments: target: Target qubit. param_name: Name of parameters. + noise: Optional noise protocols to apply. """ - super().__init__("Z", target, param_name) + super().__init__("Z", target, param_name, noise=noise) @cached_property def eigenvals_generator(self) -> Tensor: @@ -328,15 +157,17 @@ class PHASE(Parametric): def __init__( self, target: int, - param_name: str = "", + param_name: str | int | float | torch.Tensor = "", + noise: NoiseProtocol | dict[str, NoiseProtocol] | None = None, ): """Initializes PHASE. Arguments: target: Target qubit. param_name: Name of parameters. + noise: Optional noise protocols to apply. """ - super().__init__("I", target, param_name) + super().__init__("I", target, param_name, noise=noise) @cached_property def eigenvals_generator(self) -> Tensor: @@ -350,9 +181,9 @@ def eigenvals_generator(self) -> Tensor: """ return torch.tensor([[0.0], [2.0]], dtype=self.dtype, device=self.device) - def unitary( + def _construct_parametric_base_op( self, values: dict[str, Tensor] = dict(), embedding: Embedding | None = None - ) -> Operator: + ) -> Tensor: """ Get the corresponding unitary. @@ -371,7 +202,7 @@ def unitary( def jacobian( self, values: dict[str, Tensor] = dict(), embedding: Embedding | None = None - ) -> Operator: + ) -> Tensor: """ Get the corresponding unitary of the jacobian. @@ -391,116 +222,6 @@ def jacobian( return batch_mat -class ControlledRotationGate(Parametric): - """ - Primitives for controlled rotation operations. - - Attributes: - control: Control qubit(s). - qubit_support: Qubits acted on. - """ - - n_params = 1 - - def __init__( - self, - gate: str, - control: int | Tuple[int, ...], - target: int, - param_name: str = "", - ): - """Initializes a ControlledRotationGate. - - Arguments: - gate: Rotation gate. - control: Control qubit(s). - target: Target qubit. - param_name: Name of parameters. - """ - self.control = control if isinstance(control, tuple) else (control,) - super().__init__(gate, target, param_name) - self.qubit_support = self.control + (self.target,) # type: ignore[operator] - # In this class, target is always an int but herit from Parametric and Primitive that: - # target : int | tuple[int,...] - - def extra_repr(self) -> str: - """String representation of the operation. - - Returns: - String with information on operation. - """ - return ( - f"control: {self.control}, target:{(self.target,)}, param:{self.param_name}" - ) - - @cached_property - def eigenvals_generator(self) -> Tensor: - """Get eigenvalues of the underlying generator. - - Arguments: - values: Parameter values. - - Returns: - Eigenvalues of the generator operator. - """ - return torch.cat( - ( - torch.zeros( - 2 ** len(self.qubit_support) - 2, - device=self.device, - dtype=self.dtype, - ), - pauli_singleq_eigenvalues.flatten().to( - device=self.device, dtype=self.dtype - ), - ) - ).reshape(-1, 1) - - def unitary( - self, values: dict[str, Tensor] = dict(), embedding: Embedding | None = None - ) -> Operator: - """ - Get the corresponding unitary. - - Arguments: - values: A dict containing a Parameter name and value. - embedding: An optional embedding for parameters. - - Returns: - The unitary representation. - """ - thetas = self.parse_values(values, embedding) - batch_size = len(thetas) - mat = _unitary(thetas, self.pauli, self.identity, batch_size) - return _controlled(mat, batch_size, len(self.control)) - - def jacobian( - self, values: dict[str, Tensor] = dict(), embedding: Embedding | None = None - ) -> Operator: - """ - Get the corresponding unitary of the jacobian. - - Arguments: - values: Parameter value. - - Returns: - The unitary representation of the jacobian. - """ - thetas = self.parse_values(values, embedding) - batch_size = len(thetas) - n_control = len(self.control) - jU = _jacobian(thetas, self.pauli, self.identity, batch_size) - n_dim = 2 ** (n_control + 1) - jC = ( - torch.zeros((n_dim, n_dim), dtype=self.identity.dtype) - .unsqueeze(2) - .repeat(1, 1, batch_size) - ) - unitary_idx = 2 ** (n_control + 1) - 2 - jC[unitary_idx:, unitary_idx:, :] = jU - return jC - - class CRX(ControlledRotationGate): """ Primitive for the controlled RX gate. @@ -508,9 +229,10 @@ class CRX(ControlledRotationGate): def __init__( self, - control: int | Tuple[int, ...], + control: int | tuple[int, ...], target: int, - param_name: str = "", + param_name: str | int | float | Tensor = "", + noise: NoiseProtocol | dict[str, NoiseProtocol] | None = None, ): """Initializes controlled RX. @@ -518,8 +240,9 @@ def __init__( control: Control qubit(s). target: Target qubit. param_name: Name of parameters. + noise: Optional noise protocols to apply. """ - super().__init__("X", control, target, param_name) + super().__init__("X", control, target, param_name, noise=noise) @cached_property def eigenvals_generator(self) -> Tensor: @@ -538,7 +261,9 @@ def eigenvals_generator(self) -> Tensor: device=self.device, dtype=self.dtype, ), - torch.linalg.eigvalsh(self.pauli), + pauli_singleq_eigenvalues.flatten().to( + device=self.device, dtype=self.dtype + ), ) ).reshape(-1, 1) @@ -550,9 +275,10 @@ class CRY(ControlledRotationGate): def __init__( self, - control: int | Tuple[int, ...], + control: int | tuple[int, ...], target: int, - param_name: str = "", + param_name: str | int | float | Tensor = "", + noise: NoiseProtocol | dict[str, NoiseProtocol] | None = None, ): """Initializes controlled RY. @@ -560,8 +286,9 @@ def __init__( control: Control qubit(s). target: Target qubit. param_name: Name of parameters. + noise: Optional noise protocols to apply. """ - super().__init__("Y", control, target, param_name) + super().__init__("Y", control, target, param_name, noise=noise) @cached_property def eigenvals_generator(self) -> Tensor: @@ -594,9 +321,10 @@ class CRZ(ControlledRotationGate): def __init__( self, - control: int | Tuple[int, ...], + control: int | tuple[int, ...], target: int, - param_name: str = "", + param_name: str | int | float | Tensor = "", + noise: NoiseProtocol | dict[str, NoiseProtocol] | None = None, ): """Initializes controlled RZ. @@ -604,8 +332,9 @@ def __init__( control: Control qubit(s). target: Target qubit. param_name: Name of parameters. + noise: Optional noise protocols to apply. """ - super().__init__("Z", control, target, param_name) + super().__init__("Z", control, target, param_name, noise=noise) @cached_property def eigenvals_generator(self) -> Tensor: @@ -640,9 +369,10 @@ class CPHASE(ControlledRotationGate): def __init__( self, - control: int | Tuple[int, ...], + control: int | tuple[int, ...], target: int, - param_name: str = "", + param_name: str | int | float | Tensor = "", + noise: NoiseProtocol | dict[str, NoiseProtocol] | None = None, ): """Initializes controlled PHASE. @@ -650,8 +380,9 @@ def __init__( control: Control qubit(s). target: Target qubit. param_name: Name of parameters. + noise: Optional noise protocols to apply. """ - super().__init__("I", control, target, param_name) + super().__init__("I", control, target, param_name, noise=noise) @cached_property def eigenvals_generator(self) -> Tensor: @@ -674,9 +405,9 @@ def eigenvals_generator(self) -> Tensor: ) ).reshape(-1, 1) - def unitary( + def _construct_parametric_base_op( self, values: dict[str, Tensor] = dict(), embedding: Embedding | None = None - ) -> Operator: + ) -> Tensor: """ Get the corresponding unitary. @@ -691,11 +422,11 @@ def unitary( batch_size = len(thetas) mat = self.identity.unsqueeze(2).repeat(1, 1, batch_size) mat[1, 1, :] = torch.exp(1.0j * thetas).unsqueeze(0).unsqueeze(1) - return _controlled(mat, batch_size, len(self.control)) + return controlled(mat, batch_size, len(self.control)) def jacobian( self, values: dict[str, Tensor] = dict(), embedding: Embedding | None = None - ) -> Operator: + ) -> Tensor: """ Get the corresponding unitary of the jacobian. @@ -724,7 +455,7 @@ def jacobian( jC[unitary_idx:, unitary_idx:, :] = jU return jC - def to(self, *args: Any, **kwargs: Any) -> Primitive: + def to(self, *args: Any, **kwargs: Any) -> Parametric: """Set device of primitive. Returns: @@ -756,18 +487,27 @@ class U(Parametric): n_params = 3 - def __init__(self, target: int, phi: str, theta: str, omega: str): + def __init__( + self, + target: int, + phi: str, + theta: str, + omega: str, + noise: NoiseProtocol | dict[str, NoiseProtocol] | None = None, + ): """Initializes U gate. Arguments: + target: Target qubit. phi: Phi parameter. theta: Theta parameter. omega: Omega parameter. + noise: Optional noise protocols to apply. """ self.phi = phi self.theta = theta self.omega = omega - super().__init__("X", target, param_name="") + super().__init__("X", target, param_name="", noise=noise) self.register_buffer( "a", torch.tensor([[1, 0], [0, 0]], dtype=DEFAULT_MATRIX_DTYPE).unsqueeze(2) @@ -794,9 +534,9 @@ def eigenvals_generator(self) -> Tensor: """ return pauli_singleq_eigenvalues.to(device=self.device, dtype=self.dtype) - def unitary( + def _construct_parametric_base_op( self, values: dict[str, Tensor] = dict(), embedding: Embedding | None = None - ) -> Operator: + ) -> Tensor: """ Get the corresponding unitary. @@ -830,7 +570,7 @@ def unitary( def jacobian( self, values: dict[str, Tensor] = {}, embedding: Embedding | None = None - ) -> Operator: + ) -> Tensor: """ Get the corresponding unitary of the jacobian. @@ -852,13 +592,14 @@ def digital_decomposition(self) -> list[Parametric]: Returns: The digital decomposition. """ + target = self.target[0] return [ - RZ(self.qubit_support[0], self.phi), - RY(self.qubit_support[0], self.theta), - RZ(self.qubit_support[0], self.omega), + RZ(target, self.phi), + RY(target, self.theta), + RZ(target, self.omega), ] - def jacobian_decomposed(self, values: dict[str, Tensor] = dict()) -> list[Operator]: + def jacobian_decomposed(self, values: dict[str, Tensor] = dict()) -> list[Tensor]: """ Get the corresponding unitary decomposition of the jacobian. @@ -872,3 +613,8 @@ def jacobian_decomposed(self, values: dict[str, Tensor] = dict()) -> list[Operat NotImplementedError """ return [op.jacobian(values) for op in self.digital_decomposition()] + + +OPS_PARAM_1Q = {PHASE, RX, RY, RZ, U} +OPS_PARAM_2Q = {CPHASE, CRX, CRY, CRZ} +OPS_PARAM = OPS_PARAM_1Q.union(OPS_PARAM_2Q) diff --git a/pyqtorch/primitives/primitive.py b/pyqtorch/primitives/primitive.py new file mode 100644 index 00000000..e6ce3cbc --- /dev/null +++ b/pyqtorch/primitives/primitive.py @@ -0,0 +1,92 @@ +from __future__ import annotations + +from functools import cached_property +from typing import Any + +import torch +from torch import Tensor + +from pyqtorch.matrices import OPERATIONS_DICT, controlled +from pyqtorch.noise import NoiseProtocol, _repr_noise +from pyqtorch.quantum_operation import QuantumOperation, Support + + +class Primitive(QuantumOperation): + """Primitive operators based on a fixed matrix U. + + + Attributes: + operation (Tensor): Matrix U. + qubit_support: List of qubits the QuantumOperation acts on. + generator (Tensor): A tensor G s.t. U = exp(-iG). + """ + + def __init__( + self, + operation: Tensor, + qubit_support: int | tuple[int, ...] | Support, + generator: Tensor | None = None, + noise: NoiseProtocol | dict[str, NoiseProtocol] | None = None, + ) -> None: + super().__init__(operation, qubit_support, noise=noise) + self.generator = generator + + def to(self, *args: Any, **kwargs: Any) -> Primitive: + """Do device or dtype conversions. + + Returns: + Primitive: Converted instance. + """ + super().to(*args, **kwargs) + if self.generator is not None: + self.generator.to(*args, **kwargs) + return self + + @cached_property + def eigenvals_generator(self) -> Tensor: + """Get eigenvalues of the underlying generator. + + Note that for a primitive, the generator is unclear + so we execute pass. + + Arguments: + values: Parameter values. + + Returns: + Eigenvalues of the generator operator. + """ + if self.generator is not None: + return torch.linalg.eigvalsh(self.generator).reshape(-1, 1) + pass + + +class ControlledPrimitive(Primitive): + """Primitive applied depending on control qubits. + + Attributes: + operation (Tensor): Unitary tensor U. + control (int | tuple[int, ...]): List of qubits acting as controls. + target (int | tuple[int, ...]): List of qubits operations acts on. + """ + + def __init__( + self, + operation: str | Tensor, + control: int | tuple[int, ...], + target: int | tuple[int, ...], + noise: NoiseProtocol | dict[str, NoiseProtocol] | None = None, + ): + support = Support(target, control) + if isinstance(operation, str): + operation = OPERATIONS_DICT[operation] + operation = controlled( + operation=operation.unsqueeze(2), + batch_size=1, + n_control_qubits=len(support.control), + ).squeeze(2) + super().__init__(operation, support, noise=noise) + + def extra_repr(self) -> str: + return f"control: {self.control}, target: {self.target}" + _repr_noise( + self.noise + ) diff --git a/pyqtorch/primitives/primitive_gates.py b/pyqtorch/primitives/primitive_gates.py new file mode 100644 index 00000000..1b0aa43b --- /dev/null +++ b/pyqtorch/primitives/primitive_gates.py @@ -0,0 +1,192 @@ +from __future__ import annotations + +from pyqtorch.matrices import OPERATIONS_DICT +from pyqtorch.noise import NoiseProtocol +from pyqtorch.quantum_operation import Support +from pyqtorch.utils import ( + product_state, + qubit_support_as_tuple, +) + +from .primitive import ControlledPrimitive, Primitive + + +class X(Primitive): + def __init__( + self, + target: int, + noise: NoiseProtocol | dict[str, NoiseProtocol] | None = None, + ): + super().__init__(OPERATIONS_DICT["X"], target, noise=noise) + + +class Y(Primitive): + def __init__( + self, + target: int, + noise: NoiseProtocol | dict[str, NoiseProtocol] | None = None, + ): + super().__init__(OPERATIONS_DICT["Y"], target, noise=noise) + + +class Z(Primitive): + def __init__( + self, + target: int, + noise: NoiseProtocol | dict[str, NoiseProtocol] | None = None, + ): + super().__init__(OPERATIONS_DICT["Z"], target, noise=noise) + + +class I(Primitive): # noqa: E742 + def __init__( + self, + target: int, + noise: NoiseProtocol | dict[str, NoiseProtocol] | None = None, + ): + super().__init__(OPERATIONS_DICT["I"], target, noise=noise) + + +class H(Primitive): + def __init__( + self, + target: int, + noise: NoiseProtocol | dict[str, NoiseProtocol] | None = None, + ): + super().__init__(OPERATIONS_DICT["H"], target, noise=noise) + + +class T(Primitive): + def __init__( + self, + target: int, + noise: NoiseProtocol | dict[str, NoiseProtocol] | None = None, + ): + super().__init__(OPERATIONS_DICT["T"], target, noise=noise) + + +class S(Primitive): + def __init__( + self, + target: int, + noise: NoiseProtocol | dict[str, NoiseProtocol] | None = None, + ): + super().__init__( + OPERATIONS_DICT["S"], target, 0.5 * OPERATIONS_DICT["Z"], noise=noise + ) + + +class SDagger(Primitive): + def __init__( + self, + target: int, + noise: NoiseProtocol | dict[str, NoiseProtocol] | None = None, + ): + super().__init__( + OPERATIONS_DICT["SDAGGER"], target, -0.5 * OPERATIONS_DICT["Z"], noise=noise + ) + + +class Projector(Primitive): + def __init__( + self, + qubit_support: int | tuple[int, ...], + ket: str, + bra: str, + noise: NoiseProtocol | dict[str, NoiseProtocol] | None = None, + ): + + qubit_support = qubit_support_as_tuple(qubit_support) + if len(ket) != len(bra): + raise ValueError("Input ket and bra bitstrings must be of same length.") + if len(qubit_support) != len(ket): + raise ValueError( + "Qubit support must have the same number of qubits of ket and bra states." + ) + ket_state = product_state(ket).flatten() + bra_state = product_state(bra).flatten() + super().__init__( + OPERATIONS_DICT["PROJ"](ket_state, bra_state), qubit_support, noise=noise + ) + + +class N(Primitive): + def __init__( + self, + target: int, + noise: NoiseProtocol | dict[str, NoiseProtocol] | None = None, + ): + super().__init__(OPERATIONS_DICT["N"], target, noise=noise) + + +class SWAP(Primitive): + def __init__( + self, + i: int, + j: int, + noise: NoiseProtocol | dict[str, NoiseProtocol] | None = None, + ): + super().__init__(OPERATIONS_DICT["SWAP"], (i, j), noise=noise) + + +class CSWAP(Primitive): + def __init__( + self, + control: int, + target: tuple[int, ...], + noise: NoiseProtocol | dict[str, NoiseProtocol] | None = None, + ): + if not isinstance(target, tuple) or len(target) != 2: + raise ValueError("Target qubits must be a tuple with two qubits.") + support = Support(target=qubit_support_as_tuple(control) + target) + super().__init__(OPERATIONS_DICT["CSWAP"], support) + + +class CNOT(ControlledPrimitive): + def __init__( + self, + control: int | tuple[int, ...], + target: int, + noise: NoiseProtocol | dict[str, NoiseProtocol] | None = None, + ): + super().__init__("X", control, target, noise=noise) + + +CX = CNOT + + +class CY(ControlledPrimitive): + def __init__( + self, + control: int | tuple[int, ...], + target: int, + noise: NoiseProtocol | dict[str, NoiseProtocol] | None = None, + ): + super().__init__("Y", control, target, noise=noise) + + +class CZ(ControlledPrimitive): + def __init__( + self, + control: int | tuple[int, ...], + target: int, + noise: NoiseProtocol | dict[str, NoiseProtocol] | None = None, + ): + super().__init__("Z", control, target, noise=noise) + + +class Toffoli(ControlledPrimitive): + def __init__( + self, + control: int | tuple[int, ...], + target: int, + noise: NoiseProtocol | dict[str, NoiseProtocol] | None = None, + ): + super().__init__("X", control, target, noise=noise) + + +OPS_PAULI = {X, Y, Z, I} +OPS_1Q = OPS_PAULI.union({H, S, T}) +OPS_2Q = {CNOT, CY, CZ, SWAP} +OPS_3Q = {Toffoli, CSWAP} +OPS_DIGITAL = OPS_1Q.union(OPS_2Q, OPS_3Q) diff --git a/pyqtorch/quantum_operation.py b/pyqtorch/quantum_operation.py new file mode 100644 index 00000000..0b408968 --- /dev/null +++ b/pyqtorch/quantum_operation.py @@ -0,0 +1,396 @@ +from __future__ import annotations + +import logging +from functools import cached_property +from logging import getLogger +from math import log2 +from typing import Any, Callable + +import torch +from torch import Tensor + +from pyqtorch.apply import apply_density_mat, apply_operator +from pyqtorch.embed import Embedding +from pyqtorch.matrices import _dagger +from pyqtorch.noise import NoiseProtocol, _repr_noise +from pyqtorch.utils import ( + DensityMatrix, + density_mat, + expand_operator, + permute_basis, + qubit_support_as_tuple, +) + +logger = getLogger(__name__) + + +def forward_hook(*args, **kwargs) -> None: # type: ignore[no-untyped-def] + torch.cuda.nvtx.range_pop() + + +def pre_forward_hook(*args, **kwargs) -> None: # type: ignore[no-untyped-def] + torch.cuda.nvtx.range_push("QuantumOperation.forward") + + +def backward_hook(*args, **kwargs) -> None: # type: ignore[no-untyped-def] + torch.cuda.nvtx.range_pop() + + +def pre_backward_hook(*args, **kwargs) -> None: # type: ignore[no-untyped-def] + torch.cuda.nvtx.range_push("QuantumOperation.backward") + + +class Support: + """ + Generic representation of the qubit support. For single qubit operations, + a multiple index support indicates apply the operation for each index in the + support. + + Both target and control lists must be ordered! + + Attributes: + target = Index or indices where the operation is applied. + control = Index or indices to which the operation is conditioned to. + """ + + def __init__( + self, + target: int | tuple[int, ...], + control: int | tuple[int, ...] | None = None, + ) -> None: + self.target = qubit_support_as_tuple(target) + self.control = qubit_support_as_tuple(control) if control is not None else () + # if self.qubits != tuple(set(self.qubits)): + # raise ValueError("One or more qubits are defined both as control and target.") + + @classmethod + def target_all(cls) -> Support: + return Support(target=()) + + def __eq__(self, other: object) -> bool: + if not isinstance(other, Support): + return NotImplemented + + return self.target == other.target and self.control == other.control + + def __len__(self): + return len(self.qubits) + + @cached_property + def qubits(self) -> tuple[int, ...]: + return self.control + self.target + + @cached_property + def sorted_qubits(self) -> tuple[int, ...]: + return tuple(sorted(self.qubits)) + + def __repr__(self) -> str: + if not self.target: + return f"{self.__class__.__name__}.target_all()" + + subspace = f"target: {self.target}" + if self.control: + subspace += f", control: {self.control}" + + return f"{self.__class__.__name__}({subspace})" + + +class QuantumOperation(torch.nn.Module): + """Generic QuantumOperation class storing a tensor operation to represent either + a quantum operator or a tensor generator inferring the QuantumOperation. + + Attributes: + operation (Tensor): Tensor used to infer the QuantumOperation + directly or indirectly. + qubit_support (int | tuple[int, ...]): Tuple of qubits + the QuantumOperation acts on. + operator_function (Callable | None, optional): Function to generate the base operator + from operation. If None, we consider returning operation itself. + + """ + + def __init__( + self, + operation: Tensor, + qubit_support: int | tuple[int, ...] | Support, + operator_function: Callable | None = None, + noise: NoiseProtocol | dict[str, NoiseProtocol] | None = None, + ) -> None: + """Initializes QuantumOperation + + Args: + operation (Tensor): Tensor used to infer the QuantumOperation. + qubit_support (int | tuple[int, ...]): List of qubits + the QuantumOperation acts on. + + Raises: + ValueError: When operation has incompatible shape + with qubit_support + """ + super().__init__() + + self._qubit_support = ( + qubit_support + if isinstance(qubit_support, Support) + else Support(target=qubit_support) + ) + + self.register_buffer("operation", operation) + self._device = self.operation.device + self._dtype = self.operation.dtype + + is_primitive = operator_function is None + dim_nomatch = len(self.qubit_support) != int(log2(operation.shape[0])) + if is_primitive and dim_nomatch: + raise ValueError( + "The operation shape should match the length of the qubit_support." + ) + + if operator_function is None: + self._operator_function = self._default_operator_function + else: + self._operator_function = operator_function + + self.noise: NoiseProtocol | dict[str, NoiseProtocol] | None = noise + + if logger.isEnabledFor(logging.DEBUG): + # When Debugging let's add logging and NVTX markers + # WARNING: incurs performance penalty + self.register_forward_hook(forward_hook, always_call=True) + self.register_full_backward_hook(backward_hook) + self.register_forward_pre_hook(pre_forward_hook) + self.register_full_backward_pre_hook(pre_backward_hook) + + def to(self, *args: Any, **kwargs: Any) -> QuantumOperation: + """Do device or dtype conversions. + + Returns: + QuantumOperation: Converted instance. + """ + super().to(*args, **kwargs) + self._device = self.operation.device + self._dtype = self.operation.dtype + return self + + @property + def qubit_support(self) -> tuple[int, ...]: + """Getter qubit_support. + + Returns: + Support: Tuple of sorted qubits. + """ + return self._qubit_support.sorted_qubits + + @property + def target(self) -> tuple[int, ...]: + """Get target qubits. + + Returns: + tuple[int, ...]: The target qubits + """ + return self._qubit_support.target + + @property + def control(self) -> tuple[int, ...]: + """Get control qubits. + + Returns: + tuple[int, ...]: The control qubits + """ + return self._qubit_support.control + + @property + def operator_function(self) -> Callable[..., Any]: + """Getter operator_function. + + Returns: + Callable[..., Any]: Function for + getting base operator. + """ + return self._operator_function + + def __hash__(self) -> int: + return hash(self.qubit_support) + + def extra_repr(self) -> str: + return f"target: {self.qubit_support}" + _repr_noise(self.noise) + + @property + def device(self) -> torch.device: + """Returns device. + + Returns: + torch.device: Device. + """ + return self._device + + @property + def dtype(self) -> torch.dtype: + """Returns dtype. + + Returns: + torch.dtype: Dtype. + """ + return self._dtype + + @cached_property + def eigenvals_generator(self) -> Tensor: + """Get eigenvalues of the underlying operation. + + Arguments: + values: Parameter values. + + Returns: + Eigenvalues of the operation. + """ + return torch.linalg.eigvalsh(self.operation).reshape(-1, 1) + + @cached_property + def spectral_gap(self) -> Tensor: + """Difference between the moduli of the two largest eigenvalues of the generator. + + Returns: + Tensor: Spectral gap value. + """ + spectrum = self.eigenvals_generator + spectral_gap = torch.unique(torch.abs(torch.tril(spectrum - spectrum.T))) + return spectral_gap[spectral_gap.nonzero()] + + def _default_operator_function( + self, + values: dict[str, Tensor] | Tensor = dict(), + embedding: Embedding | None = None, + ) -> Tensor: + """Default operator_function simply returns the operation. + + Args: + values (dict[str, Tensor] | Tensor, optional): Parameter values. Defaults to dict(). + embedding (Embedding | None, optional): Optional embedding. Defaults to None. + + Returns: + Tensor: Base operator. + """ + operation = ( + self.operation.unsqueeze(2) + if len(self.operation.shape) == 2 + else self.operation + ) + return operation + + def dagger( + self, + values: dict[str, Tensor] | Tensor = dict(), + embedding: Embedding | None = None, + ) -> Tensor: + """Apply the dagger to operator. + + Args: + values (dict[str, Tensor] | Tensor, optional): Parameter values. Defaults to dict(). + embedding (Embedding | None, optional): Optional embedding. Defaults to None. + + Returns: + Tensor: conjugate transpose operator. + """ + return _dagger(self.operator_function(values, embedding)) + + def tensor( + self, + values: dict[str, Tensor] = dict(), + embedding: Embedding | None = None, + full_support: tuple[int, ...] | None = None, + ) -> Tensor: + """Get tensor of the QuantumOperation. + + Args: + values (dict[str, Tensor], optional): Parameter values. Defaults to dict(). + embedding (Embedding | None, optional): Optional embedding. Defaults to None. + full_support (tuple[int, ...] | None, optional): The qubits the returned tensor + will be defined over. Defaults to None for only using the qubit_support. + + Returns: + Tensor: Tensor representation of QuantumOperation. + """ + blockmat = self.operator_function(values, embedding) + if self._qubit_support.qubits != self.qubit_support: + blockmat = permute_basis(blockmat, self._qubit_support.qubits) + if full_support is None: + return blockmat + else: + return expand_operator(blockmat, self.qubit_support, full_support) + + def _forward( + self, + state: Tensor, + values: dict[str, Tensor] | Tensor = dict(), + embedding: Embedding | None = None, + ) -> Tensor: + if isinstance(state, DensityMatrix): + n_qubits = int(log2(state.size(1))) + full_support = tuple(range(n_qubits)) + return apply_density_mat( + self.tensor(values, embedding, full_support=full_support), state + ) + else: + return apply_operator( + state, + self.tensor(values, embedding), + self.qubit_support, + len(state.size()) - 1, + ) + + def _noise_forward( + self, + state: Tensor, + values: dict[str, Tensor] | Tensor = dict(), + embedding: Embedding | None = None, + ) -> Tensor: + if not isinstance(state, DensityMatrix): + state = density_mat(state) + n_qubits = int(log2(state.size(1))) + full_support = tuple(range(n_qubits)) + state = apply_density_mat( + self.tensor(values, embedding, full_support=full_support), state + ) + if isinstance(self.noise, dict): + for noise_instance in self.noise.values(): + protocol = noise_instance.protocol_to_gate() + noise_gate = protocol( + target=( + noise_instance.target + if noise_instance.target is not None + else self.target + ), + error_probability=noise_instance.error_probability, + ) + state = noise_gate(state, values) + return state + elif isinstance(self.noise, NoiseProtocol): + protocol = self.noise.protocol_to_gate() + noise_gate = protocol( + target=( + self.noise.target if self.noise.target is not None else self.target + ), + error_probability=self.noise.error_probability, + ) + return noise_gate(state, values) + + def forward( + self, + state: Tensor, + values: dict[str, Tensor] | Tensor = dict(), + embedding: Embedding | None = None, + ) -> Tensor: + """Apply the operation on input state or density matrix. + + Args: + state (Tensor): Input state. + values (dict[str, Tensor], optional): Parameter values. Defaults to dict(). + embedding (Embedding | None, optional): Optional embedding. Defaults to None. + + Returns: + Tensor: _description_ + """ + if self.noise: + return self._noise_forward(state, values, embedding) + else: + return self._forward(state, values, embedding) diff --git a/pyqtorch/utils.py b/pyqtorch/utils.py index 9d95a009..dc7b18cd 100644 --- a/pyqtorch/utils.py +++ b/pyqtorch/utils.py @@ -1,6 +1,7 @@ from __future__ import annotations import logging +from collections import Counter from dataclasses import dataclass from enum import Enum from functools import lru_cache, partial, wraps @@ -9,12 +10,13 @@ from string import ascii_uppercase as ABC from typing import Any, Callable, Sequence +import numpy as np import torch -from numpy import arange, array, delete, log2 +from numpy import arange, argsort, array, delete, log2 from numpy import ndarray as NDArray from torch import Tensor, moveaxis -from pyqtorch.matrices import DEFAULT_MATRIX_DTYPE, DEFAULT_REAL_DTYPE +from pyqtorch.matrices import DEFAULT_MATRIX_DTYPE, DEFAULT_REAL_DTYPE, IMAT State = Tensor Operator = Tensor @@ -22,12 +24,30 @@ ATOL = 1e-06 ATOL_embedding = 1e-03 RTOL = 0.0 -GRADCHECK_ATOL = 1e-06 +GRADCHECK_ATOL = 1e-05 +GRADCHECK_sampling_ATOL = 1e-01 +PSR_ACCEPTANCE = 1e-05 +GPSR_ACCEPTANCE = 1e-05 ABC_ARRAY: NDArray = array(list(ABC)) logger = getLogger(__name__) +def qubit_support_as_tuple(support: int | tuple[int, ...]) -> tuple[int, ...]: + """Make sure support returned is a tuple of integers. + + Args: + support (int | tuple[int, ...]): Qubit support. + + Returns: + tuple[int, ...]: Qubit support as tuple. + """ + if isinstance(support, np.integer): + return (support.item(),) + qubit_support = (support,) if isinstance(support, int) else tuple(support) + return qubit_support + + def inner_prod(bra: Tensor, ket: Tensor) -> Tensor: """ Compute the inner product :math:`\\langle\\bra|\\ket\\rangle` @@ -67,6 +87,44 @@ def overlap(bra: Tensor, ket: Tensor) -> Tensor: return torch.pow(inner_prod(bra, ket).real, 2) +def sample_multinomial( + probs: Tensor, + length_bitstring: int, + n_samples: int, + return_counter: bool = True, + minlength: int = 0, +) -> Counter | Tensor: + """Sample bitstrings from a probability distribution. + + Args: + probs (Tensor): Probability distribution + length_bitstring (int): Maximal length of bitstring. + n_samples (int): Number of samples to extract. + instead of ratios. + return_counter (bool): If True, return Counter object. + Otherwise, the result of torch.bincount is returned. + minlength (int): minimum number of bins. Should be non-negative. + + Returns: + Counter: Sampled bitstrings with their frequencies or probabilities. + """ + + bincount_output = torch.bincount( + torch.multinomial(input=probs, num_samples=n_samples, replacement=True), + minlength=minlength, + ) + + if return_counter: + return Counter( + { + format(k, "0{}b".format(length_bitstring)): count.item() + for k, count in enumerate(bincount_output) + if count > 0 + } + ) + return bincount_output + + class StrEnum(str, Enum): def __str__(self) -> str: """Used when dumping enum fields in a schema.""" @@ -94,6 +152,21 @@ class DiffMode(StrEnum): """The generalized parameter-shift rule""" +class DropoutMode(StrEnum): + """ + Which Dropout mode to use, using the methods stated in https://arxiv.org/abs/2310.04120. + + Options: rotational - Randomly drops entangling rotational gates. + entangling - Randomly drops entangling gates. + canonical_fwd - Randomly drops rotational gates and next immediate entangling + gates whose target bit is located on dropped rotational gates. + """ + + ROTATIONAL = "rotational_dropout" + ENTANGLING = "entangling_dropout" + CANONICAL_FWD = "canonical_fwd_dropout" + + def is_normalized(state: Tensor, atol: float = ATOL) -> bool: """ Function to test if the probabilities from the state sum up to 1. @@ -298,10 +371,34 @@ def operator_kron(op1: Tensor, op2: Tensor) -> Tensor: ) +def expand_operator( + operator: Tensor, qubit_support: tuple[int, ...], full_support: tuple[int, ...] +) -> Tensor: + """ + Expands an operator acting on a given qubit_support to act on a larger full_support + by explicitly filling in identity matrices on all remaining qubits. + """ + full_support = tuple(sorted(full_support)) + if not set(qubit_support).issubset(set(full_support)): + raise ValueError( + "Expanding tensor operation requires a `full_support` argument " + "larger than or equal to the `qubit_support`." + ) + device, dtype = operator.device, operator.dtype + for i in set(full_support) - set(qubit_support): + qubit_support += (i,) + other = IMAT.clone().to(device=device, dtype=dtype).unsqueeze(2) + operator = torch.kron(operator.contiguous(), other) + operator = permute_basis(operator, qubit_support) + return operator + + def promote_operator(operator: Tensor, target: int, n_qubits: int) -> Tensor: - from pyqtorch.primitive import I + from pyqtorch.primitives import I """ + FIXME: Remove and replace usage with the `expand_operator` above. + Promotes `operator` to the size of the circuit (number of qubits and batch). Targeting the first qubit implies target = 0, so target > n_qubits - 1. @@ -326,12 +423,36 @@ def promote_operator(operator: Tensor, target: int, n_qubits: int) -> Tensor: for qubit in qubits: operator = torch.where( target > qubit, - operator_kron(I(target).unitary(), operator), - operator_kron(operator, I(target).unitary()), + operator_kron(I(target).tensor(), operator), + operator_kron(operator, I(target).tensor()), ) return operator +def permute_basis(operator: Tensor, qubit_support: tuple) -> Tensor: + """Takes an operator tensor and permutes the rows and + columns according to the order of the qubit support. + + Args: + operator (Tensor): Operator to permute over. + qubit_support (tuple): Qubit support. + + Returns: + Tensor: Permuted operator. + """ + ordered_support = argsort(qubit_support) + n_qubits = len(qubit_support) + if all(a == b for a, b in zip(ordered_support, list(range(n_qubits)))): + return operator + batch_size = operator.size(-1) + operator = operator.view([2] * 2 * n_qubits + [batch_size]) + + perm = list( + tuple(ordered_support) + tuple(ordered_support + n_qubits) + (2 * n_qubits,) + ) + return operator.permute(perm).reshape([2**n_qubits, 2**n_qubits, batch_size]) + + def random_dm_promotion( target: int, dm_input: DensityMatrix, n_qubits: int ) -> DensityMatrix: diff --git a/tests/conftest.py b/tests/conftest.py index a6fd1195..fcd0b7e3 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -9,19 +9,41 @@ from pytest import FixtureRequest from torch import Tensor -from pyqtorch.apply import apply_operator from pyqtorch.noise import ( AmplitudeDamping, BitFlip, Depolarizing, GeneralizedAmplitudeDamping, Noise, + NoiseProtocol, PauliChannel, PhaseDamping, PhaseFlip, ) -from pyqtorch.parametric import PHASE, RX, RY, RZ -from pyqtorch.primitive import H, I, Primitive, S, T, X, Y, Z +from pyqtorch.primitives import ( + CNOT, + CPHASE, + CRX, + CRY, + CRZ, + CY, + CZ, + PHASE, + RX, + RY, + RZ, + ControlledPrimitive, + ControlledRotationGate, + H, + I, + Parametric, + Primitive, + S, + T, + X, + Y, + Z, +) from pyqtorch.utils import ( DensityMatrix, density_mat, @@ -31,27 +53,6 @@ ) -def _calc_mat_vec_wavefunction( - block: Primitive, n_qubits: int, init_state: torch.Tensor, values: dict = {} -) -> torch.Tensor: - """Get the result of applying the matrix representation of a block to an initial state. - - Args: - block: The black operator to apply. - n_qubits: Number of qubits in the circuit. - init_state: Initial state to apply block on. - values: Values of parameter if block is parametric. - - Returns: - Tensor: The new wavefunction after applying the block. - - """ - mat = block.tensor(n_qubits=len(block.qubit_support), values=values) - return apply_operator( - init_state, mat, qubits=block.qubit_support, n_qubits=n_qubits - ) - - @pytest.fixture(params=[I, X, Y, Z, H, T, S]) def gate(request: Primitive) -> Any: return request.param @@ -85,6 +86,11 @@ def random_input_state(n_qubits: int, batch_size: int) -> Any: return random_state(n_qubits, batch_size) +@pytest.fixture +def random_input_dm(random_input_state: Tensor) -> Any: + return density_mat(random_input_state) + + @pytest.fixture def rho_input(batch_size: int, target: int, n_qubits: int) -> Any: rho_0: DensityMatrix = density_mat(product_state("0", batch_size)) @@ -92,8 +98,14 @@ def rho_input(batch_size: int, target: int, n_qubits: int) -> Any: return rho_input +# TODO: create random noisy protocols +@pytest.fixture +def random_noise(): + pass + + @pytest.fixture -def random_gate(target: int) -> Any: +def random_single_qubit_gate(target: int) -> Any: GATES = [X, Y, Z, I, H] gate = random.choice(GATES) return gate(target) @@ -106,6 +118,43 @@ def random_rotation_gate(target: int) -> Any: return rotation_gate(target, "theta") +@pytest.fixture +def random_controlled_gate(n_qubits: int, target: int) -> Any: + if n_qubits < 2: + raise ValueError("The controlled gates are defined on 2 qubits minimum") + CONTROLLED_GATES = [CNOT, CY, CZ] + controlled_gate = random.choice(CONTROLLED_GATES) + control = random.choice([i for i in range(n_qubits) if i != target]) + return controlled_gate(control, target) + + +@pytest.fixture +def random_rotation_control_gate(n_qubits: int, target: int) -> Any: + if n_qubits < 2: + raise ValueError("The controlled gates are defined on 2 qubits minimum") + ROTATION_CONTROL_GATES = [CRX, CRY, CRZ, CPHASE] + controlled_gate = random.choice(ROTATION_CONTROL_GATES) + control = random.choice([i for i in range(n_qubits) if i != target]) + return controlled_gate(control, target, "theta") + + +@pytest.fixture +def random_unitary_gate( + random_single_qubit_gate: Primitive, + random_rotation_gate: Parametric, + random_controlled_gate: ControlledPrimitive, + random_rotation_control_gate: ControlledRotationGate, +) -> Any: + UNITARY_GATES = [ + random_single_qubit_gate, + random_controlled_gate, + random_rotation_gate, + random_rotation_control_gate, + ] + unitary_gate = random.choice(UNITARY_GATES) + return unitary_gate + + @pytest.fixture def random_noise_gate(n_qubits: int) -> Any: NOISE_GATES = [ @@ -126,13 +175,90 @@ def random_noise_gate(n_qubits: int) -> Any: return noise_gate(torch.randint(0, n_qubits, (1,)).item(), noise_prob) elif noise_gate == GeneralizedAmplitudeDamping: noise_prob = torch.rand(size=(2,)) - p, r = noise_prob[0], noise_prob[1] - return noise_gate(torch.randint(0, n_qubits, (1,)).item(), p, r) + return noise_gate(torch.randint(0, n_qubits, (1,)).item(), noise_prob) else: noise_prob = torch.rand(size=(1,)).item() return noise_gate(torch.randint(0, n_qubits, (1,)).item(), noise_prob) +@pytest.fixture +def random_noisy_protocol() -> Any: + NOISE_PROTOCOLS = [ + NoiseProtocol.BITFLIP, + NoiseProtocol.PHASEFLIP, + NoiseProtocol.PAULI_CHANNEL, + NoiseProtocol.DEPOLARIZING, + NoiseProtocol.AMPLITUDE_DAMPING, + NoiseProtocol.PHASE_DAMPING, + NoiseProtocol.GENERALIZED_AMPLITUDE_DAMPING, + ] + noise_protocol = random.choice(NOISE_PROTOCOLS) + if noise_protocol == NoiseProtocol.PAULI_CHANNEL: + noise_prob = torch.rand(size=(3,)) + noise_prob = noise_prob / ( + noise_prob.sum(dim=0, keepdim=True) + torch.rand((1,)).item() + ) + elif noise_protocol == NoiseProtocol.GENERALIZED_AMPLITUDE_DAMPING: + noise_prob = torch.rand(size=(2,)) + else: + noise_prob = torch.rand(size=(1,)).item() + return NoiseProtocol( + protocol=noise_protocol, + options={"error_probability": noise_prob}, + ) + + +@pytest.fixture +def random_noisy_unitary_gate( + n_qubits: int, target: int, random_noisy_protocol: NoiseProtocol +) -> Any: + if n_qubits < 2: + raise ValueError("The controlled gates are defined on 2 qubits minimum") + SINGLE_GATES = [X, Y, Z, I, H] + ROTATION_GATES = [RX, RY, RZ, PHASE] + CONTROLLED_GATES = [CNOT, CY, CZ] + ROTATION_CONTROL_GATES = [CRX, CRY, CRZ, CPHASE] + UNITARY_GATES = ( + SINGLE_GATES + ROTATION_GATES + CONTROLLED_GATES + ROTATION_CONTROL_GATES + ) + unitary_gate = random.choice(UNITARY_GATES) + protocol_gate = random_noisy_protocol.protocol_to_gate() + noise_gate = protocol_gate(target, random_noisy_protocol.error_probability) + if unitary_gate in SINGLE_GATES: + return ( + unitary_gate(target=target, noise=random_noisy_protocol), # type: ignore[call-arg] + unitary_gate(target=target), # type: ignore[call-arg] + noise_gate, + ) # type: ignore[call-arg] + if unitary_gate in ROTATION_GATES: + return ( + unitary_gate( + target=target, param_name="theta", noise=random_noisy_protocol + ), # type: ignore[call-arg] + unitary_gate(target=target, param_name="theta"), # type: ignore[call-arg] + noise_gate, + ) + if unitary_gate in CONTROLLED_GATES: + control = random.choice([i for i in range(n_qubits) if i != target]) + return ( + unitary_gate(control=control, target=target, noise=random_noisy_protocol), # type: ignore[call-arg] + unitary_gate(control=control, target=target), # type: ignore[call-arg] + noise_gate, + ) + if unitary_gate in ROTATION_CONTROL_GATES: + control = random.choice([i for i in range(n_qubits) if i != target]) + return ( + unitary_gate( + control=control, + target=target, + param_name="theta", + noise=random_noisy_protocol, + ), # type: ignore[call-arg] + unitary_gate(control=control, target=target, param_name="theta"), # type: ignore[call-arg] + noise_gate, + ) + + @pytest.fixture def random_flip_gate() -> Any: FLIP_GATES = [BitFlip, PhaseFlip, Depolarizing, PauliChannel] @@ -195,9 +321,9 @@ def flip_expected_state( @pytest.fixture def flip_gates_prob_0(random_flip_gate: Noise, target: int) -> Any: if random_flip_gate == PauliChannel: - FlipGate_0 = random_flip_gate(target, probabilities=(0, 0, 0)) + FlipGate_0 = random_flip_gate(target, error_probability=(0, 0, 0)) else: - FlipGate_0 = random_flip_gate(target, probability=0) + FlipGate_0 = random_flip_gate(target, error_probability=0) return FlipGate_0 @@ -206,13 +332,13 @@ def flip_gates_prob_1( random_flip_gate: Noise, target: int, random_input_state: Tensor ) -> Any: if random_flip_gate == BitFlip: - FlipGate_1 = random_flip_gate(target, probability=1) + FlipGate_1 = random_flip_gate(target, error_probability=1) expected_op = density_mat(X(target)(random_input_state)) elif random_flip_gate == PhaseFlip: - FlipGate_1 = random_flip_gate(target, probability=1) + FlipGate_1 = random_flip_gate(target, error_probability=1) expected_op = density_mat(Z(target)(random_input_state)) elif random_flip_gate == Depolarizing: - FlipGate_1 = random_flip_gate(target, probability=1) + FlipGate_1 = random_flip_gate(target, error_probability=1) expected_op = ( 1 / 3 @@ -224,7 +350,7 @@ def flip_gates_prob_1( ) elif random_flip_gate == PauliChannel: px, py, pz = 1 / 3, 1 / 3, 1 / 3 - FlipGate_1 = random_flip_gate(target, probabilities=(px, py, pz)) + FlipGate_1 = random_flip_gate(target, error_probability=(px, py, pz)) expected_op = ( px * density_mat(X(target)(random_input_state)) # type: ignore + py * density_mat(Y(target)(random_input_state)) @@ -265,7 +391,7 @@ def damping_expected_state( [[[1 - (r - p * r)], [0]], [[0], [r - p * r]]], dtype=torch.cdouble ) ) - DampingGate = random_damping_gate(target, p, r) + DampingGate = random_damping_gate(target, damping_rate) expected_state = random_dm_promotion(target, expected_state, n_qubits).repeat( 1, 1, batch_size ) @@ -278,9 +404,9 @@ def damping_expected_state( @pytest.fixture def damping_gates_prob_0(random_damping_gate: Noise, target: int) -> Any: if random_damping_gate == GeneralizedAmplitudeDamping: - damping_gate_0 = random_damping_gate(target, probability=0, rate=0) + damping_gate_0 = random_damping_gate(target, (0, 0)) else: - damping_gate_0 = random_damping_gate(target, rate=0) + damping_gate_0 = random_damping_gate(target, 0) return damping_gate_0 diff --git a/tests/helpers.py b/tests/helpers.py new file mode 100644 index 00000000..3a1df9f2 --- /dev/null +++ b/tests/helpers.py @@ -0,0 +1,105 @@ +from __future__ import annotations + +import random + +import torch + +from pyqtorch.apply import apply_operator +from pyqtorch.composite import Add, Scale, Sequence +from pyqtorch.primitives import ( + OPS_1Q, + OPS_2Q, + OPS_3Q, + OPS_PARAM_1Q, + OPS_PARAM_2Q, + OPS_PAULI, + Parametric, + Primitive, + Toffoli, +) + + +def calc_mat_vec_wavefunction( + block: Primitive | Sequence, + init_state: torch.Tensor, + values: dict = dict(), + full_support: tuple | None = None, +) -> torch.Tensor: + """Get the result of applying the matrix representation of a block to an initial state. + + Args: + block: The black operator to apply. + init_state: Initial state to apply block on. + values: Values of parameter if block is parametric. + + Returns: + Tensor: The new wavefunction after applying the block. + + """ + mat = block.tensor(values=values, full_support=full_support) + qubit_support = block.qubit_support if full_support is None else full_support + return apply_operator( + init_state, + mat, + qubits=qubit_support, + ) + + +def get_op_support(op: type[Primitive] | type[Parametric], n_qubits: int) -> tuple: + """Decides a random qubit support for any gate, up to a some max n_qubits.""" + if op in OPS_1Q.union(OPS_PARAM_1Q): + supp: tuple = (random.randint(0, n_qubits - 1),) + elif op in OPS_2Q.union(OPS_PARAM_2Q): + supp = tuple(random.sample(range(n_qubits), 2)) + elif op in OPS_3Q: + i, j, k = tuple(random.sample(range(n_qubits), 3)) + supp = ((i, j), k) if op == Toffoli else (i, (j, k)) + return supp + + +def random_pauli_hamiltonian( + n_qubits: int, + k_1q: int = 5, + k_2q: int = 10, + make_param: bool = False, + default_scale_coeffs: float | None = None, +) -> tuple[Sequence, list]: + """Creates a random Pauli Hamiltonian as a sum of k_1q + k_2q terms. + + Args: + n_qubits (int): Number of qubits. + k_1q (int, optional): Number of one-qubit terms. Defaults to 5. + k_2q (int, optional): Number of two-qubit terms. Defaults to 10. + make_param (bool, optional): Coefficients as parameters. Defaults to False. + default_scale_coeffs (float | None, optional): Default value for the parameter + of Scale operations. Defaults to None. + + Returns: + tuple[Sequence, list]: Hamiltonian and list of parameters. + """ + OPS_PAULI_choice = list(OPS_PAULI) + one_q_terms: list = random.choices(OPS_PAULI_choice, k=k_1q) + two_q_terms: list = [random.choices(OPS_PAULI_choice, k=2) for _ in range(k_2q)] + terms: list = [] + for term in one_q_terms: + supp = random.sample(range(n_qubits), 1) + terms.append(term(supp)) + for term in two_q_terms: + supp = random.sample(range(n_qubits), 2) + terms.append(Sequence([term[0](supp[0]), term[1](supp[1])])) + param_list = [] + for i, t in enumerate(terms): + if random.random() > 0.5: + if make_param: + terms[i] = Scale(t, f"p_{i}") + param_list.append(f"p_{i}") + else: + terms[i] = Scale( + t, + ( + torch.rand(1) + if not default_scale_coeffs + else torch.tensor([default_scale_coeffs]) + ), + ) + return Add(terms), param_list diff --git a/tests/test_analog.py b/tests/test_analog.py index f7d60e6f..6f148ca2 100644 --- a/tests/test_analog.py +++ b/tests/test_analog.py @@ -3,13 +3,12 @@ from collections import OrderedDict from typing import Callable -import numpy as np import pytest import torch -from conftest import _calc_mat_vec_wavefunction +from helpers import calc_mat_vec_wavefunction, random_pauli_hamiltonian import pyqtorch as pyq -from pyqtorch.analog import GeneratorType +from pyqtorch.hamiltonians import GeneratorType from pyqtorch.matrices import ( DEFAULT_MATRIX_DTYPE, DEFAULT_REAL_DTYPE, @@ -160,69 +159,6 @@ def test_hamiltonianevolution_with_types( assert torch.allclose(result, target, rtol=RTOL, atol=ATOL) -@pytest.mark.parametrize("n_qubits", [2, 4, 6]) -def test_hamevo_parametric_gen(n_qubits: int) -> None: - dim = torch.randint(1, n_qubits + 1, (1,)).item() - vparam = "theta" - sup = tuple(range(dim)) - parametric = True - ops = [pyq.X, pyq.Y, pyq.Z] - qubit_targets = np.random.choice(dim, len(ops), replace=True) - generator = pyq.Add([pyq.Scale(op(q), vparam) for op, q in zip(ops, qubit_targets)]) - hamevo = pyq.HamiltonianEvolution( - generator, vparam, sup, parametric, cache_length=2 - ) - assert hamevo.generator_type == GeneratorType.PARAMETRIC_OPERATION - assert len(hamevo._cache_hamiltonian_evo) == 0 - - def apply_hamevo_and_compare_expected(psi, vals): - psi_star = hamevo(psi, vals) - psi_expected = _calc_mat_vec_wavefunction(hamevo, n_qubits, psi, vals) - assert torch.allclose(psi_star, psi_expected, rtol=RTOL, atol=ATOL) - - vals = {vparam: torch.rand(1)} - psi = random_state(n_qubits) - apply_hamevo_and_compare_expected(psi, vals) - - # test cached - assert len(hamevo._cache_hamiltonian_evo) == 1 - apply_hamevo_and_compare_expected(psi, vals) - assert len(hamevo._cache_hamiltonian_evo) == 1 - - # test caching new value - vals[vparam] += 0.1 - apply_hamevo_and_compare_expected(psi, vals) - assert len(hamevo._cache_hamiltonian_evo) == 2 - - # changing input state should not change the cache - psi = random_state(n_qubits) - apply_hamevo_and_compare_expected(psi, vals) - assert len(hamevo._cache_hamiltonian_evo) == 2 - - # test limit cache - previous_cache_keys = hamevo._cache_hamiltonian_evo.keys() - vals[vparam] += 0.1 - values_cache_key = str(OrderedDict(vals)) - assert values_cache_key not in previous_cache_keys - - apply_hamevo_and_compare_expected(psi, vals) - assert len(hamevo._cache_hamiltonian_evo) == 2 - assert values_cache_key in previous_cache_keys - - -def test_hevo_constant_gen() -> None: - sup = (0, 1) - generator = pyq.Add( - [pyq.Scale(pyq.Z(0), torch.rand(1)), pyq.Scale(pyq.Z(1), torch.rand(1))] - ) - hamevo = pyq.HamiltonianEvolution(generator, torch.rand(1), sup) - assert hamevo.generator_type == GeneratorType.OPERATION - psi = pyq.zero_state(2) - psi_star = hamevo(psi) - psi_expected = _calc_mat_vec_wavefunction(hamevo, len(sup), psi) - assert torch.allclose(psi_star, psi_expected, rtol=RTOL, atol=ATOL) - - @pytest.mark.parametrize( "H, t_evo, expected_state", [ @@ -285,95 +221,6 @@ def test_hamevo_fixed_tensor_result() -> None: assert torch.allclose(hamevo.tensor(), expected_evo_result, atol=1.0e-4) -@pytest.mark.parametrize("n_qubits", [2, 4, 6]) -def test_hamevo_tensor(n_qubits: int) -> None: - dim = torch.randint(1, n_qubits + 1, (1,)).item() - vparam = "theta" - sup = tuple(range(dim)) - parametric = True - - h = torch.rand(2**dim, 2**dim, dtype=DEFAULT_REAL_DTYPE) - hermitian_matrix = h + torch.conj(torch.transpose(h, 0, 1)) - hamevo = pyq.HamiltonianEvolution(hermitian_matrix, vparam, sup, parametric) - - psi = random_state(n_qubits) - vals = {vparam: torch.rand(1)} - psi_star = hamevo(psi, vals) - psi_expected = _calc_mat_vec_wavefunction(hamevo, n_qubits, psi, vals) - assert torch.allclose(psi_star, psi_expected, rtol=RTOL, atol=ATOL) - - -@pytest.mark.parametrize("n_qubits", [2, 4, 6]) -@pytest.mark.parametrize("dim", [1, 2, 3, 4, 5, 6]) -@pytest.mark.parametrize("same_qubit_case", [True, False]) -def test_hamevo_tensor_from_circuit( - n_qubits: int, dim: int, same_qubit_case: bool -) -> None: - dim = min(n_qubits, dim) - vparam = "theta" - parametric = True - ops = [pyq.X, pyq.Y] * 2 - if same_qubit_case: - qubit_targets = [dim] * len(ops) - else: - qubit_targets = np.random.choice(dim, len(ops), replace=True) - generator = pyq.QuantumCircuit( - n_qubits, - [ - pyq.Add([op(q) for op, q in zip(ops, qubit_targets)]), - *[op(q) for op, q in zip(ops, qubit_targets)], - ], - ) - generator = generator.tensor(n_qubits=n_qubits) - generator = generator + torch.conj(torch.transpose(generator, 0, 1)) - hamevo = pyq.HamiltonianEvolution( - generator, vparam, tuple(range(n_qubits)), parametric - ) - assert hamevo.generator_type == GeneratorType.TENSOR - vals = {vparam: torch.tensor([0.5])} - psi = random_state(n_qubits) - psi_star = hamevo(psi, vals) - psi_expected = _calc_mat_vec_wavefunction(hamevo, n_qubits, psi, vals) - assert torch.allclose(psi_star, psi_expected, rtol=RTOL, atol=ATOL) - - -@pytest.mark.parametrize("n_qubits", [2, 4, 6]) -@pytest.mark.parametrize("dim", [1, 2, 3, 4, 5, 6]) -@pytest.mark.parametrize("same_qubit_case", [True, False]) -def test_hamevo_tensor_from_paramcircuit( - n_qubits: int, dim: int, same_qubit_case: bool -) -> None: - dim = min(n_qubits, dim) - tparam = "theta" - vparam = "rtheta" - parametric = True - ops = [pyq.RX, pyq.RY] * 2 - if same_qubit_case: - qubit_targets = [dim] * len(ops) - else: - qubit_targets = np.random.choice(dim, len(ops), replace=True) - generator = pyq.QuantumCircuit( - n_qubits, - [ - pyq.Add([op(q, vparam) for op, q in zip(ops, qubit_targets)]), - *[op(q, vparam) for op, q in zip(ops, qubit_targets)], - ], - ) - generator = generator.tensor( - n_qubits=n_qubits, values={vparam: torch.tensor([0.5])} - ) - generator = generator + torch.conj(torch.transpose(generator, 0, 1)) - hamevo = pyq.HamiltonianEvolution( - generator, tparam, tuple(range(n_qubits)), parametric - ) - assert hamevo.generator_type == GeneratorType.TENSOR - vals = {tparam: torch.tensor([0.5])} - psi = random_state(n_qubits) - psi_star = hamevo(psi, vals) - psi_expected = _calc_mat_vec_wavefunction(hamevo, n_qubits, psi, vals) - assert torch.allclose(psi_star, psi_expected, rtol=RTOL, atol=ATOL) - - @pytest.mark.parametrize( "state_fn", [pyq.random_state, pyq.zero_state, pyq.uniform_state] ) @@ -461,9 +308,7 @@ def test_timedependent( *[op(q, vparam) for op, q in zip(ops, qubit_targets)], ], ) - generator = generator.tensor( - n_qubits=n_qubits, values={vparam: torch.tensor([0.5])} - ) + generator = generator.tensor(values={vparam: torch.tensor([0.5])}) generator = generator + torch.conj(torch.transpose(generator, 0, 1)) hamevo = pyq.HamiltonianEvolution( generator, tparam, tuple(range(n_qubits)), parametric, is_time_dependent=True @@ -476,6 +321,67 @@ def test_timedependent( embedding = Embedding(vparam_names=[vparam], tparam_name=tparam) psi_star = hamevo(psi, vals, embedding) # TODO vytautas FIX - psi_expected = _calc_mat_vec_wavefunction(hamevo, n_qubits, psi, vals) + psi_expected = calc_mat_vec_wavefunction(hamevo, psi, vals) assert torch.allclose(psi_star, psi_expected, rtol=RTOL, atol=ATOL) + + +@pytest.mark.parametrize("n_qubits", [2, 4, 6]) +@pytest.mark.parametrize("batch_size", [1, 2]) +def test_hamevo_parametric_gen(n_qubits: int, batch_size: int) -> None: + k_1q = 2 * n_qubits # Number of 1-qubit terms + k_2q = n_qubits**2 # Number of 2-qubit terms + generator, param_list = random_pauli_hamiltonian( + n_qubits, k_1q, k_2q, make_param=True + ) + + tparam = "t" + + param_list.append("t") + + hamevo = pyq.HamiltonianEvolution( + generator, tparam, generator_parametric=True, cache_length=2 + ) + + assert hamevo.generator_type == GeneratorType.PARAMETRIC_OPERATION + assert len(hamevo._cache_hamiltonian_evo) == 0 + + def apply_hamevo_and_compare_expected(psi, values): + psi_star = hamevo(psi, values) + psi_expected = calc_mat_vec_wavefunction(hamevo, psi, values) + assert torch.allclose(psi_star, psi_expected, rtol=RTOL, atol=ATOL) + + values = {param: torch.rand(batch_size) for param in param_list} + + psi = random_state(n_qubits) + apply_hamevo_and_compare_expected(psi, values) + + # test cached + assert len(hamevo._cache_hamiltonian_evo) == 1 + apply_hamevo_and_compare_expected(psi, values) + assert len(hamevo._cache_hamiltonian_evo) == 1 + + # test caching new value + for param in param_list: + values[param] += 0.1 + + apply_hamevo_and_compare_expected(psi, values) + assert len(hamevo._cache_hamiltonian_evo) == 2 + + # changing input state should not change the cache + psi = random_state(n_qubits) + apply_hamevo_and_compare_expected(psi, values) + assert len(hamevo._cache_hamiltonian_evo) == 2 + + # test limit cache + previous_cache_keys = hamevo._cache_hamiltonian_evo.keys() + + for param in param_list: + values[param] += 0.1 + + values_cache_key = str(OrderedDict(values)) + assert values_cache_key not in previous_cache_keys + + apply_hamevo_and_compare_expected(psi, values) + assert len(hamevo._cache_hamiltonian_evo) == 2 + assert values_cache_key in previous_cache_keys diff --git a/tests/test_circuit.py b/tests/test_circuit.py index ce2aa622..26095677 100644 --- a/tests/test_circuit.py +++ b/tests/test_circuit.py @@ -4,16 +4,10 @@ import pytest import torch -from torch import Tensor import pyqtorch as pyq from pyqtorch import run, sample -from pyqtorch.circuit import QuantumCircuit -from pyqtorch.noise import Noise -from pyqtorch.parametric import Parametric -from pyqtorch.primitive import Primitive from pyqtorch.utils import ( - DensityMatrix, product_state, ) @@ -26,7 +20,7 @@ def test_device_inference() -> None: @pytest.mark.parametrize("fn", [pyq.X, pyq.Z, pyq.Y]) -def test_scale(fn: pyq.primitive.Primitive) -> None: +def test_scale(fn: pyq.primitives.Primitive) -> None: n_qubits = torch.randint(low=1, high=4, size=(1,)).item() target = random.choice([i for i in range(n_qubits)]) state = pyq.random_state(n_qubits) @@ -99,33 +93,6 @@ def test_merge_nested_dict() -> None: mergecirc(pyq.random_state(2), vals) -@pytest.mark.parametrize("n_qubits", [{"low": 2, "high": 5}], indirect=True) -@pytest.mark.parametrize("batch_size", [{"low": 1, "high": 5}], indirect=True) -def test_noise_circ( - n_qubits: int, - batch_size: int, - random_input_state: Tensor, - random_gate: Primitive, - random_noise_gate: Noise, - random_rotation_gate: Parametric, -) -> None: - OPERATORS = [random_gate, random_noise_gate, random_rotation_gate] - random.shuffle(OPERATORS) - circ = QuantumCircuit(n_qubits, OPERATORS) - - values = {random_rotation_gate.param_name: torch.rand(1)} - output_state = circ(random_input_state, values) - assert isinstance(output_state, DensityMatrix) - assert output_state.shape == torch.Size([2**n_qubits, 2**n_qubits, batch_size]) - - diag_sums = [] - for i in range(batch_size): - diag_batch = torch.diagonal(output_state[:, :, i], dim1=0, dim2=1) - diag_sums.append(torch.sum(diag_batch)) - diag_sum = torch.stack(diag_sums) - assert torch.allclose(diag_sum, torch.ones((batch_size,), dtype=torch.cdouble)) - - def test_sample_run() -> None: ops = [pyq.X(0), pyq.X(1)] circ = pyq.QuantumCircuit(4, ops) diff --git a/tests/test_differentiation.py b/tests/test_differentiation.py index 05324751..2be244ad 100644 --- a/tests/test_differentiation.py +++ b/tests/test_differentiation.py @@ -6,13 +6,12 @@ import pyqtorch as pyq from pyqtorch import DiffMode, expectation from pyqtorch.matrices import COMPLEX_TO_REAL_DTYPES -from pyqtorch.primitive import Primitive +from pyqtorch.primitives import Primitive from pyqtorch.utils import ( GRADCHECK_ATOL, ) -# TODO add GPSR when multigap is implemented for this test @pytest.mark.parametrize("n_qubits", [3, 4, 5]) @pytest.mark.parametrize("n_layers", [1, 2, 3]) def test_adjoint_diff(n_qubits: int, n_layers: int) -> None: @@ -22,7 +21,7 @@ def test_adjoint_diff(n_qubits: int, n_layers: int) -> None: cnot = pyq.CNOT(1, 2) ops = [rx, cry, rz, cnot] * n_layers circ = pyq.QuantumCircuit(n_qubits, ops) - obs = pyq.Observable(n_qubits, [pyq.Z(0)]) + obs = pyq.Observable(pyq.Z(0)) theta_0_value = torch.pi / 2 theta_1_value = torch.pi @@ -60,77 +59,59 @@ def test_adjoint_diff(n_qubits: int, n_layers: int) -> None: for i in range(len(grad_ad)): assert torch.allclose(grad_ad[i], grad_adjoint[i], atol=GRADCHECK_ATOL) + # TODO higher order adjoint is not yet supported. + # gradgrad_adjoint = torch.autograd.grad( + # grad_adjoint, tuple(values_adjoint.values()), torch.ones_like(grad_adjoint) + # ) + +@pytest.mark.parametrize("n_qubits", [3, 5]) +@pytest.mark.parametrize("batch_size", [1, 2]) +@pytest.mark.parametrize("ops_op", [pyq.Z, pyq.Y]) @pytest.mark.parametrize("dtype", [torch.complex64, torch.complex128]) -@pytest.mark.parametrize("batch_size", [1, 5]) -@pytest.mark.parametrize("n_qubits", [3, 4]) -def test_differentiate_circuit( - dtype: torch.dtype, batch_size: int, n_qubits: int +def test_sampled_diff( + n_qubits: int, + batch_size: int, + ops_op: Primitive, + dtype: torch.dtype, ) -> None: - ops = [ - pyq.Y(1), - pyq.RX(0, "theta_0"), - pyq.PHASE(0, "theta_1"), - pyq.CSWAP(0, (1, 2)), - pyq.CRX(1, 2, "theta_2"), - pyq.CPHASE(1, 2, "theta_3"), - pyq.CNOT(0, 1), - pyq.Toffoli((2, 1), 0), - ] - theta_0_value = torch.rand(1, dtype=dtype) - theta_1_value = torch.rand(1, dtype=dtype) - theta_2_value = torch.rand(1, dtype=dtype) - theta_3_value = torch.rand(1, dtype=dtype) + rx = pyq.RX(0, param_name="theta_0") + cry = pyq.CPHASE(0, 1, param_name="theta_1") + rz = pyq.RZ(2, param_name="theta_2") + cnot = pyq.CNOT(1, 2) + ops = [rx, cry, rz, cnot] circ = pyq.QuantumCircuit(n_qubits, ops).to(dtype) - state = pyq.random_state(n_qubits, batch_size, dtype=dtype) - - theta_0_ad = torch.tensor([theta_0_value], requires_grad=True) - theta_0_adjoint = torch.tensor([theta_0_value], requires_grad=True) + obs = pyq.Observable(pyq.Add([ops_op(i) for i in range(n_qubits)])).to(dtype) - theta_1_ad = torch.tensor([theta_1_value], requires_grad=True) - theta_1_adjoint = torch.tensor([theta_1_value], requires_grad=True) - - theta_2_ad = torch.tensor([theta_2_value], requires_grad=True) - theta_2_adjoint = torch.tensor([theta_2_value], requires_grad=True) - - theta_3_ad = torch.tensor([theta_3_value], requires_grad=True) - theta_3_adjoint = torch.tensor([theta_3_value], requires_grad=True) - - values_ad = torch.nn.ParameterDict( - { - "theta_0": theta_0_ad, - "theta_1": theta_1_ad, - "theta_2": theta_2_ad, - "theta_3": theta_3_ad, - } - ).to(COMPLEX_TO_REAL_DTYPES[dtype]) - values_adjoint = torch.nn.ParameterDict( - { - "theta_0": theta_0_adjoint, - "theta_1": theta_1_adjoint, - "theta_2": theta_2_adjoint, - "theta_3": theta_3_adjoint, - } - ).to(COMPLEX_TO_REAL_DTYPES[dtype]) - - obs = pyq.QuantumCircuit(n_qubits, [pyq.Z(0)]).to(dtype) - exp_ad = expectation(circ, state, values_ad, obs, DiffMode.AD) - exp_adjoint = expectation(circ, state, values_adjoint, obs, DiffMode.ADJOINT) + theta_0_value = torch.pi / 2 + theta_1_value = torch.pi + theta_2_value = torch.pi / 4 - grad_ad = torch.autograd.grad( - exp_ad, tuple(values_ad.values()), torch.ones_like(exp_ad) + state = pyq.random_state(n_qubits, batch_size=batch_size, dtype=dtype) + theta_0_ad = torch.tensor( + [theta_0_value], requires_grad=True, dtype=COMPLEX_TO_REAL_DTYPES[dtype] ) - - grad_adjoint = torch.autograd.grad( - exp_adjoint, tuple(values_adjoint.values()), torch.ones_like(exp_adjoint) + theta_1_ad = torch.tensor( + [theta_1_value], requires_grad=True, dtype=COMPLEX_TO_REAL_DTYPES[dtype] + ) + theta_2_ad = torch.tensor( + [theta_2_value], requires_grad=True, dtype=COMPLEX_TO_REAL_DTYPES[dtype] ) + values = {"theta_0": theta_0_ad, "theta_1": theta_1_ad, "theta_2": theta_2_ad} - assert len(grad_ad) == len(grad_adjoint) - for i in range(len(grad_ad)): - assert torch.allclose(grad_ad[i], grad_adjoint[i], atol=GRADCHECK_ATOL) + exp_ad = expectation(circ, state, values, obs, DiffMode.AD) + exp_ad_sampled = expectation( + circ, + state, + values, + obs, + DiffMode.AD, + n_shots=100000, + ) + assert torch.allclose(exp_ad, exp_ad_sampled, atol=1e-01) -@pytest.mark.xfail # investigate +@pytest.mark.xfail # Adjoint Scale is currently not supported @pytest.mark.parametrize("dtype", [torch.complex64, torch.complex128]) @pytest.mark.parametrize("batch_size", [1, 5]) @pytest.mark.parametrize("n_qubits", [2]) @@ -156,7 +137,7 @@ def test_adjoint_scale(dtype: torch.dtype, batch_size: int, n_qubits: int) -> No } ).to(COMPLEX_TO_REAL_DTYPES[dtype]) - obs = pyq.Observable(n_qubits, [pyq.Z(0)]).to(dtype) + obs = pyq.Observable(pyq.Z(0)).to(dtype) exp_ad = expectation(circ, state, values_ad, obs, DiffMode.AD) exp_adjoint = expectation(circ, state, values_adjoint, obs, DiffMode.ADJOINT) @@ -193,7 +174,7 @@ def test_all_diff_singlegap( ops = [ops_rx, ops_rz, cnot] circ = pyq.QuantumCircuit(n_qubits, ops).to(dtype) - obs = pyq.Observable(n_qubits, [op_obs(i) for i in range(n_qubits)]).to(dtype) + obs = pyq.Observable([op_obs(i) for i in range(n_qubits)]).to(dtype) state = pyq.random_state(n_qubits, batch_size, dtype=dtype) dtype_val = COMPLEX_TO_REAL_DTYPES[dtype] @@ -264,43 +245,3 @@ def test_all_diff_singlegap( # check second order gradients for j in range(len(gradgrad_ad)): assert torch.allclose(gradgrad_ad[j], gradgrad_gpsr[j], atol=GRADCHECK_ATOL) - - -@pytest.mark.parametrize("gate_type", ["scale", "hamevo", "same", ""]) -def test_compatibility_gpsr(gate_type: str) -> None: - - pname = "theta_0" - if gate_type == "scale": - seq_gate = pyq.Sequence([pyq.X(0)]) - scale = pyq.Scale(seq_gate, pname) - ops = [scale] - elif gate_type == "hamevo": - hamevo = pyq.HamiltonianEvolution(pyq.Sequence([pyq.X(0)]), pname, (0,)) - ops = [hamevo] - elif gate_type == "same": - ops = [pyq.RY(0, pname), pyq.RZ(0, pname)] - else: - # check that CNOT is not tested on spectral gap call - ops = [pyq.RY(0, pname), pyq.CNOT(0, 1)] - - circ = pyq.QuantumCircuit(2, ops) - obs = pyq.Observable(2, [pyq.Z(0)]) - state = pyq.zero_state(2) - - param_value = torch.pi / 2 - values = {"theta_0": torch.tensor([param_value], requires_grad=True)} - - if gate_type != "": - with pytest.raises(ValueError): - exp_gpsr = expectation(circ, state, values, obs, DiffMode.GPSR) - - grad_gpsr = torch.autograd.grad( - exp_gpsr, tuple(values.values()), torch.ones_like(exp_gpsr) - ) - else: - exp_gpsr = expectation(circ, state, values, obs, DiffMode.GPSR) - - grad_gpsr = torch.autograd.grad( - exp_gpsr, tuple(values.values()), torch.ones_like(exp_gpsr) - ) - assert len(grad_gpsr) > 0 diff --git a/tests/test_digital.py b/tests/test_digital.py index 4e29b752..0605172b 100644 --- a/tests/test_digital.py +++ b/tests/test_digital.py @@ -6,41 +6,17 @@ import pytest import torch -from conftest import _calc_mat_vec_wavefunction -from numpy import array, rollaxis, sort -from qutip import Qobj -from torch import Tensor, moveaxis +from torch import Tensor import pyqtorch as pyq -from pyqtorch.apply import apply_operator, operator_product +from pyqtorch.apply import apply_operator from pyqtorch.matrices import ( DEFAULT_MATRIX_DTYPE, - HMAT, - IMAT, - XMAT, - YMAT, - ZMAT, - _dagger, ) -from pyqtorch.noise import ( - AmplitudeDamping, - GeneralizedAmplitudeDamping, - Noise, - PhaseDamping, -) -from pyqtorch.parametric import Parametric -from pyqtorch.primitive import H, I, Primitive, S, T, X, Y, Z +from pyqtorch.primitives import Parametric, Primitive from pyqtorch.utils import ( ATOL, - RTOL, - DensityMatrix, - density_mat, - dm_partial_trace, - generate_dm, - operator_kron, product_state, - promote_operator, - random_state, ) state_000 = product_state("000") @@ -59,7 +35,7 @@ def test_identity() -> None: assert torch.allclose(product_state("0"), pyq.I(0)(product_state("0"), None)) - assert torch.allclose(product_state("1"), pyq.I(1)(product_state("1"))) + assert torch.allclose(product_state("1"), pyq.I(0)(product_state("1"))) def test_N() -> None: @@ -68,42 +44,6 @@ def test_N() -> None: assert torch.allclose(product_state("1"), pyq.N(0)(product_state("1"), None)) -@pytest.mark.parametrize("gate", [I, X, Y, Z, H, T, S]) -@pytest.mark.parametrize("n_qubits", [1, 2, 4]) -def test_single_qubit_gates(gate: Primitive, n_qubits: int) -> None: - target = torch.randint(0, n_qubits, (1,)).item() - block = gate(target) - init_state = pyq.random_state(n_qubits) - wf_pyq = block(init_state, None) - wf_mat = _calc_mat_vec_wavefunction(block, n_qubits, init_state) - assert torch.allclose(wf_mat, wf_pyq, rtol=RTOL, atol=ATOL) - - -@pytest.mark.parametrize("batch_size", [i for i in range(2, 10)]) -@pytest.mark.parametrize("gate", [pyq.RX, pyq.RY, pyq.RZ, pyq.PHASE]) -@pytest.mark.parametrize("n_qubits", [1, 2, 4]) -def test_rotation_gates(batch_size: int, gate: Primitive, n_qubits: int) -> None: - params = [f"th{i}" for i in range(gate.n_params)] - values = {param: torch.rand(batch_size) for param in params} - target = torch.randint(0, n_qubits, (1,)).item() - - init_state = pyq.random_state(n_qubits) - block = gate(target, *params) - wf_pyq = block(init_state, values) - wf_mat = _calc_mat_vec_wavefunction(block, n_qubits, init_state, values=values) - assert torch.allclose(wf_mat, wf_pyq, rtol=RTOL, atol=ATOL) - assert block.spectral_gap == 2.0 - if not isinstance(block, pyq.PHASE): - assert torch.allclose( - block.eigenvals_generator, - torch.tensor([[-1.0], [1.0]], dtype=block.dtype), - ) - else: - assert torch.allclose( - block.eigenvals_generator, torch.tensor([[0.0], [2.0]], dtype=block.dtype) - ) - - @pytest.mark.parametrize("n_qubits", [1, 2, 3]) def test_projectors(n_qubits: int) -> None: qubit_support = tuple(range(n_qubits)) @@ -129,47 +69,6 @@ def test_projectors(n_qubits: int) -> None: ) -@pytest.mark.parametrize( - "projector, exp_projector_mat", - [ - ( - pyq.Projector(0, bra="1", ket="1"), - torch.tensor( - [[0.0 + 0.0j, 0.0 + 0.0j], [0.0 + 0.0j, 1.0 + 0.0j]], - dtype=torch.complex128, - ), - ), - ( - pyq.N(0), - (IMAT - ZMAT) / 2.0, - ), - ( - pyq.CNOT(0, 1), - torch.tensor( - [ - [ - [1.0 + 0.0j, 0.0 + 0.0j, 0.0 + 0.0j, 0.0 + 0.0j], - [0.0 + 0.0j, 1.0 + 0.0j, 0.0 + 0.0j, 0.0 + 0.0j], - [0.0 + 0.0j, 0.0 + 0.0j, 0.0 + 0.0j, 1.0 + 0.0j], - [0.0 + 0.0j, 0.0 + 0.0j, 1.0 + 0.0j, 0.0 + 0.0j], - ] - ], - dtype=torch.complex128, - ), - ), - ], -) -def test_projector_tensor( - projector: Primitive, exp_projector_mat: torch.Tensor -) -> None: - - nbqubits = int(log2(exp_projector_mat.shape[-1])) - projector_mat = projector.tensor( - n_qubits=nbqubits, values={"theta": torch.Tensor([1.0])} - ).squeeze(-1) - assert torch.allclose(projector_mat, exp_projector_mat, atol=1.0e-4) - - czop_example = pyq.CZ(control=(0, 1), target=2) crxop_example = pyq.CRX(control=(0, 1), target=2, param_name="theta") @@ -313,7 +212,7 @@ def test_multi_controlled_gates( initial_state = initial_state.to(device=device, dtype=dtype) rot_gate = getattr(pyq, gate) controlled_rot_gate = getattr(pyq, "C" + gate) - phi = torch.rand(batch_size) + phi = torch.rand(batch_size).to(device=device) n_qubits = int(log2(torch.numel(initial_state))) qubits = tuple([i for i in range(n_qubits)]) op = controlled_rot_gate(qubits[:-1], qubits[-1], "phi").to( @@ -378,7 +277,7 @@ def test_dagger_single_qubit() -> None: values = ( {param_name: torch.rand(1)} if param_name == "theta" else torch.rand(1) ) - new_state = apply_operator(state, op.unitary(values), [target]) + new_state = apply_operator(state, op.tensor(values), [target]) daggered_back = apply_operator(new_state, op.dagger(values), [target]) assert torch.allclose(daggered_back, state) @@ -411,7 +310,7 @@ def test_dagger_nqubit() -> None: values = ( {param_name: torch.rand(1)} if param_name == "theta" else torch.rand(1) ) - new_state = apply_operator(state, op.unitary(values), qubit_support) + new_state = apply_operator(state, op.tensor(values), qubit_support) daggered_back = apply_operator(new_state, op.dagger(values), qubit_support) assert torch.allclose(daggered_back, state) @@ -420,7 +319,7 @@ def test_U() -> None: n_qubits = torch.randint(low=1, high=8, size=(1,)).item() target = random.choice([i for i in range(n_qubits)]) params = ["phi", "theta", "omega"] - u = pyq.U(target, *params) + u = pyq.U(target, phi=params[0], theta=params[1], omega=params[2]) values = {param: torch.rand(1) for param in params} state = pyq.random_state(n_qubits) assert torch.allclose( @@ -430,209 +329,8 @@ def test_U() -> None: assert u.spectral_gap == 2.0 -def test_dm(n_qubits: int, batch_size: int) -> None: - state = random_state(n_qubits) - projector = torch.outer(state.flatten(), state.conj().flatten()).view( - 2**n_qubits, 2**n_qubits, 1 - ) - dm = density_mat(state) - assert dm.size() == torch.Size([2**n_qubits, 2**n_qubits, 1]) - assert torch.allclose(dm, projector) - assert torch.allclose(dm.squeeze(), dm.squeeze() @ dm.squeeze()) - states = [] - projectors = [] - for batch in range(batch_size): - state = random_state(n_qubits) - states.append(state) - projector = torch.outer(state.flatten(), state.conj().flatten()).view( - 2**n_qubits, 2**n_qubits, 1 - ) - projectors.append(projector) - dm_proj = torch.cat(projectors, dim=2) - state_cat = torch.cat(states, dim=n_qubits) - dm = density_mat(state_cat) - assert dm.size() == torch.Size([2**n_qubits, 2**n_qubits, batch_size]) - assert torch.allclose(dm, dm_proj) - - -def test_promote(random_gate: Primitive, n_qubits: int, target: int) -> None: - op_prom = promote_operator(random_gate.unitary(), target, n_qubits) - assert op_prom.size() == torch.Size([2**n_qubits, 2**n_qubits, 1]) - assert torch.allclose( - operator_product(op_prom, _dagger(op_prom), target), - torch.eye(2**n_qubits, dtype=torch.cdouble).unsqueeze(2), - ) - - -def test_operator_product(random_gate: Primitive, n_qubits: int, target: int) -> None: - op = random_gate - batch_size_1 = torch.randint(low=1, high=5, size=(1,)).item() - batch_size_2 = torch.randint(low=1, high=5, size=(1,)).item() - max_batch = max(batch_size_2, batch_size_1) - op_prom = promote_operator(op.unitary(), target, n_qubits).repeat( - 1, 1, batch_size_1 - ) - op_mul = operator_product( - op.unitary().repeat(1, 1, batch_size_2), _dagger(op_prom), target - ) - assert op_mul.size() == torch.Size([2**n_qubits, 2**n_qubits, max_batch]) - assert torch.allclose( - op_mul, - torch.eye(2**n_qubits, dtype=torch.cdouble) - .unsqueeze(2) - .repeat(1, 1, max_batch), - ) - - -@pytest.mark.parametrize( - "operator,matrix", [(I, IMAT), (X, XMAT), (Z, ZMAT), (Y, YMAT), (H, HMAT)] -) -def test_operator_kron(operator: Tensor, matrix: Tensor) -> None: - n_qubits = torch.randint(low=1, high=5, size=(1,)).item() - batch_size = torch.randint(low=1, high=5, size=(1,)).item() - states, krons = [], [] - for batch in range(batch_size): - state = random_state(n_qubits) - states.append(state) - kron = torch.kron(density_mat(state).squeeze(), matrix).unsqueeze(2) - krons.append(kron) - input_state = torch.cat(states, dim=n_qubits) - kron_out = operator_kron(density_mat(input_state), operator(0).dagger()) - assert kron_out.size() == torch.Size( - [2 ** (n_qubits + 1), 2 ** (n_qubits + 1), batch_size] - ) - kron_expect = torch.cat(krons, dim=2) - assert torch.allclose(kron_out, kron_expect) - assert torch.allclose( - torch.kron(operator(0).dagger().contiguous(), I(0).unitary()), - operator_kron(operator(0).dagger(), I(0).unitary()), - ) - - -def test_kron_batch() -> None: - n_qubits = torch.randint(low=1, high=5, size=(1,)).item() - batch_size_1 = torch.randint(low=1, high=5, size=(1,)).item() - batch_size_2 = torch.randint(low=1, high=5, size=(1,)).item() - max_batch = max(batch_size_1, batch_size_2) - dm_1 = density_mat(random_state(n_qubits, batch_size_1)) - dm_2 = density_mat(random_state(n_qubits, batch_size_2)) - dm_out = operator_kron(dm_1, dm_2) - assert dm_out.size() == torch.Size( - [2 ** (2 * n_qubits), 2 ** (2 * n_qubits), max_batch] - ) - if batch_size_1 > batch_size_2: - dm_2 = dm_2.repeat(1, 1, batch_size_1)[:, :, :batch_size_1] - elif batch_size_2 > batch_size_1: - dm_1 = dm_1.repeat(1, 1, batch_size_2)[:, :, :batch_size_2] - density_matrices = [] - for batch in range(max_batch): - density_matrice = torch.kron(dm_1[:, :, batch], dm_2[:, :, batch]).unsqueeze(2) - density_matrices.append(density_matrice) - dm_expect = torch.cat(density_matrices, dim=2) - assert torch.allclose(dm_out, dm_expect) - - -def test_circuit_tensor() -> None: - ops = [pyq.RX(0, "theta_0"), pyq.RY(0, "theta_1"), pyq.RX(1, "theta_2")] - circ = pyq.QuantumCircuit(2, ops) - values = {f"theta_{i}": torch.Tensor([float(i)]) for i in range(3)} - tensorcirc = circ.tensor(values) - assert tensorcirc.size() == (4, 4, 1) - assert torch.allclose( - tensorcirc, - torch.tensor( - [ - [ - [0.4742 + 0.0000j], - [0.0000 - 0.7385j], - [-0.2590 + 0.0000j], - [0.0000 + 0.4034j], - ], - [ - [0.0000 - 0.7385j], - [0.4742 + 0.0000j], - [0.0000 + 0.4034j], - [-0.2590 + 0.0000j], - ], - [ - [0.2590 + 0.0000j], - [0.0000 - 0.4034j], - [0.4742 + 0.0000j], - [0.0000 - 0.7385j], - ], - [ - [0.0000 - 0.4034j], - [0.2590 + 0.0000j], - [0.0000 - 0.7385j], - [0.4742 + 0.0000j], - ], - ], - dtype=torch.complex128, - ), - atol=1.0e-4, - ) - - -@pytest.mark.parametrize("n_qubits", [{"low": 2, "high": 5}], indirect=True) -def test_flip_gates( - n_qubits: int, - target: int, - batch_size: int, - rho_input: Tensor, - random_flip_gate: Noise, - flip_expected_state: DensityMatrix, - flip_probability: Tensor | float, - flip_gates_prob_0: Noise, - flip_gates_prob_1: tuple, - random_input_state: Tensor, -) -> None: - FlipGate = random_flip_gate - output_state: DensityMatrix = FlipGate(target, flip_probability)(rho_input) - assert output_state.size() == torch.Size([2**n_qubits, 2**n_qubits, batch_size]) - assert torch.allclose(output_state, flip_expected_state) - - input_state = random_input_state # fix the same random state for every call - assert torch.allclose( - flip_gates_prob_0(density_mat(input_state)), density_mat(input_state) - ) - - FlipGate_1, expected_op = flip_gates_prob_1 - assert torch.allclose(FlipGate_1(density_mat(input_state)), expected_op) - - -def test_damping_gates( - n_qubits: int, - target: int, - batch_size: int, - random_damping_gate: Noise, - damping_expected_state: tuple, - damping_rate: Tensor, - damping_gates_prob_0: Tensor, - random_input_state: Tensor, - rho_input: Tensor, -) -> None: - DampingGate, expected_state = damping_expected_state - apply_gate = DampingGate(rho_input) - assert apply_gate.size() == torch.Size([2**n_qubits, 2**n_qubits, batch_size]) - assert torch.allclose(apply_gate, expected_state) - - input_state = random_input_state - assert torch.allclose( - damping_gates_prob_0(input_state), density_mat(I(target)(input_state)) - ) - - rho_0: DensityMatrix = density_mat(product_state("0", batch_size)) - rho_1: DensityMatrix = density_mat(product_state("1", batch_size)) - if DampingGate == AmplitudeDamping: - assert torch.allclose(DampingGate(target, rate=1)(rho_1), rho_0) - elif DampingGate == PhaseDamping: - assert torch.allclose(DampingGate(target, rate=1)(rho_1), I(target)(rho_1)) - elif DampingGate == GeneralizedAmplitudeDamping: - assert torch.allclose(DampingGate(target, probability=1, rate=1)(rho_1), rho_0) - - @pytest.mark.parametrize("gate", [pyq.RX, pyq.RY, pyq.RZ]) -def test_parametric_constantparam(gate: pyq.parametric.Parametric) -> None: +def test_parametric_constantparam(gate: Parametric) -> None: n_qubits = 4 max_batch_size = 10 target = torch.randint(0, n_qubits, (1,)).item() @@ -642,39 +340,3 @@ def test_parametric_constantparam(gate: pyq.parametric.Parametric) -> None: gate(target, "theta")(state, {"theta": param_val}), gate(target, param_val)(state), ) - - -def test_dm_partial_trace() -> None: - n_qubits = torch.randint(low=1, high=7, size=(1,)).item() - batch_size = torch.randint(low=1, high=5, size=(1,)).item() - state_str = "".join(random.choice("01") for _ in range(n_qubits)) - - # testing swaps - rho = density_mat(product_state(state_str, batch_size=batch_size)) - keep_indices = random.sample(range(n_qubits), k=random.randint(1, n_qubits)) - n_keep = len(keep_indices) - state_reduce_str = "".join([state_str[i] for i in keep_indices]) - rho_reduce = dm_partial_trace(rho, keep_indices) - - assert rho_reduce.shape == torch.Size([2**n_keep, 2**n_keep, batch_size]) - assert torch.allclose( - rho_reduce, - density_mat(product_state(state_reduce_str, batch_size=batch_size)), - ) - - # testing reduced density matrix - rho_list = generate_dm(n_qubits, batch_size) - - rho_sub = torch.from_numpy( - array( - [ - Qobj( - rollaxis(rho_list.numpy(), 2, 0)[i], dims=[[2] * n_qubits] * 2 - ).ptrace(sort(keep_indices)) - for i in range(batch_size) - ] - ) - ) - rho_sub = moveaxis(rho_sub, 0, 2).type(torch.cfloat) - - assert torch.allclose(dm_partial_trace(rho_list, sort(keep_indices)), rho_sub) diff --git a/tests/test_dropout.py b/tests/test_dropout.py new file mode 100644 index 00000000..e51d31c3 --- /dev/null +++ b/tests/test_dropout.py @@ -0,0 +1,49 @@ +from __future__ import annotations + +import pytest +import torch + +import pyqtorch as pyq +from pyqtorch.utils import ATOL, DropoutMode + + +@pytest.mark.parametrize( + "dropout_mode", + [DropoutMode.ROTATIONAL, DropoutMode.ENTANGLING, DropoutMode.CANONICAL_FWD], +) +def test_dropout( + dropout_mode: DropoutMode, n_qubits: int = 4, theta: float = 0.6 +) -> None: + """Test when dropout prob = 1.0 if all gates the dropout method affects are dropped.""" + circ_ops_1 = [pyq.RY(i, f"theta_{str(i)}") for i in range(n_qubits)] + circ_ops_2 = [pyq.CNOT(i, i + 1) for i in range(n_qubits - 1)] + + drop_circ_ops = circ_ops_1 + circ_ops_2 + + operations = { + DropoutMode.ROTATIONAL: circ_ops_2, + DropoutMode.ENTANGLING: circ_ops_1, + DropoutMode.CANONICAL_FWD: [], + } + + values = { + f"theta_{str(i)}": torch.tensor([theta], requires_grad=True) + for i in range(n_qubits) + } + + state = pyq.zero_state(n_qubits=n_qubits) + circ = pyq.QuantumCircuit(n_qubits=n_qubits, operations=operations[dropout_mode]) + dropout_circ = pyq.DropoutQuantumCircuit( + n_qubits=n_qubits, + operations=drop_circ_ops, + dropout_mode=dropout_mode, + dropout_prob=1.0, + ) + obs = pyq.Observable([pyq.Z(1)]) + + exp_circ = pyq.expectation(circuit=circ, state=state, values=values, observable=obs) + exp_dropout_circ = pyq.expectation( + circuit=dropout_circ, state=state, values=values, observable=obs + ) + + assert torch.allclose(exp_circ, exp_dropout_circ, atol=ATOL) diff --git a/tests/test_embedding.py b/tests/test_embedding.py index 1c4e1f5d..f4015f0c 100644 --- a/tests/test_embedding.py +++ b/tests/test_embedding.py @@ -11,7 +11,7 @@ import pyqtorch as pyq from pyqtorch.embed import ConcretizedCallable, Embedding -from pyqtorch.primitive import Primitive +from pyqtorch.primitives import Primitive from pyqtorch.utils import ATOL_embedding @@ -145,7 +145,7 @@ def test_sample_run_expectation_grads_with_embedding(diff_mode) -> None: ops = [rx, cry, phase, ry, cnot] n_qubits = 3 circ = pyq.QuantumCircuit(n_qubits, ops) - obs = pyq.Observable(n_qubits, [pyq.Z(0)]) + obs = pyq.Observable(pyq.Z(0)) state = pyq.zero_state(n_qubits) @@ -156,10 +156,12 @@ def test_sample_run_expectation_grads_with_embedding(diff_mode) -> None: embedded_params = embedding(values_ad) wf = pyq.run(circ, state, embedded_params, embedding) samples = pyq.sample(circ, state, embedded_params, 100, embedding) - exp_ad = pyq.expectation(circ, state, embedded_params, obs, diff_mode, embedding) + exp_ad = pyq.expectation( + circ, state, embedded_params, obs, diff_mode, embedding=embedding + ) assert torch.autograd.gradcheck( lambda x, y: pyq.expectation( - circ, state, {"x": x, "y": y}, obs, diff_mode, embedding + circ, state, {"x": x, "y": y}, obs, diff_mode, embedding=embedding ), (x, y), atol=1e-1, # torch.autograd.gradcheck is very susceptible to small numerical errors diff --git a/tests/test_gpsr.py b/tests/test_gpsr.py new file mode 100644 index 00000000..e6d14e73 --- /dev/null +++ b/tests/test_gpsr.py @@ -0,0 +1,205 @@ +from __future__ import annotations + +from typing import Callable + +import pytest +import torch +from helpers import random_pauli_hamiltonian + +import pyqtorch as pyq +from pyqtorch import DiffMode, expectation +from pyqtorch.circuit import QuantumCircuit +from pyqtorch.hamiltonians import Observable +from pyqtorch.matrices import COMPLEX_TO_REAL_DTYPES +from pyqtorch.primitives import Parametric +from pyqtorch.utils import GPSR_ACCEPTANCE, PSR_ACCEPTANCE, GRADCHECK_sampling_ATOL + + +def circuit_psr(n_qubits: int) -> QuantumCircuit: + """Helper function to make an example circuit using single gap PSR.""" + + ops = [ + pyq.RX(0, "x"), + pyq.RY(1, "y"), + pyq.RX(0, "theta"), + pyq.RY(1, torch.pi / 2), + pyq.CNOT(0, 1), + ] + + circ = QuantumCircuit(n_qubits, ops) + + return circ + + +def circuit_gpsr(n_qubits: int) -> QuantumCircuit: + """Helper function to make an example circuit using multi gap GPSR.""" + + ops = [ + pyq.Y(1), + pyq.RX(0, "theta_0"), + pyq.PHASE(0, "theta_1"), + pyq.CSWAP(0, (1, 2)), + pyq.CRX(1, 2, "theta_2"), + pyq.CPHASE(1, 2, "theta_3"), + pyq.CNOT(0, 1), + pyq.Toffoli((0, 1), 2), + ] + + circ = QuantumCircuit(n_qubits, ops) + + return circ + + +def circuit_sequence(n_qubits: int) -> QuantumCircuit: + """Helper function to make an example circuit using Sequences of rotations.""" + name_angles = "theta" + + ops_rx = pyq.Sequence( + [pyq.RX(i, param_name=name_angles + "_x_" + str(i)) for i in range(n_qubits)] + ) + ops_rz = pyq.Sequence( + [pyq.RZ(i, param_name=name_angles + "_z_" + str(i)) for i in range(n_qubits)] + ) + cnot = pyq.CNOT(1, 2) + ops = [ops_rx, ops_rz, cnot] + circ = QuantumCircuit(n_qubits, ops) + return circ + + +@pytest.mark.parametrize( + ["n_qubits", "batch_size", "circuit_fn"], + [ + (2, 1, circuit_psr), + (5, 10, circuit_psr), + (3, 1, circuit_gpsr), + (5, 10, circuit_gpsr), + (3, 1, circuit_sequence), + (5, 10, circuit_sequence), + ], +) +@pytest.mark.parametrize("dtype", [torch.complex64, torch.complex128]) +def test_expectation_gpsr( + n_qubits: int, + batch_size: int, + circuit_fn: Callable, + dtype: torch.dtype, +) -> None: + torch.manual_seed(42) + circ = circuit_fn(n_qubits).to(dtype) + obs = Observable( + random_pauli_hamiltonian( + n_qubits, k_1q=n_qubits, k_2q=0, default_scale_coeffs=1.0 + )[0] + ).to(dtype) + values = { + op.param_name: torch.rand( + batch_size, requires_grad=True, dtype=COMPLEX_TO_REAL_DTYPES[dtype] + ) + for op in circ.flatten() + if isinstance(op, Parametric) and isinstance(op.param_name, str) + } + state = pyq.random_state(n_qubits, dtype=dtype) + + # Apply adjoint + exp_ad = expectation(circ, state, values, obs, DiffMode.AD) + grad_ad = torch.autograd.grad( + exp_ad, tuple(values.values()), torch.ones_like(exp_ad), create_graph=True + ) + + # Apply PSR + exp_gpsr = expectation(circ, state, values, obs, DiffMode.GPSR) + grad_gpsr = torch.autograd.grad( + exp_gpsr, tuple(values.values()), torch.ones_like(exp_gpsr), create_graph=True + ) + + exp_gpsr_sampled = expectation( + circ, + state, + values, + obs, + DiffMode.GPSR, + n_shots=100000, + ) + grad_gpsr_sampled = torch.autograd.grad( + exp_gpsr_sampled, + tuple(values.values()), + torch.ones_like(exp_gpsr_sampled), + create_graph=True, + ) + assert torch.allclose(exp_gpsr, exp_gpsr_sampled, atol=1e-01) + + atol = PSR_ACCEPTANCE if circuit_fn != circuit_gpsr else GPSR_ACCEPTANCE + + # first order checks + + for i in range(len(grad_ad)): + assert torch.allclose(grad_ad[i], grad_gpsr[i], atol=atol) + assert torch.allclose( + grad_gpsr[i], grad_gpsr_sampled[i], atol=GRADCHECK_sampling_ATOL + ) + + # second order checks + for i in range(len(grad_ad)): + gradgrad_ad = torch.autograd.grad( + grad_ad[i], + tuple(values.values()), + torch.ones_like(grad_ad[i]), + create_graph=True, + ) + + gradgrad_gpsr = torch.autograd.grad( + grad_gpsr[i], + tuple(values.values()), + torch.ones_like(grad_gpsr[i]), + create_graph=True, + ) + + assert len(gradgrad_ad) == len(gradgrad_gpsr) + + # check second order gradients + for j in range(len(gradgrad_ad)): + assert torch.allclose(gradgrad_ad[j], gradgrad_gpsr[j], atol=atol) + + +@pytest.mark.parametrize("gate_type", ["scale", "hamevo", "same", ""]) +@pytest.mark.parametrize("sequence_circuit", [True, False]) +def test_compatibility_gpsr(gate_type: str, sequence_circuit: bool) -> None: + + pname = "theta_0" + if gate_type == "scale": + seq_gate = pyq.Sequence([pyq.X(0)]) + scale = pyq.Scale(seq_gate, pname) + ops = [scale] + elif gate_type == "hamevo": + hamevo = pyq.HamiltonianEvolution(pyq.Sequence([pyq.X(0)]), pname, (0,)) + ops = [hamevo] + elif gate_type == "same": + ops = [pyq.RY(0, pname), pyq.RZ(0, pname)] + else: + # check that CNOT is not tested on spectral gap call + ops = [pyq.RY(0, pname), pyq.CNOT(0, 1)] + + if sequence_circuit: + circ = pyq.QuantumCircuit(2, pyq.Sequence(ops)) + else: + circ = pyq.QuantumCircuit(2, ops) + obs = pyq.Observable([pyq.Z(0)]) + state = pyq.zero_state(2) + + param_value = torch.pi / 2 + values = {"theta_0": torch.tensor([param_value], requires_grad=True)} + + if gate_type != "": + with pytest.raises(ValueError): + exp_gpsr = expectation(circ, state, values, obs, DiffMode.GPSR) + + grad_gpsr = torch.autograd.grad( + exp_gpsr, tuple(values.values()), torch.ones_like(exp_gpsr) + ) + else: + exp_gpsr = expectation(circ, state, values, obs, DiffMode.GPSR) + + grad_gpsr = torch.autograd.grad( + exp_gpsr, tuple(values.values()), torch.ones_like(exp_gpsr) + ) + assert len(grad_gpsr) > 0 diff --git a/tests/test_noise.py b/tests/test_noise.py new file mode 100644 index 00000000..67c4b1b8 --- /dev/null +++ b/tests/test_noise.py @@ -0,0 +1,262 @@ +from __future__ import annotations + +import random + +import pytest +import torch +from torch import Tensor + +from pyqtorch.apply import apply_density_mat, operator_product +from pyqtorch.circuit import QuantumCircuit +from pyqtorch.matrices import ( + HMAT, + IMAT, + XMAT, + YMAT, + ZMAT, + _dagger, +) +from pyqtorch.noise import ( + AmplitudeDamping, + GeneralizedAmplitudeDamping, + Noise, + PhaseDamping, +) +from pyqtorch.primitives import ( + ControlledPrimitive, + ControlledRotationGate, + H, + I, + Parametric, + Primitive, + X, + Y, + Z, +) +from pyqtorch.utils import ( + DensityMatrix, + density_mat, + operator_kron, + product_state, + random_state, +) + + +def test_dm(n_qubits: int, batch_size: int) -> None: + state = random_state(n_qubits) + projector = torch.outer(state.flatten(), state.conj().flatten()).view( + 2**n_qubits, 2**n_qubits, 1 + ) + dm = density_mat(state) + assert dm.size() == torch.Size([2**n_qubits, 2**n_qubits, 1]) + assert torch.allclose(dm, projector) + assert torch.allclose(dm.squeeze(), dm.squeeze() @ dm.squeeze()) + states = [] + projectors = [] + for batch in range(batch_size): + state = random_state(n_qubits) + states.append(state) + projector = torch.outer(state.flatten(), state.conj().flatten()).view( + 2**n_qubits, 2**n_qubits, 1 + ) + projectors.append(projector) + dm_proj = torch.cat(projectors, dim=2) + state_cat = torch.cat(states, dim=n_qubits) + dm = density_mat(state_cat) + assert dm.size() == torch.Size([2**n_qubits, 2**n_qubits, batch_size]) + assert torch.allclose(dm, dm_proj) + + +@pytest.mark.parametrize("n_qubits", [{"low": 2, "high": 5}], indirect=True) +def test_operator_product( + random_unitary_gate: Primitive | Parametric, + n_qubits: int, +) -> None: + batch_size_1 = torch.randint(low=1, high=5, size=(1,)).item() + batch_size_2 = torch.randint(low=1, high=5, size=(1,)).item() + max_batch = max(batch_size_2, batch_size_1) + values = {"theta": torch.rand(1)} + op = random_unitary_gate.tensor(values=values, full_support=tuple(range(n_qubits))) + op_mul = operator_product( + op.repeat(1, 1, batch_size_1), _dagger(op.repeat(1, 1, batch_size_2)) + ) + assert op_mul.size() == torch.Size([2**n_qubits, 2**n_qubits, max_batch]) + assert torch.allclose( + op_mul, + torch.eye(2**n_qubits, dtype=torch.cdouble) + .unsqueeze(2) + .repeat(1, 1, max_batch), + ) + + +@pytest.mark.parametrize("n_qubits", [{"low": 2, "high": 5}], indirect=True) +def test_apply_density_mat( + random_unitary_gate: Primitive | Parametric, + n_qubits: int, + batch_size: int, + random_input_dm: DensityMatrix, +) -> None: + values = {"theta": torch.rand(1)} + op = random_unitary_gate.tensor(values=values, full_support=tuple(range(n_qubits))) + rho = random_input_dm + rho_evol = apply_density_mat(op, rho) + assert rho_evol.size() == torch.Size([2**n_qubits, 2**n_qubits, batch_size]) + rho_expected = operator_product(op, operator_product(rho, _dagger(op))) + assert torch.allclose(rho_evol, rho_expected) + + +@pytest.mark.parametrize( + "operator,matrix", [(I, IMAT), (X, XMAT), (Z, ZMAT), (Y, YMAT), (H, HMAT)] +) +def test_operator_kron(operator: Tensor, matrix: Tensor) -> None: + n_qubits = torch.randint(low=1, high=5, size=(1,)).item() + batch_size = torch.randint(low=1, high=5, size=(1,)).item() + states, krons = [], [] + for batch in range(batch_size): + state = random_state(n_qubits) + states.append(state) + kron = torch.kron(density_mat(state).squeeze(), matrix).unsqueeze(2) + krons.append(kron) + input_state = torch.cat(states, dim=n_qubits) + kron_out = operator_kron(density_mat(input_state), operator(0).dagger()) + assert kron_out.size() == torch.Size( + [2 ** (n_qubits + 1), 2 ** (n_qubits + 1), batch_size] + ) + kron_expect = torch.cat(krons, dim=2) + assert torch.allclose(kron_out, kron_expect) + assert torch.allclose( + torch.kron(operator(0).dagger().contiguous(), I(0).tensor()), + operator_kron(operator(0).dagger(), I(0).tensor()), + ) + + +def test_kron_batch() -> None: + n_qubits = torch.randint(low=1, high=5, size=(1,)).item() + batch_size_1 = torch.randint(low=1, high=5, size=(1,)).item() + batch_size_2 = torch.randint(low=1, high=5, size=(1,)).item() + max_batch = max(batch_size_1, batch_size_2) + dm_1 = density_mat(random_state(n_qubits, batch_size_1)) + dm_2 = density_mat(random_state(n_qubits, batch_size_2)) + dm_out = operator_kron(dm_1, dm_2) + assert dm_out.size() == torch.Size( + [2 ** (2 * n_qubits), 2 ** (2 * n_qubits), max_batch] + ) + if batch_size_1 > batch_size_2: + dm_2 = dm_2.repeat(1, 1, batch_size_1)[:, :, :batch_size_1] + elif batch_size_2 > batch_size_1: + dm_1 = dm_1.repeat(1, 1, batch_size_2)[:, :, :batch_size_2] + density_matrices = [] + for batch in range(max_batch): + density_matrice = torch.kron(dm_1[:, :, batch], dm_2[:, :, batch]).unsqueeze(2) + density_matrices.append(density_matrice) + dm_expect = torch.cat(density_matrices, dim=2) + assert torch.allclose(dm_out, dm_expect) + + +@pytest.mark.parametrize("n_qubits", [{"low": 2, "high": 5}], indirect=True) +def test_flip_gates( + n_qubits: int, + target: int, + batch_size: int, + rho_input: Tensor, + random_flip_gate: Noise, + flip_expected_state: DensityMatrix, + flip_probability: Tensor | float, + flip_gates_prob_0: Noise, + flip_gates_prob_1: tuple, + random_input_dm: DensityMatrix, +) -> None: + FlipGate = random_flip_gate + output_state: DensityMatrix = FlipGate(target, flip_probability)(rho_input) + assert output_state.size() == torch.Size([2**n_qubits, 2**n_qubits, batch_size]) + assert torch.allclose(output_state, flip_expected_state) + + input_state = random_input_dm # fix the same random state for every call + assert torch.allclose(flip_gates_prob_0(input_state), input_state) + + FlipGate_1, expected_op = flip_gates_prob_1 + assert torch.allclose(FlipGate_1(input_state), expected_op) + + +def test_damping_gates( + n_qubits: int, + target: int, + batch_size: int, + damping_expected_state: tuple, + damping_gates_prob_0: Tensor, + random_input_dm: DensityMatrix, + rho_input: Tensor, +) -> None: + DampingGate, expected_state = damping_expected_state + apply_gate = DampingGate(rho_input) + assert apply_gate.size() == torch.Size([2**n_qubits, 2**n_qubits, batch_size]) + assert torch.allclose(apply_gate, expected_state) + + input_state = random_input_dm + assert torch.allclose( + damping_gates_prob_0(input_state), + I(target)(input_state), + ) + + rho_0: DensityMatrix = density_mat(product_state("0", batch_size)) + rho_1: DensityMatrix = density_mat(product_state("1", batch_size)) + if DampingGate == AmplitudeDamping: + assert torch.allclose(DampingGate(target, rate=1)(rho_1), rho_0) + elif DampingGate == PhaseDamping: + assert torch.allclose(DampingGate(target, rate=1)(rho_1), I(target)(rho_1)) + elif DampingGate == GeneralizedAmplitudeDamping: + assert torch.allclose( + DampingGate(target, error_probability=(1, 1))(rho_1), rho_0 + ) + + +@pytest.mark.parametrize("n_qubits", [{"low": 2, "high": 5}], indirect=True) +def test_noisy_primitive( + random_noisy_unitary_gate: tuple, + random_input_dm: DensityMatrix, + n_qubits: int, + batch_size: int, +) -> None: + noisy_primitive, primitve_gate, noise_gate = random_noisy_unitary_gate + state = random_input_dm + values = {"theta": torch.rand(1)} + rho_evol = noisy_primitive(state, values) + assert rho_evol.size() == torch.Size([2**n_qubits, 2**n_qubits, batch_size]) + rho_expected = noise_gate(primitve_gate(state, values)) + assert rho_expected.size() == torch.Size([2**n_qubits, 2**n_qubits, batch_size]) + assert torch.allclose(rho_evol, rho_expected) + + +@pytest.mark.parametrize("n_qubits", [{"low": 2, "high": 5}], indirect=True) +@pytest.mark.parametrize("batch_size", [{"low": 1, "high": 5}], indirect=True) +def test_noise_circ( + n_qubits: int, + batch_size: int, + random_input_dm: DensityMatrix, + random_single_qubit_gate: Primitive, + random_noise_gate: Noise, + random_rotation_gate: Parametric, + random_controlled_gate: ControlledPrimitive, + random_rotation_control_gate: ControlledRotationGate, +) -> None: + OPERATORS = [ + random_single_qubit_gate, + random_noise_gate, + random_rotation_gate, + random_controlled_gate, + random_rotation_control_gate, + ] + random.shuffle(OPERATORS) + circ = QuantumCircuit(n_qubits, OPERATORS) + + values = {random_rotation_gate.param_name: torch.rand(1)} + output_state = circ(random_input_dm, values) + assert isinstance(output_state, DensityMatrix) + assert output_state.shape == torch.Size([2**n_qubits, 2**n_qubits, batch_size]) + + diag_sums = [] + for i in range(batch_size): + diag_batch = torch.diagonal(output_state[:, :, i], dim1=0, dim2=1) + diag_sums.append(torch.sum(diag_batch)) + diag_sum = torch.stack(diag_sums) + assert torch.allclose(diag_sum, torch.ones((batch_size,), dtype=torch.cdouble)) diff --git a/tests/test_tensor.py b/tests/test_tensor.py new file mode 100644 index 00000000..503920d2 --- /dev/null +++ b/tests/test_tensor.py @@ -0,0 +1,235 @@ +from __future__ import annotations + +import random + +import pytest +import torch +from helpers import calc_mat_vec_wavefunction, get_op_support, random_pauli_hamiltonian + +from pyqtorch.composite import Add, Scale, Sequence +from pyqtorch.hamiltonians import GeneratorType, HamiltonianEvolution +from pyqtorch.primitives import ( + CNOT, + OPS_DIGITAL, + OPS_PARAM, + SWAP, + N, + Parametric, + Primitive, + Projector, +) +from pyqtorch.utils import ( + ATOL, + RTOL, + random_state, +) + +pi = torch.tensor(torch.pi) + + +@pytest.mark.parametrize("use_full_support", [True, False]) +@pytest.mark.parametrize("n_qubits", [4, 5]) +@pytest.mark.parametrize("batch_size", [1, 5]) +def test_digital_tensor(n_qubits: int, batch_size: int, use_full_support: bool) -> None: + """ + Goes through all non-parametric gates and tests their application to a random state + in comparison with the `tensor` method, either using just the qubit support of the gate + or expanding its matrix to the maximum qubit support of the full circuit. + """ + op: type[Primitive] + for op in OPS_DIGITAL: + supp = get_op_support(op, n_qubits) + op_concrete = op(*supp) + psi_init = random_state(n_qubits, batch_size) + psi_star = op_concrete(psi_init) + full_support = tuple(range(n_qubits)) if use_full_support else None + psi_expected = calc_mat_vec_wavefunction( + op_concrete, psi_init, full_support=full_support + ) + assert torch.allclose(psi_star, psi_expected, rtol=RTOL, atol=ATOL) + + +@pytest.mark.parametrize("use_full_support", [True, False]) +@pytest.mark.parametrize("n_qubits", [4, 5]) +@pytest.mark.parametrize("batch_size", [1, 5]) +def test_param_tensor(n_qubits: int, batch_size: int, use_full_support: bool) -> None: + """ + Goes through all parametric gates and tests their application to a random state + in comparison with the `tensor` method, either using just the qubit support of the gate + or expanding its matrix to the maximum qubit support of the full circuit. + """ + op: type[Parametric] + for op in OPS_PARAM: + supp = get_op_support(op, n_qubits) + params = [f"th{i}" for i in range(op.n_params)] + op_concrete = op(*supp, *params) + psi_init = random_state(n_qubits) + values = {param: torch.rand(batch_size) for param in params} + psi_star = op_concrete(psi_init, values) + full_support = tuple(range(n_qubits)) if use_full_support else None + psi_expected = calc_mat_vec_wavefunction( + op_concrete, psi_init, values=values, full_support=full_support + ) + assert torch.allclose(psi_star, psi_expected, rtol=RTOL, atol=ATOL) + + +@pytest.mark.parametrize("use_full_support", [True, False]) +@pytest.mark.parametrize("n_qubits", [4, 5]) +@pytest.mark.parametrize("batch_size", [1, 5]) +@pytest.mark.parametrize("compose", [Sequence, Add]) +def test_sequence_tensor( + n_qubits: int, + batch_size: int, + use_full_support: bool, + compose: type[Sequence] | type[Add], +) -> None: + op_list = [] + values = {} + op: type[Primitive] | type[Parametric] + """ + Builds a Sequence or Add composition of all possible gates on random qubit + supports. Also assigns a Scale of a random value to the non-parametric gates. + Tests the forward method (which goes through each gate individually) to the + `tensor` method, which builds the full operator matrix and applies it. + """ + for op in OPS_DIGITAL: + supp = get_op_support(op, n_qubits) + op_concrete = Scale(op(*supp), torch.rand(1)) + op_list.append(op_concrete) + for op in OPS_PARAM: + supp = get_op_support(op, n_qubits) + params = [f"{op.__name__}_th{i}" for i in range(op.n_params)] + values.update({param: torch.rand(batch_size) for param in params}) + op_concrete = op(*supp, *params) + op_list.append(op_concrete) + random.shuffle(op_list) + op_composite = compose(op_list) + psi_init = random_state(n_qubits, batch_size) + psi_star = op_composite(psi_init, values) + full_support = tuple(range(n_qubits)) if use_full_support else None + psi_expected = calc_mat_vec_wavefunction( + op_composite, psi_init, values=values, full_support=full_support + ) + assert torch.allclose(psi_star, psi_expected, rtol=RTOL, atol=ATOL) + + +@pytest.mark.parametrize("use_full_support", [True, False]) +@pytest.mark.parametrize("n_qubits", [4, 5]) +@pytest.mark.parametrize("n_proj", [1, 3]) +@pytest.mark.parametrize("batch_size", [1, 5]) +def test_projector_tensor( + n_qubits: int, n_proj: int, batch_size: int, use_full_support: bool +) -> None: + """ + Instantiates various random projectors on arbitrary qubit support + and compares the forward method with directly applying the tensor. + """ + iterations = 5 + for _ in range(iterations): + rand_int_1 = random.randint(0, 2**n_proj - 1) + rand_int_2 = random.randint(0, 2**n_proj - 1) + bitstring1 = "{0:b}".format(rand_int_1).zfill(n_proj) + bitstring2 = "{0:b}".format(rand_int_2).zfill(n_proj) + supp = tuple(random.sample(range(n_qubits), n_proj)) + op = Projector(supp, bitstring1, bitstring2) + psi_init = random_state(n_qubits, batch_size) + psi_star = op(psi_init) + full_support = tuple(range(n_qubits)) if use_full_support else None + psi_expected = calc_mat_vec_wavefunction( + op, psi_init, full_support=full_support + ) + assert torch.allclose(psi_star, psi_expected, rtol=RTOL, atol=ATOL) + + +@pytest.mark.parametrize("n_qubits", [4, 5]) +@pytest.mark.parametrize("operator", [N, CNOT, SWAP]) +@pytest.mark.parametrize("use_full_support", [True, False]) +def test_projector_vs_operator( + n_qubits: int, + operator: type[Primitive], + use_full_support: bool, +) -> None: + if operator == N: + supp: tuple = (random.randint(0, n_qubits - 1),) + op_concrete = N(*supp) + projector = Projector(supp, "1", "1") + if operator == CNOT: + supp = tuple(random.sample(range(n_qubits), 2)) + op_concrete = CNOT(*supp) + projector = Add( + [ + Projector(supp, "00", "00"), + Projector(supp, "01", "01"), + Projector(supp, "10", "11"), + Projector(supp, "11", "10"), + ] + ) + if operator == SWAP: + supp = tuple(random.sample(range(n_qubits), 2)) + op_concrete = SWAP(*supp) + projector = Add( + [ + Projector(supp, "00", "00"), + Projector(supp, "01", "10"), + Projector(supp, "10", "01"), + Projector(supp, "11", "11"), + ] + ) + full_support = tuple(range(n_qubits)) if use_full_support else None + projector_mat = projector.tensor(full_support=full_support) + operator_mat = op_concrete.tensor(full_support=full_support) + assert torch.allclose(projector_mat, operator_mat, rtol=RTOL, atol=ATOL) + + +@pytest.mark.parametrize("n_qubits", [4, 5]) +@pytest.mark.parametrize("make_param", [True, False]) +@pytest.mark.parametrize("use_full_support", [True, False]) +@pytest.mark.parametrize("batch_size", [1, 5]) +def test_hevo_pauli_tensor( + n_qubits: int, make_param: bool, use_full_support: bool, batch_size: int +) -> None: + k_1q = 2 * n_qubits # Number of 1-qubit terms + k_2q = n_qubits**2 # Number of 2-qubit terms + generator, param_list = random_pauli_hamiltonian(n_qubits, k_1q, k_2q, make_param) + values = {param: torch.rand(batch_size) for param in param_list} + psi_init = random_state(n_qubits, batch_size) + full_support = tuple(range(n_qubits)) if use_full_support else None + # Test the generator itself + psi_star = generator(psi_init, values) + psi_expected = calc_mat_vec_wavefunction(generator, psi_init, values, full_support) + assert torch.allclose(psi_star, psi_expected, rtol=RTOL, atol=ATOL) + # Test the hamiltonian evolution + tparam = "t" + operator = HamiltonianEvolution(generator, tparam, generator_parametric=make_param) + if make_param: + assert operator.generator_type == GeneratorType.PARAMETRIC_OPERATION + else: + assert operator.generator_type == GeneratorType.OPERATION + values[tparam] = torch.rand(batch_size) + psi_star = operator(psi_init, values) + psi_expected = calc_mat_vec_wavefunction(operator, psi_init, values, full_support) + assert torch.allclose(psi_star, psi_expected, rtol=RTOL, atol=ATOL) + + +@pytest.mark.parametrize("gen_qubits", [3, 4]) +@pytest.mark.parametrize("n_qubits", [4, 5]) +@pytest.mark.parametrize("use_full_support", [True, False]) +@pytest.mark.parametrize("batch_size", [1, 5]) +def test_hevo_tensor_tensor( + gen_qubits: int, n_qubits: int, use_full_support: bool, batch_size: int +) -> None: + k_1q = 2 * gen_qubits # Number of 1-qubit terms + k_2q = gen_qubits**2 # Number of 2-qubit terms + generator, _ = random_pauli_hamiltonian(gen_qubits, k_1q, k_2q) + psi_init = random_state(n_qubits, batch_size) + full_support = tuple(range(n_qubits)) if use_full_support else None + # Test the hamiltonian evolution + generator_matrix = generator.tensor() + supp = tuple(random.sample(range(n_qubits), gen_qubits)) + tparam = "t" + operator = HamiltonianEvolution(generator_matrix, tparam, supp) + assert operator.generator_type == GeneratorType.TENSOR + values = {tparam: torch.rand(batch_size)} + psi_star = operator(psi_init, values) + psi_expected = calc_mat_vec_wavefunction(operator, psi_init, values, full_support) + assert torch.allclose(psi_star, psi_expected, rtol=RTOL, atol=ATOL)