Add dense Haar random state helper

Author: rezaquant

docs/api/package.md

@@ -3,6 +3,7 @@
 Top-level convenience exports.
 
 ```{eval-rst}
+.. autofunction:: pepsy.haar_random_state
 .. autofunction:: pepsy.tn_fidelity
 .. autofunction:: pepsy.tn_norm
 ```

docs/development/package_layout.md

@@ -28,7 +28,7 @@ Use the clearer namespaces for submodule imports:
 from pepsy.boundary import BdyMPS, contract_boundary
 from pepsy.optimizers import SweepOptimizer, MpsOptimizer
 from pepsy.operators import rx, rzz, gate
-from pepsy.tensors import ps_to_peps, ps_to_mps
+from pepsy.tensors import haar_random_state, ps_to_peps, ps_to_mps
 ```
 
 When a leaf module is needed, import the new implementation path directly:

src/pepsy/__init__.py

@@ -104,6 +104,7 @@
     "build_optimizer": ".tensors",
     "contract_hypercompressed_tn": ".tensors",
     "expec_mpo": ".tensors",
+    "haar_random_state": ".tensors",
     "hrps_to_mps": ".tensors",
     "hrps_to_peps": ".tensors",
     "id_to_mpo": ".tensors",
@@ -200,6 +201,7 @@
     "FDSolver",
     "OneDMap",
     "expec_mpo",
+    "haar_random_state",
     "hrps_to_mps",
     "hrps_to_peps",
     "id_to_mpo",
@@ -297,6 +299,7 @@ def __getattr__(name):
     from .tensors import (
         OneDMap,
         expec_mpo,
+        haar_random_state,
         hrps_to_mps,
         hrps_to_peps,
         id_to_mpo,

src/pepsy/tensors/__init__.py

@@ -16,6 +16,7 @@
     expec_mpo,
     get_default_array_backend,
     get_default_grad_backend,
+    haar_random_state,
     hrps_to_mps,
     hrps_to_peps,
     id_to_mpo,
@@ -52,6 +53,7 @@
     "expec_mpo",
     "get_default_array_backend",
     "get_default_grad_backend",
+    "haar_random_state",
     "hrps_to_mps",
     "hrps_to_peps",
     "id_to_mpo",

src/pepsy/tensors/constructors.py

@@ -2,6 +2,7 @@
 
 from .core import (
     add_cycle,
+    haar_random_state,
     hrps_to_mps,
     hrps_to_peps,
     id_to_mpo,
@@ -15,6 +16,7 @@
 
 __all__ = [
     "add_cycle",
+    "haar_random_state",
     "hrps_to_mps",
     "hrps_to_peps",
     "id_to_mpo",

src/pepsy/tensors/core.py

@@ -1,6 +1,7 @@
 """Shared DMRG backend, optimizer, and fidelity helpers."""
 
 import math
+import warnings
 from numbers import Integral
 from string import Formatter
 from typing import Any
@@ -44,6 +45,7 @@
     "ps_to_mps",
     "ps_to_pepo",
     "ps_to_mpo",
+    "haar_random_state",
     "random_haar_qubit",
     "hrps_to_peps",
     "hrps_to_mps",
@@ -1847,6 +1849,74 @@ def random_haar_qubit(seed=None, perturb=0.0):
     return float(theta), float(phi)
 
 
+def haar_random_state(
+    L: int,
+    dtype: str = "complex128",
+    seed=None,
+    L_max: int = 20,
+    as_tensor: bool = False,
+):
+    """Create a dense Haar-random ``L``-qubit state.
+
+    This samples a full Hilbert-space state, so the result is generally
+    entangled. Unlike :func:`hrps_to_mps`, this is not a product-state tensor
+    network: it returns dense amplitudes with ``2**L`` entries.
+
+    Parameters
+    ----------
+    L : int
+        Number of qubits.
+    dtype : str, optional
+        Complex numpy dtype for the returned amplitudes.
+    seed : int | None, optional
+        Seed for deterministic samples.
+    L_max : int, optional
+        Maximum allowed number of qubits. Values above 20 are capped to 20
+        with a warning because this helper constructs a dense state.
+    as_tensor : bool, optional
+        If True, return the amplitudes reshaped as a dense tensor with shape
+        ``(2,) * L``. Otherwise return a dense vector with shape ``(2**L,)``.
+
+    Returns
+    -------
+    numpy.ndarray
+        Normalized dense Haar-random state amplitudes.
+    """
+    if not isinstance(L, Integral) or L < 0:
+        raise ValueError("L must be a non-negative integer.")
+    if not isinstance(L_max, Integral) or L_max < 0:
+        raise ValueError("L_max must be a non-negative integer.")
+    if L_max > 20:
+        warnings.warn(
+            "haar_random_state constructs dense entangled states and is "
+            "intended for L <= 20; capping L_max to 20.",
+            UserWarning,
+            stacklevel=2,
+        )
+        L_max = 20
+    if L > L_max:
+        raise ValueError(
+            "haar_random_state constructs a dense entangled state and only "
+            f"supports L <= L_max (got L={L}, L_max={L_max})."
+        )
+
+    dtype = np.dtype(dtype)
+    if dtype.kind != "c":
+        raise TypeError("dtype must be a complex numpy dtype.")
+
+    real_dtype = np.float32 if dtype == np.dtype("complex64") else np.float64
+    rng = np.random.default_rng(seed)
+    dim = 2 ** int(L)
+    state = rng.normal(size=dim).astype(real_dtype)
+    state = state + 1j * rng.normal(size=dim).astype(real_dtype)
+    state = state.astype(dtype, copy=False)
+    state /= np.linalg.norm(state)
+
+    if as_tensor:
+        return state.reshape((2,) * int(L))
+    return state
+
+
 def hrps_to_peps(
     Lx: int,
     Ly: int,

tests/test_package_layout.py

@@ -9,7 +9,13 @@
 from pepsy.optimizers import MpsOptimizer, SweepOptimizer
 from pepsy.sampling import MpsSampler, PepsBpSampler
 from pepsy.solvers import FDSolver
-from pepsy.tensors import OneDMap, backend_torch, ps_to_peps, reg_complex_svd_torch
+from pepsy.tensors import (
+    OneDMap,
+    backend_torch,
+    haar_random_state,
+    ps_to_peps,
+    reg_complex_svd_torch,
+)
 
 
 def test_new_namespace_imports_resolve():
@@ -25,6 +31,7 @@ def test_new_namespace_imports_resolve():
     assert FDSolver is not None
     assert OneDMap is not None
     assert callable(backend_torch)
+    assert callable(haar_random_state)
     assert callable(ps_to_peps)
     assert callable(reg_complex_svd_torch)
 

tests/test_public_api.py

@@ -21,7 +21,7 @@ def test_package_version_available():
     "pauli", "x", "y", "z", "s", "sdg", "t", "tdg", "h", "hadamard",
     "cnot", "cx", "cy", "cz", "swap", "iswap", "phase", "u1", "u2",
     "cphase", "crx", "cry", "crz", "cu1", "cu2", "cu3", "rx", "ry", "rz",
-    "rxx", "ryy", "rzz", "u3", "su4", "fsim", "fsimg", "ps_to_peps", "expec_mpo",
+    "rxx", "ryy", "rzz", "u3", "su4", "fsim", "fsimg", "haar_random_state", "ps_to_peps", "expec_mpo",
     "id_to_mpo", "id_to_pepo", "ps_to_pepo", "ps_to_mpo", "make_numpy_array_caster", "SweepOptimizer",
     "FDSolver", "MpsOptimizer", "MpoOptimizer", "PEPSSampleResult",
     "PepsBpSampler", "MpsSampler", "MpsSampleResult", "VecSampler", "gate", "tn_fidelity", "tn_norm",
@@ -64,7 +64,7 @@ def test_internal_symbol_not_exported(name):
     "x", "y", "z", "s", "sdg", "t", "tdg", "h", "hadamard",
     "cnot", "cx", "cy", "cz", "swap", "iswap", "phase", "u1", "u2",
     "cphase", "crx", "cry", "crz", "cu1", "cu2", "cu3", "rx", "ry", "rz",
-    "rxx", "ryy", "rzz", "u3", "su4", "fsim", "fsimg", "ps_to_peps", "expec_mpo",
+    "rxx", "ryy", "rzz", "u3", "su4", "fsim", "fsimg", "haar_random_state", "ps_to_peps", "expec_mpo",
     "id_to_mpo", "id_to_pepo", "ps_to_pepo", "ps_to_mpo", "SweepOptimizer",
     "FDSolver", "MpsOptimizer", "MpoOptimizer", "PEPSSampleResult", "PepsBpSampler",
     "tn_fidelity", "tn_norm",

tests/test_tensor_constructors.py

@@ -0,0 +1,57 @@
+"""Tests for tensor and dense-state constructors."""
+
+import numpy as np
+import pytest
+
+from pepsy import haar_random_state
+from pepsy.tensors.constructors import haar_random_state as constructors_haar_random_state
+
+
+def test_haar_random_state_returns_normalized_dense_vector():
+    """Dense Haar state should have full Hilbert-space shape and unit norm."""
+    state = haar_random_state(3, seed=123)
+
+    assert state.shape == (8,)
+    assert state.dtype == np.dtype("complex128")
+    assert np.isclose(np.linalg.norm(state), 1.0)
+
+
+def test_haar_random_state_tensor_shape_matches_vector_seed():
+    """The tensor form should be a reshape of the seeded dense vector."""
+    vector = haar_random_state(4, seed=7)
+    tensor = haar_random_state(4, seed=7, as_tensor=True)
+
+    assert tensor.shape == (2, 2, 2, 2)
+    assert np.allclose(tensor.reshape(-1), vector)
+
+
+def test_haar_random_state_sample_is_generally_entangled():
+    """A full Haar sample should not be restricted to product-state rank."""
+    state = haar_random_state(3, seed=11)
+
+    assert np.linalg.matrix_rank(state.reshape(2, -1), tol=1e-12) > 1
+
+
+def test_haar_random_state_rejects_large_dense_state():
+    """The default guard should keep dense state construction small."""
+    with pytest.raises(ValueError, match="L <= L_max"):
+        haar_random_state(21)
+
+
+def test_haar_random_state_caps_lmax_above_twenty():
+    """Raising L_max above the documented dense-state limit should warn."""
+    with pytest.warns(UserWarning, match="L <= 20"):
+        state = haar_random_state(2, seed=5, L_max=21)
+
+    assert state.shape == (4,)
+
+
+def test_haar_random_state_rejects_real_dtype():
+    """Haar-random qubit amplitudes require a complex dtype."""
+    with pytest.raises(TypeError, match="complex"):
+        haar_random_state(2, dtype="float64")
+
+
+def test_haar_random_state_constructor_facade_resolves():
+    """Constructor facade should expose the dense Haar helper."""
+    assert constructors_haar_random_state is haar_random_state