diff --git a/CHANGELOG.md b/CHANGELOG.md index 625cb21ea..59c380dac 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -12,6 +12,10 @@ releases may include breaking changes. ### Added +- added custom one-site matrix observables for non-qubit local dimensions + + , including named position observables ([#497]) ([**@linusschulte**]) + - added direct MPO Process Tensor construction and temporal entanglement ([#508]) ([**@aaronleesander**]) @@ -229,6 +233,7 @@ changelogs._ +[#497]: https://github.com/munich-quantum-toolkit/yaqs/pull/497 [#519]: https://github.com/munich-quantum-toolkit/yaqs/pull/519 [#518]: https://github.com/munich-quantum-toolkit/yaqs/pull/518 [#516]: https://github.com/munich-quantum-toolkit/yaqs/pull/516 diff --git a/docs/examples/simulation_parameters.md b/docs/examples/simulation_parameters.md index e1c75ad71..8df8c166d 100644 --- a/docs/examples/simulation_parameters.md +++ b/docs/examples/simulation_parameters.md @@ -33,6 +33,7 @@ do not import gate classes for standard measurements. | `"x"`, `"y"`, `"z"` | Single-qubit Pauli operators | `Observable("z", sites=0)` | | `"h"`, `"s"`, `"t"`, `"rx"`, … | Other single-qubit gates from the built-in library | `Observable("h", sites=0)` | | `"xx"`, `"yy"`, `"zz"` | Two-qubit Pauli strings | `Observable("zz", sites=[0, 1])` | +| `"position"` | Position operator for a supplied local position basis | `Observable("position", 0, positions=grid)` | | `"entropy"` | Bipartite entanglement entropy across a cut | `Observable("entropy", sites=cut)` | | `"schmidt_spectrum"` | Schmidt spectrum across a cut | `Observable("schmidt_spectrum", sites=cut)` | | bitstring / `"pvm"` | Projection-valued measurement onto a computational basis state | see {doc}`circuit_observables` | @@ -40,6 +41,10 @@ do not import gate classes for standard measurements. For custom unitaries and circuit gates, use {doc}`custom_gates` — those workflows still use `GateLibrary` or Qiskit circuits directly. +Named observables that require configuration accept keyword-only factory +arguments. Missing or unknown arguments raise `TypeError`, so misspelled +parameters are not silently ignored. + ## Start with a preset You do **not** need to tune every numerical knob before running a simulation. diff --git a/docs/examples/trapped_ion.md b/docs/examples/trapped_ion.md index deb9ed897..c4c7c71e2 100644 --- a/docs/examples/trapped_ion.md +++ b/docs/examples/trapped_ion.md @@ -30,7 +30,7 @@ reaches the opposite turning point. ```{code-cell} ipython3 import numpy as np -from mqt.yaqs import Hamiltonian, MPO, State +from mqt.yaqs import Hamiltonian, MPO, Observable, State omega = 1.0 initial_displacement = 1.0 @@ -44,6 +44,7 @@ initial_grid_state /= np.linalg.norm(initial_grid_state) hamiltonian = Hamiltonian.from_mpo(MPO.trapped_ion(positions, masses=[1.0], omega=omega)) state = State(length=1, vector=initial_grid_state, physical_dimensions=[grid_dim]) +position_observable = Observable("position", 0, positions=positions) ``` ## 2. Noiseless evolution to $T/2$ @@ -52,7 +53,7 @@ state = State(length=1, vector=initial_grid_state, physical_dimensions=[grid_dim from mqt.yaqs import AnalogSimParams, Simulator params = AnalogSimParams( - observables=[], + observables=[position_observable], elapsed_time=half_period, dt=half_period / 16, max_bond_dim=None, @@ -60,16 +61,18 @@ params = AnalogSimParams( krylov_tol=1e-12, preset="exact", get_state=True, - sample_timesteps=False, + sample_timesteps=True, ) result = Simulator(show_progress=False).run(state, hamiltonian, params) final_state = result.output_state.vector -final_x = float(np.sum(positions * np.abs(final_state) ** 2)) +position_expectation = np.real(result.expectation_values[0]) +final_x = float(position_expectation[-1]) ``` -The final $\langle x\rangle$ is close to $-x_0$ but not exact because the -simulation uses a finite grid and a finite-difference kinetic operator. +The position observable is a custom one-site matrix on the grid basis. The final +$\langle x\rangle$ is close to $-x_0$ but not exact because the simulation uses +a finite grid and a finite-difference kinetic operator. ```{code-cell} ipython3 print(f"Initial = {initial_displacement:.6f}") @@ -77,7 +80,7 @@ print(f"Final at T/2 = {final_x:.6f}") print(f"Continuum target = {-initial_displacement:.6f}") ``` -## 3. Wavepacket at $t=0$ and $t=T/2$ +## 3. Wavepacket over time ```{code-cell} ipython3 --- @@ -88,24 +91,27 @@ mystnb: --- import matplotlib.pyplot as plt -prob_initial = np.abs(initial_grid_state) ** 2 -prob_final = np.abs(final_state) ** 2 - -fig, axes = plt.subplots(1, 2, figsize=(8, 3.2), layout="constrained", sharey=True) -axes[0].fill_between(positions, prob_initial, alpha=0.35, color="tab:blue") -axes[0].plot(positions, prob_initial, color="tab:blue", lw=1.5) -axes[0].set_title(r"$t = 0$") -axes[0].set_xlabel(r"$x$") -axes[0].set_ylabel(r"$|\psi(x)|^2$") -axes[0].grid(alpha=0.3) - -axes[1].fill_between(positions, prob_final, alpha=0.35, color="tab:orange") -axes[1].plot(positions, prob_final, color="tab:orange", lw=1.5) -axes[1].set_title(rf"$t = T/2$") -axes[1].set_xlabel(r"$x$") -axes[1].grid(alpha=0.3) - -fig.suptitle("Harmonic wavepacket reflection on a position grid") +dense_hamiltonian = hamiltonian.to_matrix() +eigenvalues, eigenvectors = np.linalg.eigh(dense_hamiltonian) +coefficients = eigenvectors.conj().T @ initial_grid_state +phases = np.exp(-1j * eigenvalues[:, None] * params.times[None, :]) +states = eigenvectors @ (coefficients[:, None] * phases) +probability_density = np.abs(states) ** 2 + +fig, ax = plt.subplots(figsize=(7.2, 3.6), layout="constrained") +image = ax.imshow( + probability_density, + aspect="auto", + origin="lower", + extent=(params.times[0], params.times[-1], positions[0], positions[-1]), + cmap="viridis", +) +ax.plot(params.times, position_expectation, color="white", lw=1.4, label=r"$\langle x\rangle$") +ax.set_xlabel(r"$t$") +ax.set_ylabel(r"$x$") +ax.set_title("Position-grid wavepacket density") +ax.legend(loc="upper right") +fig.colorbar(image, ax=ax, label=r"$|\psi(x,t)|^2$") plt.show() ``` diff --git a/src/mqt/yaqs/core/data_structures/mps.py b/src/mqt/yaqs/core/data_structures/mps.py index 04b4a8481..43b831291 100644 --- a/src/mqt/yaqs/core/data_structures/mps.py +++ b/src/mqt/yaqs/core/data_structures/mps.py @@ -968,32 +968,34 @@ def local_expect(self, operator: Observable, sites: int | list[int]) -> np.compl Returns: np.complex128: The computed expectation value (typically, its real part is of interest). + Raises: + ValueError: If the observable is not supported or its matrix shape does not match the target site. + Notes: A deep copy of the state is used to prevent modifications to the original MPS. Requires :meth:`check_covers_sites` to hold for ``sites``; prefer :meth:`expect` for gauge-safe evaluation. """ temp_state = copy.deepcopy(self) - if operator.gate.matrix.shape[0] == 2: # Local observable - i = None - if isinstance(sites, list): - i = sites[0] - elif isinstance(sites, int): - i = sites - - if isinstance(operator.sites, list): - assert operator.sites[0] == i, f"Operator sites mismatch {operator.sites[0]}, {i}" - elif isinstance(operator.sites, int): - assert operator.sites == i, f"Operator sites mismatch {operator.sites}, {i}" + sites_list = [sites] if isinstance(sites, int) else list(sites) + operator_sites = [operator.sites] if isinstance(operator.sites, int) else list(operator.sites) - assert i is not None, f"Invalid type for 'sites': expected int or list[int], got {type(sites).__name__}" + if operator.gate.interaction == 1: + assert len(sites_list) == 1, f"One-site observable requires one site, got {sites_list}." + assert operator_sites == sites_list, f"Operator sites mismatch {operator_sites}, {sites_list}" + i = sites_list[0] a = temp_state.tensors[i] - temp_state.tensors[i] = oe.contract("ab, bcd->acd", operator.gate.matrix, a) + local_dim = a.shape[0] + matrix = np.asarray(operator.gate.matrix, dtype=np.complex128) + if matrix.shape != (local_dim, local_dim): + msg = f"Local observable matrix shape {matrix.shape} does not match site {i} dimension {local_dim}." + raise ValueError(msg) + temp_state.tensors[i] = oe.contract("ab, bcd->acd", matrix, a) - elif operator.gate.matrix.shape[0] == 4: # Two-site correlator + elif operator.gate.interaction == 2: assert isinstance(sites, list) assert isinstance(operator.sites, list) - i, j = sites + i, j = sites_list assert operator.sites[0] == i, "Observable sites mismatch" assert operator.sites[1] == j, "Observable sites mismatch" @@ -1005,6 +1007,10 @@ def local_expect(self, operator: Observable, sites: int | list[int]) -> np.compl b = temp_state.tensors[j] d_i, left, _ = a.shape d_j, _, right = b.shape + matrix = np.asarray(operator.gate.matrix, dtype=np.complex128) + if matrix.shape != (d_i * d_j, d_i * d_j): + msg = f"Two-site observable matrix shape {matrix.shape} does not match site dimensions {d_i} and {d_j}." + raise ValueError(msg) # 1) merge A,B into theta of shape (l, d_i*d_j, r) theta = np.tensordot(a, b, axes=(2, 1)) # (d_i, l, d_j, r) @@ -1012,7 +1018,7 @@ def local_expect(self, operator: Observable, sites: int | list[int]) -> np.compl theta = theta.reshape(left, d_i * d_j, right) # (l, d_i*d_j, r) # 2) apply operator on the combined phys index - theta = oe.contract("ab, cbd->cad", operator.gate.matrix, theta) # (l, d_i*d_j, r) + theta = oe.contract("ab, cbd->cad", matrix, theta) # (l, d_i*d_j, r) theta = theta.reshape(left, d_i, d_j, right) # back to (l, d_i, d_j, r) # 3) split via SVD @@ -1030,6 +1036,9 @@ def local_expect(self, operator: Observable, sites: int | list[int]) -> np.compl temp_state.tensors[i] = a_new temp_state.tensors[j] = b_new + else: + msg = "Local observable must be one-site or nearest-neighbor two-site." + raise ValueError(msg) return self.scalar_product(temp_state, sites) @@ -1094,26 +1103,41 @@ def bubble_swaps_backward(state: MPS) -> None: for i in reversed(range(state.length - 2)): apply_two_site_nn_inplace(state, i, sw) - sites = [observable.sites] if isinstance(observable.sites, int) else observable.sites + sites = [observable.sites] if isinstance(observable.sites, int) else list(observable.sites) - if observable.gate.matrix.shape[0] == 2: + if observable.gate.interaction == 1: + if len(sites) != 1: + msg = f"One-site local observable requires one site, got {sites}." + raise ValueError(msg) site = sites[0] - self.tensors[site] = oe.contract("ab, bcd->acd", observable.gate.matrix, self.tensors[site]) + local_dim = self.tensors[site].shape[0] + matrix = np.asarray(observable.gate.matrix, dtype=np.complex128) + if matrix.shape != (local_dim, local_dim): + msg = f"Local observable matrix shape {matrix.shape} does not match site {site} dimension {local_dim}." + raise ValueError(msg) + self.tensors[site] = oe.contract("ab, bcd->acd", matrix, self.tensors[site]) return - if observable.gate.matrix.shape[0] == 4: + if observable.gate.interaction == 2: + if len(sites) != 2: + msg = f"Two-site local observable requires two sites, got {sites}." + raise ValueError(msg) i, j = int(sites[0]), int(sites[1]) length = self.length + mat = np.asarray(observable.gate.matrix, dtype=np.complex128) + d_i = self.tensors[i].shape[0] + d_j = self.tensors[j].shape[0] + if mat.shape != (d_i * d_j, d_i * d_j): + msg = f"Two-site observable matrix shape {mat.shape} does not match site dimensions {d_i} and {d_j}." + raise ValueError(msg) if length == 2: if i == length - 1 and j == 0: - mat = np.asarray(observable.gate.matrix, dtype=np.complex128) g_merged = permuted_periodic_wrap(mat) apply_two_site_nn_inplace(self, 0, g_merged) return i, j = min(i, j), max(i, j) elif (i == length - 1 and j == 0) or (i == 0 and j == length - 1): - mat = np.asarray(observable.gate.matrix, dtype=np.complex128) bubble_swaps_forward(self) g_merged = permuted_periodic_wrap(mat) apply_two_site_nn_inplace(self, length - 2, g_merged) @@ -1124,7 +1148,7 @@ def bubble_swaps_backward(state: MPS) -> None: msg = "Only nearest-neighbor two-site observables are currently implemented." raise ValueError(msg) - apply_two_site_nn_inplace(self, i, np.asarray(observable.gate.matrix, dtype=np.complex128)) + apply_two_site_nn_inplace(self, i, mat) return msg = "Local observable must be one-site or nearest-neighbor two-site." diff --git a/src/mqt/yaqs/core/data_structures/simulation_parameters.py b/src/mqt/yaqs/core/data_structures/simulation_parameters.py index 95f9b794f..9398f978d 100644 --- a/src/mqt/yaqs/core/data_structures/simulation_parameters.py +++ b/src/mqt/yaqs/core/data_structures/simulation_parameters.py @@ -23,10 +23,10 @@ import numpy as np -from mqt.yaqs.core.libraries.gate_library import GateLibrary +from mqt.yaqs.core.libraries.gate_library import BaseGate, GateLibrary if TYPE_CHECKING: - from mqt.yaqs.core.libraries.gate_library import BaseGate + from numpy.typing import ArrayLike SimulationPreset = Literal["fast", "balanced", "accurate", "exact"] GateMode = Literal["tdvp", "full-tdvp", "swaps", "mpo"] @@ -233,31 +233,49 @@ class Observable: sites: The site or site indices on which this observable is measured. """ - def __init__(self, gate: BaseGate | str, sites: int | list[int] | None = None) -> None: + def __init__( + self, + gate: BaseGate | str | ArrayLike, + sites: int | list[int] | None = None, + **gate_kwargs: object, + ) -> None: """Initializes an Observable instance. Args: - gate: The gate that will act as the observable. + gate: The gate or one-site local matrix that will act as the observable. sites: The qubit or site indices on which this observable is measured. + **gate_kwargs: Keyword-only arguments for a named gate or observable factory. + + Raises: + TypeError: If factory arguments are missing, unexpected, or supplied for a gate instance or matrix. """ if isinstance(gate, str): - if gate == "entropy": - gate = GateLibrary.entropy() - elif gate == "schmidt_spectrum": - gate = GateLibrary.schmidt_spectrum() - elif gate == "pvm": - gate = GateLibrary.pvm(gate) + if gate == "pvm": + if gate_kwargs: + msg = "'pvm' does not accept observable parameters." + raise TypeError(msg) + resolved_gate = GateLibrary.pvm(gate) elif hasattr(GateLibrary, gate): attr = getattr(GateLibrary, gate) - try: - gate = attr() - except TypeError: - gate = GateLibrary.pvm(gate) + resolved_gate = attr(**gate_kwargs) else: - gate = GateLibrary.pvm(gate) - assert hasattr(GateLibrary, gate.name), f"Observable {gate.name} not found in GateLibrary." - self.gate = copy.deepcopy(gate) - if gate.name != "pvm": + if gate_kwargs: + msg = f"Unknown observable {gate!r} does not accept observable parameters." + raise TypeError(msg) + resolved_gate = GateLibrary.pvm(gate) + elif isinstance(gate, BaseGate): + if gate_kwargs: + msg = "Observable parameters are only supported for named observables." + raise TypeError(msg) + resolved_gate = gate + else: + if gate_kwargs: + msg = "Observable parameters are only supported for named observables." + raise TypeError(msg) + resolved_gate = GateLibrary.local(gate) + assert hasattr(GateLibrary, resolved_gate.name), f"Observable {resolved_gate.name} not found in GateLibrary." + self.gate: BaseGate = copy.deepcopy(resolved_gate) + if resolved_gate.name != "pvm": assert sites is not None self.sites = sites self.gate.set_sites(self.sites) diff --git a/src/mqt/yaqs/core/libraries/gate_library.py b/src/mqt/yaqs/core/libraries/gate_library.py index 0cba89819..ea3de5e2b 100644 --- a/src/mqt/yaqs/core/libraries/gate_library.py +++ b/src/mqt/yaqs/core/libraries/gate_library.py @@ -1650,6 +1650,65 @@ def __init__(self, bitstring: str) -> None: super().__init__(mat) +class LocalOperator(BaseGate): + """Custom one-site operator for arbitrary local Hilbert-space dimensions. + + This gate is intended for observables such as position-grid operators on + qudits or oscillator truncations. Unlike :class:`BaseGate`, it does not + interpret the matrix dimension as a qubit interaction count. + """ + + name = "local" + + def __init__(self, matrix: ArrayLike) -> None: + """Create a one-site local operator. + + Args: + matrix: Square matrix acting on one local site. + + Raises: + ValueError: If ``matrix`` is not a square two-dimensional array. + """ + mat = np.asarray(matrix, dtype=np.complex128) + if mat.ndim != 2: + msg = "Local operator matrix must be a 2-D array." + raise ValueError(msg) + if mat.shape[0] != mat.shape[1]: + msg = "Local operator matrix must be square." + raise ValueError(msg) + self.matrix = mat + self.tensor = mat + self.interaction = 1 + + +class Position(LocalOperator): + """One-site position operator for a supplied position basis.""" + + name = "position" + + def __init__(self, *, positions: ArrayLike) -> None: + """Create a position operator that is diagonal in the supplied basis. + + Args: + positions: One-dimensional position values defining the local basis. + + Raises: + ValueError: If ``positions`` is complex or not a non-empty, finite one-dimensional array. + """ + position_values = np.asarray(positions) + if np.iscomplexobj(position_values): + msg = "positions must contain only real values." + raise ValueError(msg) + position_values = np.asarray(position_values, dtype=np.float64) + if position_values.ndim != 1 or position_values.size == 0: + msg = "positions must be a non-empty one-dimensional array." + raise ValueError(msg) + if not np.all(np.isfinite(position_values)): + msg = "positions must contain only finite values." + raise ValueError(msg) + super().__init__(np.diag(position_values)) + + class Entropy(BaseGate): """Meta-observable for bipartite entanglement entropy across a cut. @@ -1761,6 +1820,8 @@ class GateLibrary: p0: Class for projector ``|0⟩⟨0|``. p1: Class for projector ``|1⟩⟨1|``. pvm: Class for projection-valued measurement onto a given bitstring. + local: Class for arbitrary one-site local operators. + position: Class for a one-site position operator in a supplied position basis. entropy: Class representing a request for bipartite entanglement entropy across a cut. schmidt_spectrum: Class representing a request for the Schmidt spectrum across a cut. @@ -1811,6 +1872,8 @@ class GateLibrary: p0 = P0 p1 = P1 pvm = PVM + local = LocalOperator + position = Position entropy = Entropy schmidt_spectrum = SchmidtSpectrum diff --git a/tests/core/data_structures/test_mpo.py b/tests/core/data_structures/test_mpo.py index 309cce894..9aea6cdc6 100644 --- a/tests/core/data_structures/test_mpo.py +++ b/tests/core/data_structures/test_mpo.py @@ -513,6 +513,47 @@ def test_trapped_ion_coulomb_truncation() -> None: np.testing.assert_allclose(truncated_coulomb, expected_rank_2, atol=1e-12) +def test_trapped_ion_one_ion_position_observable_centers_on_trap() -> None: + """The one-ion ground-state position expectation follows the static trap center.""" + positions = np.linspace(-2.0, 2.0, 9, dtype=np.float64) + trap_center = 0.4 + mpo = MPO.trapped_ion(positions, [1.0], omega=1.0, trap_center=trap_center) + _energy, eigenvectors = np.linalg.eigh(mpo.to_matrix()) + ground_state = eigenvectors[:, 0] + mps = MPS( + length=1, + tensors=[ground_state.reshape(positions.size, 1, 1)], + physical_dimensions=[positions.size], + ) + + position = Observable("position", 0, positions=positions) + + np.testing.assert_allclose(mps.expect(position), trap_center, atol=6e-2) + + +def test_trapped_ion_two_ion_coulomb_increases_ground_state_separation() -> None: + """Softened Coulomb repulsion increases the two-ion ground-state separation.""" + positions = np.linspace(-3.0, 3.0, 9, dtype=np.float64) + grid_dim = positions.size + separation = np.abs(positions[:, None] - positions[None, :]) + + def ground_state_separation(coulomb_strength: float) -> float: + mpo = MPO.trapped_ion( + positions, + [1.0, 1.0], + omega=0.6, + coulomb_strength=coulomb_strength, + ) + _energy, eigenvectors = np.linalg.eigh(mpo.to_matrix()) + probabilities = np.abs(eigenvectors[:, 0].reshape(grid_dim, grid_dim)) ** 2 + return float(np.sum(separation * probabilities)) + + uncoupled_separation = ground_state_separation(0.0) + repulsive_separation = ground_state_separation(0.8) + + assert repulsive_separation > uncoupled_separation + 0.2 + + @pytest.mark.parametrize( ("kwargs", "match"), [ diff --git a/tests/core/data_structures/test_mps.py b/tests/core/data_structures/test_mps.py index c05d82d56..1c3472567 100644 --- a/tests/core/data_structures/test_mps.py +++ b/tests/core/data_structures/test_mps.py @@ -441,6 +441,89 @@ def test_local_expect_x_on_plus_state() -> None: np.testing.assert_allclose(val, 1.0, atol=1e-12) +def test_non_qubit_local_expectation_from_matrix_observable() -> None: + """A matrix observable can be measured on a non-qubit local site.""" + amplitudes = np.sqrt(np.array([0.2, 0.3, 0.5], dtype=np.float64)).astype(np.complex128) + tensor0 = amplitudes.reshape(3, 1, 1) + tensor1 = np.array([1.0, 0.0], dtype=np.complex128).reshape(2, 1, 1) + psi_mps = MPS(length=2, tensors=[tensor0, tensor1], physical_dimensions=[3, 2]) + position = np.diag(np.array([-1.0, 0.5, 2.0], dtype=np.float64)) + + val = psi_mps.expect(Observable(position, 0)) + + expected = 0.2 * -1.0 + 0.3 * 0.5 + 0.5 * 2.0 + np.testing.assert_allclose(val, expected, atol=1e-12) + + +def test_four_level_local_observable_is_not_treated_as_two_site() -> None: + """A ``4 x 4`` matrix observable acts on one four-level site when wrapped as local.""" + amplitudes = np.array([0.5, 0.5j, -0.5, 0.5], dtype=np.complex128) + psi_mps = MPS(length=1, tensors=[amplitudes.reshape(4, 1, 1)], physical_dimensions=[4]) + matrix = np.diag(np.array([0.0, 1.0, 2.0, 3.0], dtype=np.float64)) + + val = psi_mps.expect(Observable(matrix, 0)) + + np.testing.assert_allclose(val, 1.5, atol=1e-12) + + +def test_local_observable_dimension_mismatch_raises() -> None: + """One-site matrix observables must match the measured site's local dimension.""" + psi_mps = MPS(length=1, physical_dimensions=[3], state="zeros") + + with pytest.raises(ValueError, match="does not match site 0 dimension 3"): + psi_mps.expect(Observable(np.eye(2), 0)) + + +def test_two_site_local_observable_dimension_mismatch_raises() -> None: + """Two-site observables must match the product of both local dimensions.""" + psi_mps = MPS(length=2, physical_dimensions=[3, 2], state="zeros") + observable = Observable(BaseGate(np.eye(4)), [0, 1]) + + with pytest.raises(ValueError, match="does not match site dimensions 3 and 2"): + psi_mps.local_expect(observable, [0, 1]) + + +def test_local_expect_rejects_observables_with_unsupported_interaction() -> None: + """Local expectation values support at most two-site observables.""" + psi_mps = MPS(length=3, state="zeros") + observable = Observable(BaseGate(np.eye(8)), [0, 1, 2]) + + with pytest.raises(ValueError, match="Local observable must be one-site or nearest-neighbor two-site"): + psi_mps.local_expect(observable, [0, 1, 2]) + + +def test_apply_local_rejects_one_site_observable_with_multiple_sites() -> None: + """One-site observables retain a defensive site-count check during application.""" + psi_mps = MPS(length=2, state="zeros") + observable = Observable(X(), 0) + observable.sites = [0, 1] + + with pytest.raises(ValueError, match=r"One-site local observable requires one site, got \[0, 1\]"): + psi_mps.apply_local(observable) + + +def test_apply_local_rejects_mismatched_one_site_dimension() -> None: + """One-site observables must match the local dimension during direct application.""" + psi_mps = MPS(length=1, physical_dimensions=[3], state="zeros") + + with pytest.raises(ValueError, match="does not match site 0 dimension 3"): + psi_mps.apply_local(Observable(np.eye(2), 0)) + + +def test_apply_local_rejects_invalid_two_site_observable_shape() -> None: + """Two-site observables must supply two sites and match both local dimensions.""" + psi_mps = MPS(length=2, physical_dimensions=[3, 2], state="zeros") + observable = Observable(BaseGate(np.eye(4)), [0, 1]) + + observable.sites = [0] + with pytest.raises(ValueError, match=r"requires two sites, got \[0\]"): + psi_mps.apply_local(observable) + + observable.sites = [0, 1] + with pytest.raises(ValueError, match="does not match site dimensions 3 and 2"): + psi_mps.apply_local(observable) + + def test_mps_apply_local_l2_periodic_wrap_matches_permuted_nn() -> None: """For ``L == 2``, wrap-ordered and permuted NN applications must agree.""" length = 2 diff --git a/tests/core/data_structures/test_simulation_parameters.py b/tests/core/data_structures/test_simulation_parameters.py index a049098ef..73db2262d 100644 --- a/tests/core/data_structures/test_simulation_parameters.py +++ b/tests/core/data_structures/test_simulation_parameters.py @@ -35,7 +35,7 @@ Observable, _validate_tdvp_sweeps, ) -from mqt.yaqs.core.libraries.gate_library import GateLibrary, X +from mqt.yaqs.core.libraries.gate_library import BaseGate, GateLibrary, X from mqt.yaqs.core.methods.tdvp import primitives as tdvp_primitives if TYPE_CHECKING: @@ -61,6 +61,85 @@ def test_observable_creation_valid() -> None: assert obs.sites == site +def test_observable_accepts_custom_local_matrix() -> None: + """Observable accepts square matrices as arbitrary one-site local operators.""" + matrix = np.diag(np.array([-1.0, 0.25, 2.0])) + + obs = Observable(matrix, 0) + + assert obs.gate.name == "local" + assert obs.gate.interaction == 1 + np.testing.assert_allclose(obs.gate.matrix, matrix) + assert obs.sites == 0 + + +def test_observable_accepts_named_position_operator() -> None: + """Position observables build a diagonal local operator from the supplied basis.""" + positions = np.array([-1.5, 0.0, 2.5]) + + obs = Observable("position", 1, positions=positions) + + assert obs.gate.name == "position" + assert obs.gate.interaction == 1 + np.testing.assert_allclose(obs.gate.matrix, np.diag(positions)) + assert obs.sites == 1 + + +def test_position_observable_requires_positions() -> None: + """Position observables require their basis values as a keyword argument.""" + with pytest.raises(TypeError, match="required keyword-only argument: 'positions'"): + Observable("position", 0) + + +@pytest.mark.parametrize( + ("gate", "kwargs", "match"), + [ + ("position", {"position_values": [0.0, 1.0]}, "unexpected keyword argument 'position_values'"), + ("z", {"positions": [0.0, 1.0]}, "unexpected keyword argument 'positions'"), + ], +) +def test_named_observable_rejects_unexpected_parameters( + gate: str, + kwargs: dict[str, object], + match: str, +) -> None: + """Named observable factories reject misspelled or inapplicable parameters.""" + with pytest.raises(TypeError, match=match): + Observable(gate, 0, **kwargs) + + +def test_matrix_observable_rejects_named_parameters() -> None: + """Factory parameters cannot be supplied with a matrix observable.""" + with pytest.raises(TypeError, match="only supported for named observables"): + Observable(np.eye(2), 0, positions=[0.0, 1.0]) + + +def test_observable_rejects_parameters_without_a_matching_factory() -> None: + """Only recognized named factories accept additional observable parameters.""" + with pytest.raises(TypeError, match="'pvm' does not accept observable parameters"): + Observable("pvm", bitstring="0") + + with pytest.raises(TypeError, match="Unknown observable 'unknown'"): + Observable("unknown", parameter=1) + + with pytest.raises(TypeError, match="only supported for named observables"): + Observable(BaseGate(np.eye(2)), 0, parameter=1) + + +@pytest.mark.parametrize("positions", [np.array([]), np.array([0.0, np.nan]), np.array([0.0, 1.0j])]) +def test_position_observable_rejects_invalid_positions(positions: np.ndarray) -> None: + """Position bases must be non-empty, finite, and real.""" + with pytest.raises(ValueError, match="positions must"): + Observable("position", 0, positions=positions) + + +@pytest.mark.parametrize("matrix", [np.ones(3), np.ones((2, 3))]) +def test_observable_rejects_invalid_custom_local_matrix(matrix: np.ndarray) -> None: + """Matrix observables must be two-dimensional and square.""" + with pytest.raises(ValueError, match="Local operator matrix"): + Observable(matrix, 0) + + def test_analog_simparams_basic() -> None: """Test that AnalogSimParams is initialized with correct parameters. @@ -392,6 +471,15 @@ def test_observable_from_string_falls_back_to_pvm() -> None: assert np.allclose(obs.gate.matrix, np.eye(2)) +def test_observable_from_explicit_pvm_string_falls_back_to_pvm() -> None: + """The explicit PVM name retains its historical string-resolution behavior.""" + obs = Observable("pvm") + + assert obs.gate.name == "pvm" + assert hasattr(obs.gate, "bitstring") + assert obs.gate.bitstring == "pvm" + + def test_observable_from_gate_instance_keeps_gate_and_sites_int() -> None: """Passing a concrete BaseGate instance should be preserved and sites can be an int.""" x_gate = GateLibrary.x()