diff --git a/CHANGELOG.md b/CHANGELOG.md index 5699c572..3a7ac4ff 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -9,6 +9,10 @@ This project adheres to [Semantic Versioning], with the exception that minor rel ## [Unreleased] +### Added + +- added direct MPO Process Tensor construction and temporal entanglement ([#508]) ([**@aaronleesander**]) + ### Fixed - fixed TJM jump selection to align processes with site-sweep probabilities ([#506]) ([**@Pouri96**]) @@ -154,6 +158,7 @@ _📚 Refer to the [GitHub Release Notes](https://github.com/munich-quantum-tool +[#508]: https://github.com/munich-quantum-toolkit/yaqs/pull/508 [#506]: https://github.com/munich-quantum-toolkit/yaqs/pull/506 [#288]: https://github.com/munich-quantum-toolkit/yaqs/pull/288 [#482]: https://github.com/munich-quantum-toolkit/yaqs/pull/482 diff --git a/docs/examples/characterization.md b/docs/examples/characterization.md index 2b9f8f2b..1db52929 100644 --- a/docs/examples/characterization.md +++ b/docs/examples/characterization.md @@ -17,7 +17,9 @@ mystnb: Open quantum systems in YAQS couple a **probe qubit** (site 0) to an **environment** simulated by the remaining chain. **Environmental memory** measures how long the environment keeps past control and measurement choices relevant for future probe responses, evaluated at a temporal cut $c$ in a sequence of interventions. -Use {meth}`~mqt.yaqs.memory_characterizer.MemoryCharacterizer.characterize` to probe the process: assemble the weighted **response matrix** $\widetilde{V}(c)$, then read $S_V(c)$, $R(c)=\exp(S_V(c))$, and the mode spectrum. +Use {meth}`~mqt.yaqs.memory_characterizer.MemoryCharacterizer.characterize` to probe **operational memory**: assemble the weighted **response matrix** $\widetilde{V}(c)$, then read $S_V(c)$, $R(c)=\exp(S_V(c))$, and the mode spectrum. + +Alternatively, build a process tensor (default: direct MPO) and call `compute_temporal_entropy` for **temporal entanglement** $S_{PT}(c)$ of the multi-time process itself — a distinct quantity from $S_V(c)$. For fast dynamics under control sequences, see {doc}`memory_surrogate`. ## Setup @@ -209,6 +211,50 @@ fig.tight_layout() `MemoryCharacterizer(representation="auto")` mirrors `Simulator`: `"vector"` selects MCWF, `"mps"` selects TJM for the **environment** chain. With `"auto"`, MCWF is used when `hamiltonian.length <= vector_max_qubits` (default 10). +## Temporal entanglement from a process tensor + +Operational memory ($S_V$) comes from probe responses. +**Temporal entanglement** $S_{PT}(c)$ is computed directly from a process tensor at the same causal cut. +By default, `build_process_tensor` uses direct MPO construction (`return_type="mpo"`). +Pass `return_type="dense"` for exhaustive tomography (required for `noise_model`): + +```{code-cell} ipython3 +k = 3 +cut_pt = 2 +timesteps = [0.1] * (k + 1) + +pt_mpo = mc.build_process_tensor(ham, params, timesteps=timesteps) +pt_dense = mc.build_process_tensor( + ham, + params, + timesteps=timesteps, + return_type="dense", +) + +s_mpo = pt_mpo.compute_temporal_entropy(cut_pt) +s_dense = pt_dense.compute_temporal_entropy(cut_pt) +print( + f"S_PT(c={cut_pt}): mpo={s_mpo['entropy']:.4f}, " + f"dense={s_dense['entropy']:.4f}, schmidt_rank={s_mpo['schmidt_rank']}" +) + +# Same process tensor also supports operational memory via characterize: +pt_result = mc.characterize( + pt_mpo, + cut=cut_pt, + num_interventions=k, + n_pasts=6, + n_futures=6, + rng=np.random.default_rng(7), +) +print(f"S_V(c={cut_pt}) from process-tensor probes: {pt_result.entropy(cut_pt):.4f}") +``` + +Dense and uncapped MPO construction (`max_bond_dim=None`) agree on $S_{PT}$ for small $k$. +Use `return_type="dense"` when you need noise. The default `max_bond_dim=64` keeps direct +construction scalable; pass `max_bond_dim=None` for an exact noiseless MPO. +`characterize(pt, ...)` still builds $S_V$ from probe responses (native MPO `evaluate_probes`, without densifying for the V-matrix path). + ## Related topics - {doc}`quickstart` — minimal characterize and surrogate predict snippets diff --git a/docs/examples/memory_surrogate.md b/docs/examples/memory_surrogate.md index 9e8b235c..0e7aba83 100644 --- a/docs/examples/memory_surrogate.md +++ b/docs/examples/memory_surrogate.md @@ -196,21 +196,21 @@ Extend the per-leg list when `num_interventions > 1` to probe multi-step sequenc ## Validate against exact references -Build exhaustive process tensors for the same schedule. +Build process tensors for the same schedule. +By default, `build_process_tensor` returns an MPO from direct construction (noiseless). +Pass `return_type="dense"` for exhaustive tomography. For process tensors, `rho0` in `predict` must match `pt.initial_rho` (the site-0 state after the initial leg of the reference schedule). **Dense** and **MPO** implementations should agree on identical interventions. Compare all three backends on a **stochastic sequence drawn from the training style** (`measure_prepare` here): pass a **fresh** `np.random.default_rng(seed)` to each `predict` call (reusing one RNG object advances its state between calls). ```{code-cell} ipython3 +pt_mpo = mc.build_process_tensor(ham, params, timesteps=timesteps) pt_dense = mc.build_process_tensor( ham, params, timesteps=timesteps, return_type="dense", num_trajectories=48, ) -pt_mpo = mc.build_process_tensor( - ham, params, timesteps=timesteps, return_type="mpo", num_trajectories=48, -) -rho0 = pt_dense.initial_rho +rho0 = pt_mpo.initial_rho compare_seed = 7 rho_dense = mc.predict( @@ -248,7 +248,7 @@ fig.tight_layout() ``` Dense and MPO should overlap; the surrogate approximates the same draw at the reference `rho0`. -Held-out Hamiltonian rollouts above use random probe `rho0` values from data generation — a different setup than the fixed reference state baked into process tensors. +Held-out Hamiltonian rollouts above use random probe `rho0` values from data generation — a different setup than the fixed reference state stored on process tensors. The same information functionals are available on either backend. For this short horizon, conditional mutual information is near zero while QMI grows when more past legs are included: diff --git a/docs/examples/quickstart.md b/docs/examples/quickstart.md index 32186623..04f0a254 100644 --- a/docs/examples/quickstart.md +++ b/docs/examples/quickstart.md @@ -135,7 +135,8 @@ For larger circuits, compiler passes, and OpenQASM inputs, see {doc}`equivalence ## 4. Characterize environmental memory -Probe a probe qubit coupled to a short chain at an interior temporal cut. The memory spectrum and response matrix show how many independent past branches remain visible at the cut: +Probe a probe qubit coupled to a short chain at an interior temporal cut. The memory spectrum and response matrix show how many independent past branches remain visible at the cut ($S_V$). +For temporal entanglement of a process tensor ($S_{PT}$), see {doc}`characterization`. ```{code-cell} ipython3 import numpy as np diff --git a/src/mqt/yaqs/characterization/memory/backends/__init__.py b/src/mqt/yaqs/characterization/memory/backends/__init__.py index 2983b99b..aa452534 100644 --- a/src/mqt/yaqs/characterization/memory/backends/__init__.py +++ b/src/mqt/yaqs/characterization/memory/backends/__init__.py @@ -12,7 +12,7 @@ - :mod:`.sequences` — :func:`~mqt.yaqs.characterization.memory.backends.sequences.simulate_sequences` and pool workers - :mod:`.exact` — :class:`~mqt.yaqs.characterization.memory.backends.exact.ExactBackend` and :func:`~mqt.yaqs.characterization.memory.backends.exact.simulate_exact` -- :mod:`.tomography` — reference dense/MPO process tensors via exhaustive tomography +- :mod:`.tomography` — MPO process tensors via direct construction (default); dense via exhaustive tomography - :mod:`.surrogates` — :class:`~mqt.yaqs.characterization.memory.backends.surrogates.model.ProcessTensorSurrogate` and training-data workflow """ diff --git a/src/mqt/yaqs/characterization/memory/backends/tomography/__init__.py b/src/mqt/yaqs/characterization/memory/backends/tomography/__init__.py index f8613976..f3f696c9 100644 --- a/src/mqt/yaqs/characterization/memory/backends/tomography/__init__.py +++ b/src/mqt/yaqs/characterization/memory/backends/tomography/__init__.py @@ -5,12 +5,14 @@ # # Licensed under the MIT License -"""Process-tensor tomography (exact/exhaustive). +"""Process-tensor construction (dense tomography and direct MPO). -This subpackage constructs a process tensor from exhaustive discrete intervention sequences (size -``16**num_interventions`` for ``num_interventions`` steps), optionally under MCWF/TJM noise, and returns a -:class:`~mqt.yaqs.characterization.memory.backends.tomography.process_tensors.DenseProcessTensor` or +``return_type="mpo"`` (default) builds an MPO by direct construction (noiseless) and returns :class:`~mqt.yaqs.characterization.memory.backends.tomography.process_tensors.MPOProcessTensor`. + +``return_type="dense"`` runs exhaustive discrete-basis tomography +(``16**num_interventions`` sequences) and returns +:class:`~mqt.yaqs.characterization.memory.backends.tomography.process_tensors.DenseProcessTensor`. """ from .basis import TomographyBasis as TomographyBasis @@ -21,3 +23,5 @@ from .data import SequenceData as SequenceData from .process_tensors import DenseProcessTensor as DenseProcessTensor from .process_tensors import MPOProcessTensor as MPOProcessTensor +from .process_tensors import compute_temporal_entropy as compute_temporal_entropy +from .process_tensors import evaluate_probes as evaluate_probes diff --git a/src/mqt/yaqs/characterization/memory/backends/tomography/constructor.py b/src/mqt/yaqs/characterization/memory/backends/tomography/constructor.py index 68862234..233269ee 100644 --- a/src/mqt/yaqs/characterization/memory/backends/tomography/constructor.py +++ b/src/mqt/yaqs/characterization/memory/backends/tomography/constructor.py @@ -5,18 +5,22 @@ # # Licensed under the MIT License -"""Process-tensor tomography workflow: exhaustive discrete-basis simulation. +"""Process-tensor construction: dense tomography or direct MPO growth. **Public** (see ``__all__`` in :mod:`mqt.yaqs.characterization.memory.backends.tomography`): :func:`build_process_tensor` (this module, :mod:`mqt.yaqs.characterization.memory.backends.tomography.constructor`). -:func:`build_process_tensor` is the high-level user entry point returning a process tensor directly -(:class:`~mqt.yaqs.characterization.memory.backends.tomography.process_tensors.DenseProcessTensor` or -:class:`~mqt.yaqs.characterization.memory.backends.tomography.process_tensors.MPOProcessTensor`). +:func:`build_process_tensor` dispatches on ``return_type``: + +- ``"mpo"`` (default) — direct MPO construction (noiseless); returns + :class:`~mqt.yaqs.characterization.memory.backends.tomography.process_tensors.MPOProcessTensor`. +- ``"dense"`` — exhaustive discrete-basis tomography (``16**num_interventions`` sequences), + optionally with noise; returns + :class:`~mqt.yaqs.characterization.memory.backends.tomography.process_tensors.DenseProcessTensor`. + The lower-level :func:`run_all_sequences` returns -:class:`~mqt.yaqs.characterization.memory.backends.tomography.data.SequenceData` covering all -``16**num_interventions`` Choi index sequences for ``num_interventions`` steps. +:class:`~mqt.yaqs.characterization.memory.backends.tomography.data.SequenceData` for the dense path. **Execution model** — Same pattern as :mod:`mqt.yaqs.simulator` and :mod:`mqt.yaqs.characterization.memory.backends.surrogates.workflow`: build a picklable payload, optionally @@ -383,48 +387,46 @@ def build_process_tensor( num_trajectories: int = 100, basis: TomographyBasis = "tetrahedral", basis_seed: int | None = None, - return_type: Literal["dense", "mpo"] = "dense", + return_type: Literal["dense", "mpo"] = "mpo", # Dense reconstruction check: bool = True, atol: float = 1e-8, - # MPO reconstruction - compress_every: int = 100, + # Direct MPO construction + compress_every: int = 16, tol: float = 1e-12, - max_bond_dim: int | None = None, + max_bond_dim: int | None = 64, n_sweeps: int = 2, solver: StochasticSolver | None = None, initial_rho: np.ndarray | None = None, initial_rho_atol: float = 1e-8, _execution: ExecutionConfig | None = None, ) -> DenseProcessTensor | MPOProcessTensor: - """Construct a process tensor via exhaustive discrete-basis tomography. - - This simulates **every** ``16**num_interventions`` discrete basis sequence and returns a - process tensor directly: + """Construct a process tensor as dense tomography or a direct MPO. - - ``return_type="dense"``: reconstruct and return a :class:`DenseProcessTensor`. - - ``return_type="mpo"``: build and return an :class:`MPOProcessTensor`. + - ``return_type="mpo"`` (default): direct MPO construction (noiseless only). + - ``return_type="dense"``: exhaustive discrete-basis tomography + (``16**num_interventions`` sequences; supports ``noise_model``). Args: operator: Hamiltonian MPO. sim_params: Analog simulation parameters. timesteps: Optional process-tensor schedule evolution durations (length ``num_interventions + 1``; defaults to ``[dt, dt]`` for one intervention leg). - noise_model: Optional open-system noise model. - parallel: Whether to parallelize over sequences. - num_trajectories: MCWF trajectories per sequence (forced to 1 when noiseless). - basis: Tomography basis name. + noise_model: Optional open-system noise model (dense tomography only). + parallel: Whether to parallelize dense tomography sequences or MPO construction. + num_trajectories: MCWF trajectories per sequence (dense path only). + basis: Tomography / Choi basis name. basis_seed: Optional seed when ``basis="random"``. - return_type: ``"dense"`` or ``"mpo"`` process-tensor representation. + return_type: ``"mpo"`` (direct construction, default) or ``"dense"`` (tomography). check: Run self-consistency check for dense reconstruction. atol: Absolute tolerance for the dense self-check. - compress_every: MPO rank-1 accumulation compress interval. + compress_every: Direct-MPO rank-1 accumulation compress interval. tol: MPO compression tolerance. - max_bond_dim: Optional MPO bond-dimension cap. + max_bond_dim: Cap on the branch ensemble / MPO bond dimension for direct construction. + Defaults to ``64`` for scalability; pass ``None`` for exact uncapped construction. n_sweeps: MPO compression sweeps. solver: Stochastic solver (``"MCWF"`` or ``"TJM"``). - initial_rho: Optional expected site-0 reference after ``U_0``; validated against the - computed tomography reference when provided. + initial_rho: Optional expected site-0 reference after ``U_0``. initial_rho_atol: Tolerance for optional ``initial_rho`` validation. _execution: Optional internal execution configuration. @@ -432,8 +434,36 @@ def build_process_tensor( Dense or MPO process-tensor wrapper depending on ``return_type``. Raises: - ValueError: If ``return_type`` is not ``"dense"`` or ``"mpo"``. + ValueError: If ``return_type`` is invalid, or ``noise_model`` is set with ``"mpo"``. """ + if return_type == "mpo": + if noise_model is not None: + msg = ( + "return_type='mpo' uses direct construction and does not support noise_model; use return_type='dense'." + ) + raise ValueError(msg) + from .direct import build_process_tensor_direct # noqa: PLC0415 + + return build_process_tensor_direct( + operator, + sim_params, + timesteps, + basis=basis, + basis_seed=basis_seed, + tol=tol, + max_bond_dim=max_bond_dim, + n_sweeps=n_sweeps, + compress_every=compress_every, + solver=solver, + initial_rho=initial_rho, + initial_rho_atol=initial_rho_atol, + parallel=parallel, + _execution=_execution, + ) + if return_type != "dense": + msg = f"Unknown return_type {return_type!r} (expected 'dense' or 'mpo')." + raise ValueError(msg) + data = _construct_data( operator, sim_params, @@ -450,14 +480,4 @@ def build_process_tensor( if initial_rho is not None: validate_initial_rho(coerce_rho_matrix(initial_rho), data.initial_rho, atol=initial_rho_atol) - if return_type == "dense": - return data.to_dense_process_tensor(check=check, atol=atol) - if return_type == "mpo": - return data.to_mpo_process_tensor( - compress_every=compress_every, - tol=tol, - max_bond_dim=max_bond_dim, - n_sweeps=n_sweeps, - ) - msg = f"Unknown return_type {return_type!r} (expected 'dense' or 'mpo')." - raise ValueError(msg) + return data.to_dense_process_tensor(check=check, atol=atol) diff --git a/src/mqt/yaqs/characterization/memory/backends/tomography/data.py b/src/mqt/yaqs/characterization/memory/backends/tomography/data.py index acb5bf92..2a198885 100644 --- a/src/mqt/yaqs/characterization/memory/backends/tomography/data.py +++ b/src/mqt/yaqs/characterization/memory/backends/tomography/data.py @@ -5,11 +5,12 @@ # # Licensed under the MIT License -"""Exhaustive discrete-basis process-tensor data + reconstruction helpers. +"""Exhaustive discrete-basis process-tensor data + dense reconstruction helpers. -The main product of :func:`~mqt.yaqs.characterization.memory.backends.tomography.constructor.build_process_tensor` -is :class:`SequenceData`. It can be converted to dense or MPO process-tensor representations via -:meth:`SequenceData.to_dense_process_tensor` and :meth:`SequenceData.to_mpo_process_tensor`. +:class:`SequenceData` holds exhaustive tomography outputs and converts to a dense process +tensor via :meth:`SequenceData.to_dense_process_tensor`. Rank-1 MPO helpers in this module +support direct MPO construction in +:mod:`mqt.yaqs.characterization.memory.backends.tomography.direct`. """ from __future__ import annotations @@ -21,7 +22,7 @@ from mqt.yaqs.core.data_structures.mpo import MPO -from .process_tensors import DenseProcessTensor, MPOProcessTensor +from .process_tensors import DenseProcessTensor if TYPE_CHECKING: from collections.abc import Iterable @@ -30,7 +31,14 @@ def _num_intervention_steps(timesteps: list[float]) -> int: - """Return intervention-leg count ``k`` from a process-tensor schedule of length ``k + 1``.""" + """Return intervention-leg count ``k`` from a process-tensor schedule of length ``k + 1``. + + Args: + timesteps: Process-tensor schedule durations. + + Returns: + Number of intervention legs ``max(0, len(timesteps) - 1)``. + """ return max(0, len(timesteps) - 1) @@ -123,22 +131,6 @@ def pack_sequence_outputs(data: SequenceData) -> tuple[NDArray[np.complex128], N return out_vecs, seq_weights -def _iter_rank1_terms(data: SequenceData) -> Iterable[MPO]: - """Yield rank-1 MPO terms for MPO process-tensor construction. - - Args: - data: SequenceData instance. - - Yields: - MPO rank-1 terms. - """ - for i, alpha in enumerate(data.sequences): - rho_out = data.outputs[i] - w = float(data.weights[i]) - dual_ops = [data.choi_duals[a].T for a in alpha] - yield _rank1_mpo_term(rho_out, dual_ops, weight=w) - - def assemble_upsilon( *, out_vecs: NDArray[np.complex128], @@ -256,34 +248,3 @@ def to_dense_process_tensor(self, *, check: bool = True, atol: float = 1e-8) -> atol=atol, ) return DenseProcessTensor(upsilon, list(self.timesteps), initial_rho=self.initial_rho.copy()) - - def to_mpo_process_tensor( - self, - *, - compress_every: int = 100, - tol: float = 1e-12, - max_bond_dim: int | None = None, - n_sweeps: int = 2, - ) -> MPOProcessTensor: - """Build an MPO process tensor via rank-1 accumulation. - - Args: - compress_every: Compress after this many terms. - tol: Compression tolerance. - max_bond_dim: Optional maximum bond dimension. - n_sweeps: Number of compression sweeps. - - Returns: - MPO process-tensor representation. - """ - num_steps = _num_intervention_steps(self.timesteps) - mpo = accumulate_rank1_terms( - _iter_rank1_terms(self), - num_steps=num_steps, - dims=(2, 2), - compress_every=compress_every, - tol=tol, - max_bond_dim=max_bond_dim, - n_sweeps=n_sweeps, - ) - return MPOProcessTensor(mpo, self.timesteps, initial_rho=self.initial_rho.copy()) diff --git a/src/mqt/yaqs/characterization/memory/backends/tomography/direct.py b/src/mqt/yaqs/characterization/memory/backends/tomography/direct.py new file mode 100644 index 00000000..37a1e866 --- /dev/null +++ b/src/mqt/yaqs/characterization/memory/backends/tomography/direct.py @@ -0,0 +1,410 @@ +# Copyright (c) 2025 - 2026 Chair for Design Automation, TUM +# All rights reserved. +# +# SPDX-License-Identifier: MIT +# +# Licensed under the MIT License + +"""Leg-by-leg process-tensor MPO construction without exhaustive ``16**k`` tomography.""" + +from __future__ import annotations + +import copy +from dataclasses import dataclass +from typing import TYPE_CHECKING, Any + +import numpy as np +from numpy.typing import NDArray + +from mqt.yaqs.core.data_structures.mps import MPS +from mqt.yaqs.core.parallel_utils import ExecutionConfig, merge_execution_config, resolve_worker_ctx, run_indexed_jobs + +from ...shared.encoding import normalize_backend_rho +from ...shared.intervention_steps import apply_intervention_to_backend +from ...shared.utils import ( + StochasticSolver, + _evolve_backend_state, + _initialize_backend_state, + extract_site0_rho, + make_mcwf_static_context, + resolve_stochastic_solver, +) +from ..sequences.workers import _get_times_cached +from .basis import TomographyBasis, assemble_fixed_basis, compute_dual_choi_basis +from .constructor import _reference_initial_rho +from .data import _rank1_mpo_term, accumulate_rank1_terms +from .process_tensors import MPOProcessTensor, validate_initial_rho + +if TYPE_CHECKING: + from mqt.yaqs.analog.mcwf import MCWFContext + from mqt.yaqs.core.data_structures.mpo import MPO + from mqt.yaqs.core.data_structures.simulation_parameters import AnalogSimParams + +_N_CHOI = 16 +BackendState = MPS | NDArray[np.complex128] + + +@dataclass +class _Branch: + """One definite past intervention history and backend state before the next leg.""" + + history: tuple[int, ...] + psi: BackendState + weight: float + + +def _clone_backend_state(state: BackendState) -> BackendState: + """Return an independent copy of a dense or MPS backend state. + + Args: + state: Dense vector (MCWF) or MPS (TJM). + + Returns: + Deep-copied MPS or a contiguous dense copy. + """ + if isinstance(state, MPS): + return copy.deepcopy(state) + return np.asarray(state, dtype=np.complex128).reshape(-1).copy() + + +def _compress_branches( + branches: list[_Branch], + *, + max_bond_dim: int | None, + tol: float, +) -> list[_Branch]: + """Compress the branch ensemble to at most ``max_bond_dim`` states. + + Dense (MCWF) ensembles use a weighted SVD. MPS (TJM) ensembles keep the + highest-weight branches without densifying. + + Args: + branches: Current ensemble of weighted backend states. + max_bond_dim: Optional cap on retained branches; ``None`` keeps all. + tol: Singular-value floor when selecting kept dense modes. + + Returns: + Compressed branch list (unchanged when already within the cap). + """ + if max_bond_dim is None or len(branches) <= max_bond_dim: + return branches + if len(branches) == 1: + return branches + + if isinstance(branches[0].psi, MPS): + ordered = sorted(branches, key=lambda br: br.weight, reverse=True) + return ordered[: int(max_bond_dim)] + + dim = int(np.asarray(branches[0].psi, dtype=np.complex128).reshape(-1).size) + n = len(branches) + mat = np.zeros((dim, n), dtype=np.complex128) + for col, br in enumerate(branches): + scale = float(np.sqrt(max(br.weight, 0.0))) + mat[:, col] = scale * np.asarray(br.psi, dtype=np.complex128).reshape(-1) + + _u, singular_values, vh = np.linalg.svd(mat, full_matrices=False) + keep = int(np.sum(singular_values > tol)) + keep = min(keep, int(max_bond_dim)) + keep = max(1, keep) + + out: list[_Branch] = [] + for row in range(keep): + coeffs = vh[row, :] + i_dom = int(np.argmax(np.abs(coeffs))) + psi = mat @ coeffs.conj() + norm = float(np.linalg.norm(psi)) + if norm <= 1e-15: + psi = np.asarray(branches[i_dom].psi, dtype=np.complex128).reshape(-1).copy() + norm = float(np.linalg.norm(psi)) + else: + psi /= norm + weight = float(singular_values[row] ** 2) + out.append(_Branch(history=branches[i_dom].history, psi=psi, weight=weight)) + return out + + +def _prepare_step( + operator: MPO, + sim_params: AnalogSimParams, + duration: float, + *, + solver: StochasticSolver, +) -> tuple[AnalogSimParams, MCWFContext | None]: + """Build local sim params and MCWF context for one schedule slot. + + Args: + operator: Hamiltonian MPO. + sim_params: Analog simulation parameters. + duration: Evolution duration for this slot. + solver: Stochastic solver name. + + Returns: + Tuple ``(step_params, static_ctx)``. + """ + local_params = copy.copy(sim_params) + local_params.get_state = True + local_params.num_traj = 1 + static_ctx = make_mcwf_static_context(operator, local_params, noise_model=None) if solver == "MCWF" else None + times_cache: dict[tuple[float, float], np.ndarray] = {} + step_params = copy.copy(local_params) + step_params.elapsed_time = float(duration) + step_params.times = _get_times_cached(times_cache, dt=float(step_params.dt), duration=float(duration)) + return step_params, static_ctx + + +def _evolve_initial_state( + operator: MPO, + sim_params: AnalogSimParams, + duration: float, + *, + solver: StochasticSolver, +) -> BackendState: + """Evolve from ``|0...0>`` for one schedule slot. + + Args: + operator: Hamiltonian MPO. + sim_params: Analog simulation parameters. + duration: Evolution duration for ``U_0``. + solver: Stochastic solver name. + + Returns: + Backend state after the initial evolution (dense for MCWF, MPS for TJM). + """ + step_params, static_ctx = _prepare_step(operator, sim_params, duration, solver=solver) + state = _initialize_backend_state(operator, solver) + return _evolve_backend_state( + state, + operator, + None, + step_params, + solver, + traj_idx=0, + static_ctx=static_ctx, + ) + + +def _branch_extension_worker( + job_idx: int, + job_payload: dict[str, Any] | None = None, +) -> tuple[tuple[int, ...], BackendState, float, NDArray[np.complex128]] | None: + """Extend one branch by one Choi-basis intervention and post-evolution. + + Flat index layout is ``branch_index * 16 + choi_index``. Branch states keep their + native backend type (dense for MCWF, MPS for TJM). + + Args: + job_idx: Flat index over branches and Choi-basis elements. + job_payload: Shared step payload; defaults to + :data:`~mqt.yaqs.core.parallel_utils.WORKER_CTX`. + + Returns: + ``(history, psi, weight, rho_out)`` for kept extensions, or ``None`` when the + cumulative weight falls below the discard threshold. + """ + ctx = resolve_worker_ctx(job_payload) + br_idx, choi_idx = divmod(int(job_idx), _N_CHOI) + br: _Branch = ctx["branches"][br_idx] + prep_idx, meas_idx = ctx["choi_indices"][choi_idx] + basis_set = ctx["basis_set"] + meas_psi, prep_psi = basis_set[meas_idx][1], basis_set[prep_idx][1] + + state = _clone_backend_state(br.psi) + state, step_prob = apply_intervention_to_backend( + state, + (meas_psi, prep_psi), + solver=ctx["solver"], + chain_length=int(ctx["chain_length"]), + ) + weight = float(br.weight) * float(step_prob) + if weight <= 1e-30: + return None + + state = _evolve_backend_state( + state, + ctx["operator"], + None, + ctx["sim_params"], + ctx["solver"], + traj_idx=0, + static_ctx=ctx["static_ctx"], + ) + rho_out = normalize_backend_rho(extract_site0_rho(state)) + history = (*br.history, choi_idx) + return history, state, weight, rho_out + + +def _apply_timestep( + branches: list[_Branch], + *, + operator: MPO, + sim_params: AnalogSimParams, + duration: float, + basis_set: list[tuple[str, NDArray[np.complex128], NDArray[np.complex128]]], + choi_indices: list[tuple[int, int]], + choi_duals: list[NDArray[np.complex128]], + solver: StochasticSolver, + execution: ExecutionConfig, +) -> tuple[list[_Branch], list[MPO]]: + """Extend every branch by one local CPTP leg via indexed jobs. + + Args: + branches: Ensemble before the intervention leg. + operator: Hamiltonian MPO. + sim_params: Analog simulation parameters. + duration: Evolution duration after the intervention. + basis_set: Discrete basis kets ``(label, ket, dual)``. + choi_indices: Map from Choi index to ``(prep_idx, meas_idx)``. + choi_duals: Dual Choi operators for rank-1 MPO assembly. + solver: Stochastic solver name. + execution: Parallel / serial job-dispatch configuration. + + Returns: + Expanded branches and the corresponding rank-1 MPO terms. + """ + step_params, static_ctx = _prepare_step(operator, sim_params, duration, solver=solver) + n_jobs = len(branches) * _N_CHOI + payload: dict[str, Any] = { + "branches": branches, + "operator": operator, + "sim_params": step_params, + "basis_set": basis_set, + "choi_indices": choi_indices, + "solver": solver, + "static_ctx": static_ctx, + "chain_length": int(operator.length), + } + job_results = run_indexed_jobs( + _branch_extension_worker, + payload=payload, + n_jobs=n_jobs, + config=execution, + desc=f"MPO construction ({len(branches)} branches)", + ) + + expanded: list[_Branch] = [] + terms: list[MPO] = [] + for job_idx in range(n_jobs): + out = job_results[job_idx] + if out is None: + continue + history, state, weight, rho_out = out + dual_ops = [choi_duals[idx].T for idx in history] + terms.append(_rank1_mpo_term(rho_out, dual_ops, weight=weight)) + expanded.append(_Branch(history=history, psi=state, weight=weight)) + return expanded, terms + + +def build_process_tensor_direct( + operator: MPO, + sim_params: AnalogSimParams, + timesteps: list[float] | None = None, + *, + basis: TomographyBasis = "tetrahedral", + basis_seed: int | None = None, + tol: float = 1e-12, + max_bond_dim: int | None = 64, + n_sweeps: int = 2, + compress_every: int = 16, + solver: StochasticSolver | None = None, + initial_rho: np.ndarray | None = None, + initial_rho_atol: float = 1e-8, + parallel: bool = True, + _execution: ExecutionConfig | None = None, +) -> MPOProcessTensor: + """Build a process-tensor MPO by leg-by-leg contraction. + + At each timestep only ``16 * chi`` local basis updates are simulated, where ``chi`` is the + compressed branch count from the previous step. This avoids enumerating all ``16**k`` sequences. + Construction is noiseless (site-0 interventions only). Within each intervention step, branch + extensions are dispatched with :func:`~mqt.yaqs.core.parallel_utils.run_indexed_jobs`. + + Args: + operator: Hamiltonian MPO. + sim_params: Analog simulation parameters. + timesteps: Process-tensor schedule of length ``num_interventions + 1``. + basis: Discrete Choi basis name. + basis_seed: Optional seed when ``basis="random"``. + tol: MPO compression tolerance. + max_bond_dim: Cap on the branch ensemble / MPO bond dimension. Defaults to ``64`` so + branch compression runs for scalability; pass ``None`` for exact uncapped construction. + n_sweeps: MPO compression sweeps after each step. + compress_every: Rank-1 accumulation batch size before intermediate compression. + solver: Stochastic solver (``"MCWF"`` or ``"TJM"``). + initial_rho: Optional reference site-0 state after ``U_0``. + initial_rho_atol: Tolerance for optional ``initial_rho`` validation. + parallel: Whether to parallelize branch extensions within each intervention step. + _execution: Optional internal execution configuration. + + Returns: + MPO process-tensor wrapper. + + Raises: + ValueError: If ``num_interventions`` is zero or the solver is unsupported. + """ + if timesteps is None: + dt = float(sim_params.dt) + timesteps = [dt, dt] + + stochastic_solver = resolve_stochastic_solver(sim_params, solver=solver) + if stochastic_solver not in {"MCWF", "TJM"}: + msg = f"Direct construction requires solvers MCWF or TJM, got {stochastic_solver!r}." + raise ValueError(msg) + + num_interventions = len(timesteps) - 1 + if num_interventions <= 0: + msg = "Direct construction requires at least one intervention leg." + raise ValueError(msg) + + basis_set, choi_basis, choi_indices, _choi_feat = assemble_fixed_basis(basis=basis, basis_seed=basis_seed) + choi_duals = compute_dual_choi_basis(choi_basis) + execution = merge_execution_config(_execution, parallel=parallel) + + ref_rho = _reference_initial_rho( + operator, + sim_params, + timesteps, + noise_model=None, + solver=stochastic_solver, + num_trajectories=1, + ) + if initial_rho is not None: + validate_initial_rho(np.asarray(initial_rho, dtype=np.complex128), ref_rho, atol=initial_rho_atol) + + psi0 = _evolve_initial_state( + operator, + sim_params, + float(timesteps[0]), + solver=stochastic_solver, + ) + branches = [_Branch(history=(), psi=psi0, weight=1.0)] + + comb: MPO | None = None + for step_idx in range(num_interventions): + branches, terms = _apply_timestep( + branches, + operator=operator, + sim_params=sim_params, + duration=float(timesteps[step_idx + 1]), + basis_set=basis_set, + choi_indices=choi_indices, + choi_duals=choi_duals, + solver=stochastic_solver, + execution=execution, + ) + if not terms: + msg = f"Direct construction produced no rank-1 terms at leg {step_idx + 1}." + raise ValueError(msg) + comb = accumulate_rank1_terms( + terms, + num_steps=step_idx + 1, + tol=tol, + max_bond_dim=max_bond_dim, + n_sweeps=n_sweeps, + compress_every=compress_every, + ) + branches = _compress_branches(branches, max_bond_dim=max_bond_dim, tol=tol) + + if comb is None: + comb = _rank1_mpo_term(ref_rho, [], weight=1.0) + + return MPOProcessTensor(comb, list(timesteps), initial_rho=ref_rho.copy()) diff --git a/src/mqt/yaqs/characterization/memory/backends/tomography/process_tensors.py b/src/mqt/yaqs/characterization/memory/backends/tomography/process_tensors.py index 1cd924bb..51116f4c 100644 --- a/src/mqt/yaqs/characterization/memory/backends/tomography/process_tensors.py +++ b/src/mqt/yaqs/characterization/memory/backends/tomography/process_tensors.py @@ -9,7 +9,7 @@ from __future__ import annotations -from typing import TYPE_CHECKING, cast +from typing import TYPE_CHECKING, Protocol, cast import numpy as np @@ -27,6 +27,17 @@ from ...operational_memory.samples import ProbeSet +class SupportsPredict(Protocol): + """Process-tensor backends that map intervention sequences to a final state.""" + + def predict( + self, + interventions: list[Callable[[NDArray[np.complex128]], NDArray[np.complex128]]], + ) -> NDArray[np.complex128]: + """Predict the final reduced state for a sequence of interventions.""" + ... + + def validate_initial_rho( rho0: NDArray[np.complex128], reference: NDArray[np.complex128], @@ -53,7 +64,7 @@ def validate_initial_rho( def convert_probe_callable( step: AnyInterventionStep, ) -> Callable[[NDArray[np.complex128]], NDArray[np.complex128]]: - """Convert a probe-grid step to a CPTP map callable for :meth:`DenseProcessTensor.predict`. + """Convert a probe-grid step to a CPTP map callable for :meth:`~SupportsPredict.predict`. Args: step: Structured dict step or measure/prepare ket pair. @@ -72,11 +83,14 @@ def unitary_map(rho: NDArray[np.complex128]) -> NDArray[np.complex128]: return inter -def evaluate_dense_probes(process_tensor: DenseProcessTensor, probe_set: ProbeSet) -> np.ndarray: - """Evaluate split-cut probe Pauli responses on a dense process tensor. +def evaluate_probes(process_tensor: SupportsPredict, probe_set: ProbeSet) -> np.ndarray: + """Evaluate split-cut probe Pauli responses via process-tensor :meth:`predict`. + + Shared by dense and MPO process tensors for operational-memory V-matrix assembly. + Does not densify MPO tensors. Args: - process_tensor: Dense reference process-tensor backend. + process_tensor: Backend implementing :meth:`~SupportsPredict.predict`. probe_set: Sampled split-cut probes. Returns: @@ -170,6 +184,152 @@ def compute_entropy_dense(r: NDArray[np.complex128], base: int = 2) -> float: return float(-(nz * (np.log(nz) / log_base)).sum()) +def _validate_cut(cut: int, num_interventions: int) -> None: + if cut < 1 or cut > num_interventions: + msg = f"cut must satisfy 1 <= cut <= num_interventions ({num_interventions}), got {cut}." + raise ValueError(msg) + + +def _unfuse_slot_index(fused: int, *, out_first: bool = True) -> tuple[int, int]: + """Split a fused 4-index Choi leg into ``(output, input)`` qubit indices. + + ``encode_cptp_choi`` uses ``kron(output, input)`` so ``f = 2 * out + in`` by default. + + Returns: + Tuple ``(output_index, input_index)`` each in ``{0, 1}``. + """ + if out_first: + return fused // 2, fused % 2 + return fused % 2, fused // 2 + + +def _upsilon_to_unfused_operator( + upsilon: NDArray[np.complex128], + num_interventions: int, + *, + out_first: bool = True, +) -> NDArray[np.complex128]: + """Reshape a process-tensor Choi operator into explicit ket/bra qubit axes. + + Subsystem order in ``upsilon`` is ``[final(2), slot_1(4), …, slot_k(4)]`` with + ``slot_t = output_t ⊗ input_t`` and ``f = 2 * output + input`` when ``out_first=True``. + + Returns: + Tensor with axes ``final_ket/bra`` then per-slot ``out/in`` ket/bra pairs. + + Raises: + ValueError: If ``upsilon`` shape is inconsistent with ``num_interventions``. + """ + k = num_interventions + expected = 2 * (4**k) + ups = np.asarray(upsilon, dtype=np.complex128) + if ups.shape != (expected, expected): + msg = f"Expected upsilon shape ({expected}, {expected}) for k={k}, got {ups.shape}." + raise ValueError(msg) + dims = [2] + [4] * k + mat = ups.reshape(*dims, *dims) + out = np.zeros([2, 2] + [2, 2, 2, 2] * k, dtype=np.complex128) + for idx in np.ndindex(*dims, *dims): + sub_k = idx[: k + 1] + sub_b = idx[k + 1 :] + coords: list[int] = [sub_k[0], sub_b[0]] + for t in range(k): + ok, ik = _unfuse_slot_index(sub_k[t + 1], out_first=out_first) + ob, ib = _unfuse_slot_index(sub_b[t + 1], out_first=out_first) + coords.extend([ok, ik, ob, ib]) + out[tuple(coords)] = mat[idx] + return out + + +def _block_axis_indices(num_interventions: int) -> list[list[int]]: + """Return unfused tensor axis indices for causal blocks ``B_0 … B_k``. + + Axis numbering matches :func:`_upsilon_to_unfused_operator`: + + - ``final_ket=0``, ``final_bra=1`` + - slot ``t`` (0-based): ``out_ket=2+4t``, ``in_ket=3+4t``, ``out_bra=4+4t``, ``in_bra=5+4t`` + + Args: + num_interventions: Number of intervention slots ``k``. + + Returns: + List of ``k + 1`` blocks of axis indices. + """ + k = num_interventions + blocks: list[list[int]] = [[3, 5]] + blocks.extend([2 + 4 * t, 3 + 4 * (t + 1), 4 + 4 * t, 5 + 4 * (t + 1)] for t in range(k - 1)) + blocks.append([2 + 4 * (k - 1), 0, 4 + 4 * (k - 1), 1]) + return blocks + + +def compute_temporal_entropy( + upsilon: NDArray[np.complex128], + num_interventions: int, + cut: int, + *, + rtol: float = 1e-12, + weight_tol: float = 1e-30, +) -> dict[str, NDArray[np.float64] | float | int]: + r"""Compute temporal entanglement of the process tensor at a causal cut. + + Partitions causal blocks ``B_0, \ldots, B_k`` at cut ``c`` as:: + + LEFT = B_0, …, B_{c-1} + RIGHT = B_c, …, B_k + + and computes the operator-Schmidt spectrum of the unfused Choi operator without + partial tracing or trace normalization. The result is temporal entanglement + :math:`S_{PT}(c)`, distinct from operational response entropy :math:`S_V(c)`. + + Args: + upsilon: Dense process-tensor Choi matrix. + num_interventions: Intervention count ``k``. + cut: Causal cut index ``c`` matching the response protocol. + rtol: Relative threshold ``s_i > rtol * s_0`` for resolved Schmidt rank. + weight_tol: Absolute floor on ``sum(s**2)``; below this raises ``ValueError``. + + Returns: + Dictionary with keys ``entropy`` (:math:`S_{PT}`), ``effective_rank``, + ``schmidt_rank``, ``singular_values``, and ``weights``. + + Raises: + ValueError: If ``cut`` is invalid or the squared-Schmidt weight sum is below ``weight_tol``. + """ + _validate_cut(cut, num_interventions) + op = _upsilon_to_unfused_operator(upsilon, num_interventions) + blocks = _block_axis_indices(num_interventions) + left_axes = [i for b in blocks[:cut] for i in b] + right_axes = [i for b in blocks[cut:] for i in b] + perm = left_axes + right_axes + tensor_perm = np.transpose(op, perm) + dim_left = int(np.prod([tensor_perm.shape[i] for i in range(len(left_axes))], dtype=np.int64)) + dim_right = int( + np.prod([tensor_perm.shape[i] for i in range(len(left_axes), len(left_axes) + len(right_axes))], dtype=np.int64) + ) + mat = tensor_perm.reshape(dim_left, dim_right) + singular_values = np.linalg.svd(mat, compute_uv=False).astype(np.float64) + total_weight = float(np.sum(singular_values**2)) + if total_weight < weight_tol: + msg = f"Operator-Schmidt weight sum {total_weight:.3e} below tolerance {weight_tol:.3e}." + raise ValueError(msg) + weights = singular_values**2 / total_weight + nz = weights > weight_tol + entropy = float(-np.sum(weights[nz] * np.log(weights[nz]))) if np.any(nz) else 0.0 + if singular_values.size and singular_values[0] > 0.0: + resolved = singular_values > rtol * singular_values[0] + else: + resolved = singular_values > 0.0 + schmidt_rank = int(np.sum(resolved)) + effective_rank = float(np.exp(entropy)) if entropy > 0.0 else 1.0 + return { + "entropy": entropy, + "effective_rank": effective_rank, + "schmidt_rank": schmidt_rank, + "singular_values": singular_values, + "weights": weights, + } + + class DenseProcessTensor: """Wrapper around a dense process-tensor Choi operator Upsilon.""" @@ -227,6 +387,31 @@ def _num_interventions(self) -> int: size = self.upsilon.shape[0] return int(np.round(np.log2(size / 2) / 2)) + def compute_temporal_entropy( + self, + cut: int, + *, + rtol: float = 1e-12, + weight_tol: float = 1e-30, + ) -> dict[str, NDArray[np.float64] | float | int]: + """Compute temporal entanglement :math:`S_{PT}(c)` at ``cut``. + + Args: + cut: Causal cut index ``c`` matching the response protocol. + rtol: Relative Schmidt threshold for ``schmidt_rank``. + weight_tol: Absolute floor on ``sum(s**2)``. + + Returns: + Result dictionary from :func:`compute_temporal_entropy`. + """ + return compute_temporal_entropy( + self.upsilon, + self._num_interventions(), + cut, + rtol=rtol, + weight_tol=weight_tol, + ) + def _predict_raw( self, interventions: list[Callable[[NDArray[np.complex128]], NDArray[np.complex128]]], @@ -296,12 +481,15 @@ def _num_interventions_for_probe(self) -> int: return self._num_interventions() def evaluate_probes(self, probe_set: ProbeSet) -> np.ndarray: - """Evaluate split-cut probe Pauli responses. + """Evaluate split-cut probe Pauli responses for V-matrix assembly. + + Args: + probe_set: Sampled split-cut probes. Returns: Array of shape ``(n_pasts, n_futures, 4)``. """ - return evaluate_dense_probes(self, probe_set) + return evaluate_probes(self, probe_set) def qmi( self, @@ -472,13 +660,39 @@ def to_dense(self) -> DenseProcessTensor: def _num_interventions_for_probe(self) -> int: return int(self.length) - 1 + def compute_temporal_entropy( + self, + cut: int, + *, + rtol: float = 1e-12, + weight_tol: float = 1e-30, + ) -> dict[str, NDArray[np.float64] | float | int]: + """Compute temporal entanglement :math:`S_{PT}(c)` at ``cut``. + + Delegates to the dense representation via :meth:`to_dense`. + + Args: + cut: Causal cut index ``c`` matching the response protocol. + rtol: Relative Schmidt threshold for ``schmidt_rank``. + weight_tol: Absolute floor on ``sum(s**2)``. + + Returns: + Result dictionary from :func:`compute_temporal_entropy`. + """ + return self.to_dense().compute_temporal_entropy(cut, rtol=rtol, weight_tol=weight_tol) + def evaluate_probes(self, probe_set: ProbeSet) -> np.ndarray: - """Evaluate split-cut probe Pauli responses. + """Evaluate split-cut probe Pauli responses for V-matrix assembly. + + Uses native MPO :meth:`predict` (does not densify the process tensor). + + Args: + probe_set: Sampled split-cut probes. Returns: Array of shape ``(n_pasts, n_futures, 4)`` with Pauli tomography coefficients. """ - return self.to_dense().evaluate_probes(probe_set) + return evaluate_probes(self, probe_set) def predict( self, diff --git a/src/mqt/yaqs/memory_characterizer.py b/src/mqt/yaqs/memory_characterizer.py index 2b7a2e1a..09e1251d 100644 --- a/src/mqt/yaqs/memory_characterizer.py +++ b/src/mqt/yaqs/memory_characterizer.py @@ -321,41 +321,46 @@ def build_process_tensor( num_trajectories: int = 100, basis: TomographyBasis = "tetrahedral", basis_seed: int | None = None, - return_type: Literal["dense", "mpo"] = "dense", + return_type: Literal["dense", "mpo"] = "mpo", check: bool = True, atol: float = 1e-8, - compress_every: int = 100, + compress_every: int = 16, tol: float = 1e-12, - max_bond_dim: int | None = None, + max_bond_dim: int | None = 64, n_sweeps: int = 2, parallel: bool | None = None, initial_rho: np.ndarray | None = None, initial_rho_atol: float = 1e-8, ) -> DenseProcessTensor | MPOProcessTensor: - """Build an exhaustive reference process tensor (validation only; scales as ``16**num_interventions``). + """Build a process tensor via dense tomography or direct MPO construction. + + - ``return_type="mpo"`` (default): direct MPO construction (noiseless only). + - ``return_type="dense"``: exhaustive tomography (scales as ``16**num_interventions``; + supports ``noise_model``). Args: hamiltonian: System Hamiltonian. sim_params: Analog simulation parameters. timesteps: Optional process-tensor schedule evolution durations (length ``num_interventions + 1``; defaults to ``[dt, dt]`` for one intervention leg). - noise_model: Optional noise model during tomography sequences. - num_trajectories: Monte Carlo trajectories per tomography sample. - basis: Intervention basis for process-tensor tomography. + noise_model: Optional noise model (dense tomography only). + num_trajectories: Monte Carlo trajectories per tomography sample (dense only). + basis: Intervention / Choi basis name. basis_seed: Optional RNG seed for basis construction. - return_type: ``"dense"`` or ``"mpo"`` process-tensor storage. - check: Whether to validate CPTP properties during construction. + return_type: ``"mpo"`` (direct construction, default) or ``"dense"`` (tomography). + check: Whether to validate CPTP properties during dense construction. atol: CPTP check tolerance. - compress_every: MPO compression cadence during construction. + compress_every: How often to compress while accumulating direct-MPO terms. tol: MPO compression tolerance. - max_bond_dim: Optional MPO bond-dimension cap. - n_sweeps: MPO variational refinement sweeps. - parallel: Override instance parallel setting. + max_bond_dim: Cap on the branch ensemble / MPO bond dimension for direct construction. + Defaults to ``64`` for scalability; pass ``None`` for exact uncapped construction. + n_sweeps: MPO compression sweeps. + parallel: Override instance parallel setting for dense tomography or MPO construction. initial_rho: Optional expected site-0 reference after ``U_0``; validated when provided. initial_rho_atol: Tolerance for optional ``initial_rho`` validation. Returns: - Dense or MPO reference process tensor for small-horizon validation. + Dense or MPO process tensor depending on ``return_type``. """ operator = _require_hamiltonian(hamiltonian) execution = self._execution if parallel is None else merge_execution_config(self._execution, parallel=parallel) diff --git a/tests/characterization/memory/backends/tomography/test_constructor.py b/tests/characterization/memory/backends/tomography/test_constructor.py index cf339bde..29bd1fbd 100644 --- a/tests/characterization/memory/backends/tomography/test_constructor.py +++ b/tests/characterization/memory/backends/tomography/test_constructor.py @@ -38,30 +38,35 @@ def test_build_process_tensor_returns_dense_and_mpo_smoke() -> None: dense = mc.build_process_tensor(ham, params, timesteps=[0.0, 0.0], return_type="dense") assert dense.to_matrix().shape == (8, 8) - mpo = mc.build_process_tensor(ham, params, timesteps=[0.0, 0.0], return_type="mpo", compress_every=1) + mpo = mc.build_process_tensor(ham, params, timesteps=[0.0, 0.0], compress_every=1) mat = mpo.to_matrix() assert mat.shape == (8, 8) np.testing.assert_allclose(mat, dense.to_matrix(), atol=1e-8) def test_build_process_tensor_rejects_k_zero() -> None: - """Zero-step tomography is rejected before sequence enumeration.""" + """Zero-step schedules are rejected for both MPO and dense construction.""" op = MPO.ising(length=1, J=0.0, g=0.0) params = AnalogSimParams(dt=0.1, max_bond_dim=8) - with pytest.raises(ValueError, match="No sequences for num_interventions=0"): + with pytest.raises(ValueError, match="at least one intervention"): build_process_tensor(op, params, timesteps=[]) - with pytest.raises(ValueError, match="No sequences for num_interventions=0"): + with pytest.raises(ValueError, match="at least one intervention"): build_process_tensor(op, params, timesteps=[0.1]) + with pytest.raises(ValueError, match="No sequences for num_interventions=0"): + build_process_tensor(op, params, timesteps=[], return_type="dense") + with pytest.raises(ValueError, match="No sequences for num_interventions=0"): + build_process_tensor(op, params, timesteps=[0.1], return_type="dense") def test_build_process_tensor_parallel_smoke() -> None: - """build_process_tensor runs with parallel execution enabled.""" + """build_process_tensor runs with parallel execution enabled for dense and MPO.""" ham = Hamiltonian.ising(length=1, J=0.0, g=0.0) params = AnalogSimParams(dt=0.1, max_bond_dim=8) - dense = MemoryCharacterizer(parallel=True, max_workers=2, show_progress=False).build_process_tensor( - ham, params, timesteps=[0.0, 0.0], return_type="dense" - ) + mc = MemoryCharacterizer(parallel=True, max_workers=2, show_progress=False) + dense = mc.build_process_tensor(ham, params, timesteps=[0.0, 0.0], return_type="dense") assert dense.to_matrix().shape == (8, 8) + mpo = mc.build_process_tensor(ham, params, timesteps=[0.0, 0.0], compress_every=1) + assert mpo.to_matrix().shape == (8, 8) def test_build_process_tensor_stores_reference_initial_rho() -> None: diff --git a/tests/characterization/memory/backends/tomography/test_data.py b/tests/characterization/memory/backends/tomography/test_data.py index f8dc043f..6bfac8f3 100644 --- a/tests/characterization/memory/backends/tomography/test_data.py +++ b/tests/characterization/memory/backends/tomography/test_data.py @@ -58,29 +58,3 @@ def test_to_dense_sequence_data_zero_step_weighted() -> None: atol=1e-8, ) np.testing.assert_allclose(rho_w, 0.25 * rho, atol=1e-12) - - -def test_to_mpo_sequence_data_minimal() -> None: - """Smoke test SequenceData.to_mpo_process_tensor on minimal data.""" - rho = np.eye(2, dtype=np.complex128) - seqs: list[tuple[int, ...]] = [(0,)] - outputs = [rho] - weights = [1.0] - choi_basis = [np.eye(4, dtype=np.complex128)] * 16 - choi_indices = [(0, 0)] * 16 - choi_duals = [np.eye(4, dtype=np.complex128)] * 16 - timesteps = [0.1, 0.1] - - data = SequenceData( - sequences=seqs, - outputs=outputs, - weights=weights, - choi_basis=choi_basis, - choi_indices=choi_indices, - choi_duals=choi_duals, - timesteps=timesteps, - initial_rho=_REF_RHO0, - ) - pt = data.to_mpo_process_tensor(compress_every=1) - mat = pt.to_matrix() - assert mat.shape == (2 * 4, 2 * 4) diff --git a/tests/characterization/memory/backends/tomography/test_direct.py b/tests/characterization/memory/backends/tomography/test_direct.py new file mode 100644 index 00000000..27b7b5e4 --- /dev/null +++ b/tests/characterization/memory/backends/tomography/test_direct.py @@ -0,0 +1,191 @@ +# Copyright (c) 2025 - 2026 Chair for Design Automation, TUM +# All rights reserved. +# +# SPDX-License-Identifier: MIT +# +# Licensed under the MIT License + +"""Tests for leg-by-leg direct process-tensor construction.""" + +from __future__ import annotations + +from typing import TYPE_CHECKING, cast + +import numpy as np +import pytest + +from mqt.yaqs import AnalogSimParams, Hamiltonian, MemoryCharacterizer +from mqt.yaqs.characterization.memory.backends.tomography.constructor import build_process_tensor +from mqt.yaqs.characterization.memory.backends.tomography.process_tensors import MPOProcessTensor +from mqt.yaqs.core.data_structures.noise_model import NoiseModel + +if TYPE_CHECKING: + from mqt.yaqs.characterization.memory.backends.tomography.process_tensors import DenseProcessTensor + + +def test_build_process_tensor_defaults_to_mpo() -> None: + """Default return_type is direct MPO construction.""" + ham = Hamiltonian.ising(length=1, J=0.0, g=0.0) + params = AnalogSimParams(dt=0.1, max_bond_dim=8, order=1) + pt = MemoryCharacterizer(parallel=False, show_progress=False).build_process_tensor( + ham, params, timesteps=[0.0, 0.0], compress_every=1 + ) + assert isinstance(pt, MPOProcessTensor) + + +def test_default_mpo_recreates_dense_process_tensor() -> None: + """Uncapped default MPO path matches the dense Choi matrix on a small schedule.""" + ham = Hamiltonian.ising(length=2, J=1.0, g=1.0) + params = AnalogSimParams(dt=0.1, max_bond_dim=8, order=1) + timesteps = [0.1, 0.1, 0.1] + mc = MemoryCharacterizer(parallel=False, show_progress=False) + + pt_mpo = mc.build_process_tensor(ham, params, timesteps=timesteps, max_bond_dim=None, compress_every=1) + pt_dense = mc.build_process_tensor(ham, params, timesteps=timesteps, return_type="dense") + + assert isinstance(pt_mpo, MPOProcessTensor) + np.testing.assert_allclose(pt_mpo.to_matrix(), pt_dense.to_matrix(), atol=1e-8) + np.testing.assert_allclose(pt_mpo.initial_rho, pt_dense.initial_rho, atol=1e-10) + + +@pytest.mark.parametrize("j_val", [0.0, 1.0]) +@pytest.mark.parametrize("num_interventions", [1, 2]) +def test_direct_mpo_matches_dense_temporal_entropy(j_val: float, num_interventions: int) -> None: + """Direct MPO and dense tomography agree on temporal entanglement at small k.""" + ham = Hamiltonian.ising(length=2, J=float(j_val), g=1.0) + params = AnalogSimParams(dt=0.1, max_bond_dim=8, order=1) + timesteps = [0.1] * (num_interventions + 1) + mc = MemoryCharacterizer(parallel=False, show_progress=False) + + pt_dense = cast( + "DenseProcessTensor", + mc.build_process_tensor( + ham, + params, + timesteps=timesteps, + return_type="dense", + num_trajectories=12, + ), + ) + pt_mpo = cast( + "MPOProcessTensor", + build_process_tensor(ham.mpo, params, timesteps=timesteps, max_bond_dim=None, compress_every=1), + ) + + for cut in range(1, num_interventions + 1): + dense = pt_dense.compute_temporal_entropy(cut) + mpo = pt_mpo.compute_temporal_entropy(cut) + assert float(cast("float", mpo["entropy"])) == pytest.approx( + float(cast("float", dense["entropy"])), + abs=1e-6, + ) + assert cast("int", mpo["schmidt_rank"]) == cast("int", dense["schmidt_rank"]) + + +def test_direct_j_zero_matches_dense_temporal_entropy() -> None: + """Direct MPO matches dense tomography for a memoryless chain.""" + ham = Hamiltonian.ising(length=2, J=0.0, g=0.0) + params = AnalogSimParams(dt=0.1, max_bond_dim=8, order=1) + timesteps = [0.0, 0.0, 0.0] + mc = MemoryCharacterizer(parallel=False, show_progress=False) + pt_dense = cast( + "DenseProcessTensor", + mc.build_process_tensor(ham, params, timesteps=timesteps, return_type="dense"), + ) + pt_mpo = cast( + "MPOProcessTensor", + build_process_tensor(ham.mpo, params, timesteps=timesteps, max_bond_dim=None, compress_every=1), + ) + for cut in (1, 2): + dense = pt_dense.compute_temporal_entropy(cut) + mpo = pt_mpo.compute_temporal_entropy(cut) + assert float(cast("float", mpo["entropy"])) == pytest.approx( + float(cast("float", dense["entropy"])), + abs=1e-6, + ) + + +def test_mpo_rejects_noise_model() -> None: + """Direct MPO construction rejects a noise model.""" + ham = Hamiltonian.ising(length=2, J=0.0, g=0.0) + params = AnalogSimParams(dt=0.1, max_bond_dim=8, order=1) + + with pytest.raises(ValueError, match="does not support noise_model"): + build_process_tensor(ham.mpo, params, timesteps=[0.0, 0.0], noise_model=NoiseModel([])) + + +def test_direct_parallel_matches_serial() -> None: + """Parallel and serial direct construction agree on the process-tensor matrix.""" + ham = Hamiltonian.ising(length=2, J=1.0, g=1.0) + params = AnalogSimParams(dt=0.1, max_bond_dim=8, order=1) + timesteps = [0.1, 0.1] + serial = cast( + "MPOProcessTensor", + MemoryCharacterizer(parallel=False, show_progress=False).build_process_tensor( + ham, + params, + timesteps=timesteps, + compress_every=1, + ), + ) + parallel = cast( + "MPOProcessTensor", + MemoryCharacterizer(parallel=True, max_workers=2, show_progress=False).build_process_tensor( + ham, + params, + timesteps=timesteps, + compress_every=1, + ), + ) + np.testing.assert_allclose(parallel.to_matrix(), serial.to_matrix(), atol=1e-8) + + +def test_direct_parallel_temporal_entropy_matches_dense() -> None: + """Parallel direct MPO agrees with dense tomography on temporal entanglement.""" + ham = Hamiltonian.ising(length=2, J=0.0, g=1.0) + params = AnalogSimParams(dt=0.1, max_bond_dim=8, order=1) + timesteps = [0.1, 0.1, 0.1] + mc = MemoryCharacterizer(parallel=True, max_workers=2, show_progress=False) + pt_dense = cast( + "DenseProcessTensor", + mc.build_process_tensor(ham, params, timesteps=timesteps, return_type="dense"), + ) + pt_mpo = cast( + "MPOProcessTensor", + mc.build_process_tensor(ham, params, timesteps=timesteps, max_bond_dim=None, compress_every=1), + ) + for cut in (1, 2): + dense = pt_dense.compute_temporal_entropy(cut) + mpo = pt_mpo.compute_temporal_entropy(cut) + assert float(cast("float", mpo["entropy"])) == pytest.approx( + float(cast("float", dense["entropy"])), + abs=1e-6, + ) + + +def test_direct_tjm_matches_mcwf() -> None: + """Direct MPO construction preserves MPS states under TJM and matches MCWF.""" + ham = Hamiltonian.ising(length=2, J=1.0, g=1.0) + params = AnalogSimParams(dt=0.1, max_bond_dim=8, order=1) + timesteps = [0.1, 0.1] + mcwf = cast( + "MPOProcessTensor", + MemoryCharacterizer(representation="vector", parallel=False, show_progress=False).build_process_tensor( + ham, + params, + timesteps=timesteps, + max_bond_dim=4, + compress_every=1, + ), + ) + tjm = cast( + "MPOProcessTensor", + MemoryCharacterizer(representation="mps", parallel=False, show_progress=False).build_process_tensor( + ham, + params, + timesteps=timesteps, + max_bond_dim=4, + compress_every=1, + ), + ) + np.testing.assert_allclose(tjm.to_matrix(), mcwf.to_matrix(), atol=1e-6) diff --git a/tests/characterization/memory/backends/tomography/test_process_tensors.py b/tests/characterization/memory/backends/tomography/test_process_tensors.py index 87f7e3da..6b27794d 100644 --- a/tests/characterization/memory/backends/tomography/test_process_tensors.py +++ b/tests/characterization/memory/backends/tomography/test_process_tensors.py @@ -17,14 +17,15 @@ import pytest from mqt.yaqs import AnalogSimParams, Hamiltonian, MemoryCharacterizer -from mqt.yaqs.characterization.memory.backends.tomography.data import SequenceData +from mqt.yaqs.characterization.memory.backends.tomography.constructor import build_process_tensor from mqt.yaqs.characterization.memory.backends.tomography.process_tensors import ( DenseProcessTensor, MPOProcessTensor, compute_entropy_dense, + compute_temporal_entropy, convert_probe_callable, encode_cptp_choi, - evaluate_dense_probes, + evaluate_probes, trace_partial_dense, ) from mqt.yaqs.characterization.memory.operational_memory.samples import sample_probes @@ -35,6 +36,40 @@ _REF_RHO0 = np.array([[1.0, 0.0], [0.0, 0.0]], dtype=np.complex128) +def _tiny_mpo_process_tensor(*, num_interventions: int = 1) -> MPOProcessTensor: + """Build a noiseless direct MPO process tensor for wrapper unit tests. + + Returns: + MPO process tensor with the requested number of intervention legs. + """ + ham = Hamiltonian.ising(length=1, J=0.0, g=0.0) + params = AnalogSimParams(dt=0.1, max_bond_dim=8) + timesteps = [0.0] * (num_interventions + 1) + return cast( + "MPOProcessTensor", + build_process_tensor( + ham.mpo, + params, + timesteps=timesteps, + return_type="mpo", + max_bond_dim=None, + compress_every=1, + ), + ) + + +def _single_site_mpo_pt(rho: np.ndarray) -> MPOProcessTensor: + """Wrap a single-site output density matrix as a zero-intervention MPO PT. + + Returns: + MPO process tensor whose only site encodes ``rho``. + """ + tensors = [np.asarray(rho, dtype=np.complex128).reshape(2, 2, 1, 1)] + mpo = MPO() + mpo.custom(tensors, transpose=False) + return MPOProcessTensor(mpo, [], initial_rho=_REF_RHO0.copy()) + + def test_dense_process_tensor_predict_matches_helper() -> None: """DenseProcessTensor._predict_raw matches the Choi contraction; predict physicalizes.""" ups = np.eye(2 * 4, dtype=np.complex128) @@ -99,18 +134,7 @@ def test_mpo_process_tensor_qmi_fallback_to_dense() -> None: def test_mpo_process_tensor_predict_smoke_identity_map() -> None: """MPOProcessTensor.predict returns a physical density matrix for a trivial intervention.""" - rho = np.array([[1.0, 0.0], [0.0, 0.0]], dtype=np.complex128) - data = SequenceData( - sequences=[(0,)], - outputs=[rho], - weights=[1.0], - choi_basis=[np.eye(4, dtype=np.complex128)] * 16, - choi_indices=[(0, 0)] * 16, - choi_duals=[np.eye(4, dtype=np.complex128)] * 16, - timesteps=[0.1], - initial_rho=_REF_RHO0, - ) - pt = data.to_mpo_process_tensor(compress_every=1) + pt = _tiny_mpo_process_tensor(num_interventions=1) def id_map(x: np.ndarray) -> np.ndarray: return x @@ -123,17 +147,7 @@ def id_map(x: np.ndarray) -> np.ndarray: def test_mpo_process_tensor_predict_raises_on_empty_interventions() -> None: """Predict rejects empty interventions when num_interventions>0.""" - data = SequenceData( - sequences=[(0,)], - outputs=[np.eye(2, dtype=np.complex128)], - weights=[1.0], - choi_basis=[np.eye(4, dtype=np.complex128)] * 16, - choi_indices=[(0, 0)] * 16, - choi_duals=[np.eye(4, dtype=np.complex128)] * 16, - timesteps=[0.1], - initial_rho=_REF_RHO0, - ) - pt = data.to_mpo_process_tensor(compress_every=1) + pt = _tiny_mpo_process_tensor(num_interventions=1) with pytest.raises(ValueError, match="interventions list must be non-empty"): pt.predict([]) @@ -141,34 +155,14 @@ def test_mpo_process_tensor_predict_raises_on_empty_interventions() -> None: def test_mpo_process_tensor_predict_zero_steps() -> None: """MPOProcessTensor.predict([]) returns the stored output when num_interventions=0.""" rho = np.array([[0.6, 0.1 + 0.0j], [0.1 - 0.0j, 0.4]], dtype=np.complex128) - data = SequenceData( - sequences=[()], - outputs=[rho], - weights=[1.0], - choi_basis=[], - choi_indices=[], - choi_duals=[], - timesteps=[], - initial_rho=_REF_RHO0, - ) - pt = data.to_mpo_process_tensor(compress_every=1) + pt = _single_site_mpo_pt(rho) rho_out = pt.predict([]) np.testing.assert_allclose(rho_out, rho, atol=1e-10) def test_mpo_process_tensor_predict_raises_on_length_mismatch() -> None: """Predict rejects intervention lists whose length mismatches the process tensor.""" - data = SequenceData( - sequences=[(0,)], - outputs=[np.eye(2, dtype=np.complex128)], - weights=[1.0], - choi_basis=[np.eye(4, dtype=np.complex128)] * 16, - choi_indices=[(0, 0)] * 16, - choi_duals=[np.eye(4, dtype=np.complex128)] * 16, - timesteps=[0.1], - initial_rho=_REF_RHO0, - ) - pt = data.to_mpo_process_tensor(compress_every=1) + pt = _tiny_mpo_process_tensor(num_interventions=1) def id_map(x: np.ndarray) -> np.ndarray: return x @@ -314,14 +308,16 @@ def test_dense_process_tensor_evaluate_probes_smoke() -> None: """Dense process-tensor probe evaluation returns Pauli tomography coefficients.""" pt = _tiny_process_tensor(num_interventions=1) probe_set = sample_probes(cut=1, num_interventions=1, n_pasts=2, n_futures=2, rng=np.random.default_rng(0)) - pauli = evaluate_dense_probes(pt, probe_set) + pauli = evaluate_probes(pt, probe_set) assert pauli.shape == (2, 2, 4) wrapped = pt.evaluate_probes(probe_set) np.testing.assert_allclose(wrapped, pauli) -def test_mpo_process_tensor_evaluate_probes_and_cmi_delegates() -> None: - """MPOProcessTensor wrappers delegate probe and information metrics to dense.""" +def test_mpo_process_tensor_evaluate_probes_matches_dense_without_densifying( + monkeypatch: pytest.MonkeyPatch, +) -> None: + """MPO probe evaluation matches dense and does not call :meth:`to_dense`.""" ham = Hamiltonian.ising(length=1, J=0.0, g=0.0) params = AnalogSimParams(dt=0.1, max_bond_dim=8) mc = MemoryCharacterizer(parallel=False, show_progress=False) @@ -332,9 +328,66 @@ def test_mpo_process_tensor_evaluate_probes_and_cmi_delegates() -> None: dense_pt = mpo_pt.to_dense() probe_set = sample_probes(cut=1, num_interventions=2, n_pasts=2, n_futures=2, rng=np.random.default_rng(1)) + + def _fail_to_dense(self: MPOProcessTensor) -> DenseProcessTensor: + _ = self + msg = "evaluate_probes must not densify the MPO process tensor" + raise AssertionError(msg) + + monkeypatch.setattr(MPOProcessTensor, "to_dense", _fail_to_dense) mpo_pauli = mpo_pt.evaluate_probes(probe_set) + dense_pauli = dense_pt.evaluate_probes(probe_set) assert mpo_pauli.shape == dense_pauli.shape == (2, 2, 4) + np.testing.assert_allclose(mpo_pauli, dense_pauli, atol=1e-6) + # Restore before methods that still densify (qmi/cmi). + monkeypatch.undo() assert isinstance(mpo_pt.cmi(), float) assert mpo_pt._num_interventions_for_probe() == 2 + + +def test_compute_temporal_entropy_markov_j0() -> None: + """Uncoupled Ising process has vanishing temporal entanglement at every cut.""" + ham = Hamiltonian.ising(length=6, J=0.0, g=1.0) + params = AnalogSimParams(dt=0.1, max_bond_dim=64, order=1) + pt = cast( + "DenseProcessTensor", + MemoryCharacterizer(parallel=False, show_progress=False).build_process_tensor( + ham, + params, + timesteps=[0.1] * 4, + return_type="dense", + ), + ) + for cut in (1, 2, 3): + result = pt.compute_temporal_entropy(cut) + assert cast("int", result["schmidt_rank"]) == 1 + assert float(cast("float", result["entropy"])) == pytest.approx(0.0, abs=1e-10) + + +def test_compute_temporal_entropy_correlated_j1() -> None: + """Correlated process has positive temporal entanglement at the center cut.""" + ham = Hamiltonian.ising(length=6, J=1.0, g=1.0) + params = AnalogSimParams(dt=0.1, max_bond_dim=64, order=1) + pt = cast( + "DenseProcessTensor", + MemoryCharacterizer(parallel=False, show_progress=False).build_process_tensor( + ham, + params, + timesteps=[0.1] * 4, + return_type="dense", + ), + ) + result = pt.compute_temporal_entropy(2) + assert cast("int", result["schmidt_rank"]) > 1 + assert float(cast("float", result["entropy"])) > 0.0 + + +def test_compute_temporal_entropy_scale_invariant() -> None: + """Overall scaling of upsilon does not change temporal entanglement.""" + k = 2 + ups = np.eye(2 * 4**k, dtype=np.complex128) + base = float(cast("float", compute_temporal_entropy(ups, k, 1)["entropy"])) + scaled = float(cast("float", compute_temporal_entropy(2.5 * ups, k, 1)["entropy"])) + assert base == pytest.approx(scaled, abs=1e-12) diff --git a/tests/characterization/noise/backends/test_cma.py b/tests/characterization/noise/backends/test_cma.py index 44f6cef9..0ab2897f 100644 --- a/tests/characterization/noise/backends/test_cma.py +++ b/tests/characterization/noise/backends/test_cma.py @@ -125,6 +125,8 @@ def test_backend_exports_cma_opt() -> None: assert callable(cma_opt) +@pytest.mark.filterwarnings("ignore:sigma change np.exp:UserWarning") +@pytest.mark.filterwarnings("ignore:Initial solution argument x0.*:UserWarning") def test_cma_opt_integration_smoke() -> None: """Real CMA-ES backend minimizes a simple quadratic objective.""" pytest.importorskip("cma") @@ -133,16 +135,18 @@ class Objective: def __call__(self, x: np.ndarray) -> float: return float(np.sum(x**2)) + # Mild step size / short run avoids CMA sigma-clip advisories that vary by + # cma/numpy version under pytest's warnings-as-errors policy. xbest, fbest, loss_history, param_history = cma_backend.cma_opt( Objective(), - np.array([1.0, 1.0]), - sigma0=0.2, - max_iter=3, + np.array([0.5, 0.5]), + sigma0=0.05, + max_iter=2, popsize=4, seed=0, ) - assert fbest < 2.0 + assert fbest < 1.0 assert len(loss_history) >= 4 assert len(param_history) == len(loss_history) assert xbest.shape == (2,)