From db1c8d19a6022ad70f2092bd9ba6df530a345024 Mon Sep 17 00:00:00 2001 From: mateuszkupper Date: Thu, 9 Oct 2025 13:44:09 +0100 Subject: [PATCH 01/19] Tidying up backends.py --- optyx/core/backends.py | 94 ++++++++++++++++++++++++++---------------- 1 file changed, 58 insertions(+), 36 deletions(-) diff --git a/optyx/core/backends.py b/optyx/core/backends.py index 552933f..d70b7dc 100644 --- a/optyx/core/backends.py +++ b/optyx/core/backends.py @@ -66,7 +66,7 @@ from abc import ABC, abstractmethod from dataclasses import dataclass -from typing import Any, Union +from typing import Any from collections import defaultdict from enum import Enum from cotengra import ( @@ -79,8 +79,7 @@ import numpy as np import perceval as pcvl from quimb.tensor import TensorNetwork -from optyx.core.channel import Diagram -from optyx.core.channel import Ty, mode, bit +from optyx.core.channel import Diagram, Ty, mode, bit from optyx.utils.utils import preprocess_quimb_tensors_safe @@ -113,16 +112,16 @@ def tensor(self) -> discopy_tensor.Box: return self._tensor @property - def density_matrix(self) -> discopy_tensor.Box: + def density_matrix(self) -> np.ndarray: """ Get the density matrix from the result tensor. Returns: - tensor.Box: The density matrix. + np.ndarray: The density matrix. """ if len(self.tensor.dom) != 0: raise ValueError( - "Result tensor must represent a state with inputs." + "Result tensor must represent a state with no inputs." ) if self.state_type not in {StateType.AMP, StateType.DM}: raise TypeError( @@ -130,10 +129,10 @@ def density_matrix(self) -> discopy_tensor.Box: ) if self.state_type is StateType.AMP: density_matrix = self.tensor.dagger() >> self.tensor - return density_matrix + return density_matrix.array return self.tensor.array - def amplitudes(self, normalise=True) -> dict[tuple[int, ...], float]: + def amplitudes(self, normalise=True) -> dict[tuple[int, ...], complex]: """ Get the amplitudes from the result tensor. Returns: @@ -147,7 +146,7 @@ def amplitudes(self, normalise=True) -> dict[tuple[int, ...], float]: ) if len(self.tensor.dom) != 0: raise ValueError( - "Result tensor must represent a state with inputs." + "Result tensor must represent a state with no inputs." ) dic = self._convert_array_to_dict(self.tensor.array) @@ -167,7 +166,7 @@ def prob_dist(self, round_digits: int = None) -> dict: """ if len(self.tensor.dom) != 0: raise ValueError( - "Result tensor must represent a state with inputs." + "Result tensor must represent a state with no inputs." ) if self.state_type is StateType.AMP: return self._prob_dist_pure(round_digits) @@ -231,6 +230,8 @@ def _prob_dist_pure(self, round_digits: int = None) -> dict: round_digits=round_digits ) sum_ = np.sum(np.abs(list(values.values())) ** 2) + if sum_ == 0: + raise ValueError("The probability distribution sums to zero.") return {key: (abs(value) ** 2)/sum_ for key, value in values.items()} def _prob_dist_mixed( @@ -241,6 +242,10 @@ def _prob_dist_mixed( This method computes the probability distribution by aggregating occupation configurations based on the output types of the tensor. + Assumes the output types contain at least one 'bit' or 'mode' + These will be treated as measured registers while 'qubit' and 'qmode' + are treated as unmeasured and traced out. + Args: round_digits: Optional number of digits to round the probabilities. @@ -251,10 +256,11 @@ def _prob_dist_mixed( if not any(t in {bit, mode} for t in self.output_types): raise ValueError( - "Output types must contain at least one 'bit' or 'mode'." + "Output types must contain at least one 'bit' or 'mode'." + + "These will be treated as measured registers." ) - values = self._convert_array_to_dict(self.tensor.array, round_digits) + values = self._convert_array_to_dict(self.tensor.array, None) mask_flat = np.concatenate( [[1] if t in {bit, mode} else [0, 0] for t in self.output_types] ) @@ -270,15 +276,25 @@ def _prob_dist_mixed( zip(key, mask_flat) if not m) if all(occs_unmeasured[i] == occs_unmeasured[i + 1] for i in range(0, len(occs_unmeasured) - 1, 2)): - probs[occ_measured] += amp + + val = float(np.real_if_close(amp)) + if val < 0 and abs(val) < 1e-12: + val = 0.0 + probs[occ_measured] += val for occ in all_measured: probs.setdefault(occ, 0.0) sum_ = np.sum(list(probs.values())) + if sum_ == 0: + raise ValueError("The probability distribution sums to zero.") prob = { key: value / sum_ for key, value in probs.items() } + + if round_digits is not None: + prob = {k: round(v, round_digits) for k, v in prob.items()} + return prob @@ -361,12 +377,13 @@ class QuimbBackend(AbstractBackend): def __init__( self, - hyperoptimiser: Union[ - HyperOptimizer, - ReusableHyperOptimizer, - HyperCompressedOptimizer, - ReusableHyperCompressedOptimizer - ] = None, + hyperoptimiser: ( + HyperOptimizer | + ReusableHyperOptimizer | + HyperCompressedOptimizer | + ReusableHyperCompressedOptimizer | + None + ) = None, contraction_params: dict = None): """ Initialize the Quimb backend. @@ -395,7 +412,7 @@ def eval( tensor_diagram = self._get_discopy_tensor(diagram) - if hasattr(diagram, 'terms'): + if hasattr(tensor_diagram, 'terms'): results = sum( self._process_term(term) for term in tensor_diagram.terms ) @@ -418,12 +435,12 @@ def eval( state_type=state_type ) - def _process_term(self, term: Diagram) -> np.ndarray: + def _process_term(self, term: discopy_tensor.Diagram) -> np.ndarray: """ Process a term in a sum of diagrams. Args: - term (Diagram): The term to process. + term (discopy.tensor.Diagram): The term to process. Returns: np.ndarray: The processed term as a numpy array. @@ -540,29 +557,34 @@ def eval( The result of the evaluation. """ - if extra: + perceval_state = extra.get( + "perceval_state", + pcvl.StateVector([1] * len(diagram.dom)) + ) + + if not isinstance(perceval_state, pcvl.StateVector): try: - perceval_state: pcvl.StateVector = extra["perceval_state"] - except KeyError as error: + perceval_state = pcvl.StateVector(list(perceval_state)) + except Exception as e: raise TypeError( - "PercevalBackend.eval requires " + - "a 'perceval_state=' keyword." - ) from error - else: - perceval_state = pcvl.StateVector( - [1] * len(diagram.dom) - ) + "perceval_state must be a perceval.StateVector" + " or a sequence of non-negative " + + "integers (occupation numbers)." + ) from e + tensor_diagram = self._get_discopy_tensor(diagram) sim = pcvl.Simulator(self.perceval_backend) matrix = self._get_matrix(diagram) - if not np.allclose( - np.eye(matrix.shape[0]), - matrix.dot(matrix.conj().T) + if ( + matrix.shape[0] != matrix.shape[1] or + not np.allclose(np.eye(matrix.shape[0]), + matrix.dot(matrix.conj().T)) ): raise ValueError( - "The provided diagram does not represent a unitary operation." + "The provided diagram does not represent a" + + " unitary operation or the matrix is not square." ) perceval_circuit = self._umatrix_to_perceval_circuit(matrix) From e0dfb4a7a96d42273a68b00280ec254499455fa3 Mon Sep 17 00:00:00 2001 From: mateuszkupper Date: Thu, 9 Oct 2025 14:40:31 +0100 Subject: [PATCH 02/19] Add extra tests to for backends --- test/test_backends.py | 171 +++++++++++++++++++++++++++++++++++++++++- 1 file changed, 168 insertions(+), 3 deletions(-) diff --git a/test/test_backends.py b/test/test_backends.py index 02ef527..32ec44c 100644 --- a/test/test_backends.py +++ b/test/test_backends.py @@ -1,11 +1,12 @@ import pytest -from optyx import photonic, qubits, classical +from optyx import photonic, qubits, classical, mode, qubit, qmode, bit from cotengra import ReusableHyperCompressedOptimizer -from optyx.core.backends import QuimbBackend, PercevalBackend, DiscopyBackend +from optyx.core.backends import QuimbBackend, PercevalBackend, DiscopyBackend, EvalResult, StateType import numpy as np import math from itertools import chain import perceval as pcvl +import discopy.tensor as discopy_tensor @pytest.mark.skip(reason="Helper function for testing") def chip_mzi(w, l): @@ -221,6 +222,159 @@ def test_probs_sum_to_one(self, circuit): total_prob = sum(prob_dist.values()) assert math.isclose(total_prob, 1.0) + def test_prob_lookup_and_missing_key(self): + v = np.array([1/np.sqrt(2), 1j/np.sqrt(2)], dtype=complex) + box = discopy_tensor.Tensor(v, discopy_tensor.Dim(1), discopy_tensor.Dim(2)) + ev = EvalResult(_tensor=box, output_types=(mode,), state_type=StateType.AMP) + + p = ev.prob_dist() + assert ev.prob((0,)) == pytest.approx(p[(0,)], rel=1e-12) + assert ev.prob((1,)) == pytest.approx(p[(1,)], rel=1e-12) + assert ev.prob((2,)) == 0.0 + + def test_density_matrix_from_amp_is_outer_product(self): + v = np.array([1/np.sqrt(3), np.sqrt(2/3)], dtype=complex) + box = discopy_tensor.Tensor(v, discopy_tensor.Dim(1), discopy_tensor.Dim(2)) + ev = EvalResult(_tensor=box, output_types=(mode,), state_type=StateType.AMP) + dm = ev.density_matrix + assert dm.shape == (2, 2) + assert np.allclose(dm, np.outer(v.conj(), v), atol=1e-12) + + def test_mixed_interleaved_partial_trace_numeric_hygiene(self): + dm = np.zeros((2, 2, 2), dtype=complex) + dm[0, 0, 0] = 0.499999999999 + 1e-14j + dm[0, 1, 1] = 0.100000000001 - 1e-14j + dm[1, 0, 0] = 0.2000000000005 + 1e-14j + dm[1, 1, 1] = 0.2000000000005 - 1e-14j + #off-diagonals that should be ignored by the diagonal check + dm[0, 0, 1] = -1e-13 + dm[1, 1, 0] = -1e-13 + + box = discopy_tensor.Tensor(dm, discopy_tensor.Dim(1), discopy_tensor.Dim(*dm.shape)) + #1st axis measured (mode), 2ns is quantum (unmeasured, bra/ket) + ev = EvalResult(_tensor=box, output_types=(mode, object()), state_type=StateType.DM) + probs = ev.prob_dist() + + assert pytest.approx(probs[(0,)], rel=1e-12) == 0.6 + assert pytest.approx(probs[(1,)], rel=1e-12) == 0.4 + assert pytest.approx(sum(probs.values()), rel=1e-12) == 1.0 + + def test_prob_branch_rounding_and_sum(self): + P = np.array([[1/3, 1/6], + [1/6, 1/3]], dtype=float) + box = discopy_tensor.Tensor(P, discopy_tensor.Dim(1), discopy_tensor.Dim(*P.shape)) + ev = EvalResult(_tensor=box, output_types=(mode, mode), state_type=StateType.PROB) + + d = ev.prob_dist(round_digits=4) + assert pytest.approx(sum(d.values()), rel=1e-12) == 1.0 + print(d) + assert d[(0, 0)] == round(float(P[0, 0]), 4) + assert d[(1, 1)] == round(float(P[1, 1]), 4) + + def test_density_matrix_from_amp_complex_phase(self): + v = np.array([np.exp(1j*0.2)/np.sqrt(3), + np.exp(-1j*0.3)*np.sqrt(2/3)], dtype=complex) + box = discopy_tensor.Tensor(v, discopy_tensor.Dim(1), discopy_tensor.Dim(2)) + ev = EvalResult(_tensor=box, output_types=(mode,), state_type=StateType.AMP) + dm = ev.density_matrix + assert dm.shape == (2, 2) + assert np.allclose(dm, dm.conj().T, atol=1e-12) + assert np.allclose(dm, np.outer(v.conj(), v), atol=1e-12) + eig = np.linalg.eigvalsh(dm) + assert np.all(eig >= -1e-12) + + def test_density_matrix_passthrough_for_DM(self): + dm = np.array([[0.7, 0.1j], + [-0.1j, 0.3]], dtype=complex) + box = discopy_tensor.Tensor(dm, discopy_tensor.Dim(1), discopy_tensor.Dim(*dm.shape)) + ev = EvalResult(_tensor=box, output_types=(mode,), state_type=StateType.DM) + out = ev.density_matrix + assert out.shape == (2, 2) + assert np.allclose(out, dm, atol=1e-12) + + def test_mixed_simple_single_measured_single_unmeasured(self): + dm = np.zeros((2, 2, 2), dtype=complex) + dm[0, 0, 0] = 0.6 + 1e-14j + dm[1, 1, 1] = 0.4 - 1e-14j + dm[0, 1, 0] = -1e-13 + dm[1, 0, 1] = -1e-13 + + box = discopy_tensor.Tensor(dm, discopy_tensor.Dim(1), discopy_tensor.Dim(*dm.shape)) + ev = EvalResult(_tensor=box, output_types=(mode, qubit), state_type=StateType.DM) + probs = ev.prob_dist() + assert pytest.approx(probs[(0,)], rel=1e-12) == 0.6 + assert pytest.approx(probs[(1,)], rel=1e-12) == 0.4 + assert pytest.approx(sum(probs.values()), rel=1e-12) == 1.0 + + def test_mixed_two_measured_one_unmeasured_matches_manual_trace(self): + P = np.array([[0.1, 0.2], + [0.3, 0.4]], dtype=float) + W = np.array([0.6, 0.4], dtype=float) + dm = np.zeros((2, 2, 2, 2), dtype=complex) + + for k in range(2): + for l in range(2): + dm[k, l, 0, 0] = P[k, l] * W[0] + dm[k, l, 1, 1] = P[k, l] * W[1] + dm[0, 0, 0, 1] = 1e-13 + dm[1, 1, 1, 0] = -1e-13 + + box = discopy_tensor.Tensor(dm, discopy_tensor.Dim(1), discopy_tensor.Dim(*dm.shape)) + ev = EvalResult(_tensor=box, output_types=(mode, bit, qubit), state_type=StateType.DM) + + probs = ev.prob_dist() + manual = {(k, l): (dm[k, l, 0, 0].real + dm[k, l, 1, 1].real) for k in range(2) for l in range(2)} + Z = sum(manual.values()) + manual = {k: v / Z for k, v in manual.items()} + + for k in range(2): + for l in range(2): + assert pytest.approx(probs[(k, l)], rel=1e-12) == manual[(k, l)] + assert pytest.approx(sum(probs.values()), rel=1e-12) == 1.0 + + def test_mixed_round_digits_applied_before_normalization(self): + dm = np.zeros((2, 2, 2), dtype=complex) + dm[0, 0, 0] = 0.333333333333 + dm[1, 1, 1] = 0.666666666667 + box = discopy_tensor.Tensor(dm, discopy_tensor.Dim(1), discopy_tensor.Dim(*dm.shape)) + ev = EvalResult(_tensor=box, output_types=(mode, qubit), state_type=StateType.DM) + d = ev.prob_dist(round_digits=6) + assert pytest.approx(sum(d.values()), rel=1e-12) == 1.0 + assert d[(0,)] == pytest.approx(round(0.333333333333, 6) / (round(0.333333333333, 6) + round(0.666666666667, 6)), rel=1e-12) + assert d[(1,)] == pytest.approx(round(0.666666666667, 6) / (round(0.333333333333, 6) + round(0.666666666667, 6)), rel=1e-12) + + def test_mixed_requires_at_least_one_measured_type(self): + dm = np.eye(4, dtype=complex).reshape(2, 2, 2, 2) / 4.0 + box = discopy_tensor.Tensor(dm, discopy_tensor.Dim(1), discopy_tensor.Dim(*dm.shape)) + ev = EvalResult(_tensor=box, output_types=(qubit, qmode), state_type=StateType.DM) + with pytest.raises(ValueError): + _ = ev.prob_dist() + + def test_mixed_two_unmeasured_pairs_partial_trace(self): + + dm = np.zeros((2, 2, 2, 2, 2), dtype=float) + + w_m = {0: 0.55, 1: 0.45} + w_u1 = np.array([0.6, 0.4]) + w_u2 = np.array([0.7, 0.3]) + + for m in (0, 1): + for r1 in (0, 1): + for r2 in (0, 1): + dm[m, r1, r1, r2, r2] = w_m[m] * w_u1[r1] * w_u2[r2] + + dm[0, 0, 1, 0, 0] = 1e-13 + dm[1, 1, 0, 1, 0] = -1e-13 + + box = discopy_tensor.Tensor(dm, discopy_tensor.Dim(1), discopy_tensor.Dim(*dm.shape)) + ev = EvalResult(_tensor=box, output_types=(mode, object(), object()), state_type=StateType.DM) + probs = ev.prob_dist() + + assert pytest.approx(sum(probs.values()), rel=1e-12) == 1.0 + assert pytest.approx(probs[(0,)], rel=1e-12) == 0.55 + assert pytest.approx(probs[(1,)], rel=1e-12) == 0.45 + + class TestExceptions: # EvalResult # density matrix/amps/prob from not a state @@ -265,4 +419,15 @@ def test_abstract_backend_not_lo(self, circuit): with pytest.raises(AssertionError): backend = PercevalBackend() perceval_state = pcvl.BasicState(get_state(circuit)) - result_perceval = circuit.eval(backend, perceval_state=perceval_state) \ No newline at end of file + result_perceval = circuit.eval(backend, perceval_state=perceval_state) + + def test_mixed_all_off_diagonal_mass_raises_zero_total(self): + dm = np.zeros((2, 2, 2), dtype=float) + dm[0, 0, 1] = 0.6 + dm[1, 1, 0] = 0.4 + + box = discopy_tensor.Tensor(dm, discopy_tensor.Dim(1), discopy_tensor.Dim(*dm.shape)) + ev = EvalResult(_tensor=box, output_types=(mode, object()), state_type=StateType.DM) + + with pytest.raises(ValueError): + _ = ev.prob_dist() \ No newline at end of file From 9e5181e7463f97142824462f4dde5afeaa80996e Mon Sep 17 00:00:00 2001 From: mateuszkupper Date: Thu, 9 Oct 2025 16:57:29 +0100 Subject: [PATCH 03/19] Path permanent eval as a backend --- optyx/core/backends.py | 57 ++++++++++++++++++++++++++++++++++++++ test/test_backends.py | 62 ++++++++++++++++++++++++++++++++++++++++-- 2 files changed, 117 insertions(+), 2 deletions(-) diff --git a/optyx/core/backends.py b/optyx/core/backends.py index d70b7dc..c289057 100644 --- a/optyx/core/backends.py +++ b/optyx/core/backends.py @@ -619,3 +619,60 @@ def eval( output_types=diagram.cod, state_type=StateType.PROB ) + + +class PermanentBackend(AbstractBackend): + """ + Backend implementation using optyx' Path module to compute matrix + permanents. + """ + + def eval( + self, + diagram: Diagram, + **extra: Any + ): + """ + Evaluate the diagram using the Permanent/Path backend. + Works only for LO circuits. + + Args: + diagram (Diagram): The diagram to evaluate. + **extra: Additional arguments for the evaluation, + including 'n_photons' for optyx.core.path.Matrix.eval. + """ + try: + matrix = diagram.to_path() + except AssertionError as error: + raise ValueError( + "The diagram cannot be converted to a matrix. " + + "It is not linear optical." + ) from error + + n_photons = extra.get( + "n_photons", + 0 + ) + + if ( + len(matrix.creations) == 0 and + n_photons == 0 + ): + raise ValueError( + "The diagram does not include any photon creations. " + + "n_photons must be greater than 0." + ) + + tensor_diagram = self._get_discopy_tensor(diagram) + result = matrix.eval(n_photons=n_photons, as_tensor=True) + print(result) + return EvalResult( + discopy_tensor.Box( + "Result", + tensor_diagram.dom, + tensor_diagram.cod, + result.array + ), + output_types=diagram.cod, + state_type=StateType.AMP + ) diff --git a/test/test_backends.py b/test/test_backends.py index 32ec44c..038d99b 100644 --- a/test/test_backends.py +++ b/test/test_backends.py @@ -1,7 +1,14 @@ import pytest from optyx import photonic, qubits, classical, mode, qubit, qmode, bit from cotengra import ReusableHyperCompressedOptimizer -from optyx.core.backends import QuimbBackend, PercevalBackend, DiscopyBackend, EvalResult, StateType +from optyx.core.backends import ( + QuimbBackend, + PercevalBackend, + DiscopyBackend, + EvalResult, + StateType, + PermanentBackend +) import numpy as np import math from itertools import chain @@ -430,4 +437,55 @@ def test_mixed_all_off_diagonal_mass_raises_zero_total(self): ev = EvalResult(_tensor=box, output_types=(mode, object()), state_type=StateType.DM) with pytest.raises(ValueError): - _ = ev.prob_dist() \ No newline at end of file + _ = ev.prob_dist() + + @pytest.mark.parametrize( + "circuit", + MIXED_CIRCUITS_TO_TEST + CIRCUITS_WITH_DISCARDS_TO_TEST, + ) + def test_non_lo_circuits_raise_notimplemented(self, circuit): + backend = PermanentBackend() + with pytest.raises(ValueError): + _ = circuit.eval(backend) + + @pytest.mark.parametrize("circuit", PURE_CIRCUITS_TO_TEST) + def test_zero_photons_no_creations_raises_valueerror(self, circuit): + backend = PermanentBackend() + with pytest.raises(ValueError): + _ = circuit.eval(backend, n_photons=0) + +class TestPermanentBackendVsQuimb: + @pytest.mark.parametrize("circuit", PURE_CIRCUITS_TO_TEST) + def test_permanent_amp_matches_quimb_with_create(self, circuit): + """ + Diagram has explicit Create(...), so PermanentBackend.eval returns a **state** (AMP). + Compare amplitudes, probs, and raw arrays with Quimb. + """ + state_occ = get_state(circuit) + state = photonic.Create(*state_occ) + diagram = state >> circuit + + result_quimb = diagram.eval() + + backend = PermanentBackend() + result_perm = diagram.eval(backend) + + assert dict_allclose(result_perm.amplitudes(), result_quimb.amplitudes()) + assert dict_allclose(result_perm.prob_dist(), result_quimb.prob_dist()) + assert np.allclose(result_perm.tensor.array, result_quimb.tensor.array, atol=1e-12) + + @pytest.mark.parametrize("circuit", PURE_CIRCUITS_TO_TEST) + def test_permanent_prob_matches_quimb_with_create(self, circuit): + """ + Same as above but request PROB directly from PermanentBackend. + """ + state_occ = get_state(circuit) + state = photonic.Create(*state_occ) + diagram = state >> circuit + + result_quimb = diagram.eval() + + backend = PermanentBackend() + result_perm = diagram.eval(backend, return_kind="prob") + + assert dict_allclose(result_perm.prob_dist(), result_quimb.prob_dist()) \ No newline at end of file From 0f8b5850095e298978457c9d5def7e5fcd330848 Mon Sep 17 00:00:00 2001 From: mateuszkupper Date: Thu, 9 Oct 2025 17:24:53 +0100 Subject: [PATCH 04/19] Added support for sums for the permanent backend --- docs/notebooks/bosonic-vqe-2.ipynb | 809 ++++++++++++++++------------- optyx/core/backends.py | 52 +- test/test_backends.py | 2 +- 3 files changed, 470 insertions(+), 393 deletions(-) diff --git a/docs/notebooks/bosonic-vqe-2.ipynb b/docs/notebooks/bosonic-vqe-2.ipynb index d028b85..dc860e5 100644 --- a/docs/notebooks/bosonic-vqe-2.ipynb +++ b/docs/notebooks/bosonic-vqe-2.ipynb @@ -17,7 +17,7 @@ " \n", " \n", " \n", - " 2025-09-05T20:09:47.097250\n", + " 2025-10-09T17:20:26.561086\n", " image/svg+xml\n", " \n", " \n", @@ -46,62 +46,62 @@ "L 223.2 7.2 \n", "L 7.2 7.2 \n", "z\n", - "\" clip-path=\"url(#pd1650474ac)\" style=\"fill: #ffffff; stroke: #ffffff; stroke-linejoin: miter\"/>\n", + "\" clip-path=\"url(#p605fbb3927)\" style=\"fill: #ffffff; stroke: #ffffff; stroke-linejoin: miter\"/>\n", " \n", " \n", " \n", + "\" clip-path=\"url(#p605fbb3927)\" style=\"fill: none; stroke: #000000; stroke-linejoin: miter\"/>\n", " \n", " \n", " \n", + "\" clip-path=\"url(#p605fbb3927)\" style=\"fill: none; stroke: #000000; stroke-linejoin: miter\"/>\n", " \n", " \n", " \n", + "\" clip-path=\"url(#p605fbb3927)\" style=\"fill: none; stroke: #000000; stroke-linejoin: miter\"/>\n", " \n", " \n", " \n", + "\" clip-path=\"url(#p605fbb3927)\" style=\"fill: none; stroke: #000000; stroke-linejoin: miter\"/>\n", " \n", " \n", " \n", + "\" clip-path=\"url(#p605fbb3927)\" style=\"fill: none; stroke: #000000; stroke-linejoin: miter\"/>\n", " \n", " \n", " \n", + "\" clip-path=\"url(#p605fbb3927)\" style=\"fill: none; stroke: #000000; stroke-linejoin: miter\"/>\n", " \n", " \n", " \n", + "\" clip-path=\"url(#p605fbb3927)\" style=\"fill: none; stroke: #000000; stroke-linejoin: miter\"/>\n", " \n", " \n", " \n", + "\" clip-path=\"url(#p605fbb3927)\" style=\"fill: none; stroke: #000000; stroke-linejoin: miter\"/>\n", " \n", " \n", " \n", + "\" clip-path=\"url(#p605fbb3927)\" style=\"fill: none; stroke: #000000; stroke-linejoin: miter\"/>\n", " \n", " \n", " \n", + "\" clip-path=\"url(#p605fbb3927)\" style=\"fill: none; stroke: #000000; stroke-linejoin: miter\"/>\n", " \n", " \n", " \n", + "\" clip-path=\"url(#p605fbb3927)\" style=\"fill: none; stroke: #000000; stroke-linejoin: miter\"/>\n", " \n", " \n", " \n", + "\" clip-path=\"url(#p605fbb3927)\" style=\"fill: #ffffff; stroke: #000000; stroke-linejoin: miter\"/>\n", " \n", " \n", " \n", + "\" clip-path=\"url(#p605fbb3927)\" style=\"fill: #ffffff; stroke: #000000; stroke-linejoin: miter\"/>\n", " \n", " \n", " \n", + "\" clip-path=\"url(#p605fbb3927)\" style=\"fill: #ffffff; stroke: #000000; stroke-linejoin: miter\"/>\n", " \n", " \n", " \n", + "\" clip-path=\"url(#p605fbb3927)\" style=\"fill: #ffffff; stroke: #000000; stroke-linejoin: miter\"/>\n", " \n", " \n", " \n", + "\" clip-path=\"url(#p605fbb3927)\" style=\"fill: #ffffff; stroke: #000000; stroke-linejoin: miter\"/>\n", " \n", " \n", " \n", + "\" clip-path=\"url(#p605fbb3927)\" style=\"fill: #ffffff; stroke: #000000; stroke-linejoin: miter\"/>\n", " \n", " \n", " \n", @@ -878,7 +878,7 @@ " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", @@ -911,12 +911,12 @@ "\n", "\n", - "\n", + "\n", " \n", " \n", " \n", " \n", - " 2025-09-05T20:09:47.853222\n", + " 2025-10-09T17:20:27.356367\n", " image/svg+xml\n", " \n", " \n", @@ -931,161 +931,209 @@ " \n", " \n", " \n", - " \n", " \n", " \n", " \n", - " \n", + "\" clip-path=\"url(#pe741f2f76b)\" style=\"fill: #ffffff; stroke: #ffffff; stroke-linejoin: miter\"/>\n", " \n", " \n", - " \n", + " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", + "\" clip-path=\"url(#pe741f2f76b)\" style=\"fill: #ffffff; stroke: #000000; stroke-linejoin: miter\"/>\n", " \n", " \n", - " \n", + "\" clip-path=\"url(#pe741f2f76b)\" style=\"fill: #ffffff; stroke: #000000; stroke-linejoin: miter\"/>\n", " \n", " \n", - " \n", + "\" clip-path=\"url(#pe741f2f76b)\" style=\"fill: #ffffff; stroke: #000000; stroke-linejoin: miter\"/>\n", " \n", " \n", - " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "\" clip-path=\"url(#pe741f2f76b)\" style=\"fill: #ffffff; stroke: #000000; stroke-linejoin: miter\"/>\n", + " \n", + " \n", + " \n", " \n", " \n", " \n", - " \n", " \n", - " \n", - " \n", + " \n", + " \n", " \n", - " \n", - " \n", + " \n", + " \n", " \n", - " \n", - " \n", + " \n", + " \n", " \n", " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", @@ -1235,7 +1283,7 @@ " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", @@ -1245,7 +1293,7 @@ " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", @@ -1255,7 +1303,7 @@ " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", @@ -1265,7 +1313,7 @@ " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", @@ -1275,7 +1323,7 @@ " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", @@ -1285,7 +1333,7 @@ " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", @@ -1295,7 +1343,7 @@ " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", @@ -1305,7 +1353,7 @@ " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", @@ -1315,7 +1363,7 @@ " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", @@ -1325,7 +1373,7 @@ " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", @@ -1335,7 +1383,7 @@ " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", @@ -1345,7 +1393,7 @@ " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", @@ -1355,7 +1403,7 @@ " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", @@ -1365,7 +1413,7 @@ " \n", " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", @@ -1374,70 +1422,69 @@ " \n", " \n", " \n", - " \n", - " \n", + " \n", + " \n", " \n", - " \n", - " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", " \n", - " \n", - " \n", - " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", " \n", - " \n", - " \n", " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", " \n", " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", " \n", " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", " \n", " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", " \n", " \n", - " \n", - " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", " \n", - " \n", - " \n", + " \n", + " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", " \n", - " \n", - " \n", + " \n", + " \n", " \n", - " \n", + " \n", " \n", " \n", " \n", @@ -1717,14 +1776,14 @@ " \n", " \n", " \n", - " \n", - " \n", + " \n", + " \n", " \n", " \n", "\n" ], "text/plain": [ - "
" + "
" ] }, "metadata": {}, @@ -1770,12 +1829,12 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 4, "id": "f4e5462d-464d-4b09-8574-6867522b61ea", "metadata": {}, "outputs": [], "source": [ - "from optyx.core.backends import DiscopyBackend\n", + "from optyx.core.backends import PermanentBackend\n", "\n", "def to_float(x):\n", " if isinstance(x, complex):\n", @@ -1785,18 +1844,26 @@ "\n", "free_syms = list(expectation.free_symbols)\n", "\n", - "f_exp = lambda xs: to_float(expectation.lambdify(*free_syms)(*xs).eval().tensor.array)\n", + "f_exp = lambda xs: to_float(\n", + " expectation.lambdify(*free_syms)(*xs)\n", + " .eval(PermanentBackend())\n", + " .tensor\n", + " .array\n", + ")\n", "\n", "def d_f_exp(xs):\n", " return [\n", - " expectation.grad(s).lambdify(*free_syms)(*xs).eval().tensor.array\n", + " expectation.grad(s).lambdify(*free_syms)(*xs)\n", + " .eval(PermanentBackend())\n", + " .tensor\n", + " .array\n", " for s in free_syms\n", " ]" ] }, { "cell_type": "code", - "execution_count": 7, + "execution_count": 5, "id": "24355282-4c9c-40b3-b643-461981f9475f", "metadata": {}, "outputs": [], @@ -1827,7 +1894,7 @@ }, { "cell_type": "code", - "execution_count": 8, + "execution_count": 6, "id": "e96227b0-db8e-4871-8c12-b1d034a4bc86", "metadata": {}, "outputs": [ @@ -1835,7 +1902,7 @@ "name": "stderr", "output_type": "stream", "text": [ - "100%|██████████| 10/10 [04:56<00:00, 29.62s/it]\n" + "100%|██████████| 10/10 [03:00<00:00, 18.01s/it]\n" ] } ], diff --git a/optyx/core/backends.py b/optyx/core/backends.py index c289057..1d6d7a7 100644 --- a/optyx/core/backends.py +++ b/optyx/core/backends.py @@ -80,6 +80,7 @@ import perceval as pcvl from quimb.tensor import TensorNetwork from optyx.core.channel import Diagram, Ty, mode, bit +from optyx.core.path import Matrix from optyx.utils.utils import preprocess_quimb_tensors_safe @@ -308,7 +309,7 @@ class AbstractBackend(ABC): def _get_matrix( self, diagram: Diagram - ) -> np.ndarray: + ) -> Matrix: """ Get the matrix representation of the diagram. @@ -319,7 +320,7 @@ def _get_matrix( np.ndarray: The matrix representation of the diagram. """ try: - return diagram.to_path().array + return diagram.to_path() except NotImplementedError as error: raise NotImplementedError( "The diagram cannot be converted to a matrix. " + @@ -575,7 +576,7 @@ def eval( tensor_diagram = self._get_discopy_tensor(diagram) sim = pcvl.Simulator(self.perceval_backend) - matrix = self._get_matrix(diagram) + matrix = self._get_matrix(diagram).array if ( matrix.shape[0] != matrix.shape[1] or @@ -641,37 +642,46 @@ def eval( **extra: Additional arguments for the evaluation, including 'n_photons' for optyx.core.path.Matrix.eval. """ - try: - matrix = diagram.to_path() - except AssertionError as error: - raise ValueError( - "The diagram cannot be converted to a matrix. " + - "It is not linear optical." - ) from error n_photons = extra.get( "n_photons", 0 ) - if ( - len(matrix.creations) == 0 and - n_photons == 0 - ): - raise ValueError( - "The diagram does not include any photon creations. " + - "n_photons must be greater than 0." - ) + def check_creations(matrix): + if ( + len(matrix.creations) == 0 and + n_photons == 0 + ): + raise ValueError( + "The diagram does not include any photon creations. " + + "n_photons must be greater than 0." + ) tensor_diagram = self._get_discopy_tensor(diagram) - result = matrix.eval(n_photons=n_photons, as_tensor=True) - print(result) + if hasattr( + diagram, + "terms" + ): + result = 0 + for term in diagram.terms: + matrix = self._get_matrix(term) + check_creations(matrix) + result += matrix.eval( + n_photons=n_photons, + as_tensor=True + ).array + else: + matrix = self._get_matrix(diagram) + check_creations(matrix) + result = matrix.eval(n_photons=n_photons, as_tensor=True).array + return EvalResult( discopy_tensor.Box( "Result", tensor_diagram.dom, tensor_diagram.cod, - result.array + result ), output_types=diagram.cod, state_type=StateType.AMP diff --git a/test/test_backends.py b/test/test_backends.py index 038d99b..46b3fb8 100644 --- a/test/test_backends.py +++ b/test/test_backends.py @@ -445,7 +445,7 @@ def test_mixed_all_off_diagonal_mass_raises_zero_total(self): ) def test_non_lo_circuits_raise_notimplemented(self, circuit): backend = PermanentBackend() - with pytest.raises(ValueError): + with pytest.raises(AssertionError): _ = circuit.eval(backend) @pytest.mark.parametrize("circuit", PURE_CIRCUITS_TO_TEST) From 1b1fdd5422f136393db5de5935c8b140d84dc923 Mon Sep 17 00:00:00 2001 From: mateuszkupper Date: Fri, 10 Oct 2025 18:14:17 +0100 Subject: [PATCH 05/19] Added support for producing amplitudes with Perceval --- optyx/core/backends.py | 258 +++++++++++++++++++++++++++++++++-------- optyx/core/path.py | 4 +- test/test_backends.py | 20 ++-- 3 files changed, 224 insertions(+), 58 deletions(-) diff --git a/optyx/core/backends.py b/optyx/core/backends.py index 1d6d7a7..4cd84bf 100644 --- a/optyx/core/backends.py +++ b/optyx/core/backends.py @@ -43,7 +43,7 @@ >>> diag = Create(1, 1) >> BS >>> backend = QuimbBackend() >>> result = diag.eval(backend) ->>> np.round(result.prob((2, 0)), 1) +>>> np.round(result.single_prob((2, 0)), 1) 0.5 **Compressed contraction (hyper-optimiser reused across calls)** @@ -52,7 +52,7 @@ >>> opt = ReusableHyperCompressedOptimizer(max_repeats=32) >>> backend = QuimbBackend(hyperoptimiser=opt) >>> result = diag.eval(backend) ->>> np.round(result.prob((2, 0)), 1) +>>> np.round(result.single_prob((2, 0)), 1) 0.5 **Unitary circuit simulation with Perceval** @@ -60,10 +60,11 @@ >>> from optyx.core.backends import PercevalBackend >>> backend = PercevalBackend() >>> result = BS.eval(backend) ->>> np.round(result.prob((2, 0)), 1) +>>> np.round(result.single_prob((2, 0)), 1) 0.5 """ +import warnings from abc import ABC, abstractmethod from dataclasses import dataclass from typing import Any @@ -88,9 +89,11 @@ class StateType(Enum): """ Enum to represent the type of state represented by the result tensor. """ - AMP = "amp" # pure-state amplitudes - DM = "dm" # density matrix - PROB = "prob" # classical probability distribution + AMP = "amp" # pure-state amplitudes + DM = "dm" # density matrix + PROB = "prob" # classical probability distribution + SINGLE_PROB = "SINGLE_PROB" # single ||^2 probability + SINGLE_AMP = "SINGLE_AMP" # single amplitude @dataclass(frozen=True) @@ -126,7 +129,7 @@ def density_matrix(self) -> np.ndarray: ) if self.state_type not in {StateType.AMP, StateType.DM}: raise TypeError( - "Cannot get density matrix from probability distribution." + f"Cannot get density matrix from {self.state_type}." ) if self.state_type is StateType.AMP: density_matrix = self.tensor.dagger() >> self.tensor @@ -142,8 +145,7 @@ def amplitudes(self, normalise=True) -> dict[tuple[int, ...], complex]: """ if self.state_type != StateType.AMP: raise TypeError( - ("Cannot get amplitudes from density " + - "matrix or probability distribution.") + (f"Cannot get amplitudes from {self.state_type}.") ) if len(self.tensor.dom) != 0: raise ValueError( @@ -182,7 +184,7 @@ def prob_dist(self, round_digits: int = None) -> dict: "Must be StateType.AMP, StateType.DM, " + "or StateType.PROB.") - def prob(self, occupation: tuple) -> float: + def single_prob(self, occupation: tuple) -> float: """ Get the probability of a specific occupation configuration. @@ -192,9 +194,28 @@ def prob(self, occupation: tuple) -> float: Returns: float: The probability of the specified occupation configuration. """ + if self.state_type == StateType.SINGLE_PROB: + return float(self.tensor.array) + prob_dist = self.prob_dist() return prob_dist.get(occupation, 0.0) + def single_amplitude(self, occupation: tuple) -> complex: + """ + Get the amplitude of a specific occupation configuration. + + Args: + occupation: The occupation configuration to query. + + Returns: + complex: The amplitude of the specified occupation configuration. + """ + if self.state_type == StateType.SINGLE_AMP: + return complex(self.tensor.array) + + dic = self.amplitudes(normalise=False) + return dic.get(occupation, 0.0) + def _convert_array_to_dict( self, array: np.ndarray, @@ -539,6 +560,7 @@ def __init__(self, perceval_backend: pcvl.ABackend = None): else: self.perceval_backend = perceval_backend + # pylint: disable=too-many-locals def eval( self, diagram: Diagram, @@ -551,51 +573,93 @@ def eval( Args: diagram (Diagram): The diagram to evaluate. - **extra: Additional arguments for the evaluation, - including 'perceval_state'. - + **extra: Additional arguments for the evaluation: + - perceval_state: A `perceval.BasicState` or a sequence of + non-negative integers (occupation numbers). Defaults + to a bosonic product state |11...1> if the + diagram does not include any photon creations. + Either the creations for all input ports are specified + by the diagram (`Create(...)`) or the user must provide + the `perceval_state` argument covering all input ports. + - task: The Perceval task to perform. Allowed values are + "probs" (default), "amps", "single_amp", "single_prob". + - out: Required if task is "single_amp" or "single_prob". + A sequence of non-negative integers (occupation numbers) + specifying the output configuration for which to compute + the amplitude or probability. Returns: - The result of the evaluation. + The result of the evaluation (EvalResult). """ + task = extra.get("task", "probs") + tensor_diagram = self._get_discopy_tensor(diagram) + matrix = self._get_matrix(diagram) + perceval_state = extra.get( "perceval_state", - pcvl.StateVector([1] * len(diagram.dom)) + self._get_state_from_creations(matrix, diagram) ) + perceval_state = self._process_state(perceval_state) + perceval_effect = None + if task in ["single_amp", "single_prob"]: + perceval_effect = self._get_effect(extra) + + # pylint: disable=protected-access + if not matrix._umatrix_is_unitary(): + matrix, perceval_effect, perceval_state = self._dilate( + matrix, task, perceval_effect, perceval_state + ) - if not isinstance(perceval_state, pcvl.StateVector): - try: - perceval_state = pcvl.StateVector(list(perceval_state)) - except Exception as e: - raise TypeError( - "perceval_state must be a perceval.StateVector" - " or a sequence of non-negative " + - "integers (occupation numbers)." - ) from e + sim = pcvl.Simulator(self.perceval_backend) + perceval_circuit = self._umatrix_to_perceval_circuit(matrix.array) + sim.set_circuit(perceval_circuit) - tensor_diagram = self._get_discopy_tensor(diagram) + m_orig = len(diagram.dom) + k_extra = matrix.dom - m_orig + result = None + p = None + if task in ("probs", "amps"): + if task == "probs": + result = sim.probs(perceval_state) + result = {tuple(k): float(v) for k, v in result.items()} + return_type = StateType.PROB + else: + sv = sim.evolve(perceval_state) + result = {tuple(k): complex(v) for k, v in sv} + return_type = StateType.AMP + + result = self._post_select_vacuum(result, m_orig, k_extra) + shape = tensor_diagram.cod.inside + output_types = diagram.cod + output_tensor_cod = tensor_diagram.cod + + elif task in ("single_prob", "single_amp"): + if task == "single_prob": + p = float(sim.probability(perceval_state, perceval_effect)) + return_type = StateType.SINGLE_PROB + + else: + p = complex( + sim.prob_amplitude(perceval_state, perceval_effect) + ) + return_type = StateType.SINGLE_AMP - sim = pcvl.Simulator(self.perceval_backend) - matrix = self._get_matrix(diagram).array + shape = (1,) + output_types = None + output_tensor_cod = discopy_tensor.Dim(1) - if ( - matrix.shape[0] != matrix.shape[1] or - not np.allclose(np.eye(matrix.shape[0]), - matrix.dot(matrix.conj().T)) - ): + else: raise ValueError( - "The provided diagram does not represent a" + - " unitary operation or the matrix is not square." + "Invalid task. Allowed values are" + + " 'probs', 'amps', 'single_amp', 'single_prob'." ) - perceval_circuit = self._umatrix_to_perceval_circuit(matrix) - sim.set_circuit(perceval_circuit) - result = sim.probs(perceval_state) - result = {tuple(k): v for k, v in result.items()} - - array = np.zeros(tensor_diagram.cod.inside) + array = np.zeros( + shape, + dtype=float if task in ("single_prob", "probs") else complex + ) - if result: + if result and task in ["probs", "amps"]: configs = np.fromiter( (i for key in result for i in key), dtype=int, @@ -604,22 +668,124 @@ def eval( coeffs = np.fromiter( result.values(), - dtype=float, + dtype=float if task == "probs" else complex, count=len(result) ) array[tuple(configs.T)] = coeffs + elif task in ["single_prob", "single_amp"]: + array[0] = p + else: + pass + return EvalResult( discopy_tensor.Box( "Result", - tensor_diagram.dom**0, - tensor_diagram.cod, + discopy_tensor.Dim(1), + output_tensor_cod, array ), - output_types=diagram.cod, - state_type=StateType.PROB + output_types=output_types, + state_type=return_type + ) + + def _get_state_from_creations( + self, + matrix, + diagram + ): + if len(matrix.creations) == 0: + warnings.warn( + "The diagram does not include any photon creations. " + + "The default perceval_state will be " + + "the bosonic product state.", + UserWarning + ) + + return pcvl.BasicState([1] * len(diagram.dom)) + if len(matrix.creations) != len(diagram.dom): + raise ValueError( + "The number of photon creations in the diagram " + + "does not match the number of input modes." + ) + return pcvl.BasicState(matrix.creations) + + def _post_select_vacuum( + self, + dist, + m_orig, + k_extra + ): + """Keep only entries where extra (ancilla) + modes are all 0, then drop them.""" + if k_extra <= 0: + return dist + return { + k[:m_orig]: v + for k, v in dist.items() + if all(x == 0 for x in k[m_orig:]) + } + + def _get_effect(self, extra): + if "out" not in extra: + raise ValueError( + "The 'out' argument must be provided for " + + "task 'single_amp' or 'single_prob'." + ) + perceval_effect = extra["out"] + if not isinstance(perceval_effect, pcvl.BasicState): + try: + perceval_effect = pcvl.BasicState(list(perceval_effect)) + except Exception as e: + raise TypeError( + "perceval_effect must be a perceval.BasicState" + " or a sequence of non-negative " + + "integers (occupation numbers)." + ) from e + return perceval_effect + + def _process_state(self, perceval_state): + if not isinstance(perceval_state, pcvl.BasicState): + try: + perceval_state = pcvl.BasicState(list(perceval_state)) + except Exception as e: + raise TypeError( + "perceval_state must be a perceval.BasicState" + " or a sequence of non-negative " + + "integers (occupation numbers)." + ) from e + return perceval_state + + def _dilate( + self, + matrix, + task, + perceval_effect, + perceval_state + ): + warnings.warn( + "The provided matrix is not unitary. " + "PercevalBackend expects a unitary matrix. " + "Dilation will be used. " + "This can impact performance.", + UserWarning, + stacklevel=2 ) + input_state_len = matrix.dom + matrix = matrix.dilate() + if len(perceval_state) < matrix.dom: + perceval_state = pcvl.BasicState( + list(perceval_state) + + [0] * (matrix.dom - input_state_len) + ) + if task in ["single_amp", "single_prob"]: + if len(perceval_effect) < matrix.cod: + perceval_effect = pcvl.BasicState( + list(perceval_effect) + + [0] * (matrix.cod - len(perceval_effect)) + ) + return matrix, perceval_effect, perceval_state class PermanentBackend(AbstractBackend): diff --git a/optyx/core/path.py b/optyx/core/path.py index 2fdb68b..4f27061 100644 --- a/optyx/core/path.py +++ b/optyx/core/path.py @@ -434,7 +434,7 @@ def prob_with_perceval( ... (n_photons=1).round(1)== Probabilities[complex]( ... [0.9+0.j, 0.1+0.j, 0.1+0.j, 0.9+0.j], dom=2, cod=2) """ - if not self._umatrix_is_is_unitary(): + if not self._umatrix_is_unitary(): self = self.dilate() circ = self._umatrix_to_perceval_circuit() @@ -491,7 +491,7 @@ def _to_perceval_post_select(self) -> pcvl.PostSelect: ] return pcvl.PostSelect(" & ".join(post_str)) - def _umatrix_is_is_unitary(self) -> bool: + def _umatrix_is_unitary(self) -> bool: m = self.umatrix.array return np.allclose(np.eye(m.shape[0]), m.dot(m.conj().T)) diff --git a/test/test_backends.py b/test/test_backends.py index 46b3fb8..4b4f9d7 100644 --- a/test/test_backends.py +++ b/test/test_backends.py @@ -133,6 +133,13 @@ def test_pure_circuit(self, circuit): result_perceval.prob_dist(), ) + result_perceval = circuit.eval(backend, perceval_state=perceval_state, task="amps") + + assert dict_allclose( + result_quimb.amplitudes(), + result_perceval.amplitudes(), + ) + class TestDiscopyBackend: # compare with matrix.probs @pytest.mark.parametrize("circuit", PURE_CIRCUITS_TO_TEST) @@ -235,9 +242,9 @@ def test_prob_lookup_and_missing_key(self): ev = EvalResult(_tensor=box, output_types=(mode,), state_type=StateType.AMP) p = ev.prob_dist() - assert ev.prob((0,)) == pytest.approx(p[(0,)], rel=1e-12) - assert ev.prob((1,)) == pytest.approx(p[(1,)], rel=1e-12) - assert ev.prob((2,)) == 0.0 + assert ev.single_prob((0,)) == pytest.approx(p[(0,)], rel=1e-12) + assert ev.single_prob((1,)) == pytest.approx(p[(1,)], rel=1e-12) + assert ev.single_prob((2,)) == 0.0 def test_density_matrix_from_amp_is_outer_product(self): v = np.array([1/np.sqrt(3), np.sqrt(2/3)], dtype=complex) @@ -457,10 +464,6 @@ def test_zero_photons_no_creations_raises_valueerror(self, circuit): class TestPermanentBackendVsQuimb: @pytest.mark.parametrize("circuit", PURE_CIRCUITS_TO_TEST) def test_permanent_amp_matches_quimb_with_create(self, circuit): - """ - Diagram has explicit Create(...), so PermanentBackend.eval returns a **state** (AMP). - Compare amplitudes, probs, and raw arrays with Quimb. - """ state_occ = get_state(circuit) state = photonic.Create(*state_occ) diagram = state >> circuit @@ -476,9 +479,6 @@ def test_permanent_amp_matches_quimb_with_create(self, circuit): @pytest.mark.parametrize("circuit", PURE_CIRCUITS_TO_TEST) def test_permanent_prob_matches_quimb_with_create(self, circuit): - """ - Same as above but request PROB directly from PermanentBackend. - """ state_occ = get_state(circuit) state = photonic.Create(*state_occ) diagram = state >> circuit From 8d7c162a19c586de0ed27746ab24af5c75f8080a Mon Sep 17 00:00:00 2001 From: mateuszkupper Date: Tue, 14 Oct 2025 17:55:39 +0100 Subject: [PATCH 06/19] bose-hubbard added as an example --- docs/notebooks/bose_hubbard.ipynb | 3353 +++++++++++++++++++++++++++++ optyx/core/backends.py | 190 +- 2 files changed, 3461 insertions(+), 82 deletions(-) create mode 100644 docs/notebooks/bose_hubbard.ipynb diff --git a/docs/notebooks/bose_hubbard.ipynb b/docs/notebooks/bose_hubbard.ipynb new file mode 100644 index 0000000..808ec29 --- /dev/null +++ b/docs/notebooks/bose_hubbard.ipynb @@ -0,0 +1,3353 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "6c56b043", + "metadata": {}, + "source": [ + "# Bose-Hubbard model" + ] + }, + { + "cell_type": "markdown", + "id": "22cefce0", + "metadata": {}, + "source": [ + "This notebook builds a **Bose–Hubbard** Hamiltonian over a graph and minimizes its energy with a **photonic variational circuit** in Optyx. We’ll:\n", + " \n", + "1) prepare a photonic ansatz, \n", + "2) define creation/annihilation and number operators on a single photonic mode, \n", + "3) assemble the Hamiltonian \\(H\\) from a NetworkX graph using *function syntax*, \n", + "4) evaluate $E(\\boldsymbol\\theta)=\\langle\\psi|H|\\psi\\rangle$ and its gradients, \n", + "5) run an optimiser with a decaying learning rate, and \n", + "6) plot convergence" + ] + }, + { + "cell_type": "markdown", + "id": "11e739e8", + "metadata": {}, + "source": [ + "## Ansatz\n", + "Define the variational ansatz:" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "12d61180-27de-4d3a-bad6-73a0e94c8b16", + "metadata": {}, + "outputs": [ + { + "data": { + "image/svg+xml": [ + "\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " 2025-10-14T17:04:54.650833\n", + " image/svg+xml\n", + " \n", + " \n", + " Matplotlib v3.10.5, https://matplotlib.org/\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "\n" + ], + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "from optyx import photonic, classical\n", + "\n", + "circuit = photonic.Create(1, 1, 1) >> photonic.ansatz(3, 4) >> photonic.Id(2) @ classical.Select(1)\n", + "circuit.draw()" + ] + }, + { + "cell_type": "markdown", + "id": "685b16c4", + "metadata": {}, + "source": [ + "## The model" + ] + }, + { + "cell_type": "markdown", + "id": "a49ddabf", + "metadata": {}, + "source": [ + "Define the creation/annihilation channels for a single photonic mode; we’ll reuse these to place operators on specific lattice sites.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5736ad24", + "metadata": {}, + "outputs": [], + "source": [ + "from optyx import Channel\n", + "from optyx.core.diagram import mode\n", + "from optyx.core.zw import W, Create, Select\n", + "\n", + "creation_op = Channel(\n", + " \"a†\",\n", + " Create(1) @ mode >> W(2).dagger()\n", + ")\n", + "\n", + "annihilation_op = Channel(\n", + " \"a\",\n", + " W(2) >> Select(1) @ mode\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "9e720322", + "metadata": {}, + "source": [ + "We’ll construct the Bose–Hubbard Hamiltonian $H(t,U,\\mu)$ on an arbitrary graph using function syntax so each term edits only its intended wire(s).\n" + ] + }, + { + "cell_type": "markdown", + "id": "77676e33", + "metadata": {}, + "source": [ + "![Bose Hubbard model](./bose_hubbard.png \"Bose Hubbard model\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1fbf6443", + "metadata": {}, + "outputs": [], + "source": [ + "import networkx as nx\n", + "from optyx import Diagram, qmode\n", + "from optyx.photonic import NumOp, Scalar\n", + "\n", + "def bose_hubbard_from_graph(\n", + " graph: nx.Graph,\n", + " t: float,\n", + " mu: float,\n", + " U: float\n", + "):\n", + " nodes = sorted(graph.nodes())\n", + " idx = {u: i for i, u in enumerate(nodes)}\n", + " N = len(nodes)\n", + "\n", + " H = None\n", + "\n", + " # hopping: -t (a_i^dagger a_j + a_j^dagger a_i)\n", + " for u, v in graph.edges():\n", + " i, j = idx[u], idx[v]\n", + "\n", + " @Diagram.from_callable(dom=qmode**N, cod=qmode**N)\n", + " def hop_ij(*in_wires):\n", + " # a_i^dagger a_j\n", + " out = list(in_wires)\n", + " out[i] = creation_op(out[i])\n", + " out[j] = annihilation_op(out[j])\n", + " Scalar(-t)()\n", + " return tuple(out)\n", + "\n", + " @Diagram.from_callable(dom=qmode**N, cod=qmode**N)\n", + " def hop_ji(*in_wires):\n", + " # a_j^dagger a_i\n", + " out = list(in_wires)\n", + " out[j] = creation_op(out[j])\n", + " out[i] = annihilation_op(out[i])\n", + " Scalar(-t)()\n", + " return tuple(out)\n", + "\n", + " H = (hop_ij + hop_ji) if H is None else (H + hop_ij + hop_ji)\n", + "\n", + " # onsite interaction: (U/2) a_i^dagger a_i^dagger a_i a_i\n", + " for u in nodes:\n", + " i = idx[u]\n", + "\n", + " @Diagram.from_callable(dom=qmode**N, cod=qmode**N)\n", + " def quartic_i(*in_wires, i=i):\n", + " out = list(in_wires)\n", + " w = out[i]\n", + " w = creation_op(w)\n", + " w = creation_op(w)\n", + " w = annihilation_op(w)\n", + " w = annihilation_op(w)\n", + " out[i] = w\n", + " Scalar(U/2)()\n", + " return tuple(out)\n", + "\n", + " H = quartic_i if H is None else (H + quartic_i)\n", + "\n", + " # -mu n_i\n", + " for u in nodes:\n", + " i = idx[u]\n", + "\n", + " @Diagram.from_callable(dom=qmode**N, cod=qmode**N)\n", + " def n_i(*in_wires, i=i):\n", + " out = list(in_wires)\n", + " out[i] = NumOp()(out[i])\n", + " Scalar(-mu)()\n", + " return tuple(out)\n", + "\n", + " H = n_i if H is None else (H + n_i)\n", + "\n", + " return H\n" + ] + }, + { + "cell_type": "markdown", + "id": "770d29fb", + "metadata": {}, + "source": [ + "Assemble hopping, on-site interaction, and chemical-potential terms into one Diagram representing the full Hamiltonian on $qmode^{\\otimes N}$. Start with a 2-site chain and some parameters:\n" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "900d0c48", + "metadata": {}, + "outputs": [], + "source": [ + "import networkx as nx\n", + "\n", + "graph = nx.path_graph(2) # 2 sites\n", + "\n", + "t, U, mu = 0.10, 4.0, 2.0\n", + "\n", + "hamiltonian = bose_hubbard_from_graph(graph, t, mu, U)" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "bb26b7d4", + "metadata": {}, + "outputs": [ + { + "data": { + "image/svg+xml": [ + "\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " 2025-10-14T17:04:54.764399\n", + " image/svg+xml\n", + " \n", + " \n", + " Matplotlib v3.10.5, https://matplotlib.org/\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "\n" + ], + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "hamiltonian.draw()" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "8a02308c-aba4-4a39-9325-faee4c1f0c2d", + "metadata": {}, + "outputs": [], + "source": [ + "expectation = circuit >> hamiltonian >> circuit.dagger()" + ] + }, + { + "cell_type": "markdown", + "id": "c76edfb7", + "metadata": {}, + "source": [ + "## Optimisation" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "id": "f4e5462d-464d-4b09-8574-6867522b61ea", + "metadata": {}, + "outputs": [], + "source": [ + "from optyx.core.backends import PermanentBackend\n", + "\n", + "\n", + "def to_float(x):\n", + " if isinstance(x, complex):\n", + " assert x.imag < 1e-8, x\n", + " return x.real\n", + " return x\n", + "\n", + "free_syms = list(expectation.free_symbols)\n", + "\n", + "f_exp = lambda xs: to_float(\n", + " expectation.lambdify(*free_syms)(*xs)\n", + " .eval(PermanentBackend())\n", + " .tensor\n", + " .array\n", + ")\n", + "\n", + "def d_f_exp(xs):\n", + " return [\n", + " expectation.grad(s).lambdify(*free_syms)(*xs)\n", + " .eval(PermanentBackend())\n", + " .tensor\n", + " .array\n", + " for s in free_syms\n", + " ]" + ] + }, + { + "cell_type": "markdown", + "id": "5adb53ee", + "metadata": {}, + "source": [ + "Run a short gradient-descent loop with a decaying learning rate to drive the energy down." + ] + }, + { + "cell_type": "code", + "execution_count": 38, + "id": "24355282-4c9c-40b3-b643-461981f9475f", + "metadata": {}, + "outputs": [], + "source": [ + "from tqdm import tqdm\n", + "\n", + "xs = []\n", + "fxs = []\n", + "dfxs = []\n", + "\n", + "def optimize(x0):\n", + " x = x0\n", + " lr = 5\n", + " steps = 10\n", + " for _ in tqdm(range(steps)):\n", + " fx = f_exp(x)\n", + " dfx = d_f_exp(x)\n", + "\n", + " xs.append(x[::])\n", + " fxs.append(fx)\n", + " dfxs.append(dfx)\n", + " for i, dfxx in enumerate(dfx):\n", + " x[i] = to_float(x[i] - lr * dfxx)\n", + "\n", + " lr *= 0.2*(i**(1/6)) # make lr smaller with each step\n", + "\n", + " xs.append(x[::])\n", + " fxs.append(f_exp(x))\n", + " dfxs.append(d_f_exp(x))" + ] + }, + { + "cell_type": "code", + "execution_count": 39, + "id": "cde44850", + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "100%|██████████| 10/10 [05:31<00:00, 33.16s/it]\n" + ] + } + ], + "source": [ + "optimize([2]*len(free_syms))" + ] + }, + { + "cell_type": "code", + "execution_count": 40, + "id": "088c8fa8-76b3-45bd-a5dd-b337be0d33b1", + "metadata": {}, + "outputs": [ + { + "data": { + "image/svg+xml": [ + "\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " 2025-10-14T17:53:51.051493\n", + " image/svg+xml\n", + " \n", + " \n", + " Matplotlib v3.10.5, https://matplotlib.org/\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "\n" + ], + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "import matplotlib.pyplot as plt\n", + "import seaborn as sns\n", + "sns.set_theme()\n", + "\n", + "fig, axs = plt.subplots(1, 1, figsize=(10, 7))\n", + "\n", + "axs.plot(range(len(xs)), fxs, c=\"#0072B2\", marker='o')\n", + "axs.set_xlabel('Iteration', fontsize=18)\n", + "axs.set_ylabel('Expected Energy', fontsize=18)\n", + "axs.grid(True)\n", + "axs.tick_params(axis='both', which='major', labelsize=16)\n", + "axs.tick_params(axis='both', which='minor', labelsize=16)\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "92dfc684", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": ".venv", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.12.3" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/optyx/core/backends.py b/optyx/core/backends.py index 4cd84bf..82f323c 100644 --- a/optyx/core/backends.py +++ b/optyx/core/backends.py @@ -593,91 +593,24 @@ def eval( task = extra.get("task", "probs") tensor_diagram = self._get_discopy_tensor(diagram) - matrix = self._get_matrix(diagram) - - perceval_state = extra.get( - "perceval_state", - self._get_state_from_creations(matrix, diagram) - ) - perceval_state = self._process_state(perceval_state) - perceval_effect = None - if task in ["single_amp", "single_prob"]: - perceval_effect = self._get_effect(extra) - - # pylint: disable=protected-access - if not matrix._umatrix_is_unitary(): - matrix, perceval_effect, perceval_state = self._dilate( - matrix, task, perceval_effect, perceval_state - ) - - sim = pcvl.Simulator(self.perceval_backend) - perceval_circuit = self._umatrix_to_perceval_circuit(matrix.array) - sim.set_circuit(perceval_circuit) - - m_orig = len(diagram.dom) - k_extra = matrix.dom - m_orig - result = None - p = None - if task in ("probs", "amps"): - if task == "probs": - result = sim.probs(perceval_state) - result = {tuple(k): float(v) for k, v in result.items()} - return_type = StateType.PROB - else: - sv = sim.evolve(perceval_state) - result = {tuple(k): complex(v) for k, v in sv} - return_type = StateType.AMP - - result = self._post_select_vacuum(result, m_orig, k_extra) - shape = tensor_diagram.cod.inside - output_types = diagram.cod - output_tensor_cod = tensor_diagram.cod - - elif task in ("single_prob", "single_amp"): - if task == "single_prob": - p = float(sim.probability(perceval_state, perceval_effect)) - return_type = StateType.SINGLE_PROB - - else: - p = complex( - sim.prob_amplitude(perceval_state, perceval_effect) - ) - return_type = StateType.SINGLE_AMP - - shape = (1,) - output_types = None - output_tensor_cod = discopy_tensor.Dim(1) + if hasattr( + diagram, + "terms" + ): + array = 0 + for term in diagram.terms: + arr, output_types, output_tensor_cod, return_type = \ + self._process_term( + term, diagram, task, tensor_diagram, extra + ) + array += arr else: - raise ValueError( - "Invalid task. Allowed values are" + - " 'probs', 'amps', 'single_amp', 'single_prob'." - ) - - array = np.zeros( - shape, - dtype=float if task in ("single_prob", "probs") else complex - ) - - if result and task in ["probs", "amps"]: - configs = np.fromiter( - (i for key in result for i in key), - dtype=int, - count=len(result) * array.ndim - ).reshape(len(result), array.ndim) - - coeffs = np.fromiter( - result.values(), - dtype=float if task == "probs" else complex, - count=len(result) - ) + array, output_types, output_tensor_cod, return_type = \ + self._process_term( + diagram, diagram, task, tensor_diagram, extra + ) - array[tuple(configs.T)] = coeffs - - elif task in ["single_prob", "single_amp"]: - array[0] = p - else: - pass return EvalResult( discopy_tensor.Box( @@ -787,6 +720,99 @@ def _dilate( ) return matrix, perceval_effect, perceval_state + def _process_term(self, term, diagram, task, tensor_diagram, extra): + """ + Process a term in a sum of diagrams. + + Args: + term (discopy.tensor.Diagram): The term to process. + """ + matrix = self._get_matrix(term) + + perceval_state = extra.get( + "perceval_state", + self._get_state_from_creations(matrix, diagram) + ) + perceval_state = self._process_state(perceval_state) + perceval_effect = None + if task in ["single_amp", "single_prob"]: + perceval_effect = self._get_effect(extra) + + # pylint: disable=protected-access + if not matrix._umatrix_is_unitary(): + matrix, perceval_effect, perceval_state = self._dilate( + matrix, task, perceval_effect, perceval_state + ) + + sim = pcvl.Simulator(self.perceval_backend) + perceval_circuit = self._umatrix_to_perceval_circuit(matrix.array) + sim.set_circuit(perceval_circuit) + + m_orig = len(diagram.dom) + k_extra = matrix.dom - m_orig + result = None + p = None + if task in ("probs", "amps"): + if task == "probs": + result = sim.probs(perceval_state) + result = {tuple(k): float(v) for k, v in result.items()} + return_type = StateType.PROB + else: + sv = sim.evolve(perceval_state) + result = {tuple(k): complex(v) for k, v in sv} + return_type = StateType.AMP + + result = self._post_select_vacuum(result, m_orig, k_extra) + shape = tensor_diagram.cod.inside + output_types = diagram.cod + output_tensor_cod = tensor_diagram.cod + + elif task in ("single_prob", "single_amp"): + if task == "single_prob": + p = float(sim.probability(perceval_state, perceval_effect)) + return_type = StateType.SINGLE_PROB + + else: + p = complex( + sim.prob_amplitude(perceval_state, perceval_effect) + ) + return_type = StateType.SINGLE_AMP + + shape = (1,) + output_types = None + output_tensor_cod = discopy_tensor.Dim(1) + + else: + raise ValueError( + "Invalid task. Allowed values are" + + " 'probs', 'amps', 'single_amp', 'single_prob'." + ) + + array = np.zeros( + shape, + dtype=float if task in ("single_prob", "probs") else complex + ) + + if result and task in ["probs", "amps"]: + configs = np.fromiter( + (i for key in result for i in key), + dtype=int, + count=len(result) * array.ndim + ).reshape(len(result), array.ndim) + + coeffs = np.fromiter( + result.values(), + dtype=float if task == "probs" else complex, + count=len(result) + ) + + array[tuple(configs.T)] = coeffs + + elif task in ["single_prob", "single_amp"]: + array[0] = p + else: + pass + return array, output_types, output_tensor_cod, return_type class PermanentBackend(AbstractBackend): """ From 59c1f311cf99ddc0f4f1111060cf2348d3e957d6 Mon Sep 17 00:00:00 2001 From: mateuszkupper Date: Tue, 14 Oct 2025 17:55:53 +0100 Subject: [PATCH 07/19] add bose hubbard --- docs/notebooks/bose_hubbard.png | Bin 0 -> 34908 bytes 1 file changed, 0 insertions(+), 0 deletions(-) create mode 100644 docs/notebooks/bose_hubbard.png diff --git a/docs/notebooks/bose_hubbard.png b/docs/notebooks/bose_hubbard.png new file mode 100644 index 0000000000000000000000000000000000000000..ea43a2e5dc3c64fd8ac36131d31c071db4585c70 GIT binary patch literal 34908 zcmd?Rhg*(+|37|35|UI>X^9frON)e3DO9Ar_tf6y9YULyhN7iZno?3mOIy2364Fu$ zO)cNY={`Qk@%tZsj{AP!Z{1hdb)M(z^?E)Z>v`R|c>e6p?X25LB+^b5k@3&l~?k>2^xx5;gw#Q(N4|-&x!hb=&0&eE6mX=O# zHqP#|n_r#4iv)-lop!Y}b+>hP;<{w(Xi2){VaX*T&UL}uol8_$RDw%bLI(d47vZ|7 z&UH%b(g~|h8WM?%q;f{#l6T7Vke30&(%-E=ubd)X>6MwPdy%A=&5`A(uO(xD>)M=Pu*)Fed=gD zu&a8+b*bjt@5-vBfRIp;hLhC$Z_u6k-#>~IKaeT@?=SfAy(BK|*rzxn_3zuKg* zlm0y|4UUevfhsSd=D4G`=y_BhI*?0BUWKJolP>1<(orfTdAhk>fBg7Rk6m|%#nO+V zFt1NOTQBEnZ!xF&@#9C=qB!eLN`-65$@|V|+~5oSIX{0x>2yL0u?$V}m%>HHrn&yD z$8A@s$l3~*5n;<_#;{@iEY&@&_eC}vzHk>=8JrvM(8$vkkn{O{@a)v@-@o6cD~7!< zchnCX;rjPJ1vxB-iAz=HI6xhI`Hm!s@A4DIPeFeJVuED-79AoP97z{ndb{*Ix!m$Z z>ffii9(_VqBL3~()~u{7g)GheJUm;xetdlMrS_JH>o;YxXz!EDwKujgTueW+ttVG& zYnEzUuwFq(RaI4!RiM=%$cv^2}13#oNU681&b9UPOh92^`I-H&(S?;$^bUXA6~ zrN>V$=3ep~HvYHet|+kE{M%pO6t?VOq9if$T_k0lOK9kj|Lahu*8kM7v)XTo?NF_p zoE$0h(4?plOOvFKkPyE?1>MYtXIeJ9{&(ev?37ghU3ufS;A=65F0i)k+OeZFw@lcG)C=f+d)HfL&-cKiaT5+McR#uM90P!zJ1RTooYeiZlYh_SaHC| zFfg0bQBoX<+wSD#)J}Wr8$~)_un1M0s%6QD&#CFjx;qR?b8p|ibvQCisSs_~@@I8b zOIL-stl#qP&n{5!=Uw~VBJ>kW+Q=D1qrv!qowjN8jK&%4BaKQ|%kX{2+R0y{{goIM zQ)E1vZE6bIPNd2Ca2NCs52q=y@emjOhf}qD_bG!AzR=H7w7%7!CSCXqUobW`H4T*6 zGSIR~eeUsA{&9$UQ(vKR)VFUIY|^fjavid=2dmyK?-XlDVc*_V*YlZrv?bx}GrL{~ zz1hx;vlN+2{Y<-en|8nW-vgxo_X)vV)H=`X_9Thh2H(C-{rFRU-lZenK0YeliGMHL zrvGgwrm4Xw;_~$j{+Nt5$I(%~A@9HYSK8-yrbh5FKE8L^R%ziPxttZ(u8|lIYu05m zZzJ~XZ~E`05p-YO^$M6;;st6yzqw993R+)C=uSS+o1?Lr^m%gf&XLPcq$HT`t>)qb zeSJ?o?tatV-A(>7a8s{9=@u<}kzMbR`I*7J&Yxbs8BYs}4BoOsMz+jfD1%DRs+@Na z4?Zbc(&^*<75jhtqUqr5alJz~z8wkXGWqnBnLEN_H8BE3>59pR!iKpYuB#F^#`aoNkzuk zeSW$6ql$D5_40I<=6@nospGjJv`r!!jQgpBV=YA0%C`lRQpe{k=v6a1O+J*_Y6_eP z3JSup^#O3&7alvjW&HqrrMYBDN4sucEltg< zEAzG_qSmzZDDx7>)Md|+24*);x*O4J*T;R%C@GmWgtL&YHN_kHH-W@mqxe{F*QdwC)w>Sv$6^A~(tPk1b`ZR$u8 z2Y7Bs(U{ZhlXQrtcYGl0)o^?F0mgG>x`S*NGCChR57&^oWLG%EtY|!ce7ut&WNbFx z{y?Y1Vw>;MkDaK?;tm6J3JMCLbvrVHcS%PBduDZ=oWSv_mjWA)^}EUIp4T-OZK>x9~$WW69F$SB|A>_73fr4-Os*6`+H@< zUN?Q`wYJf{<}^6mo1%@3#-6bz{;Scd56`w7G&T!Yc-$@hQZkN<`h2B`M(!oS&goz6 z4`jQpQ4_5QGy=ge;?{U^+PWkU-?zkG&_}@WUG)J{y;(Uy*Hh|T@ zxfBr#+VCfxEYN#KbYtNLF6vsj8{FG(WMXFPN$9> zVXHB$Hkys~6T#{xRC2UNO~f4f(-pDr*zev5qf3_1cB#%2(=fzwohRwq14E3n%=SFMcm9Fsk;Y4YJAVpW4mD6jEy4A>&@vPESuCym>3_ z`!@RhH~36|QZ()3*GxK6kK2$-vPx^;J^$Tku5gM5n_otj@>bF@=TxoKuNI!)dk*L8 z6mL4|Yit8drE$=w-i1vmocTZ9QY)3YmoHyNM@O4zX+4|jFWgKLG^lv!9doSOD2i$) z>*wF2Eu?FM&z&g(69?_c`3O)U84^&bx9eriy; ziB#dc@PymN%5$uhiU<**`zlwvUpMl1oHnH;Img z<<4QPY!YcZ6aVd?jXztRs~6(A2f5?Vsi~d0dQ~`8K8TjYeB4(&`}~9Lq@QzhwZBFh z3%=7B74~PGxkK+gHx^;mC?9PU5B#>(@k1E}>E-(3AVGryS7(lqCfB()R;CICJHyw2 zV2EEY(0@)M1*|UC-{U!#(-8T9xJR247gl zKg7s&jI;>7ifuonyN#Oq^?ecZe2$m;aKD`l^2%_x+VPsF99w0`d|>>ydJCdqmUGAMR_`L!+h z^%)I|U@qav5#~mh>9*u!{D%+UnHheSrv5-C?wRv&Dq6gnRHnymZip+)?A2m(xFM8==@ryZoj-Rjrf)YZYqNLib7>hjjcODBBIpC#R{~K#ub)YwYfG|-W6#` zmGgPKw!Yq$x`b|=@y@x8m1TT3+^ypLV{_0J)5b`4^Drg>vufW3GxQ&6nV`)k*6k@S z&4FM=My1%XV`F3Sj*0u2ncs292TsQ>s|{^F;kOu$VsK}MTzX4UNvZjVR%`VO&zWXQ zg{#dq-@IZS4*dCmHEw0P2rwZj>+1l~O<35YCul|-`AxcjKtl)%%T2=||VB){( z9qwD5x-c_%CG?SYfnNBCRB~T#_oNQaWu}MIrL7Zyt#;bu5MCQtQF8_HgqjrnMpb5)zjmDHoLbEjf+- z_~Ep;7BD}d@wg#zj{@C&qe5#mWHazA^WMBG9iGMn6ne~aY}@{WLSKA7O?*kqy)s-A zVEFYzQ#4;F9)Rc455sbcCi=yV)!F2sgc~;~t#3UL(kO}eW-r%KT5R5!B=8X6a_rts zN2_$ck}9v+5*9J5BK#=nW7Ox*Cdtmf1uS*)9G8E7kt~QBk?N~-v$hR5mmnzNQ~I>m z&a$Pzuu8IVkdE@dXSUr+zJ&!%eap+<+Mvq1t=ad2t~@ok%kvRW)2iE31#rntX(jN< z<-7(!Y|FB6F6!yEKXWCS#XOe5xUF@V-HD{V6@6F=m=nCZI8{yHr{dp@{iy}k@TGRP&)0Xf>b09q0uYi`UJ9P{Y}l#Xcr3s$dJU4!RQsa8-< z!`4^~B%XXc2%Ax^KlZWAgUAR@`PAx#K7FOt-7GAV(?3C=I*1;cfBAZD0+f}ourN2z z9BFRqO|t&W!uCTio^=Ms9@YvOs`l$A?;n%JLA@aoA07Srh<3isVsHS6SUasmBg!V0 z&MYu6uq0r4zWV~ZXKldRvJ1I<^QKLVDmuUAIy`4S#AWbG;GPX%1rDVzgr&bb@psMb z@Q+u2HwHgkv~)?4^ASDU0$7+(sFK*P^E=TPulxEf8fR_gy4I79a$vPEGtanNTY|ZY zuwiw{SMdqfx8lZ$wv&@MT4LFfuXp9tso+afScpnNH-3?Z0qh)0iQrOS;F}LT_u8aA zz!mz+Y_-TImRD9Z&64O&F>;4;>sR|eBhQ`1lJA*%Mc!Rv$14l8X?-hxcHjgf>wce; zsC9G<4CbDdI7q4)$}eQV9$&vcnIh$~Y4`5kU28Lg&qshEJcn!eLYF5Wr&-Jql8~a3 zmYbW~X%H1x@l(-NRkC({`L@3#MWd+o%4{^khMkW@g4m7Q=w7&RA%Xe0sC@A^u|nO_ z?8_O_Ob~2Hue&(ib=Ly5rb-W5U14duSfHIjY;_;`_4%0I0}}Mzp?igB<>lqca`w#| zVS9zr{g!?xtEtriR{R`IFn-ZmOi7A6dZpo8Q%sxSF{h6e25uWv*Rbusynl3hf?WDo zd=^!3pxFEtLk2fw7Es7Xc#CCp`BK9lz)#DIgU)YUQm%8CMnF`nS^Xh9{i%8%^r?Ge z;4DKn-yF=0u_GuaZ#n|O2;01W=#Wqquqp|<6=T0$pY-$R&yxlQY*f2AXh`0R(=3wC zpB$EEvNJPH(3VPk@YJ5KddA(8els>XIa#NyYy>cMvB*S$bn?-un-8_eXS!R4y=9DD&ij4(+?BC2`q+6tYj??4e7k;)mR*KRMn-1P(rLIR z0&4~RLYnW$pshG-zTUI#K!_l)a(B?50qN(0^Ef}Q_O%* z0C@ww;vXDEdB2`0W z+qXA+Rzm4A`QDRT`J|}VM*O~nJ=3Pmn<+@Lp3_BsZ-1g_lq8Qs)gk!o)4@xcV@ZxT z`1T(lYiBC&K6Li0e*~+6m!12Qst<%QV}pI)+bKfcg$c^ea$OqugA z$KAVkNnQ55a^62*H))JIxwyF4bf#~iqG7&!~Iw1eV?3GX_?P-wi!_DTj{q@amVh(sr93 zBqZ9YzJSzOs={T);^pO)Z~V$ga`qNy;7P!gtABo5{*>FLK*vSBNmEPfZgliDX}7V5 z!0XJ%oY(-kJ70L%^BdO~^fA`&sjR9pZ;m~jqnV{z-V(ts%M`Ts>lU|qC3}5a+fES? z5#r`SheAN0qz-N@FE7g_49I$XXCxj-PtT=r7BQ;X>ZMfmDF(&eSi)=Cd6)T4eN`nM z5HhxBvM`PG;D zu)Utc*?ghQ;QJ`@L(fOM*g8d0Q0CbhRsLE90!fHOjt~&Q*47p8h%3Com5gXURh3tdf!@%jg-F@p$Ww)MESor^NvOD#LAP`vVu`*DdCjYUkwa0rfY( zrlMuLpYi56{~aHwH=3>Xefg$n&8~miKks6fIe7GdYAoN~JQ=Qb^{%e2%+51nMqv7_ zEiGYNsHjApot?8gZ7hOz(b8JsmzpIqysdOfts=;C&kGBMs@%pSOHV%=z}@O}rYj~E z78VX|?Yt8el|YpWzIE5|`6%Dv!-w21dHLSIzUBjl@{5aKTHH=c8^6IxPtU!+4@f<5ec=jug-OVWUtjW2 zVYp>~p|Nv^A&*ZuK8`bJt(g3_Jm~L6<({T6j@(l> zw}%Ft*DHvZuzAllI6I1VA?a+rx;PWcGYINz;#nrgW5(Mj4I80+SUp!W>%Mm0GuwT* zcA#ESlEb%mcO@I;fX*GvX}4{2`Lj8LI^tLftlX;!Gk)_^Ax zJ|7v6YKKce!WA&6h_-#EmmKr9{q>yO>R7UKYw9Z$z8}9|?LLj|h6h*@9Zk155+(n^ zQ8v3>7c8qSuJLpTZnE^lGkY;ozDf))KMzlR_5zqy-&4b*tKoXj?Bb>(bJw^txNkjc z{JFlmlr=Jr;+=UzpM!v~+KtPO z3{IZvPZ%hCT-uJ*>rM~R|sltw6~XSfBXxdIo>>tj1aQ9{bW}bo^YE%)g@aX zU-%*KPGmrgI=aYW_|jWsa5e)^S@Mt)s2olP4GDE9y6weM7^CGQrew-73$Sk_eyo4;2KB1@EN4TUnYF3adhg6y)@B zg7&FvXLsbBva$(UnAkoEJI2*N<7wndw{ac(98zTEhAeS732@jkTwGjS`QyX;)N$vA za0W+~ZB!(JU{?MfcyIS?TK?~!s4`Bp#8TS;tN?yPWh(DVVe0%}U~SE}ZS_!suESag zjT1CXf|1kmeF61fcO~^~6vu-pj`En4mYm|Gp*)MGd(#l1+G%^agSnIiw$DK4kc2`M zV`+}PHD^iA`H#J_vhoo~icYcF79a--lJ~;*J>aJ%We$)MPMtbM%3Mldq^CFQNjrz$ ze-S1cA>+cG2^UcVEVnTP2ByYg=-b1}8VZvtf%QqM&aDxtoL)N!Y$)kfPIP1E*Lr?#A=$c8z≦M3Ce=`Y(&1Q|gcrI=jBIpz$gYRIRgiG-}T|A8&7>@mF~N;s)7G8Fd{FSYr1N2xz0t|ExOwXYh7Z z)LAHIS#EaFvx%m0;lv+GboJB9%G8k@$3wBV$CAupI@Q(HDQ6Td!CT?<|M{Q#l>STa zVsYO9JZ!M=>kuBza6c6uQ%rKlk3Z_Qvpwi`h+c7b6sLT{%EDAJ7>Gf&?;etaj1!vF z?RW2#2^m=H@p(l^>7CD=hYjZRxupAZwGK^IwKqfG6ZV|eaTqLb<2gGisaT*rSn$^Y zS|(6$IIG0=lfGmruO9}W928pFYHt>228oUe>?hmrZxnjf-E$V_TPg7i?4|hlc&lG8 zP{qUhMnQp$UIjYddUvPPu;QwM@7wP=<$(UByALSc4|vh}==3fB#DZgrdA+O>@+jEP zTgp*HfzT5i?ePdh<$S=lnKr+MV)AgJI2-52#wLe<3<`h|A&B0WDdhY0O%sRl>gSVG z*Bl3m6#i5GI}8oQ;#a(&q2ApO48x@(|CJMIZ`b3AbEKD)H*4FR2^Yde{Ruc4?XuT?kJ z*OSW%j$V0Mv}wzB<)y-T_HjdyF)-n?t_Lu;=fSk;$7khwN+{B zBnnc&3+WRc-(x#pyl|hqXD>)nZU;l{g9c9aa#RoQ``~!?%$fTF4vi7}5@vVkH!xR$ z%}u-g`uVd+pZng3ly_h(m@>H%ZlygxN`qiRvyB>B#THuO` zUO!j}xg|l!vgImUBQ%2p0t}3d--_&e*Hjr57x5-Z>k&siBKShNr-rJYm({>rs92l5 zrZY>}KT9mxMg+0)%sFrV#KC{zdz=Fcu+bD}C9ugK-Y-^ZzBkAf{(nk3U!q3jJ^9Qb z2$0ReQ^F3 z>eeF=7iCpdS5GM=xWin@PJQ`_W{u8Y-rn4!8}$gyuO%CBDb7BL2L=F#VkqO?wy>xu zFHP|=V8BNJSog)OX{93FJf-Srz)B$Q5+Ii>_}pno#{cj75}jhF!re;R2IfU%A?oVR zfT>df`%HXzme=X=N8h>nCN7=02-3A3dj%`j#y#K~G|$)C;c{h z|7H^O_ZxWgM=1CumIzWn0d(0|^XYkNxDPnzHHu8ZC*CHOTJvXsJIkN1?-EACsB$>! zwj_!Ry=hc<#b35gf^e!Bc-4f7hZsKn;e`BdN(vMt5{ZUGju||p5f{ay-!L@9{^re_ zffpV}ugwpu2Uc2p+GwlNQ4)fNO&~Ou0=<4HWf5c6KWan1u|u~_Q18#z?S2$+kdlgu zn;;mjot>S8vB6Gv*PdRnlW*IPsVV#Cl*@2)ExIAk@)^CX;GVJyry!v5 z*;i3nX*xI7Hc(=@1Iozw^i_*98sbOnp@L-({DPA@yNRTlAXxBjg>(&;v-uMM=F03{ zEMQKT3XKdUXCccA`Y{x0<`PZ~%BMiTf5RU`6a~WF&B%QE2nJVOeLXm5|0B&@VY51F z!hM8dN)9{^gsTIYhx@<*!nsfQp|8QHNK{yeGnCPvqAWqN@{Q7<5)kzs;BYsXPIF1e zen<_lmds!!0lD$O5p-OTRB68fq5){FtgQSR3e_((L_(YLtUBZX{EhpeB1yS^qec)x z$WZLwOT{eBFM#o9|8A^<*>`1~yB!|hny~8&w!JV?D}+(>@6#@L;L|{&P0*n#fS_c& zeoTK>eo*_RsXXIsRP&gRkIy44l6Jw=KMS_DwzgnrR0W>`9cbGY^2_({FWxfxhUfrc z@-4%`!iV9J{V6Zn(k0}?Ku>Qm2VQKRkgVBoU&BlYT6l5YZKaI=j0zhGLsYg;?f5~P zuFP${4IaDYI`-}SVfs`loVewkI8SIfN6Xv6S-yaGmiR#B3oRJ_bm>T(#wl45kxt7P z$AeyWAIc7!nVXjyqM`aDo}%f>E8htY4sB`ErcJGaQ=pnw-$IhJ+PrB4yUl9o!EN4QSTSJ=#R7YXQ58|aaXRxNi0P_C21EJ}+ zp>Cg(?C9_B7aXLgq=3B|-)0X_F8$6?04qiN2Z^k?5{Hf*s}q4q7My&{S+$#x)$Mw6 z9u2HWuoSvmjrxZfgS4W17V8PEsf7V-K0 zjjcEOQPiL9!u)`$i^w{T4!T~vzQyA%;7%Tyax^HOp>5Z z5m)f_>ksJ9HQnQJGZWxdoWynorkMMu<`D#pTKl5=rWS5TM#c|4 zzqi(Ans^xd|L13gj6Zz`)hQK(rKL_`%j9p1sHlLV(v%I)?`schI5sRS7_;&gn`e4w z&Ym4Xfujv=maurbC-#j|lVok3GT@GQ&VPBLYT?U&{*2~5mVFAUfl$WU+M2-%MI;xs z6o7z(Q4FuT>yU`$eFp|b5o5{G1B+KmESlJshONFQ#eL&QoMi|O4o=xi-4xSfCqcMR z>d=)_lkSbVz?Z*HwNrF3>pebO1Qzv^{putsVX`MO7#V@ju|57nuMF$|*XhC@;=pOS z<;%gC%b<)*OH0#Q-L8?>%Pe)3i(p_guAq%dmo8nZIpXklZSF#b4|G9If7x)6j5=kY zR>0IZ%$qyrl-1Q0_4HWL#qK1OSgpPZujI{f$p#D9DDi z#0Hk$UfTbgWl4aG|$m&XRU?) z9eSSWF89+QMRGoxl0=jNw{-i?oi|XEc5`w<47*rpbexC}W++9R!U22`pA)dYB0+$h zjmk!YgkkTQ79Igpot$*9_9%=Q_|r}z5GGjT&`kISp8BhY+_+!+oP9VcfL`a}jO zZw{OvGc&W453M~{Uq1q2IR$cs)Z73h6{5ZX*FFT>$kBf;2ys-{^ettAJ*)l;cUvNL zs8^t!bhb&-Oo%mt>hvwzdK|eTuzUUlczl$221HGh^9t}m1#tLSu^VoHEJU*}Gcyse z1I_Ag@8-Ms@LCyL#~kuTwGgwf&U}zaIqt^TQuH4*TKzPo$TVNYiu0Bqa)~DedwOu zbadCiRo*r>9x^uAI(7&Fu&%$x$fW?kbe{~6j>cYt`0*7FFGnMT3YdKtfLFuxDmIGt zQ-eDw;2F48JddZF(}RQs~nTK(HXB+4*9~Yp3lM>;sDhobQ^L}`xq$SH2E$3^-i_s3|xN1XsoP!yM8=Cfo=d@ClU78 z=GRXy<#HDkp8RhUL5{-5eH!NE{NM@@3}t=^LK4e^7VN%yoc-B*mU>1=09eoaMTUiG zIjlo;1iw(8N|o_A63}S%Q?4U~mi@@;d7OgQ8U%q#@ge7CZBabBpN_Mf8;$}AgL#LH zr?RD;5i~gcVzZlNLD|hp3Uqk>gj&VV$EU0<8FJ&syAx@yqs_&j6taHLrd=3a#FwAw(b*3$skKej+gINF&?Jc-`3X>O* zVjELAfJ>9s&4R-an;^8kTdabj=UKUb9RgPATK-mAws(StFY2Ts%|U$$hkS7fitu5M z3+F97GzE=bN}-ps5|3v@D)KmLdAS3Y61}G*g}(IoRkm?g;xZ#5hr;L7z<0^J3b1E)Nb#%c5yd zj)tJe=Cu6+0ulEdeUJfn*jW1mRSkQ$RjmI$R-o|JMiRy3u3SVl27ygFoo7)S^1)ze zNccP+eI(<~sj8ZRor7kVU;T25>}S=HIwPt=Rc+J()03|rm|8x}oQNp654*y#=cNvw zHZ?=f&~bGY9fXfGx3rYd+)G>|vHRd5c*3vfOIO@6d7OfRVjB3<^Mk!vXmt!B+FP|H zsZ7DGI{sj9>$|nr5+iXj_HPindiSZj)cV3>-a*KcM4-h&Q@ap4iqvuKEO;m3aRU1B zf0mbvaHfg$#?gJ%)zxH3fdi;=nP=R=49IBTStY+v!<8&YP`Msd+k1{g`wkOqHDI;8 zbFKA3;$9SJ6ui^$2{iI`M06!QPlu>q2=Mp62w#76lcWNpB8;W4s5PEL=e3#5;g#If zT0K$a_ z4nN?chQgNK>5ps?zkWH5dXo6vF+2k%B47vk<|-r!gO}dBQ5MdCc%GB!Z8`@e_FW4y zUB<-6N-;cUp1iL{5D|g&_9A3_?fm1@{{C_(li|oHNtNAs4;tly=92M;$Y&x!d=YeB zw-uXJ3BK^7?nx{atN+h8TM2eu`&HsG9!AL|c>7rLhTLtp9FU$rqX~qq2&FrD)O8n! z+(~zLaUu*oRrm&n;Wh+;{c>__X+ay(;Mf03ADf&v9f7D3s1gAW1L%G?@o}#6ET<`2 zbJsk9`5|^CEc-)O6R;0)vYk^$R1Vd61ezQXTgd$HqiEVz2(sC)*Xk4a>BBO+!7Oc3?qnj|N8MpBU6E<4 z931lVuhNbfR0yYMWOSkrJ0v8Fm{S3UzDJ|ks%XUtP(}2vjX&d@L@t;JADXN6MRX&n z^Q-rYu_=PT1dj{Z#p!fOa%E+On~UqDw&VuvWS{2_?q`&ovE!^tYsz2NEPvfmWgU*w zXg^q9V7OHKaVpul#wGF4t;mP{Q;7))31xrno+Ujv?$<{)O)Nva72zVQMcI$#j>&L= z_Q5AZfXX2V$(=&)8OF$waJx4o6NKNa<7Y+ktBY`e70pG@KJI3bh{*$vslpH9XD2^Z zSCu({a1Y=xvNIR2SD3YxV2%twTeG zeG#WEEHXQ-FG2HhSn|VuFB!#hCg(*x`q=PE0_7neMnl&NfUs%NSDUXq!!9K+L6O)| zIYZ>`&E2n;Q@tQ|*sR=hM54ty-Z`zGj^k z5O3|qMw`p)f+S~#4~T#H3#OAwJL71eLYMn_4nl*+LMM1UbUYlVx+OyCfFH~Pg^bSE zxL!%eWiM0*|M||df1bGYf_r|2($ZB1=MKPRw(fQbj)yFshyL+vdHtTr6%_ zzM)9*va;MB;sX()LxO^y&b8eRW!=WY9bQ1!*QxRp04q@!gnn}o ziQL=g$;rw5W8rTM?fRcOdiq65g>2s| z$e$ee>+8*-az$k;c-r-_1>d7YI3$EcL?~XoNQV&YkB9GKRHjii zw7Q7GMAXWysT?TK6qyu|A&Hg6nJ`pW{k>%@gd(f1?>2Y-_A?xl5Ofz4=+wG`8Qqga z>|5ONL(jBx+a4evdZ>73PM#!zL+@`J1`cA`Vf5;-7#9?diokIvZCA{G5fdFuXAiHWWDyN;(X^*aP|u z;quwJgb2?I<$yFS3znm^vr)u=)=NPJI=a1k4<%nI(38;JN#|FxaBC)1p5OXot8hbW zAX7S@8GvMo573coQ-|B8RTDni6Cu&L{Yeic>_uL?{uu1 zaA>;CWo`+##M+JiV;K^$s$aJUo6KiQ@bU2tF5OPInP30-(x=;W+sCcDIG7i%+X-4> zCzwgT91vMZ#XGYDOl^(ccwGxM=4`)1>0*V}*X6b*>G$;IJr!O|wcnnU=JwjDC5dl? z#USybxMdm!1|f&Qd%QQb^4@>>=YG;d`1xH6`M@jiVH_%PZTG8{UV#reynyQ}eIi@T zPT3+QnU?cJ`2G9L^-ig2$Nj7p{8-EfFhT~6GPCBxwRye6iG>F+h7^KSMbEU)yPlnv z=j9K+ zg8jn6@26L33Revuy^=hQy-^Rg`W{XomVmcq6zz!^`g({HKpWof87=}lmKu4(67UL3 zBIpYh`?fSV0=eiEa}(@$ ztCDq4NK^SVl}s^J2Tc5j$yGL0FmfD{M8=ACRAaB3LpCh&8H==piY!{+^(|gCix|_) zoJ2ELE)F}CI=-8o-LcpDe&wBdfSd@sBV&~2G%ySZ`L=#mVt)CTIa%fJkk_yx{Gif+ z4njfy-g`9TKAVRc3-4yA7_AL!{l-Cw;nK3fT<`3daqd=`l#<*+n)dq9MneEdnz@;% zn$_$&m>4xYF)_i`tKISPY4{{i?omnQ76_!Sy*F1OB3QIB8!E!@A?7dK(M|_HXu$D0 zOD?f=8@hHDP6;q#*>nYR{X~K``}pzWP`XS#4MX}PP}|@A`lHCGsN3C(cRHfy$;vy( zulxMgt*r=Xk>1&%|WZJzjQA)-I4>b;-=^><>nkWUddz zIEJ*jksji3Yw!7<{ltoD%a+k6VeD7Zd+i9L1O%q_%OC;z%Oyk*RN>9LNleIi`}xUo zizh3dIm5-vO9dGVQ97hGJ|iB|s?}I%8r)e_9$EmG-9da)My1ua>?Z}XzjZ{X_4D9n zTAjTx=EH$$mlB`E`POTGLtu(GHQ$kJifi#o$N#m_>APZjo?L#Fl2Rl%rCLnM0` z!q0p)PHxQ!LTQJZPx|o=f@SkwNOf>;_o8%LvW*vAko~>1bZl^K7BWOh-zn&dXAyy5 z3_GVn=*C-jDIxSxzW#jZr4Iw>UPGoz3>)(UN5X$0^HB$P)(#^LP%6!#WoZs8$(8LK zU?xFcpmt-rx;860dXwjQBzBy8jQ2%vjFiOY)(n#>zIAtgeB0g40E|<-60(1eVKsll z*-jTU4^u4`RZg`B9i5)M$-(%f`@}`abmlm)ya_(z54?ok###+ZFZ8LZsA#wc_$^M~ zRR>CB-Wg~zfdx>KFvCTrW`y}E&}Guf7F4q=1JhY_OS8PMFXH8hPQ6+#9s=f{aist3 zoBzuKpf+NLI`qUd90o-T3ybm=i?y+xEmRte#BdNX%}Z>`U(Tf;lBxclq?8^sr#8EU$ zdHni(@p(#Zg45*jHsnIxdekK5L}=E6w5H+2?mu{t63AFPVFy&5=idW#iMwVrG!E>` zIJ@g(O@M5RXg1=O0~m>VYOiD$nwvbF5*rk9>@O%xxlPrVFJHPAs=6%0x;9_j7HrZJ zF5t${ozEQKjZ^%5;y#poogxzo(rWxGRH333D~YOA5=kkjMb3Cq^t~=k2ahU85P}r9 zLqlDg1N3m?4!1wv&t*C{_G1ghhys&g0h$RcnEM!L`Z_YAZq%}QvLHFOBE&;dtEoJN z=9`s{WFXiwQ@~?v9~A7EW`--y_ni8VZqbrZg=}keIfxi`BIW{okf0RxD?M3~t)17I zt-cRqvJK=%@6fSh9o*^{|BziE%mNM(iMJkwr?5y7-w#87lsPQZsPRJR5ybfBY(&jl zJuf}*q0Au0BBfc^Z7|C~bX8b3#Q0R0bH4h*IfS$}mrjXX-jNISRC{P(Fk^Tp25Lb_ z$WNAj%Lz$nh2Fj`Jb35Hagl~_!Hw4AH1w2jawH9a9ej((_`={z zBQIrTW-{|!Z9BeO{Icb%OSzhuiIFL(^srh~E8oMf^ZxtDmAqavlpL%!5TD1L8|R?H zW^Pp6p85FkBQYyTxPTylvFmNMwVQD_<=-C_=ZW+hTb{S$VeYfrK})Lyz%1~orn;G_ z*5Xz{TW1VF;q?|E@JB=hoUS9b<01^E(`KZ;U z*_N?Ui)*KA&v-f6bZTQJ8qoVQ5h?2&My0q8)_=G!;)WPYd;~Fs4TZQQU{|Ay(vl@k za>=r6n|eU%phOp@J~C(%uLlR`%lLe(^6Ho{>+0OboI3tsqj{Gak8PXbKt(#1H_6}i z7q+zjU7O9(j)xw_f0b&_;beD{&}T8{O)l~J-qSR_h(dEqJz0`M+1bvfZ5Bf(8IhT7 z67^_?pz{fB_U|83Ktl0>?$~cW`zk^tQ{a!{!3(ufJqv|8)$}i_! zq;MQ8N2C18I-Tgya*k$eWQ)qtO@(7lP0$&NgT4&%8(t)l#bW*PHhCK{fR+FDGR@B8Wzzr zziI>wI*CeurxU&+?i^tbZ3fekvjy9@(duq7npizdLT5kw?B4TK;b zjJMF(v+$l8hDcG`W?Sw+HtRWN3h6<>DI^0~|=Aa&CnxL~?}&movc_~EICFU$_V63VmHd}epKbk`>5j}@ku z_7PTT*I%F11f>;8eq(%#2M{zs?DxKc9ayX?|9HY(bTCZkV~tB2Td^PU;Xn*Brj3|T z$!9tN?j%Da58#JAw;2h;TXC=Id=q0FA#S#a-C^XxS-xi_ZCesJ`{_7&J&2g)NqTOz z7H?2czRz6EGh!Tm6A>YYrguJERXW>FX)6_#m7h?%ONxeNF#1GGpx9V*1I|!wTbsa@ zxCqa%d*eQ2J}Z{|+y@x7Z#utr`qZiX&^A^~r-K#Ha_;6W1wnczX!rX_?V$qQXC32G zv}mWq7dupW7LdnoM3(!^S&_N65dx!*&ixjiEi)=jNscV*JeK%MkyvuwC~TH+`Xe~g z?=;LZ5#_so%$WnD^8`S={75h~eOm0_3-^E7DbyuqCgQb!e<*jnfPM%HLX({ccHJ)M z$gN_>ZmX^!Phi=7`N<`b)zsFE35v1@raVB_8iC|ucXI~Q?7 z8g9H`TrnbI4bf0zb8`4{b7AaM%s*^QrD+L8itSro$-RnrxNY#%E+zr$AA?0`@FS$1 zDO->3{yAaENL<`c|pNc3YIFA`3j8A5`ZLAn^Ce~nN{*6N|N|5lmdcYJtr zavC&?OSZg+);U-gw^|JZEx~9M#+O03PxbA@W4Oetvjy|N&NY9ih&WPo;BhhNt)&?a zfiNG}}#1OJ54lE9 zz+(7$b!8~>Mqw*A2@U9Ewv2O_>QHZ3BZMMozi`)yK^SN;yOEtk%_hc^h~SlC@ffUF6v}&T z8A*#fBShsp(^luUg|scxGcXkVnEMV%8$qRPaYK9>zrLOYLA$4U z%Lup6f~8~bzI@OLO-)T7J+|Qoi$ikcINp9&R@YH^o;I+?VxwD(;Lavr6v67I*O2$X?%P zW*o~oKNB%2`W(1fF_*a!&B3thYXDJzRW;iU6x#PNk6`aIAuCTz$Q0Xl^8&U%ggVTW zmMk4v3r4anbzGy=O1-LI>6v~>+oi3;gzEZzTOWYO-INsm;$MBFnud zXMFHYG&yrFTBoxh<)BGyBC>qnH21#2BRK`{B5M*%>kG8Q4vAjFIIxMM$+CyqAEJg! zM2RZ7wtjbuh>B{2;FxKPuRVd2U;n5v6ai@5M3#lV(#*?zCmWQ&Lsb>jJO9n~w8eW6 z-$sNaMX~Fs>imzPeSml{`d3yVD>I(+9pk|lg(yHIw78I6T8d~g;uW^On0%V($_hIz z-;YVXX0vvuOX>~G_~Hh>FdrOjUS}sojXG?->!>MTvGY-#2*H;08MUP%aNRdjDY1mg zyw&F0S9mjBzf{>GN2Fj5eTIThx<&-#kQicGYMS`&oqFuWU5WnR-|%HevBznC&(Fs$ z1RbH?G}e)}i?obMu`IVabf;gKTll*CYXl=(_acr2E_}aSAL9tYPXURaOx90elGX{AI@r^`FxsWSOH~zzPqCA~VSIfKvELK-eN&B|5T7fF0D= zz3|yelH%(mgc~fzKLZeHwn;LmEj)1e@j_BSpe@Hrf9%6^wnkI5ZGujB=43~f1gv`E zD?-$gFn}$J;mA2TujA%AEz(G94bH}nNI{FcEZOnXbWZ5iS}-=yBW!X_FHjJ3*yV(j zqMffpLDD39Lr!H|@JJ!X4o)n8-kcj^0R8PJrta_0To`2XJh~f2;#oP>JyS_G;lEk; zKB_P(lC`+kjQy6_nd%o7DAuZ~uGyH0vzkYp$xxYQE#?X~Jo9AF^ibWvIkS2WhU>?R z>(6YFRXll-f52~+wJPUfgOQ2XXZ?*DlCVN_swF;a}(mNh5 zLVP6zB2XARv`#giwo)qWZ%MjTQBgsBuY^?`jrg88)wt}j?=TF)lN`?Jh8?QBkUNi( z4c_`LX&VOWkDfh{JsN^-j{})C*=lwkla5FYTFk+jjZP-qIu-N@-RKnMH1_PhWGLFBO)eKHaU+kryjvQg9` zG4Pa~{P0t;`BwN9JCLf=(bfGAbH72?55H<`Xkb3qo%!7V&u{n6Q@!)`Y z^~XL)P~uRNeun_T-n{l2@qF~g*+yDWa zG23T3r?n9eMdRvtJ3D5Z_Nzw(Wr2(6#-C)xa9TqYXWDEH!=%S-FXT0Td!@%tzb2-p zjt*Ku>Ap7KLQHon6lmi+c5F>b?RuZExJM-VYe`#1en;fOrvEAXyl+7Hb926opp|q! zZp?!WamJUHmMS=JH>bh*nsjPuZ>L3G%=UC=3grIU9dyQcPF9?(T0?nY?cTXpBQGJj(F4DydB^_F5BGPyqn!aJFoXG ztgj@PgWHG^*Q?K|F{1(kabMJOJ263uMt@yslUAcYl)(R2*_ns+ymjyYlae^0h!cuZ z(j?NHPDp7|6fzT)L?JTEkRg(!LLs3-DpF)t8A_8>rc4#$^cnVh@3q&u*L~kB$g8`ZB2?YzBSNS_IWy!*O}ciNNHZitIw*;f?dRlLLmyGw ze%dha2`ft-_`Uw7W^N9@>EQ_s34N(V!oq|ASR~hm56N+O5~4OFe%H|}iVrh=oc%nz zM!Zb=E19;vvzfZ!7eUGpY;wv^O!Q1|Pb~6~Eo6k7(I4`l7nTn(+=%0!7O`yD}LYW^EE$Aeas9w2(!z5 zm!a7H@axRz1nXA2FQ3#SUyJf<_UZIEm7VQc*(H&@sNxSi5*t|gkZ;;>Mb$1p1=?(G zY3Vdg3a>UN-Cj^EXN~DP=D-^tlho+T%fGaNoKyULseJw-9Ye#N^cZ*YVj!`ws|Oe&E#IkD!;P5s9FrJiTOr zAjZ`HCSu%})L4rMuM;OuxLXNE6m8#K)5~CoO(dmCdXqXXrtB7vi z#SK^T0r>?5Z{C?GcP$8XKYs;_j^t}Ip}4@H@K||r_rMvqK$N}oYlG!GcJ!ec z6my{N_91E5U#Ie#%HMzUXIsikMEjhbsM3Us-Xw%t*+buZ3%oTJf1NxiK!{V}9I>{O z0E?2|u)A-PqZvIupeE;!bjhDuy5bV~&7{d|adnW*nLYxhe6o(2$I0md6~V8^&6?F$ zI(Ps6{h8xWC^Ym(ln>kKL(^IM*;y7jgG}l4gTj9JM&6k)=>$oT`jr0W=CNq}5st-S z9VVSintmNYt$#eZzf&d34$ph{?y-nOqL(Y4{u~%-Cg6H_Ny zntI!!Dx#rRwG!DSoZD<}@~TyPnp>O<*20?{7@*3Ae_@>b?Q%Rj3?m522+ICy>bWo0%ds#Exy4qsSJvTH@`xLrQXe?3GeVoA2-cj#&Z<#P4MEL zj~1k&=y;38TKE4JO5*5#btO7Fp7RS$^Oc97Z_sm3+Va7BHTwy`WHvF(x^ZI*0^Scm zrm_nXUYg6UGYNPYBNSQxrgQl{f}Ux%d>@>7*RBnG!RsGBe)KU!M|S{5*@o=W1}h;O zJsmHsqEr&psZC>!JT_nWY@OqexUn(*#cKtlw&_!8oKb=Qt%0M4r`|1!`o=FR_J~sZ zR*E^_BHNb@ky(v#cBYK2cyV#RLzex4b*O zx|!tJyKVtV`$Ybi4Kug-kMatf~=P)2ftt=^9@+2P~?x>#ccKIYqSJe0@dnsG5KIh*EyX7F)kNcUg zU+)#q*~e9S@c6~Xjq&lSloX*oWoA_tTr z+js%_J^yV(h?~w)_Lkpk-b#+z)##ykttg~~RG;l9T7IW!T}xUW@XiH#3BO7LX7SQ* zFAg456Iz$`>(>hw*B#%^>x^4Y*4V-X5+;H|Qv-8GaJX=~8Ypqp@H7p0s5G?svoy_o z-x%v$bXIyTqg@=bgontb)s0;Uy*4C60?;$L)&e0*IFuo*Qlg3&e9m#5c9`&KB?RN1 zpkcDiq3ZcBV4aGVkKXy^<%+_<0q?8~eviBt?~$F?vTvlV?K7T+1QQS0*h@N3fa@zu zfQ$D!{&19&i0&^d5`DxTLp8cuIy3gf0XJ~F?T*+-RJ!Sr6XfCH5zC1Y#1?!~P;E|` z8&yf&Ne14x(rZ3WN5*fR_??ifY(G^(Tt& z&xnS6Q=89bWr@mkE3EV9?2?LME)*13GiF|U_Wu1ozqm^jE-B3)c{qjkMa7*l-%U)^ z7^mA!7`zCaU{b_WsAu!6R(+S}9*mZKJlOW#??xQ|_Ma@?zTY!Ex%SaKYH3lEV-mLh z?}5{Ww5h^k9-ksqW;hGF|7u3#8O9-}En+)42c|8S(Yed2gOAQx)|LDjT;Dyqh|2QC zj-kJA?rB+lBxcd|Hx$CHb_|wf!u!<+HcP@)S65fjwBS05u3e!0>U>IQsC~qVJwjzP zWy%z;srKBgn`plY_S(bh`SY`uU&M(L2~2JR~joY3m6@3H>z=U(Km~T;zd)A z91nAiNxVC^3VcvFB|ISd#CfEeTr68bd6)U&Ck>4fgZxWLm9xxBNNv7`Nxm?OZdk?xiTsNAEb>NMmUitEgPrjgIC4xcVH_C6u?)80NL)pa`pDsl+HyMlC*+oBqX)J z!UY3HWkCWaTfbXYEp7T9u#_&zcf9H_msq5d1$h<{qPst2iIPD7s$CcSL`h$gy?NLg z8X*lGvaO$%*6Vjo-nVEgLN==Sae@R!w$D)gUC^yT3%vOrC-bVk zND+WY9?&6wfPKaTNOxIF!jY6Ef1OU||D>img+1%HR&<1ZE`bmJsXC_&opcc{q%%MN z<`OEPZ(Kyl@_;s1=C>zZWYo)wtDDtcElJi6^XuKCbo3g7%kwJI_77R+y7aNY^|wxs ze^v4Rn|r{SngZet3DT&em8ZB(vk387kT6<3#^CQHHCZ4;?QOmKCeT+tVg3d1z*vR7QeKgY7FL~QEC_CoW z0+IIpuseBwUs`iN3jwdzUc8GbEHQIP@`%AhD7~F3y$238-qzi2XImMWo73+T`|6kc zkailf{PV%D)6Y`i>GoP;sR~Im?6+#ws^)bbzT@Hj z@md+R6w51zTpV0mh&d%6g0^KYm&zCs5?A~KAO?O_;5#jI2 zK2cxFM1Ebf?85`aVHIDBsDe?!82++OPfLSMINEO|+2hFS&{N6!YsnCh{z~U7;}04T zRllbOGd1_Dk!pUY>td#jh}6WvCeyE7yU^8v{-~w2iDCm>M#HHjaZ8?kUF)jy2w;2T zkI$b=@h+G~14N2@(yB~PdFJc4Q-*I@&$bXyL=3NqJ9c2BvD zX8?fnWYe`|z!iUP-MTefNlzVlvGjvG^C=9NP)zH(^ z8^v!3{`J-Hk#g31fC}N#1BZJP58xVf0T6WjANq>h0?JJy|D~dy&%Ys%G4cuxV`HUo z!?uZKJ!m!{!mZuxMD}nBhfNsi_RnwDpXI`I`m<`~%BxiJ7*>FilfF;gVB7aUTma+3 zrjvA|F^)60u#j?iR}27`SXVG{`GGyVa)_&#t#X|SW};)!j3^^CLgoNIl7nRT;@sWn z(>Z**;bty^;@5XkF2B$dA|vPEyVAMzVs>zBCqF(D(oa!l_yzUo%$zbJq$m8#3d0-z zI+IUEDExITP-XFwCA*Gi-a9^i8J+NM=Y*vO)mBv0nPOE)a1G;_l)h~W-?3wK zWl0TC#mMR%hLAJ9l%8ilu)KYwqej|y>d?XIhY5Adx+NbT7?%2u!{q}q|IuMJ7X9Ty z(pK9ZzoN9Y7Q2t}E6aY=RhOYvfH|4Q2l<0RY7OZP4)R`aKAcIdsvU+&`fH>2rC|If zXeXKVdz4+l9kr8Q2g#fDM?Vq_eLk~(DJdNRQ4cuzUa1Z7^siAv6P#YGn7*>?A1kY* z#Tg>GruCl*OsuOq>xg#MNC#7VU_gMiocogK+ty~tYlZdQ77;O9RB8}I41axR0~Y%S z1o*UjK*6sZJoJbN89J}Y6sJu6W1G^nuefl`_v@z-OHohnwoarz(^ZYy2(YERw15Bp zT(woC!OD_s@b2j5-V_CzdU`_2Y*@I$IUDO4*>&A?+&9u`*01!TDtk{eWbejDk?JG# zU)#VHhhOS3saek>E#qcVEu5=~u)gfRtvEmPnrv-w4tWor3?+qyTUqAvnw%d<`DZU! zP$tQsy4wHbT;R2&reYl;zV(9m{riP66rI|NWoOi?4-W26nHwqV5x8W@nb^UDPV%gV zc{~MGe{i_z#vP0v#L;+o zQ_k<%E{wc-u@v+!B2~?=Ub{A+-nKh3ytPMT#@?Nn5}1JHPz#_)_sReTa^Msn$NXo?r?%X%zNg0 z9d9iu&j_~Gvfpt#Yu;ykEgI?UuykzZ_DFraKp~oVx7*HL{OE&GZLOl+PGRjR?v1LY zbs}XeLrOS*_KzMY^?2nZQaK@0f^lwA)jm!;EK{+AR0go~KJOVde}aF0ULRav*gx~f zll%84=r&Yn5O6;je_Q!_&DGzuekPT0Q<*fy8ohY)ra+QG0UU$CTTeA&I1moNxetaf zg$+9cK$_@bMn2Iuu`hg8M`2b)g9w+rGX3ut_x4H0j_p^c()||E*8AkQ9MsXkR4u8j_ER%hur!w;%^i`e9 z!n{pewUc83#b{F5&FbnIwF{NFemQWlX=#Igs6TGlKukM2#d(05np)*IYR%oV<7uB5 zrhDwO-`UF*_<`^2{DtUszwf7wNBbKXxLz1T8Mf!ffl%#NgmHfg40%b`*|TTQdG~X> zMkZz^VNiQyPqL~*6W7PBl!uPIDM{^I( z43I0lPCOZd9LaAb%`at>xiNMTmWkNpvl7w{J^x7|o78%=hbjar#sC-c-oGh}jH`Z! z3ul^&KS016VZG8?RKQ*dVVF?hb{N^eAW#K7Z9B3gCyh4sRy)`@j>Vcgaz~b6uqx!^ z#)K2Ewux>}xi`O{eD=FQwjz8UyM-U~JV#W#(10aaTzVAFQnP^!4rU0`}TV1eLnyn~YEoJ%9Of0})vw zdHm#wTc)>T^v62Sj40aH3|kNYr5Dsa(JU2hL7Q~b9WUG&eL;zi zVBUt@fnO<6!@uU&YEgJY9ne0`@3m0){AXiBMB$_v_1il>k;K2iJmAL224}{{y06hMeWX_#q zij%YbQx=Zfor13h4a$(LETioCPr|@8ZItYiul1b`NEzN|Xm2`IBe^gGQBMoy25Wf_ zoPKsm?0U2c&*;L)$;+QUS5rf}LwVoMZ>0cH?xpXmn);XW3#fN=F7wU&W@L?1N55g|FO8)z@kS72_d!Z zFC+R=WYLg%FqE)6aC34`^1yQ!FBZQ>>G%v;SYgm*iu(r#8=Y{K!G}lq$4CTcN3MPQ z;g~QfMU+s$6?O@b+`baL>cPW@y2`tJ^!_blgER?Oqb;;Zr%>CH2jumqmVcRrQc^h4 zbl+p_9(Lx8iJ%0GOP_q#|JI=_XOO+Qi$+5e0|44^BR`h zOEI;iZ+*+_&Ht)OG0_;5?lEt7OUVuF!0P;f4<^wcF)rGN)}7TJYiS#`9*Z0Igq442 zVmnV-D~zHN6@S!6(#}CXA&_?cL4yWOscYd};_WjnnVW#Y@s4y3Wjuz3cDqFX{#*-ElYd_no#FGQ_#%V+$yh zwtJEs1mr=4f1ydCR2y$@6V)4bZ_lkDpN@O?CbTW3{MqLo>S`{)RF*9^*}H4dEaixox7qPnn8*xs4`?EEP8nB9 zLrmyN>FCQ-jP05{r=Rn*)cDxm@l{>whlU%HkIW)5uDC=I31!qS?^xz&1(^@kIen-l z9zgFE9iLB&X!oq>vV%|$Oe&_Jq$I2ZC!PhFf^Oc5E+f_?mCpe`w1kY&4x6`;@KoIn z*|ZfVj(eZQN%0+!D1Uf(<}a|c5UvwoIbhvwTrW}*V9C!_z6m-iDT|m1!KIkw2NuYY zh`*9EoBWC|GqiLywg6dv9^GL0uz?&aK|Dja1?^N8@Hb=~SkPnd_gTK3(p*%`Ew0Qv zxe0+0_;MkiEE&fBRP5Hx$I#IOSq5{R2qEgs{Ft0>S~|FTb%V?Lobpr|h91=qwU1C4 zPB$W{7L+#Y@`e)S4Qg^ffvOfl+fjiw#f^C@Gw~f=p4{?s!ySi5_5604{-v-im1y=G z`0~YzF;`YMT5VGwnwgmy5+5J$FX=q1H$b<+cA0s;zU>$Sg0YQ2uPcsiC%lZJDGc9?>lA}H4lAYh?b z7L0Q0x0^x#pqR`Gc4?fH_s^sLnwx|lYnvVt@VRA>gB%>w!*NVJGr1#OT4OuU@WM~W z{tA&08~x(L2jwS;LKu`x5PV;~@mdM@RR`avY=WBjl!#f~woZuQpyr-qV8O-;K41w}e8u#gM+BI1gZ(~g zQM2qiMiW}Smkm*{=@0^9c`Yw*_+Ha|6VH|RS52BG{VvGz(%L`0MQ^NkT)!rHj+-9G zJ#|>ZTHQ@i0ZGGf))9hYgE)KA=+Q0dLR$~;$ z4gcs<*e=tVr6E*Ki#PnsUX=nyn=5V)T(ysGC~3+did< zDg}8@o?ITL`B36X(>!e2+X=_VIo@0L2teogyLV21chaoJ%j+kX0YukzNf+a2bhNcD zl-2BK@7^g~2vP;wl()rUHV}r8*Kf@x;}(;I!uCLr$gQr<5pRBw6Bw5BhDkp&(58$7 zX|EkrfBk5}nlSESJ=Kxx?R@_v4I@v3)+M0cU7h-kzC$ux-svXRS;8idE<+m$B73!q z>ohcn@|`lWvp17HaCd2!e;*3V$s7Wy?m8!k|J6tj@@dLSC9Dk0|WQzMRl}ZdzpdkkK8I$o2&@E_1kI zrXIf2gwQ8Bg=Z!WhK^ksOlUhc?znNo;K`=ntT>0y>&4$*8l6PZ+lCI@4O&N15`=wo zs;^vr5vXW0YR8)uB?&`k?Cm|+TMl|=x~HwCMxpA+{%fa0gUcf=W=P<&EC&^%{0tBs z71bP|S<-~$9GxB2R%O2D*A#gLy?Yr@zj(9C*{sDMPEyWJ(^jw}ioe6_q4`O?1V)p*@*0?h1Q@F^sW!G3LbcUfkf*>kJ0S>`*1hbpf3 zx_0d<^2d~=#Ym6`p%>ePJALd!@qZZDC>e%8TR~Cr(DH;Cn*6ShBdD>2dXHz+nZRar z@%0sA(#~z6ALN5#y6#@sti>EKo2~PPMZ6s4k*AjB<*>*mLVE+V)grCB_3qsX?>EUX zIvc_}R6OVB$K=$HpquO_1aP85hai44hN>7TIMil1SYi#i6<@&TWY=Qzm0IU9|{`xeLipty9f zx2tIQT08dKg$qtwvq{%y$Ig+K(2o*d0t~D0ns)05KYeL2IdO%U=0?Yiw4|Lxrd>NP zPCf4+XSfLmJ`#&I5+P#dmT@1TZIpC(t{Wai5+;vG39 zw)7S4!k#8`@jVpdB*CxjDOzJXZBv+~2}&z83{;FlM}VWDlSY@W4PPJ=tZx2NAd!$E z@S!gfcUbf8L_B9{C1h(o-)?Rln%7N?zjU+ZD<;3T7wbup{`vdIDdMRLgD`aINyf@W zBrrPVvo|x2Y8}y4+QI2`0OXfq+g_1(&i0b9_**Qew@!C+Y$ui=jPz&WqD5D^@uVbF z4Pu@rbC|B~RJ2(H;HE$Mv>esHM8X$4OxeNexW6=SHLc)>wnqqpzD(lKvv%ys`@TOh zyH%m3#FW~e#$&w4jcI*`i|tDt+84^^D=T}=4tCqQKWXqA`6NWC(!zhT|-gVxhFv~*O z>+_Ev_JA}B%F0{FRCDj$^Wbpka+hw#k=-qOJo=(U|AGqkx+q1d+e9p-f!r}=a~n=b zl<1|r#3=R}7rk)P#MWh0ZRD1%6Lpvc`+7bnIU_UEhO?b@=8|@pXL#IcQ~RRHGYU&5 zz!4IoQa?OC-<1Q_{XV>n^vxt7%oEae>1v}N3U}GlQ%%_`BlB?OUyGN}D@_a#UehBV z?d~@-{5d^*3otZApNF%(>s38R z0B}l-rp4jQFX}H>eB;=|qLl&?Isp8j9fW;2i-*W>an~@D8J_WA+q4Uw>em;!MY|Q_ zinh`JNu=X>rU26SD4f5KAZ-oGea(kwI=StMq9*JQphSDa=hI!!S2veyac1SWN5E z{9vzMOv-C3$QERxeVG}EsBdJo8IPcNc1h)xkOQ8S-MYJu=jL=K9>_;q{p_teqTcMd zjg5n>vKNeRJ2uSoSFd)Zzl!b}Fw8y6b!;p{D!=9c65w`W;Iax0n$RI;2lVPVh^fb# z+iSY6^0Mk+d1kUpoiv}b0YgyB%}Nr1D@&84F%wfvjs(Qt{i*Rp(P z94N3-kIrX|+(soCX> z+~g4&S|b6qDgqG`ivj^HV|TL)2iW46;e72WZkMRVOgV#&W{Zb7gw_m>Dg)Pqk|@75 zU)e!)=|N9bO(ZKUb6<0FW%iN{u8dr~9V1@GblMTCKooVfXIUA=u%7SWa~{(ISz?5hz4&p)jL~Q(STbS+EdpOLdg1G1c51iq za6aCm>GN}G#=vYKH5L);^Y`!hG*u)LlG-EhPp^*P`b|1qS-EDtrGuwqCY-yQ|xl`B`4eMk;oan3H-QSvF$o)Tu%i>tB6E^|1D9 z%ib(w(!oYda0e}UH%{PP!zwzSuW?O(bY@;C+ex=aZ#JP3Gi?PV$RiTJLdlLV;8+Gy zUs1E+ouW_q@~n+XpV0ye17uu$;yQ(YhGrtL8?IwJI85qtZCC}K4aMEmf52& z6wZ#4o5y9$&??Y2-xAI7u{drhNE|l2n>1yLUiV&xs+M|cHEVG{gHIO< zYNBj>uyyt(^P&3_Wh=ZM^mnDvbQ zOx5I2L#>OII0+x{Nk6y$7gbr3_dYTMK-5JBeXdF*?KU~u4j&i&$`Um9Dgj)~ZU$N_ z(WHKoQ*WHMaAC&R92hh-Fp?BYP;iHH+_jZec`S?RKmw`BO2iUS(8;R#*5;CBJqWEH19hxc)9pqJ)jX`vf?ysopni?`=_65ht}zo=x6ubdKx2*s<|gJf&U)o ztnI`4CLc=X(MPIyhCdgXd{F)QtDu5ZMeBK~oBe%m@`_#cJ=1_L338xL{4;L_Du3T}ZPM>Y4!=im+qpaR9wWx$vX2Ubs$#~6Nae!vAW^a^69a;IgnnJ-(z3P%9~50P zeDel9`q(hmY(lAPJ+x<8*Z|y*e~)(&CBh;H3~B{;t@>!8*?qBhnmLa_hiut9BS6Ct z#Mppj8T#H*5gUZFoll=DW%S2?T}u@xBV+U?bRue(EjqZd=AjSld_$-_5>1ogmk}k4 zr>U3vr_X1z^G5{LT4c4TQ`&%dOvs*l4p$FF>}jPW#`7O5qn^0AD4vUxE3P*mT(oEs zJ7=SJ8X-r>3dIZ%kTZ2PH3@t8ux&3M>nOrkGw$j5Tz4XW+eUr(9DkxsgzY86BmZ?p zIqlXL7pK_CR9rm_3Y15HH%V!l=wO{F53^82Os;h?(9Li6=xI|pt1eYN8;yOsx%=9VNIb#R;0D;n;Hg^Yod!t+R?ENKl%767Tz+(Wlg^J7cKFON63i;4ZQ5Z} zT>FhbzIYMkvo-*{n>S$k_p)qRcP2K9^&xTw9*Z>$b=K5OAR`gll=h%6zf0YS2vWJEd)eGp(C6xvuW^zWn>J(AwtSwuI{@xZ$bg|#I#iN z;9tg)XRIXS1d>lrC*w78IoLjA9WSQY_LAO}>D-e4yt!OqPx8dVbxz44Q1FD01XovL zz@C@pRQ`QFfBSX=bVOE7$I+C-t*}G z5&vF=joM-Zd0)MM|8Md}du$HV{Cr*cCQyxgh?x@@EGOsvMW?R+fUzWS#`EsUiyRj# zdp@z;6weqgvcdTx>;cI^RMEJ9>VIpzP-IG7MVE=WHB|SXc}`Oo)Ul1k_;Z1EQ(!qZ zY{O;?vtKF6XW$ny#@w>9&NN+UDdrKOGSgDN(`O^MLA={Zds3$N2`}>Zf3N)War&tT zs@v;4qcjnYeX}61h|KTAS(IAibjRKt`C<;+$TPt$JR@|sxSmj=4>P)@x)7D{(U?4m!l!!X2o&L3xNOsi3nL;eqq7+7uq literal 0 HcmV?d00001 From b629ae0a79ab5f6e8c5edfcee0caa1a1c693704d Mon Sep 17 00:00:00 2001 From: mateuszkupper Date: Thu, 16 Oct 2025 16:34:50 +0100 Subject: [PATCH 08/19] perceval backend refactor --- optyx/core/backends.py | 193 ++++++++++++++++++++++++++++++----------- 1 file changed, 140 insertions(+), 53 deletions(-) diff --git a/optyx/core/backends.py b/optyx/core/backends.py index 82f323c..45f46d4 100644 --- a/optyx/core/backends.py +++ b/optyx/core/backends.py @@ -611,7 +611,6 @@ def eval( diagram, diagram, task, tensor_diagram, extra ) - return EvalResult( discopy_tensor.Box( "Result", @@ -626,24 +625,40 @@ def eval( def _get_state_from_creations( self, matrix, - diagram + is_dom_closed, + has_diagram_creations, ): - if len(matrix.creations) == 0: - warnings.warn( - "The diagram does not include any photon creations. " + - "The default perceval_state will be " + - "the bosonic product state.", - UserWarning + if not has_diagram_creations and not is_dom_closed: + raise ValueError( + "No external 'perceval_state' provided and " + + "the diagram does not contain photon creations. " + + "The input state is not fully specified." ) - return pcvl.BasicState([1] * len(diagram.dom)) - if len(matrix.creations) != len(diagram.dom): + if has_diagram_creations and not is_dom_closed: raise ValueError( "The number of photon creations in the diagram " + - "does not match the number of input modes." + "does not match the total number of " + + "modes (partially open inputs?)." ) return pcvl.BasicState(matrix.creations) + def _get_effect_from_selections( + self, + matrix, + is_cod_closed, + has_diagram_selections, + ): + if not has_diagram_selections: + return None + if has_diagram_selections and not is_cod_closed: + raise ValueError( + "The number of photon selections in the diagram " + + "does not match the total number of " + + "modes (partially open outputs?)." + ) + return pcvl.BasicState(matrix.selections) + def _post_select_vacuum( self, dist, @@ -660,40 +675,54 @@ def _post_select_vacuum( if all(x == 0 for x in k[m_orig:]) } - def _get_effect(self, extra): - if "out" not in extra: - raise ValueError( - "The 'out' argument must be provided for " + - "task 'single_amp' or 'single_prob'." - ) - perceval_effect = extra["out"] - if not isinstance(perceval_effect, pcvl.BasicState): + # def _get_effect(self, extra): + # if "out" not in extra: + # raise ValueError( + # "The 'out' argument must be provided for " + + # "task 'single_amp' or 'single_prob'." + # ) + # perceval_effect = extra["out"] + # if not isinstance(perceval_effect, pcvl.BasicState): + # try: + # perceval_effect = pcvl.BasicState(list(perceval_effect)) + # except Exception as e: + # raise TypeError( + # "perceval_effect must be a perceval.BasicState" + # " or a sequence of non-negative " + + # "integers (occupation numbers)." + # ) from e + # return perceval_effect + + def _process_state(self, perceval_state): + if not isinstance(perceval_state, pcvl.BasicState): try: - perceval_effect = pcvl.BasicState(list(perceval_effect)) + perceval_state = pcvl.BasicState(list(perceval_state)) except Exception as e: raise TypeError( - "perceval_effect must be a perceval.BasicState" + "perceval_state must be a perceval.BasicState" " or a sequence of non-negative " + "integers (occupation numbers)." ) from e - return perceval_effect + return perceval_state - def _process_state(self, perceval_state): - if not isinstance(perceval_state, pcvl.BasicState): + def _process_effect(self, perceval_effect): + if perceval_effect is None: + return None + if not isinstance(perceval_effect, pcvl.BasicState): try: - perceval_state = pcvl.BasicState(list(perceval_state)) + perceval_effect = pcvl.BasicState(list(perceval_effect)) except Exception as e: raise TypeError( - "perceval_state must be a perceval.BasicState" + "perceval_effect must be a perceval.BasicState" " or a sequence of non-negative " + "integers (occupation numbers)." ) from e - return perceval_state + return perceval_effect def _dilate( self, matrix, - task, + single_output_task, perceval_effect, perceval_state ): @@ -705,22 +734,32 @@ def _dilate( UserWarning, stacklevel=2 ) - input_state_len = matrix.dom + input_state_len = len(perceval_state) + + if len(matrix.creations) == input_state_len: + pad_zeros = len(matrix.creations) - input_state_len + else: + pad_zeros = len(matrix.creations) + matrix = matrix.dilate() - if len(perceval_state) < matrix.dom: - perceval_state = pcvl.BasicState( - list(perceval_state) + - [0] * (matrix.dom - input_state_len) + + perceval_state = perceval_state * pcvl.BasicState( + [0] * pad_zeros + ) + if single_output_task: + perceval_effect = perceval_state * pcvl.BasicState( + [0] * pad_zeros ) - if task in ["single_amp", "single_prob"]: - if len(perceval_effect) < matrix.cod: - perceval_effect = pcvl.BasicState( - list(perceval_effect) + - [0] * (matrix.cod - len(perceval_effect)) - ) return matrix, perceval_effect, perceval_state - def _process_term(self, term, diagram, task, tensor_diagram, extra): + def _process_term( + self, + term, + diagram, + task, + tensor_diagram, + extra + ): """ Process a term in a sum of diagrams. @@ -729,19 +768,66 @@ def _process_term(self, term, diagram, task, tensor_diagram, extra): """ matrix = self._get_matrix(term) - perceval_state = extra.get( - "perceval_state", - self._get_state_from_creations(matrix, diagram) - ) - perceval_state = self._process_state(perceval_state) - perceval_effect = None - if task in ["single_amp", "single_prob"]: - perceval_effect = self._get_effect(extra) + state_provided = "perceval_state" in extra + effect_provided = "perceval_effect" in extra + has_diagram_creations = len(matrix.creations) != 0 + has_diagram_selections = len(matrix.selections) != 0 + is_dom_closed = len(diagram.dom) == 0 + is_cod_closed = len(diagram.cod) == 0 + single_output_task = task in ("single_amp", "single_prob") + + if state_provided and has_diagram_creations: + raise ValueError( + "External 'perceval_state' provided but the diagram contains " + "photon creations. Remove creations from the diagram or " + + "drop the external state." + ) + if effect_provided and has_diagram_selections: + raise ValueError( + "External 'perceval_effect' provided but the diagram contains " + "photon selections. Remove selections from the diagram or " + + "drop the external effect." + ) + + if state_provided: + perceval_state = self._process_state(extra["perceval_state"]) + else: + perceval_state = self._process_state( + self._get_state_from_creations( + matrix, + is_dom_closed, + has_diagram_creations + ) + ) + + if effect_provided: + perceval_effect = self._process_effect(extra["perceval_effect"]) + else: + perceval_effect = self._process_effect( + self._get_effect_from_selections( + matrix, + is_cod_closed, + has_diagram_selections + ) + ) + + if perceval_effect is not None: + if task == "amps": + task = "single_amp" + if task == "probs": + task = "single_prob" + + if single_output_task: + if perceval_effect is None: + raise ValueError( + "The 'perceval_effect' argument must be provided for " + + "task 'single_amp' or 'single_prob'." + ) # pylint: disable=protected-access if not matrix._umatrix_is_unitary(): matrix, perceval_effect, perceval_state = self._dilate( - matrix, task, perceval_effect, perceval_state + matrix, single_output_task, perceval_effect, perceval_state ) sim = pcvl.Simulator(self.perceval_backend) @@ -752,7 +838,7 @@ def _process_term(self, term, diagram, task, tensor_diagram, extra): k_extra = matrix.dom - m_orig result = None p = None - if task in ("probs", "amps"): + if not single_output_task: if task == "probs": result = sim.probs(perceval_state) result = {tuple(k): float(v) for k, v in result.items()} @@ -767,7 +853,7 @@ def _process_term(self, term, diagram, task, tensor_diagram, extra): output_types = diagram.cod output_tensor_cod = tensor_diagram.cod - elif task in ("single_prob", "single_amp"): + elif single_output_task: if task == "single_prob": p = float(sim.probability(perceval_state, perceval_effect)) return_type = StateType.SINGLE_PROB @@ -793,7 +879,7 @@ def _process_term(self, term, diagram, task, tensor_diagram, extra): dtype=float if task in ("single_prob", "probs") else complex ) - if result and task in ["probs", "amps"]: + if result and not single_output_task: configs = np.fromiter( (i for key in result for i in key), dtype=int, @@ -808,12 +894,13 @@ def _process_term(self, term, diagram, task, tensor_diagram, extra): array[tuple(configs.T)] = coeffs - elif task in ["single_prob", "single_amp"]: + elif single_output_task: array[0] = p else: pass return array, output_types, output_tensor_cod, return_type + class PermanentBackend(AbstractBackend): """ Backend implementation using optyx' Path module to compute matrix From cf9526b1f61d0a2c7506d6194e7324a52f572c66 Mon Sep 17 00:00:00 2001 From: mateuszkupper Date: Thu, 16 Oct 2025 16:35:24 +0100 Subject: [PATCH 09/19] Remove unused code in backends --- optyx/core/backends.py | 18 ------------------ 1 file changed, 18 deletions(-) diff --git a/optyx/core/backends.py b/optyx/core/backends.py index 45f46d4..811756e 100644 --- a/optyx/core/backends.py +++ b/optyx/core/backends.py @@ -675,24 +675,6 @@ def _post_select_vacuum( if all(x == 0 for x in k[m_orig:]) } - # def _get_effect(self, extra): - # if "out" not in extra: - # raise ValueError( - # "The 'out' argument must be provided for " + - # "task 'single_amp' or 'single_prob'." - # ) - # perceval_effect = extra["out"] - # if not isinstance(perceval_effect, pcvl.BasicState): - # try: - # perceval_effect = pcvl.BasicState(list(perceval_effect)) - # except Exception as e: - # raise TypeError( - # "perceval_effect must be a perceval.BasicState" - # " or a sequence of non-negative " + - # "integers (occupation numbers)." - # ) from e - # return perceval_effect - def _process_state(self, perceval_state): if not isinstance(perceval_state, pcvl.BasicState): try: From 1bee1748c2c9ddeed6063e61a327fc32db746ce4 Mon Sep 17 00:00:00 2001 From: mateuszkupper Date: Fri, 17 Oct 2025 21:44:42 +0100 Subject: [PATCH 10/19] Fix an issue with backends states and effects --- optyx/core/backends.py | 131 ++++++++++++++++++----------------------- 1 file changed, 57 insertions(+), 74 deletions(-) diff --git a/optyx/core/backends.py b/optyx/core/backends.py index 811756e..7ea718a 100644 --- a/optyx/core/backends.py +++ b/optyx/core/backends.py @@ -59,7 +59,7 @@ >>> from optyx.core.backends import PercevalBackend >>> backend = PercevalBackend() ->>> result = BS.eval(backend) +>>> result = (Create(1, 1) >> BS).eval(backend) >>> np.round(result.single_prob((2, 0)), 1) 0.5 """ @@ -624,40 +624,23 @@ def eval( def _get_state_from_creations( self, - matrix, - is_dom_closed, - has_diagram_creations, + creations, + external_perceval_state ): - if not has_diagram_creations and not is_dom_closed: - raise ValueError( - "No external 'perceval_state' provided and " + - "the diagram does not contain photon creations. " + - "The input state is not fully specified." - ) - - if has_diagram_creations and not is_dom_closed: - raise ValueError( - "The number of photon creations in the diagram " + - "does not match the total number of " + - "modes (partially open inputs?)." - ) - return pcvl.BasicState(matrix.creations) + return ( + external_perceval_state * + pcvl.BasicState(creations) + ) def _get_effect_from_selections( self, - matrix, - is_cod_closed, - has_diagram_selections, + selections, + external_perceval_effect ): - if not has_diagram_selections: - return None - if has_diagram_selections and not is_cod_closed: - raise ValueError( - "The number of photon selections in the diagram " + - "does not match the total number of " + - "modes (partially open outputs?)." - ) - return pcvl.BasicState(matrix.selections) + return ( + external_perceval_effect * + pcvl.BasicState(selections) + ) def _post_select_vacuum( self, @@ -704,8 +687,6 @@ def _process_effect(self, perceval_effect): def _dilate( self, matrix, - single_output_task, - perceval_effect, perceval_state ): warnings.warn( @@ -716,23 +697,15 @@ def _dilate( UserWarning, stacklevel=2 ) - input_state_len = len(perceval_state) - - if len(matrix.creations) == input_state_len: - pad_zeros = len(matrix.creations) - input_state_len - else: - pad_zeros = len(matrix.creations) - + current_n_create = len(matrix.creations) matrix = matrix.dilate() + pad_zeros = len(matrix.creations) - current_n_create perceval_state = perceval_state * pcvl.BasicState( [0] * pad_zeros ) - if single_output_task: - perceval_effect = perceval_state * pcvl.BasicState( - [0] * pad_zeros - ) - return matrix, perceval_effect, perceval_state + + return matrix, perceval_state def _process_term( self, @@ -752,47 +725,47 @@ def _process_term( state_provided = "perceval_state" in extra effect_provided = "perceval_effect" in extra - has_diagram_creations = len(matrix.creations) != 0 - has_diagram_selections = len(matrix.selections) != 0 is_dom_closed = len(diagram.dom) == 0 - is_cod_closed = len(diagram.cod) == 0 single_output_task = task in ("single_amp", "single_prob") - if state_provided and has_diagram_creations: - raise ValueError( - "External 'perceval_state' provided but the diagram contains " - "photon creations. Remove creations from the diagram or " + - "drop the external state." - ) - if effect_provided and has_diagram_selections: + if not state_provided and not is_dom_closed: raise ValueError( - "External 'perceval_effect' provided but the diagram contains " - "photon selections. Remove selections from the diagram or " + - "drop the external effect." + "External 'perceval_state' not provided but the diagram " + "has open input modes. Provide a 'perceval_state' " + + "or close all input modes with a state." ) + external_perceval_state = pcvl.BasicState([]) if state_provided: - perceval_state = self._process_state(extra["perceval_state"]) - else: - perceval_state = self._process_state( - self._get_state_from_creations( - matrix, - is_dom_closed, - has_diagram_creations + external_perceval_state = self._process_state( + extra["perceval_state"] + ) + + if external_perceval_state.m != matrix.dom: + raise ValueError( + "The provided 'perceval_state' does not match " + "the number of input modes of the diagram." ) + + perceval_state = self._process_state( + self._get_state_from_creations( + matrix.creations, + external_perceval_state ) + ) + perceval_effect = None if effect_provided: - perceval_effect = self._process_effect(extra["perceval_effect"]) - else: - perceval_effect = self._process_effect( - self._get_effect_from_selections( - matrix, - is_cod_closed, - has_diagram_selections - ) + external_perceval_effect = self._process_effect( + extra["perceval_effect"] ) + if external_perceval_effect.m != matrix.cod: + raise ValueError( + "The provided 'perceval_effect' does not match " + "the number of output modes of the diagram." + ) + if perceval_effect is not None: if task == "amps": task = "single_amp" @@ -808,11 +781,20 @@ def _process_term( # pylint: disable=protected-access if not matrix._umatrix_is_unitary(): - matrix, perceval_effect, perceval_state = self._dilate( - matrix, single_output_task, perceval_effect, perceval_state + matrix, perceval_state = self._dilate( + matrix, perceval_state ) + selections = matrix.selections + sim = pcvl.Simulator(self.perceval_backend) + postselect_conditions = [ + f"{str([i + matrix.cod])} == {s}" + for i, s in enumerate(selections) + ] + sim.set_postselection( + pcvl.PostSelect(str.join(" & ", postselect_conditions)) + ) perceval_circuit = self._umatrix_to_perceval_circuit(matrix.array) sim.set_circuit(perceval_circuit) @@ -820,6 +802,7 @@ def _process_term( k_extra = matrix.dom - m_orig result = None p = None + if not single_output_task: if task == "probs": result = sim.probs(perceval_state) From 4af9dcb98c0f88fd33ac568c42c7063125aed7e9 Mon Sep 17 00:00:00 2001 From: mateuszkupper Date: Fri, 17 Oct 2025 22:22:18 +0100 Subject: [PATCH 11/19] Added tests for perceval backend --- optyx/core/backends.py | 16 ++- test/test_backends.py | 280 ++++++++++++++++++++++++++++++++++++++++- 2 files changed, 291 insertions(+), 5 deletions(-) diff --git a/optyx/core/backends.py b/optyx/core/backends.py index 7ea718a..8eef2cd 100644 --- a/optyx/core/backends.py +++ b/optyx/core/backends.py @@ -726,7 +726,6 @@ def _process_term( state_provided = "perceval_state" in extra effect_provided = "perceval_effect" in extra is_dom_closed = len(diagram.dom) == 0 - single_output_task = task in ("single_amp", "single_prob") if not state_provided and not is_dom_closed: raise ValueError( @@ -756,22 +755,31 @@ def _process_term( perceval_effect = None if effect_provided: - external_perceval_effect = self._process_effect( + perceval_effect = self._process_effect( extra["perceval_effect"] ) - if external_perceval_effect.m != matrix.cod: + if perceval_effect.m != matrix.cod: raise ValueError( "The provided 'perceval_effect' does not match " "the number of output modes of the diagram." ) - if perceval_effect is not None: + if matrix.cod == 0: + perceval_effect = pcvl.BasicState([matrix.selections[0]]) + matrix.selections = matrix.selections[1:] + matrix.cod = 1 + + if ( + perceval_effect is not None + ): if task == "amps": task = "single_amp" if task == "probs": task = "single_prob" + single_output_task = task in ("single_amp", "single_prob") + if single_output_task: if perceval_effect is None: raise ValueError( diff --git a/test/test_backends.py b/test/test_backends.py index 4b4f9d7..eaf1700 100644 --- a/test/test_backends.py +++ b/test/test_backends.py @@ -15,6 +15,19 @@ import perceval as pcvl import discopy.tensor as discopy_tensor +unitary_circuit = photonic.BS +non_unitary_circuit = ( + photonic.MZI(0.23, 0.51) >> + photonic.NumOp() @ photonic.NumOp() +) + +bs_state = lambda nmodes: tuple(1 for _ in range(nmodes)) + +def _compare_prob_for_outcome(diagram, outcome): + res_q = diagram.eval() + d = res_q.prob_dist() + return float(d.get(outcome, 0.0)) + @pytest.mark.skip(reason="Helper function for testing") def chip_mzi(w, l): ansatz = photonic.ansatz(w, l) @@ -488,4 +501,269 @@ def test_permanent_prob_matches_quimb_with_create(self, circuit): backend = PermanentBackend() result_perm = diagram.eval(backend, return_kind="prob") - assert dict_allclose(result_perm.prob_dist(), result_quimb.prob_dist()) \ No newline at end of file + assert dict_allclose(result_perm.prob_dist(), result_quimb.prob_dist()) + +@pytest.mark.parametrize("circuit", [unitary_circuit, non_unitary_circuit]) +def test_prob_nothing_raises_no_state(circuit): + backend = PercevalBackend() + with pytest.raises(ValueError): + _ = circuit.eval(backend, task="probs") + +@pytest.mark.parametrize("circuit", [unitary_circuit, non_unitary_circuit]) +def test_prob_state_compare_with_quimb(circuit): + print(circuit) + nmodes = len(circuit.dom) + perceval_state = pcvl.BasicState(bs_state(nmodes)) + backend = PercevalBackend() + res_pcvl = circuit.eval(backend, perceval_state=perceval_state, task="probs") + + ref = photonic.Create(*bs_state(nmodes)) >> circuit + d_ref = ref.eval().prob_dist() + d_pcvl = res_pcvl.prob_dist() + + keys = set(d_ref) | set(d_pcvl) + for k in keys: + assert math.isclose(d_ref.get(k, 0.0), d_pcvl.get(k, 0.0), rel_tol=1e-9, abs_tol=1e-12) + +@pytest.mark.parametrize("circuit", [unitary_circuit, non_unitary_circuit]) +def test_prob_effect_only_errors_no_state(circuit): + nmodes = len(circuit.dom) + backend = PercevalBackend() + with pytest.raises(ValueError): + _ = circuit.eval( + backend, + perceval_effect=pcvl.BasicState([2] + [0]*(nmodes-1)), + task="probs" + ) + +@pytest.mark.parametrize("circuit", [unitary_circuit, non_unitary_circuit]) +def test_prob_state_plus_effect_coerces_to_single_prob_and_matches_quimb(circuit): + + nmodes = len(circuit.dom) + state_occ = bs_state(nmodes) + outcome = (2,) + (0,)*(nmodes-1) + + backend = PercevalBackend() + perceval_state = pcvl.BasicState(state_occ) + perceval_effect = pcvl.BasicState(outcome) + + res = circuit.eval( + backend, + perceval_state=perceval_state, + perceval_effect=perceval_effect, + task="probs" + ) + assert res.state_type is StateType.SINGLE_PROB + p_pcvl = float(res.single_prob(outcome)) + + ref = photonic.Create(*state_occ) >> circuit + p_ref = _compare_prob_for_outcome(ref, outcome) + assert math.isclose(p_pcvl, p_ref, rel_tol=1e-9, abs_tol=1e-12) + +@pytest.mark.parametrize("circuit", [unitary_circuit, non_unitary_circuit]) +def test_prob_diagram_only_effect_errors_no_state(circuit): + nmodes = len(circuit.cod) + diagram = circuit >> classical.Select(*([0]*(nmodes-1) + [2])) + backend = PercevalBackend() + with pytest.raises(ValueError): + _ = diagram.eval(backend, task="probs") + +@pytest.mark.parametrize("circuit", [unitary_circuit, non_unitary_circuit]) +def test_prob_diagram_state_and_effect_coerces_to_single_prob_and_matches_quimb(circuit): + nmodes = len(circuit.dom) + outcome = (0,) + (0,)*(nmodes-1) + + diagram = photonic.Create(*bs_state(nmodes)) >> circuit >> classical.Select(*outcome) + backend = PercevalBackend() + res = diagram.eval(backend, task="probs") + assert res.state_type is StateType.SINGLE_PROB + + p_pcvl = float(res.single_prob(outcome)) + + ref = photonic.Create(*bs_state(nmodes)) >> circuit + p_ref = _compare_prob_for_outcome(ref, outcome) + assert math.isclose(p_pcvl, p_ref, rel_tol=1e-9, abs_tol=1e-12) + +@pytest.mark.parametrize("circuit", [unitary_circuit, non_unitary_circuit]) +def test_single_prob_errors_without_state_and_effect(circuit): + backend = PercevalBackend() + with pytest.raises(ValueError): + _ = circuit.eval(backend, task="single_prob") + +@pytest.mark.parametrize("circuit", [unitary_circuit, non_unitary_circuit]) +def test_single_prob_errors_with_state_only(circuit): + backend = PercevalBackend() + nmodes = len(circuit.dom) + with pytest.raises(ValueError): + _ = circuit.eval( + backend, + perceval_state=pcvl.BasicState(bs_state(nmodes)), + task="single_prob" + ) + +@pytest.mark.parametrize("circuit", [unitary_circuit, non_unitary_circuit]) +def test_single_prob_errors_with_effect_only(circuit): + backend = PercevalBackend() + nmodes = len(circuit.dom) + with pytest.raises(ValueError): + _ = circuit.eval( + backend, + perceval_effect=pcvl.BasicState([2] + [0]*(nmodes-1)), + task="single_prob" + ) + + backend = PercevalBackend() + nmodes = len(circuit.dom) + with pytest.raises(ValueError): + _ = circuit.eval( + backend, + perceval_effect=pcvl.BasicState([2] + [0]*(nmodes-1)), + task="single_prob" + ) + +@pytest.mark.parametrize("circuit", [unitary_circuit, non_unitary_circuit]) +def test_single_prob_matches_quimb_selected_outcome(circuit): + nmodes = len(circuit.dom) + state_occ = bs_state(nmodes) + outcome = (2,) + (0,)*(nmodes-1) if nmodes >= 1 else tuple() + + backend = PercevalBackend() + res = circuit.eval( + backend, + perceval_state=pcvl.BasicState(state_occ), + perceval_effect=pcvl.BasicState(outcome), + task="single_prob" + ) + assert res.state_type is StateType.SINGLE_PROB + p_pcvl = float(res.single_prob(outcome)) + + ref = photonic.Create(*state_occ) >> circuit + p_ref = _compare_prob_for_outcome(ref, outcome) + assert math.isclose(p_pcvl, p_ref, rel_tol=1e-9, abs_tol=1e-12) + +@pytest.mark.parametrize("circuit", [unitary_circuit, non_unitary_circuit]) +def test_amp_nothing_raises_no_state(circuit): + backend = PercevalBackend() + with pytest.raises(ValueError): + _ = circuit.eval(backend, task="amps") # no state + +@pytest.mark.parametrize("circuit", [unitary_circuit, non_unitary_circuit]) +def test_amp_state_compare_with_quimb(circuit): + + nmodes = len(circuit.dom) + state_occ = bs_state(nmodes) + + backend = PercevalBackend() + res_pcvl = circuit.eval( + backend, + perceval_state=pcvl.BasicState(state_occ), + task="amps" + ) + + # reference + ref = photonic.Create(*state_occ) >> circuit + d_ref = ref.eval().amplitudes() + d_pcvl = res_pcvl.amplitudes() + + keys = set(d_ref) | set(d_pcvl) + for k in keys: + assert complex(d_ref.get(k, 0.0)) == pytest.approx(d_pcvl.get(k, 0.0), rel=1e-9, abs=1e-12) + +@pytest.mark.parametrize("circuit", [unitary_circuit, non_unitary_circuit]) +def test_amp_effect_only_errors_no_state(circuit): + backend = PercevalBackend() + nmodes = len(circuit.dom) + with pytest.raises(ValueError): + _ = circuit.eval( + backend, + perceval_effect=pcvl.BasicState([2] + [0]*(nmodes-1)), + task="amps" + ) + +@pytest.mark.parametrize("circuit", [unitary_circuit, non_unitary_circuit]) +def test_amp_state_plus_effect_coerces_to_single_amp_and_matches_quimb(circuit): + + nmodes = len(circuit.dom) + state_occ = bs_state(nmodes) + outcome = (2,) + (0,)*(nmodes-1) + + backend = PercevalBackend() + res = circuit.eval( + backend, + perceval_state=pcvl.BasicState(state_occ), + perceval_effect=pcvl.BasicState(outcome), + task="amps" + ) + assert res.state_type is StateType.SINGLE_AMP + a_pcvl = res.single_amplitude(outcome) + + ref = photonic.Create(*state_occ) >> circuit + a_ref = ref.eval().amplitudes().get(outcome, 0.0) + assert a_pcvl == pytest.approx(a_ref, rel=1e-9, abs=1e-12) + +@pytest.mark.parametrize("circuit", [unitary_circuit, non_unitary_circuit]) +def test_amp_diagram_state_and_effect_coerces_to_single_amp_and_matches_quimb(circuit): + + nmodes = len(circuit.dom) + state_occ = bs_state(nmodes) + outcome = (0,) + (0,)*(nmodes-1) + + diagram = photonic.Create(*state_occ) >> circuit >> classical.Select(*outcome) + backend = PercevalBackend() + res = diagram.eval(backend, task="amps") + assert res.state_type is StateType.SINGLE_AMP + a_pcvl = res.single_amplitude(outcome) + + ref = photonic.Create(*state_occ) >> circuit + a_ref = ref.eval().amplitudes().get(outcome, 0.0) + assert a_pcvl == pytest.approx(a_ref, rel=1e-9, abs=1e-12) + +@pytest.mark.parametrize("circuit", [unitary_circuit, non_unitary_circuit]) +def test_single_amp_errors_without_state_and_effect(circuit): + backend = PercevalBackend() + with pytest.raises(ValueError): + _ = circuit.eval(backend, task="single_amp") + +@pytest.mark.parametrize("circuit", [unitary_circuit, non_unitary_circuit]) +def test_single_amp_errors_with_state_only(circuit): + backend = PercevalBackend() + nmodes = len(circuit.dom) + with pytest.raises(ValueError): + _ = circuit.eval( + backend, + perceval_state=pcvl.BasicState(bs_state(nmodes)), + task="single_amp" + ) + + +@pytest.mark.parametrize("circuit", [unitary_circuit, non_unitary_circuit]) +def test_single_amp_errors_with_effect_only(circuit): + backend = PercevalBackend() + nmodes = len(circuit.dom) + with pytest.raises(ValueError): + _ = circuit.eval( + backend, + perceval_effect=pcvl.BasicState([2] + [0]*(nmodes-1)), + task="single_amp" + ) + +@pytest.mark.parametrize("circuit", [unitary_circuit, non_unitary_circuit]) +def test_single_amp_matches_quimb_selected_outcome(circuit): + + nmodes = len(circuit.dom) + state_occ = bs_state(nmodes) + outcome = (2,) + (0,)*(nmodes-1) if nmodes >= 1 else tuple() + + backend = PercevalBackend() + res = circuit.eval( + backend, + perceval_state=pcvl.BasicState(state_occ), + perceval_effect=pcvl.BasicState(outcome), + task="single_amp" + ) + assert res.state_type is StateType.SINGLE_AMP + a_pcvl = res.single_amplitude(outcome) + + ref = photonic.Create(*state_occ) >> circuit + a_ref = ref.eval().amplitudes().get(outcome, 0.0) + assert a_pcvl == pytest.approx(a_ref, rel=1e-9, abs=1e-12) \ No newline at end of file From 146840634a4049ce2cdcc3345634e8b06e4a62c7 Mon Sep 17 00:00:00 2001 From: mateuszkupper Date: Mon, 20 Oct 2025 09:11:49 +0100 Subject: [PATCH 12/19] Add tests for perceval_backend --- test/test_backends.py | 25 ++++++++++++++++++++++++- 1 file changed, 24 insertions(+), 1 deletion(-) diff --git a/test/test_backends.py b/test/test_backends.py index eaf1700..039a87b 100644 --- a/test/test_backends.py +++ b/test/test_backends.py @@ -766,4 +766,27 @@ def test_single_amp_matches_quimb_selected_outcome(circuit): ref = photonic.Create(*state_occ) >> circuit a_ref = ref.eval().amplitudes().get(outcome, 0.0) - assert a_pcvl == pytest.approx(a_ref, rel=1e-9, abs=1e-12) \ No newline at end of file + assert a_pcvl == pytest.approx(a_ref, rel=1e-9, abs=1e-12) + +def _test_sums(): + diagram = ( + photonic.Create(1, 1, 1, 1) >> + photonic.MZI(0.5, 0.5) @ photonic.MZI(0.5, 0.5) >> + photonic.qmode @ photonic.MZI(0.5, 0.5) @ photonic.qmode + ) + ( + photonic.Create(0, 2, 0, 2) >> + photonic.MZI(0.5, 0.5) @ photonic.MZI(0.5, 0.5) >> + photonic.qmode @ photonic.MZI(0.5, 0.5) @ photonic.qmode + ) + + res_perceval = diagram.eval(PercevalBackend()).prob_dist() + res_quimb = diagram.eval().prob_dist() + + keys = set(res_perceval) | set(res_quimb) + for k in keys: + assert math.isclose( + res_perceval.get(k, 0.0), + res_quimb.get(k, 0.0), + rel_tol=1e-9, + abs_tol=1e-12 + ) From fd4ad67e44e3ea6ee12ec92105ab9c8dab95f8d0 Mon Sep 17 00:00:00 2001 From: mateuszkupper Date: Mon, 20 Oct 2025 09:25:17 +0100 Subject: [PATCH 13/19] Test permanent and perceval sum eval --- test/test_backends.py | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/test/test_backends.py b/test/test_backends.py index 039a87b..9ab261e 100644 --- a/test/test_backends.py +++ b/test/test_backends.py @@ -780,6 +780,7 @@ def _test_sums(): ) res_perceval = diagram.eval(PercevalBackend()).prob_dist() + res_permanent = diagram.eval(PermanentBackend()).prob_dist() res_quimb = diagram.eval().prob_dist() keys = set(res_perceval) | set(res_quimb) @@ -790,3 +791,12 @@ def _test_sums(): rel_tol=1e-9, abs_tol=1e-12 ) + + keys = set(res_permanent) | set(res_quimb) + for k in keys: + assert math.isclose( + res_permanent.get(k, 0.0), + res_quimb.get(k, 0.0), + rel_tol=1e-9, + abs_tol=1e-12 + ) From fdd0113ac4e085c60e62f721e18caa8198af6732 Mon Sep 17 00:00:00 2001 From: mateuszkupper Date: Mon, 20 Oct 2025 10:41:15 +0100 Subject: [PATCH 14/19] Finalise tests for backends --- optyx/core/backends.py | 54 +++++++++++++++++++++++++++++++++--------- test/test_backends.py | 5 ++-- 2 files changed, 46 insertions(+), 13 deletions(-) diff --git a/optyx/core/backends.py b/optyx/core/backends.py index 8eef2cd..17b12b4 100644 --- a/optyx/core/backends.py +++ b/optyx/core/backends.py @@ -176,10 +176,28 @@ def prob_dist(self, round_digits: int = None) -> dict: if self.state_type is StateType.DM: return self._prob_dist_mixed(round_digits) if self.state_type is StateType.PROB: - return self._convert_array_to_dict( + values = self._convert_array_to_dict( self.tensor.array, - round_digits=round_digits + round_digits ) + + if np.allclose(np.sum(list(values.values())), 1): + return values + + probs = {} + for k, v in values.items(): + val = float(np.real_if_close(v)) + if val < 0 and abs(val) < 1e-12: + val = 0.0 + probs[k] = val + + total = float(np.sum(list(probs.values()))) + if total == 0: + raise ValueError("The probability distribution sums to zero.") + + norm = {k: v / total for k, v in probs.items()} + return norm + raise ValueError("Unsupported state_type type. " + "Must be StateType.AMP, StateType.DM, " + "or StateType.PROB.") @@ -560,7 +578,6 @@ def __init__(self, perceval_backend: pcvl.ABackend = None): else: self.perceval_backend = perceval_backend - # pylint: disable=too-many-locals def eval( self, diagram: Diagram, @@ -602,7 +619,7 @@ def eval( for term in diagram.terms: arr, output_types, output_tensor_cod, return_type = \ self._process_term( - term, diagram, task, tensor_diagram, extra + term, term, task, tensor_diagram, extra ) array += arr else: @@ -610,7 +627,6 @@ def eval( self._process_term( diagram, diagram, task, tensor_diagram, extra ) - return EvalResult( discopy_tensor.Box( "Result", @@ -910,29 +926,45 @@ def check_creations(matrix): "n_photons must be greater than 0." ) - tensor_diagram = self._get_discopy_tensor(diagram) if hasattr( diagram, "terms" ): result = 0 + dims = [] for term in diagram.terms: matrix = self._get_matrix(term) check_creations(matrix) - result += matrix.eval( + result_matrix = matrix.eval( n_photons=n_photons, as_tensor=True - ).array + ) + result += result_matrix.array + dims.append( + ( + result_matrix.dom.inside, # list + result_matrix.cod.inside # list + ) + ) else: matrix = self._get_matrix(diagram) check_creations(matrix) - result = matrix.eval(n_photons=n_photons, as_tensor=True).array + result_matrix = matrix.eval(n_photons=n_photons, as_tensor=True) + result = result_matrix.array + dims = [(result_matrix.dom.inside, result_matrix.cod.inside)] + + norm = lambda x: list(x) if isinstance( # noqa: E731 + x, (list, tuple) + ) else [x] + dom_lists, cod_lists = zip(*((norm(d), norm(c)) for d, c in dims)) + max_dom_dims = [max(vals) for vals in zip(*dom_lists)] + max_cod_dims = [max(vals) for vals in zip(*cod_lists)] return EvalResult( discopy_tensor.Box( "Result", - tensor_diagram.dom, - tensor_diagram.cod, + discopy_tensor.Dim(*tuple(max_dom_dims)), + discopy_tensor.Dim(*tuple(max_cod_dims)), result ), output_types=diagram.cod, diff --git a/test/test_backends.py b/test/test_backends.py index 9ab261e..a28764a 100644 --- a/test/test_backends.py +++ b/test/test_backends.py @@ -768,7 +768,7 @@ def test_single_amp_matches_quimb_selected_outcome(circuit): a_ref = ref.eval().amplitudes().get(outcome, 0.0) assert a_pcvl == pytest.approx(a_ref, rel=1e-9, abs=1e-12) -def _test_sums(): +def test_sums(): diagram = ( photonic.Create(1, 1, 1, 1) >> photonic.MZI(0.5, 0.5) @ photonic.MZI(0.5, 0.5) >> @@ -782,7 +782,8 @@ def _test_sums(): res_perceval = diagram.eval(PercevalBackend()).prob_dist() res_permanent = diagram.eval(PermanentBackend()).prob_dist() res_quimb = diagram.eval().prob_dist() - + print(res_perceval) + print(res_quimb) keys = set(res_perceval) | set(res_quimb) for k in keys: assert math.isclose( From 02247eecf3ff99726801d5cee9242eda57b092cc Mon Sep 17 00:00:00 2001 From: mateuszkupper Date: Mon, 20 Oct 2025 10:54:37 +0100 Subject: [PATCH 15/19] Fix an issue in zx --- optyx/core/zx.py | 39 +++++++++++++++------------------------ 1 file changed, 15 insertions(+), 24 deletions(-) diff --git a/optyx/core/zx.py b/optyx/core/zx.py index 4834aed..00686b0 100644 --- a/optyx/core/zx.py +++ b/optyx/core/zx.py @@ -128,15 +128,15 @@ def move(scan, source, target): swaps = Id(diagram.Bit(len(scan))) return scan, swaps - def make_wires_adjacent(scan, diagram, inputs): + def make_wires_adjacent(scan, dgrm, inputs): if not inputs: - return scan, diagram, len(scan) + return scan, dgrm, len(scan) offset = scan.index(inputs[0]) for i, _ in enumerate(inputs[1:]): source, target = scan.index(inputs[i + 1]), offset + i + 1 scan, swaps = move(scan, source, target) - diagram = diagram >> swaps - return scan, diagram, offset + dgrm = dgrm >> swaps + return scan, dgrm, offset missing_boundary = any( graph.type(node) == VertexType.BOUNDARY # noqa: E721 @@ -238,11 +238,17 @@ def to_pyzx(self): graph.set_inputs(graph.inputs() + (node,)) graph.set_position(node, i, 0) for row, (box, offset) in enumerate(zip(self.boxes, self.offsets)): - if isinstance(box, Spider): - node = graph.add_vertex( - (VertexType.Z if isinstance(box, Z) else VertexType.X), - phase=box.phase * 2 if box.phase else None, - ) + if isinstance(box, diagram.Spider): + if isinstance(box, Spider): + node = graph.add_vertex( + (VertexType.Z if isinstance(box, Z) else VertexType.X), + phase=box.phase * 2 if box.phase else None, + ) + else: + node = graph.add_vertex( + VertexType.Z, + phase=box.phase * 2 if box.phase else None, + ) graph.set_position(node, offset, row + 1) for i, _ in enumerate(box.dom): source, hadamard = scan[offset + i] @@ -264,21 +270,6 @@ def to_pyzx(self): elif box == H: node, hadamard = scan[offset] scan[offset] = (node, not hadamard) - elif isinstance(box, diagram.Spider): - node = graph.add_vertex( - VertexType.Z, - phase=box.phase * 2 if box.phase else None, - ) - graph.set_position(node, offset, row + 1) - for i, _ in enumerate(box.dom): - source, hadamard = scan[offset + i] - etype = EdgeType.HADAMARD if hadamard else EdgeType.SIMPLE - graph.add_edge((source, node), etype) - scan = ( - scan[:offset] - + len(box.cod) * [(node, False)] - + scan[offset + len(box.dom):] - ) else: raise NotImplementedError for i, _ in enumerate(self.cod): From 6c0d38837694c9a9a047c4b1bce2549c0329ca92 Mon Sep 17 00:00:00 2001 From: mateuszkupper Date: Mon, 20 Oct 2025 11:23:15 +0100 Subject: [PATCH 16/19] Minor refactor of backends --- optyx/core/backends.py | 39 +++++++++++++++++++-------------------- 1 file changed, 19 insertions(+), 20 deletions(-) diff --git a/optyx/core/backends.py b/optyx/core/backends.py index 17b12b4..167ea58 100644 --- a/optyx/core/backends.py +++ b/optyx/core/backends.py @@ -608,30 +608,33 @@ def eval( The result of the evaluation (EvalResult). """ - task = extra.get("task", "probs") - tensor_diagram = self._get_discopy_tensor(diagram) - if hasattr( diagram, "terms" ): array = 0 for term in diagram.terms: - arr, output_types, output_tensor_cod, return_type = \ + arr, output_types, return_type = \ self._process_term( - term, term, task, tensor_diagram, extra + term, extra ) array += arr else: - array, output_types, output_tensor_cod, return_type = \ + array, output_types, return_type = \ self._process_term( - diagram, diagram, task, tensor_diagram, extra + diagram, extra ) + + if array.shape == (1,): + cod = discopy_tensor.Dim(1) + else: + cod = discopy_tensor.Dim(*array.shape) + return EvalResult( discopy_tensor.Box( "Result", discopy_tensor.Dim(1), - output_tensor_cod, + cod, array ), output_types=output_types, @@ -726,9 +729,6 @@ def _dilate( def _process_term( self, term, - diagram, - task, - tensor_diagram, extra ): """ @@ -739,9 +739,10 @@ def _process_term( """ matrix = self._get_matrix(term) + task = extra.get("task", "probs") state_provided = "perceval_state" in extra effect_provided = "perceval_effect" in extra - is_dom_closed = len(diagram.dom) == 0 + is_dom_closed = len(term.dom) == 0 if not state_provided and not is_dom_closed: raise ValueError( @@ -822,7 +823,7 @@ def _process_term( perceval_circuit = self._umatrix_to_perceval_circuit(matrix.array) sim.set_circuit(perceval_circuit) - m_orig = len(diagram.dom) + m_orig = len(term.dom) k_extra = matrix.dom - m_orig result = None p = None @@ -838,9 +839,8 @@ def _process_term( return_type = StateType.AMP result = self._post_select_vacuum(result, m_orig, k_extra) - shape = tensor_diagram.cod.inside - output_types = diagram.cod - output_tensor_cod = tensor_diagram.cod + output_types = term.cod + output_shape = self._get_discopy_tensor(term).cod.inside elif single_output_task: if task == "single_prob": @@ -853,9 +853,8 @@ def _process_term( ) return_type = StateType.SINGLE_AMP - shape = (1,) + output_shape = (1,) output_types = None - output_tensor_cod = discopy_tensor.Dim(1) else: raise ValueError( @@ -864,7 +863,7 @@ def _process_term( ) array = np.zeros( - shape, + output_shape, dtype=float if task in ("single_prob", "probs") else complex ) @@ -887,7 +886,7 @@ def _process_term( array[0] = p else: pass - return array, output_types, output_tensor_cod, return_type + return array, output_types, return_type class PermanentBackend(AbstractBackend): From e15feaa3c619c6758c99e92157f3ccf6f643009b Mon Sep 17 00:00:00 2001 From: mateuszkupper Date: Mon, 20 Oct 2025 13:23:54 +0100 Subject: [PATCH 17/19] Change perceval default eval behaviour for closed diagrams --- optyx/core/backends.py | 241 ++++++++++++++++++++++------------------- test/test_backends.py | 185 ++++++++++--------------------- 2 files changed, 192 insertions(+), 234 deletions(-) diff --git a/optyx/core/backends.py b/optyx/core/backends.py index 167ea58..6c1f557 100644 --- a/optyx/core/backends.py +++ b/optyx/core/backends.py @@ -67,7 +67,7 @@ import warnings from abc import ABC, abstractmethod from dataclasses import dataclass -from typing import Any +from typing import Any, Literal, Sequence from collections import defaultdict from enum import Enum from cotengra import ( @@ -85,6 +85,35 @@ from optyx.utils.utils import preprocess_quimb_tensors_safe +@dataclass +class PercevalEvalConfig: + """ + Configuration for Perceval backend evaluation. + - task: The Perceval task to perform. Allowed values are + "probs" (default), "amps", "single_amp", "single_prob". + - state: A `perceval.BasicState` or a sequence of + non-negative integers (occupation numbers). Defaults + to a bosonic product state |11...1> if the + diagram does not include any photon creations. + Either the creations for all input ports are specified + by the diagram (`Create(...)`) or the user must provide + the `state` argument covering all input ports. + - effect: Required if task is "single_amp" or "single_prob". + A sequence of non-negative integers (occupation numbers) + specifying the output configuration for which to compute + the amplitude or probability. + """ + task: Literal[ + "probs", "amps", "single_amp", "single_prob" + ] = "probs" + state: pcvl.BasicState | Sequence[int] | None = None + effect: pcvl.BasicState | Sequence[int] | None = None + + +ProbDist = dict[tuple[int, ...], float] +Amps = dict[tuple[int, ...], complex] + + class StateType(Enum): """ Enum to represent the type of state represented by the result tensor. @@ -102,7 +131,7 @@ class EvalResult: Class to encapsulate the result of an evaluation of a diagram. """ _tensor: discopy_tensor.Box - output_types: Ty + output_types: Ty | None state_type: StateType @property @@ -125,7 +154,8 @@ def density_matrix(self) -> np.ndarray: """ if len(self.tensor.dom) != 0: raise ValueError( - "Result tensor must represent a state with no inputs." + "Result tensor must represent a state with no inputs. " + + f"Current domain: {self.tensor.dom}" ) if self.state_type not in {StateType.AMP, StateType.DM}: raise TypeError( @@ -136,7 +166,10 @@ def density_matrix(self) -> np.ndarray: return density_matrix.array return self.tensor.array - def amplitudes(self, normalise=True) -> dict[tuple[int, ...], complex]: + def amplitudes( + self, + normalise=True + ) -> Amps: """ Get the amplitudes from the result tensor. Returns: @@ -149,7 +182,8 @@ def amplitudes(self, normalise=True) -> dict[tuple[int, ...], complex]: ) if len(self.tensor.dom) != 0: raise ValueError( - "Result tensor must represent a state with no inputs." + "Result tensor must represent a state with no inputs. " + + f"Current domain: {self.tensor.dom}" ) dic = self._convert_array_to_dict(self.tensor.array) @@ -159,7 +193,10 @@ def amplitudes(self, normalise=True) -> dict[tuple[int, ...], complex]: for key, value in dic.items()} return dic - def prob_dist(self, round_digits: int = None) -> dict: + def prob_dist( + self, + round_digits: int = None + ) -> ProbDist: """ Get the probability distribution from the result tensor. @@ -169,38 +206,38 @@ def prob_dist(self, round_digits: int = None) -> dict: """ if len(self.tensor.dom) != 0: raise ValueError( - "Result tensor must represent a state with no inputs." + "Result tensor must represent a state with no inputs. " + + f"Current domain: {self.tensor.dom}" ) if self.state_type is StateType.AMP: - return self._prob_dist_pure(round_digits) - if self.state_type is StateType.DM: - return self._prob_dist_mixed(round_digits) - if self.state_type is StateType.PROB: + values = self._prob_dist_pure() + elif self.state_type is StateType.DM: + values = self._prob_dist_mixed() + elif self.state_type is StateType.PROB: values = self._convert_array_to_dict( - self.tensor.array, - round_digits + self.tensor.array + ) + else: + raise ValueError( + f"Unsupported state_type type: {self.state_type}. " + + "Must be StateType.AMP, StateType.DM, " + + "or StateType.PROB." ) - if np.allclose(np.sum(list(values.values())), 1): - return values - - probs = {} - for k, v in values.items(): - val = float(np.real_if_close(v)) - if val < 0 and abs(val) < 1e-12: - val = 0.0 - probs[k] = val - - total = float(np.sum(list(probs.values()))) + if not np.allclose( + np.sum(list(values.values())), 1, atol=1e-12 + ): + total = float(np.sum(list(values.values()))) if total == 0: raise ValueError("The probability distribution sums to zero.") - norm = {k: v / total for k, v in probs.items()} - return norm + prob = {k: v / total for k, v in values.items()} + else: + prob = values - raise ValueError("Unsupported state_type type. " + - "Must be StateType.AMP, StateType.DM, " + - "or StateType.PROB.") + if round_digits is None: + return prob + return {k: round(v, round_digits) for k, v in prob.items()} def single_prob(self, occupation: tuple) -> float: """ @@ -236,8 +273,7 @@ def single_amplitude(self, occupation: tuple) -> complex: def _convert_array_to_dict( self, - array: np.ndarray, - round_digits: int = None) -> dict: + array: np.ndarray) -> dict: """ Return a dict that maps multi-indices - values for all non-zero entries of an array. @@ -250,13 +286,9 @@ def _convert_array_to_dict( nz_vals = array.flat[nz_flat] nz_multi = np.vstack(np.unravel_index(nz_flat, array.shape)).T - if round_digits is not None: - return {tuple(idx): np.round(val, round_digits) for - idx, val in zip(nz_multi, nz_vals)} - return {tuple(idx): val for idx, val in zip(nz_multi, nz_vals)} - def _prob_dist_pure(self, round_digits: int = None) -> dict: + def _prob_dist_pure(self) -> ProbDist: """ Get the probability distribution for a pure state. @@ -266,17 +298,12 @@ def _prob_dist_pure(self, round_digits: int = None) -> dict: """ values = self._convert_array_to_dict( - self.tensor.array, - round_digits=round_digits + self.tensor.array ) - sum_ = np.sum(np.abs(list(values.values())) ** 2) - if sum_ == 0: - raise ValueError("The probability distribution sums to zero.") - return {key: (abs(value) ** 2)/sum_ for key, value in values.items()} - def _prob_dist_mixed( - self, - round_digits: int | None = None) -> dict[tuple[int, ...], float]: + return {k: abs(v) ** 2 for k, v in values.items()} + + def _prob_dist_mixed(self) -> ProbDist: """ Get the probability distribution from a mixed state. This method computes the probability distribution by aggregating @@ -297,10 +324,11 @@ def _prob_dist_mixed( if not any(t in {bit, mode} for t in self.output_types): raise ValueError( "Output types must contain at least one 'bit' or 'mode'." + - "These will be treated as measured registers." + "These will be treated as measured registers. " + + f"Current output types: {self.output_types}" ) - values = self._convert_array_to_dict(self.tensor.array, None) + values = self._convert_array_to_dict(self.tensor.array) mask_flat = np.concatenate( [[1] if t in {bit, mode} else [0, 0] for t in self.output_types] ) @@ -324,18 +352,8 @@ def _prob_dist_mixed( for occ in all_measured: probs.setdefault(occ, 0.0) - sum_ = np.sum(list(probs.values())) - if sum_ == 0: - raise ValueError("The probability distribution sums to zero.") - prob = { - key: value / sum_ - for key, value in probs.items() - } - if round_digits is not None: - prob = {k: round(v, round_digits) for k, v in prob.items()} - - return prob + return probs # pylint: disable=too-few-public-methods @@ -510,7 +528,9 @@ def _process_term(self, term: discopy_tensor.Diagram) -> np.ndarray: "Unsupported hyperoptimiser type. " + "Use ReusableHyperOptimizer, HyperOptimizer, " + "ReusableHyperCompressedOptimizer, or " + - "HyperCompressedOptimizer." + "HyperCompressedOptimizer. " + + f"Got: {type(self.hyperoptimiser)}" + ) if is_approx: @@ -585,25 +605,13 @@ def eval( """ Evaluate the diagram using Perceval. Works only for unitary operations. - If no `perceval_state` is provided in `extra`, - it defaults to a bosonic product state. Args: diagram (Diagram): The diagram to evaluate. **extra: Additional arguments for the evaluation: - - perceval_state: A `perceval.BasicState` or a sequence of - non-negative integers (occupation numbers). Defaults - to a bosonic product state |11...1> if the - diagram does not include any photon creations. - Either the creations for all input ports are specified - by the diagram (`Create(...)`) or the user must provide - the `perceval_state` argument covering all input ports. - - task: The Perceval task to perform. Allowed values are - "probs" (default), "amps", "single_amp", "single_prob". - - out: Required if task is "single_amp" or "single_prob". - A sequence of non-negative integers (occupation numbers) - specifying the output configuration for which to compute - the amplitude or probability. + - config: Configuration for Perceval evaluation, + see :class:`PercevalEvalConfig`. + Returns: The result of the evaluation (EvalResult). """ @@ -685,7 +693,8 @@ def _process_state(self, perceval_state): raise TypeError( "perceval_state must be a perceval.BasicState" " or a sequence of non-negative " + - "integers (occupation numbers)." + "integers (occupation numbers). " + + f"Got: {type(perceval_state)}" ) from e return perceval_state @@ -695,11 +704,12 @@ def _process_effect(self, perceval_effect): if not isinstance(perceval_effect, pcvl.BasicState): try: perceval_effect = pcvl.BasicState(list(perceval_effect)) - except Exception as e: - raise TypeError( + except ValueError as e: + raise ValueError( "perceval_effect must be a perceval.BasicState" " or a sequence of non-negative " + - "integers (occupation numbers)." + "integers (occupation numbers). " + + f"Got: {type(perceval_effect)}" ) from e return perceval_effect @@ -739,28 +749,30 @@ def _process_term( """ matrix = self._get_matrix(term) - task = extra.get("task", "probs") - state_provided = "perceval_state" in extra - effect_provided = "perceval_effect" in extra + cfg: PercevalEvalConfig = extra.get("config", PercevalEvalConfig()) + task = cfg.task + state_provided = cfg.state is not None + effect_provided = cfg.effect is not None is_dom_closed = len(term.dom) == 0 if not state_provided and not is_dom_closed: raise ValueError( - "External 'perceval_state' not provided but the diagram " - "has open input modes. Provide a 'perceval_state' " + - "or close all input modes with a state." + "External 'state' not provided but " + + "the diagram has open input modes. " + + "Provide a 'state' or close all input modes with a state. " + + f"Open input modes: {term.dom}" ) external_perceval_state = pcvl.BasicState([]) if state_provided: - external_perceval_state = self._process_state( - extra["perceval_state"] - ) + external_perceval_state = self._process_state(cfg.state) if external_perceval_state.m != matrix.dom: raise ValueError( - "The provided 'perceval_state' does not match " - "the number of input modes of the diagram." + "The provided 'state' does not match the " + + "number of input modes of the diagram. " + + f"Provided state has {external_perceval_state.m} " + + f"modes but diagram has {matrix.dom} input modes." ) perceval_state = self._process_state( @@ -772,32 +784,42 @@ def _process_term( perceval_effect = None if effect_provided: - perceval_effect = self._process_effect( - extra["perceval_effect"] - ) + perceval_effect = self._process_effect(cfg.effect) if perceval_effect.m != matrix.cod: raise ValueError( - "The provided 'perceval_effect' does not match " - "the number of output modes of the diagram." + "The provided 'effect' does not match the number " + + "of output modes of the diagram. " + + f"Provided effect has {perceval_effect.m} " + + f"modes but diagram has {matrix.cod} output modes." ) + # convert post-selections to pcvl effects and post-selections if matrix.cod == 0: - perceval_effect = pcvl.BasicState([matrix.selections[0]]) - matrix.selections = matrix.selections[1:] - matrix.cod = 1 + sel0 = matrix.selections[0] + m = Matrix( + matrix.array.copy(), + matrix.dom, + 1, + creations=list(matrix.creations), + selections=list(matrix.selections[1:]) + ) + perceval_effect = pcvl.BasicState([sel0]) + matrix = m if ( - perceval_effect is not None + perceval_effect is not None and + task in ("amps", "probs") ): - if task == "amps": - task = "single_amp" - if task == "probs": - task = "single_prob" + raise ValueError( + f"An 'effect' was provided but task='{task}'. " + "Use task='single_amp' or task='single_prob' " + + "when conditioning on an effect." + ) - single_output_task = task in ("single_amp", "single_prob") + is_single_output_task = task in ("single_amp", "single_prob") - if single_output_task: + if is_single_output_task: if perceval_effect is None: raise ValueError( "The 'perceval_effect' argument must be provided for " + @@ -828,7 +850,7 @@ def _process_term( result = None p = None - if not single_output_task: + if not is_single_output_task: if task == "probs": result = sim.probs(perceval_state) result = {tuple(k): float(v) for k, v in result.items()} @@ -842,7 +864,7 @@ def _process_term( output_types = term.cod output_shape = self._get_discopy_tensor(term).cod.inside - elif single_output_task: + elif is_single_output_task: if task == "single_prob": p = float(sim.probability(perceval_state, perceval_effect)) return_type = StateType.SINGLE_PROB @@ -859,7 +881,8 @@ def _process_term( else: raise ValueError( "Invalid task. Allowed values are" + - " 'probs', 'amps', 'single_amp', 'single_prob'." + " 'probs', 'amps', 'single_amp', 'single_prob'. " + + f"Got: {task}" ) array = np.zeros( @@ -867,7 +890,7 @@ def _process_term( dtype=float if task in ("single_prob", "probs") else complex ) - if result and not single_output_task: + if result and not is_single_output_task: configs = np.fromiter( (i for key in result for i in key), dtype=int, @@ -882,7 +905,7 @@ def _process_term( array[tuple(configs.T)] = coeffs - elif single_output_task: + elif is_single_output_task: array[0] = p else: pass diff --git a/test/test_backends.py b/test/test_backends.py index a28764a..040abe5 100644 --- a/test/test_backends.py +++ b/test/test_backends.py @@ -7,7 +7,8 @@ DiscopyBackend, EvalResult, StateType, - PermanentBackend + PermanentBackend, + PercevalEvalConfig ) import numpy as np import math @@ -135,7 +136,8 @@ class TestPercevalBackend: def test_pure_circuit(self, circuit): backend = PercevalBackend() perceval_state = pcvl.BasicState(get_state(circuit)) - result_perceval = circuit.eval(backend, perceval_state=perceval_state) + config = PercevalEvalConfig(state=perceval_state) + result_perceval = circuit.eval(backend, config=config) state = photonic.Create(*get_state(circuit)) diagram = state >> circuit @@ -146,7 +148,8 @@ def test_pure_circuit(self, circuit): result_perceval.prob_dist(), ) - result_perceval = circuit.eval(backend, perceval_state=perceval_state, task="amps") + config = PercevalEvalConfig(state=perceval_state, task="amps") + result_perceval = circuit.eval(backend, config=config) assert dict_allclose( result_quimb.amplitudes(), @@ -284,7 +287,7 @@ def test_mixed_interleaved_partial_trace_numeric_hygiene(self): assert pytest.approx(probs[(0,)], rel=1e-12) == 0.6 assert pytest.approx(probs[(1,)], rel=1e-12) == 0.4 - assert pytest.approx(sum(probs.values()), rel=1e-12) == 1.0 + assert pytest.approx(sum(probs.values()), rel=1e-6) == 1.0 def test_prob_branch_rounding_and_sum(self): P = np.array([[1/3, 1/6], @@ -421,7 +424,8 @@ def test_eval_result_not_density_matrix(self, circuit): with pytest.raises(TypeError): backend = PercevalBackend() perceval_state = pcvl.BasicState(get_state(circuit)) - result_perceval = circuit.eval(backend, perceval_state=perceval_state) + config = PercevalEvalConfig(state=perceval_state) + result_perceval = circuit.eval(backend, config=config) dm = result_perceval.density_matrix @@ -446,7 +450,8 @@ def test_abstract_backend_not_lo(self, circuit): with pytest.raises(AssertionError): backend = PercevalBackend() perceval_state = pcvl.BasicState(get_state(circuit)) - result_perceval = circuit.eval(backend, perceval_state=perceval_state) + config = PercevalEvalConfig(state=perceval_state) + result_perceval = circuit.eval(backend, config=config) def test_mixed_all_off_diagonal_mass_raises_zero_total(self): dm = np.zeros((2, 2, 2), dtype=float) @@ -515,7 +520,8 @@ def test_prob_state_compare_with_quimb(circuit): nmodes = len(circuit.dom) perceval_state = pcvl.BasicState(bs_state(nmodes)) backend = PercevalBackend() - res_pcvl = circuit.eval(backend, perceval_state=perceval_state, task="probs") + config = PercevalEvalConfig(state=perceval_state) + res_pcvl = circuit.eval(backend, config=config) ref = photonic.Create(*bs_state(nmodes)) >> circuit d_ref = ref.eval().prob_dist() @@ -529,96 +535,50 @@ def test_prob_state_compare_with_quimb(circuit): def test_prob_effect_only_errors_no_state(circuit): nmodes = len(circuit.dom) backend = PercevalBackend() + config = PercevalEvalConfig(effect=pcvl.BasicState([2] + [0]*(nmodes-1))) with pytest.raises(ValueError): _ = circuit.eval( backend, - perceval_effect=pcvl.BasicState([2] + [0]*(nmodes-1)), - task="probs" + config=config ) -@pytest.mark.parametrize("circuit", [unitary_circuit, non_unitary_circuit]) -def test_prob_state_plus_effect_coerces_to_single_prob_and_matches_quimb(circuit): - - nmodes = len(circuit.dom) - state_occ = bs_state(nmodes) - outcome = (2,) + (0,)*(nmodes-1) - - backend = PercevalBackend() - perceval_state = pcvl.BasicState(state_occ) - perceval_effect = pcvl.BasicState(outcome) - - res = circuit.eval( - backend, - perceval_state=perceval_state, - perceval_effect=perceval_effect, - task="probs" - ) - assert res.state_type is StateType.SINGLE_PROB - p_pcvl = float(res.single_prob(outcome)) - - ref = photonic.Create(*state_occ) >> circuit - p_ref = _compare_prob_for_outcome(ref, outcome) - assert math.isclose(p_pcvl, p_ref, rel_tol=1e-9, abs_tol=1e-12) - @pytest.mark.parametrize("circuit", [unitary_circuit, non_unitary_circuit]) def test_prob_diagram_only_effect_errors_no_state(circuit): nmodes = len(circuit.cod) diagram = circuit >> classical.Select(*([0]*(nmodes-1) + [2])) backend = PercevalBackend() with pytest.raises(ValueError): - _ = diagram.eval(backend, task="probs") - -@pytest.mark.parametrize("circuit", [unitary_circuit, non_unitary_circuit]) -def test_prob_diagram_state_and_effect_coerces_to_single_prob_and_matches_quimb(circuit): - nmodes = len(circuit.dom) - outcome = (0,) + (0,)*(nmodes-1) - - diagram = photonic.Create(*bs_state(nmodes)) >> circuit >> classical.Select(*outcome) - backend = PercevalBackend() - res = diagram.eval(backend, task="probs") - assert res.state_type is StateType.SINGLE_PROB - - p_pcvl = float(res.single_prob(outcome)) - - ref = photonic.Create(*bs_state(nmodes)) >> circuit - p_ref = _compare_prob_for_outcome(ref, outcome) - assert math.isclose(p_pcvl, p_ref, rel_tol=1e-9, abs_tol=1e-12) + _ = diagram.eval(backend) @pytest.mark.parametrize("circuit", [unitary_circuit, non_unitary_circuit]) def test_single_prob_errors_without_state_and_effect(circuit): backend = PercevalBackend() + config = PercevalEvalConfig(task="single_prob") with pytest.raises(ValueError): - _ = circuit.eval(backend, task="single_prob") + _ = circuit.eval(backend, config=config) @pytest.mark.parametrize("circuit", [unitary_circuit, non_unitary_circuit]) def test_single_prob_errors_with_state_only(circuit): backend = PercevalBackend() nmodes = len(circuit.dom) + config = PercevalEvalConfig(state=pcvl.BasicState(bs_state(nmodes)), task="single_prob") with pytest.raises(ValueError): - _ = circuit.eval( - backend, - perceval_state=pcvl.BasicState(bs_state(nmodes)), - task="single_prob" - ) + _ = circuit.eval(backend, config=config) + @pytest.mark.parametrize("circuit", [unitary_circuit, non_unitary_circuit]) def test_single_prob_errors_with_effect_only(circuit): backend = PercevalBackend() nmodes = len(circuit.dom) + config = PercevalEvalConfig(effect=pcvl.BasicState([2] + [0]*(nmodes-1)), task="single_prob") with pytest.raises(ValueError): - _ = circuit.eval( - backend, - perceval_effect=pcvl.BasicState([2] + [0]*(nmodes-1)), - task="single_prob" - ) + _ = circuit.eval(backend, config=config) backend = PercevalBackend() nmodes = len(circuit.dom) + config = PercevalEvalConfig(effect=pcvl.BasicState([2] + [0]*(nmodes-1)), task="single_prob") with pytest.raises(ValueError): - _ = circuit.eval( - backend, - perceval_effect=pcvl.BasicState([2] + [0]*(nmodes-1)), - task="single_prob" + _ = circuit.eval(backend, config=config ) @pytest.mark.parametrize("circuit", [unitary_circuit, non_unitary_circuit]) @@ -628,12 +588,13 @@ def test_single_prob_matches_quimb_selected_outcome(circuit): outcome = (2,) + (0,)*(nmodes-1) if nmodes >= 1 else tuple() backend = PercevalBackend() - res = circuit.eval( - backend, - perceval_state=pcvl.BasicState(state_occ), - perceval_effect=pcvl.BasicState(outcome), + config = PercevalEvalConfig( + state=pcvl.BasicState(state_occ), + effect=pcvl.BasicState(outcome), task="single_prob" ) + res = circuit.eval(backend, config=config) + assert res.state_type is StateType.SINGLE_PROB p_pcvl = float(res.single_prob(outcome)) @@ -644,8 +605,9 @@ def test_single_prob_matches_quimb_selected_outcome(circuit): @pytest.mark.parametrize("circuit", [unitary_circuit, non_unitary_circuit]) def test_amp_nothing_raises_no_state(circuit): backend = PercevalBackend() + config = PercevalEvalConfig(task="amps") with pytest.raises(ValueError): - _ = circuit.eval(backend, task="amps") # no state + _ = circuit.eval(backend, config=config) @pytest.mark.parametrize("circuit", [unitary_circuit, non_unitary_circuit]) def test_amp_state_compare_with_quimb(circuit): @@ -654,10 +616,13 @@ def test_amp_state_compare_with_quimb(circuit): state_occ = bs_state(nmodes) backend = PercevalBackend() + config = PercevalEvalConfig( + state=pcvl.BasicState(state_occ), + task="amps" + ) res_pcvl = circuit.eval( backend, - perceval_state=pcvl.BasicState(state_occ), - task="amps" + config=config ) # reference @@ -673,78 +638,45 @@ def test_amp_state_compare_with_quimb(circuit): def test_amp_effect_only_errors_no_state(circuit): backend = PercevalBackend() nmodes = len(circuit.dom) + config = PercevalEvalConfig( + effect=pcvl.BasicState([2] + [0]*(nmodes-1)), + task="amps" + ) with pytest.raises(ValueError): _ = circuit.eval( backend, - perceval_effect=pcvl.BasicState([2] + [0]*(nmodes-1)), - task="amps" + config=config ) -@pytest.mark.parametrize("circuit", [unitary_circuit, non_unitary_circuit]) -def test_amp_state_plus_effect_coerces_to_single_amp_and_matches_quimb(circuit): - - nmodes = len(circuit.dom) - state_occ = bs_state(nmodes) - outcome = (2,) + (0,)*(nmodes-1) - - backend = PercevalBackend() - res = circuit.eval( - backend, - perceval_state=pcvl.BasicState(state_occ), - perceval_effect=pcvl.BasicState(outcome), - task="amps" - ) - assert res.state_type is StateType.SINGLE_AMP - a_pcvl = res.single_amplitude(outcome) - - ref = photonic.Create(*state_occ) >> circuit - a_ref = ref.eval().amplitudes().get(outcome, 0.0) - assert a_pcvl == pytest.approx(a_ref, rel=1e-9, abs=1e-12) - -@pytest.mark.parametrize("circuit", [unitary_circuit, non_unitary_circuit]) -def test_amp_diagram_state_and_effect_coerces_to_single_amp_and_matches_quimb(circuit): - - nmodes = len(circuit.dom) - state_occ = bs_state(nmodes) - outcome = (0,) + (0,)*(nmodes-1) - - diagram = photonic.Create(*state_occ) >> circuit >> classical.Select(*outcome) - backend = PercevalBackend() - res = diagram.eval(backend, task="amps") - assert res.state_type is StateType.SINGLE_AMP - a_pcvl = res.single_amplitude(outcome) - - ref = photonic.Create(*state_occ) >> circuit - a_ref = ref.eval().amplitudes().get(outcome, 0.0) - assert a_pcvl == pytest.approx(a_ref, rel=1e-9, abs=1e-12) - @pytest.mark.parametrize("circuit", [unitary_circuit, non_unitary_circuit]) def test_single_amp_errors_without_state_and_effect(circuit): backend = PercevalBackend() + config = PercevalEvalConfig(task="single_amp") with pytest.raises(ValueError): - _ = circuit.eval(backend, task="single_amp") + _ = circuit.eval(backend, config=config) @pytest.mark.parametrize("circuit", [unitary_circuit, non_unitary_circuit]) def test_single_amp_errors_with_state_only(circuit): backend = PercevalBackend() nmodes = len(circuit.dom) + config = PercevalEvalConfig( + state=pcvl.BasicState(bs_state(nmodes)), + task="single_amp" + ) with pytest.raises(ValueError): - _ = circuit.eval( - backend, - perceval_state=pcvl.BasicState(bs_state(nmodes)), - task="single_amp" - ) + _ = circuit.eval(backend, config=config) @pytest.mark.parametrize("circuit", [unitary_circuit, non_unitary_circuit]) def test_single_amp_errors_with_effect_only(circuit): backend = PercevalBackend() nmodes = len(circuit.dom) + config = PercevalEvalConfig( + effect=pcvl.BasicState([2] + [0]*(nmodes-1)), + task="single_amp" + ) with pytest.raises(ValueError): - _ = circuit.eval( - backend, - perceval_effect=pcvl.BasicState([2] + [0]*(nmodes-1)), - task="single_amp" + _ = circuit.eval(backend, config=config ) @pytest.mark.parametrize("circuit", [unitary_circuit, non_unitary_circuit]) @@ -755,11 +687,14 @@ def test_single_amp_matches_quimb_selected_outcome(circuit): outcome = (2,) + (0,)*(nmodes-1) if nmodes >= 1 else tuple() backend = PercevalBackend() + config = PercevalEvalConfig( + state=pcvl.BasicState(state_occ), + effect=pcvl.BasicState(outcome), + task="single_amp" + ) res = circuit.eval( backend, - perceval_state=pcvl.BasicState(state_occ), - perceval_effect=pcvl.BasicState(outcome), - task="single_amp" + config=config ) assert res.state_type is StateType.SINGLE_AMP a_pcvl = res.single_amplitude(outcome) From 23330cf04e4dda1ddeb2de3abff6aa14116b2bf8 Mon Sep 17 00:00:00 2001 From: mateuszkupper Date: Tue, 21 Oct 2025 09:25:28 +0100 Subject: [PATCH 18/19] Split process_term in perceval backend into multiple methods --- optyx/core/backends.py | 206 +++++++++++++++++++++++++++++++---------- 1 file changed, 155 insertions(+), 51 deletions(-) diff --git a/optyx/core/backends.py b/optyx/core/backends.py index 6c1f557..f3bd54d 100644 --- a/optyx/core/backends.py +++ b/optyx/core/backends.py @@ -673,17 +673,20 @@ def _post_select_vacuum( self, dist, m_orig, - k_extra + k_extra, + task ): """Keep only entries where extra (ancilla) modes are all 0, then drop them.""" - if k_extra <= 0: - return dist - return { - k[:m_orig]: v - for k, v in dist.items() - if all(x == 0 for x in k[m_orig:]) - } + if task in ("probs", "amps"): + if k_extra <= 0: + return dist + return { + k[:m_orig]: v + for k, v in dist.items() + if all(x == 0 for x in k[m_orig:]) + } + return dist def _process_state(self, perceval_state): if not isinstance(perceval_state, pcvl.BasicState): @@ -736,19 +739,25 @@ def _dilate( return matrix, perceval_state - def _process_term( + def _process_io( self, term, + matrix, extra ): """ - Process a term in a sum of diagrams. + Process the input and output states/effects for the diagram. Args: term (discopy.tensor.Diagram): The term to process. + matrix (Matrix): The matrix representation of the diagram. + extra (dict): Additional arguments for the evaluation. + Returns: + perceval_state (pcvl.BasicState): The processed input state. + perceval_effect (pcvl.BasicState | None): + The processed output effect. + task (str): The Perceval task to perform. """ - matrix = self._get_matrix(term) - cfg: PercevalEvalConfig = extra.get("config", PercevalEvalConfig()) task = cfg.task state_provided = cfg.state is not None @@ -817,21 +826,22 @@ def _process_term( "when conditioning on an effect." ) - is_single_output_task = task in ("single_amp", "single_prob") - - if is_single_output_task: + if task in ("single_amp", "single_prob"): if perceval_effect is None: raise ValueError( "The 'perceval_effect' argument must be provided for " + "task 'single_amp' or 'single_prob'." ) - # pylint: disable=protected-access - if not matrix._umatrix_is_unitary(): - matrix, perceval_state = self._dilate( - matrix, perceval_state - ) + return perceval_state, perceval_effect, task + def _prepare_simulation( + self, + matrix + ): + """ + Prepare the Perceval simulator with the given matrix. + """ selections = matrix.selections sim = pcvl.Simulator(self.perceval_backend) @@ -844,53 +854,78 @@ def _process_term( ) perceval_circuit = self._umatrix_to_perceval_circuit(matrix.array) sim.set_circuit(perceval_circuit) + return sim - m_orig = len(term.dom) - k_extra = matrix.dom - m_orig - result = None - p = None + def _simulate( + self, + sim, + perceval_state, + perceval_effect, + task + ): - if not is_single_output_task: + if task in ("single_amp", "single_prob"): + if task == "single_prob": + result = float( + sim.probability(perceval_state, perceval_effect) + ) + + else: + result = complex( + sim.prob_amplitude(perceval_state, perceval_effect) + ) + else: if task == "probs": result = sim.probs(perceval_state) result = {tuple(k): float(v) for k, v in result.items()} - return_type = StateType.PROB else: sv = sim.evolve(perceval_state) result = {tuple(k): complex(v) for k, v in sv} - return_type = StateType.AMP - result = self._post_select_vacuum(result, m_orig, k_extra) - output_types = term.cod - output_shape = self._get_discopy_tensor(term).cod.inside + return result - elif is_single_output_task: - if task == "single_prob": - p = float(sim.probability(perceval_state, perceval_effect)) - return_type = StateType.SINGLE_PROB + def _get_output_params( + self, + term, + task + ): + """ + Get the output parameters for the diagram. - else: - p = complex( - sim.prob_amplitude(perceval_state, perceval_effect) - ) - return_type = StateType.SINGLE_AMP + Args: + term (discopy.tensor.Diagram): The term to process. + task (str): The Perceval task to perform. + Returns: + output_shape (tuple): The shape of the output array. + output_types (Ty | None): The output types of the diagram. + """ - output_shape = (1,) - output_types = None + if task in ("single_amp", "single_prob"): + return (1,), None + return self._get_discopy_tensor(term).cod.inside, term.cod - else: - raise ValueError( - "Invalid task. Allowed values are" + - " 'probs', 'amps', 'single_amp', 'single_prob'. " + - f"Got: {task}" - ) + def _get_array_from_result( + self, + result, + output_shape, + task + ): + """ + Get the output array from the simulation result. + Args: + result (dict | float | complex): The simulation result. + output_shape (tuple): The shape of the output array. + task (str): The Perceval task to perform. + Returns: + np.ndarray: The output array. + """ array = np.zeros( output_shape, dtype=float if task in ("single_prob", "probs") else complex ) - if result and not is_single_output_task: + if task in ("probs", "amps"): configs = np.fromiter( (i for key in result for i in key), dtype=int, @@ -905,10 +940,79 @@ def _process_term( array[tuple(configs.T)] = coeffs - elif is_single_output_task: - array[0] = p else: - pass + array[0] = result + return array + + def _process_term( + self, + term, + extra + ): + """ + Process a term in a sum of diagrams. + + Args: + term (discopy.tensor.Diagram): The term to process. + """ + matrix = self._get_matrix(term) + + perceval_state, perceval_effect, task = self._process_io( + term, + matrix, + extra + ) + + if task == "amps": + return_type = StateType.AMP + elif task == "probs": + return_type = StateType.PROB + elif task == "single_amp": + return_type = StateType.SINGLE_AMP + elif task == "single_prob": + return_type = StateType.SINGLE_PROB + else: + raise ValueError( + "Invalid task. Allowed values are" + + " 'probs', 'amps', 'single_amp', 'single_prob'. " + + f"Got: {task}" + ) + + # pylint: disable=protected-access + if not matrix._umatrix_is_unitary(): + matrix, perceval_state = self._dilate( + matrix, perceval_state + ) + + sim = self._prepare_simulation( + matrix + ) + + result = self._simulate( + sim, + perceval_state, + perceval_effect, + task + ) + + result = self._post_select_vacuum( + result, + len(term.dom), + matrix.dom - len(term.dom), + task + ) + + output_shape, output_types = self._get_output_params( + term, + task + ) + + array = self._get_array_from_result( + result, + output_shape, + task + ) + return array, output_types, return_type From 4fc7dbd9795f5ccee6402ecd12e1a4fbad42a0e6 Mon Sep 17 00:00:00 2001 From: mateuszkupper Date: Mon, 27 Oct 2025 13:44:28 +0000 Subject: [PATCH 19/19] Fix import issues --- optyx/classical.py | 2 +- optyx/photonic.py | 2 +- optyx/qubits.py | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/optyx/classical.py b/optyx/classical.py index 920f102..8955d15 100644 --- a/optyx/classical.py +++ b/optyx/classical.py @@ -141,7 +141,7 @@ diagram, path ) -from optyx import ( +from optyx.core.channel import ( bit, mode, qmode, diff --git a/optyx/photonic.py b/optyx/photonic.py index 8f3e17c..375900f 100644 --- a/optyx/photonic.py +++ b/optyx/photonic.py @@ -262,7 +262,7 @@ from optyx.classical import ClassicalFunction, DiscardMode from optyx.utils.utils import matrix_to_zw -from optyx import ( +from optyx.core.channel import ( bit, mode, qmode, diff --git a/optyx/qubits.py b/optyx/qubits.py index 1700f68..1a7deef 100644 --- a/optyx/qubits.py +++ b/optyx/qubits.py @@ -311,7 +311,7 @@ diagram, zx ) -from optyx import ( +from optyx.core.channel import ( bit, qubit, Measure as MeasureChannel,