staging > main: refactored entry point for application

justfile

@@ -103,8 +103,8 @@ build-requirements-dependencies:
 # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
 [group("exec")]
-run *args:
-    @{{PYVENV_ON}} && {{PYVENV}} -m src.main {{args}}
+run-examples *args:
+    @{{PYVENV_ON}} && {{PYVENV}} -m src.examples {{args}}
 
 # --------------------------------
 # TARGETS: development + tools

src/__paths__.py

@@ -0,0 +1,45 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+
+# ----------------------------------------------------------------
+# IMPORTS
+# ----------------------------------------------------------------
+
+from pathlib import Path
+
+# ----------------------------------------------------------------
+# EXPORTS
+# ----------------------------------------------------------------
+
+__all__ = [
+    "get_root_path",
+    "get_source_path",
+]
+
+# ----------------------------------------------------------------
+# CONSTANTS
+# ----------------------------------------------------------------
+
+_source = Path(__file__).parent.as_posix()
+_root = Path(__file__).parent.parent.as_posix()
+
+# ----------------------------------------------------------------
+# METHODS
+# ----------------------------------------------------------------
+
+
+def get_root_path(*parts: str) -> str:
+    return get_path(_root, *parts)
+
+
+def get_source_path(*parts: str) -> str:
+    return get_path(_source, *parts)
+
+
+# ----------------------------------------------------------------
+# AUXILIARY METHODS
+# ----------------------------------------------------------------
+
+
+def get_path(root: str, *parts: str) -> str:
+    return root if len(parts) == 0 else Path(root, *parts).as_posix()

src/examples.py

@@ -0,0 +1,63 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+
+"""
+- author:      Raj Dahya
+- institute:   Universität Leipzig
+- department:  Fakultät für Mathematik und Informatik
+- created:     2025-12-27
+- updated:     2026-03-29
+- description:
+    Example script to numerically verify methods presented in paper
+"""
+
+# ----------------------------------------------------------------
+# IMPORTS
+# ----------------------------------------------------------------
+
+import os
+import sys
+from pathlib import Path
+
+os.chdir(Path(__file__).parent.parent.as_posix())
+sys.path.insert(0, os.getcwd())
+
+
+import torch
+
+from ._core.logging import *
+from .experiments.choi_cholesky import *
+from .queries._console.prog_examples import *
+from .setup import *
+
+# ----------------------------------------------------------------
+# EXECUTION
+# ----------------------------------------------------------------
+
+if __name__ == "__main__":
+    # set sessions settings
+    sys.tracebacklimit = 0
+
+    # parse cli flags
+    args = CliArguments(info=INFO).parse(*sys.argv[1:])
+
+    # early terminate if only version requested
+    if args.mode == "version":
+        print(VERSION)
+        exit(0)
+
+    # set up logging
+    configure_logging()
+
+    # optionally set seed
+    if args.seed is not None:
+        torch.manual_seed(args.seed)
+
+    # execute main method
+    verify_choi_cholesky(
+        d1=args.dim1,
+        d2=args.dim2,
+        num_experiments=args.num,
+        map_kind=args.map,
+        algorithm_choice=args.algorithm,
+    )

src/main.py

@@ -1,461 +0,0 @@
-#!/usr/bin/env python3
-# -*- coding: utf-8 -*-
-
-"""
-- author:      Raj Dahya
-- institute:   Universität Leipzi
-- department:  Fakultät für Mathematik und Informatik
-- created:     2025-12-27
-- updated:     2026-03-29
-- title:       python script for experimental verification
-- description:
-    Verifies the Choi-Cholesky decomposition of a CP/CPTP map numerically,
-    as presented in the paper _The Choi-Cholesky algorithm for completely positive maps_
-    (see <https://arxiv.org/abs/2603.19444>).
-
-    Verification occurs by running the experiment multiple times without failure.
-
-    Failure criteria:
-
-    - [critical] algorithm fails to produce positive entries for the diagonal operator (within tolerance)
-    - in final result L·L^* is not close to Choi matrix
-
-    NOTE: Due to the discontinuous nature of pseudo-inverses (cf Remark 3.3 in §3),
-    the algorithm can get close to a critical error under low numerical precision (e.g. Float32).
-    For this reason we use double precision (Float64) numbers.
-
-    NOTE: The option `MAP_TYPE = "CP"` establishes a proper confirmation.
-    For completeness we have included the option `"CPTP"`,
-    but this provides a circular test, as the generation of CPTP-maps here assumes results of paper in advance.
-
-    NOTE: use `MODE="BASIC"` to run the algorithm exactly as presented in the paper (cf Algorithm B.1)
-    and use `MODE="SIMPLIFIED"` to run an equivalent which has been modified for numerical stability.
-"""
-
-# ----------------------------------------------------------------
-# USER ADJUSTABLE SETTINGS
-# ----------------------------------------------------------------
-
-MAP_TYPE = "CP"  # non-circular test
-# MAP_TYPE = "CPTP" # NOTE: circular test - see above
-NUM_EXPERIMENTS = 100  # number of randomly generated CP-maps on which to verify algorithm
-N1 = 10  # dimension of H1
-N2 = 20  # dimension of H2
-
-# choice of algorithm
-MODE="BASIC"      # as in paper
-# MODE="SIMPLIFIED" # algorithm equivalent to original but modified for numerical stability
-
-# ----------------------------------------------------------------
-# IMPORTS
-# ----------------------------------------------------------------
-
-import os
-import sys
-from pathlib import Path
-
-os.chdir(Path(__file__).parent.as_posix())
-sys.path.insert(0, os.getcwd())
-
-import logging
-import math
-from logging import getLogger
-from typing import Literal
-
-import scipy.linalg
-import torch
-from tabulate import tabulate
-
-# ----------------------------------------------------------------
-# SETTINGS, CONSTANTS
-# ----------------------------------------------------------------
-
-TOL = 1e-6
-# PRECISION = torch.float32 # works but to lower precision
-PRECISION = torch.float64
-# PRECISION = torch.double # same as f64
-
-# ----------------------------------------------------------------
-# METHODS
-# ----------------------------------------------------------------
-
-
-def main(
-    *,
-    map_type: Literal["CP", "CPTP"] = "CP",
-    num_experiments: int = 1,
-    mode: Literal["BASIC", "SIMPLIFIED"] = "BASIC",
-):
-    diffs: list[float] = []
-    for index in range(1, 1 + num_experiments):
-        logging.info(f"run experiment {index}")
-
-        """
-        Create a random CP/CPTP-map
-        """
-        match map_type:
-            case "CPTP":
-                logging.info("generate Choi matrix C of a random CPTP map Φ")
-                C = random_cptp_choicholesky(N1, N2, precision=PRECISION)
-                diff = torch.norm(trace_2(C) - torch.eye(N1)).item()
-                logging.info(f"generated positive C with ‖tr₂(C) – I‖ = {diff:.4g}")
-
-            case _:
-                logging.info("generate Choi matrix C of a random CP map Φ")
-                C = random_cp(N1, N2, precision=PRECISION)
-
-        """
-        Perform Choi-Cholesky decomposition
-        """
-        logging.info("compute Choi-Cholesky decomposition for Φ")
-        match mode:
-            case "SIMPLIFIED":
-                # method equivalent to algorithm in paper, modified for numerical stability
-                L = algorithm_bipartite_cholesky_modified(C)
-
-            # case "BASIC":
-            case _:
-                # method as in paper
-                L = algorithm_bipartite_cholesky(C)
-
-        """
-        Verify validity of Choi-Cholesky decomposition
-        """
-        C_recovered = tensor_multiply(L, tensor_conj(L))
-        diff = tensor_diff(C, C_recovered, rel=True)
-        diffs.append(diff)
-
-    """
-    Perform statistical evaluation of success rate
-    """
-
-    diffs = torch.asarray(diffs)
-    levels = [0.25, 0.5, 0.75, 0.975]
-    quantiles = torch.quantile(diffs, q=torch.asarray(levels)).tolist()
-    table = tabulate(
-        zip(levels, quantiles),
-        headers=("quantile", "Diff"),
-        tablefmt="rst",
-        colalign=("right", "right"),
-        floatfmt=(".2%", ".3e"),
-    )
-    logging.info(f"quantiles for rel. differences ‖L·L^* – C(Φ)‖/‖C(Φ)‖:\n{table}")
-    diff = diffs[-1]  # get 97.5% quantile of rel differences
-    assert diff < 0.5*TOL, \
-        f"Algorithm failed to produce decomposition that recovers original Choi matrix within tolerance {TOL:.4g}."  # fmt: skip
-    logging.info(f"Algorithm succeeded in producing decomposition that recovers original Choi matrix within tolerance {TOL:.4g}.")  # fmt: skip
-    pass
-
-
-# ----------------------------------------------------------------
-# ALGORITHMS
-# ----------------------------------------------------------------
-
-
-def algorithm_bipartite_cholesky(
-    C: torch.Tensor,
-    /,
-) -> tuple[
-    torch.Tensor,
-    torch.Tensor,
-]:
-    """
-    Bi-partite Cholesky decomposition
-    """
-    n1, _, n2, _ = C.shape
-    Eye = tensor_identity(n1, n2).type_as(C)
-    Zero = torch.zeros(C.shape).type_as(C)
-    unitL = Eye.clone()
-    unitR = Eye.clone()  # NOTE: unitR will compute unitL^-1
-    L = Zero.clone()  # NOTE: L = L · √D
-    D = Zero.clone()
-    D_pinv = Zero.clone()  # NOTE: will compute D^†
-
-    for i in range(n1):
-        # update row i of unitL, L, D
-        D[i, i, :, :] = C[i, i, :, :]
-        for j in range(i):
-            for k in range(i):
-                unitL[i, j, :, :] += C[i, k, :, :] @ matrix_conj(unitR[j, k, :, :])
-
-            unitL[i, j, :, :] = unitL[i, j, :, :] @ D_pinv[j, j, :, :]
-            L[i, j, :, :] = unitL[i, j, :, :] @ L[j, j, :, :]
-            D[i, i, :, :] = D[i, i, :, :] - L[i, j, :, :] @ matrix_conj(L[i, j, :, :])
-
-        # NOTE: due to numerical imperfections, raise error if not positive within tolerance, otherwise coerce to positive
-        D[i, i, :, :] = assert_matrix_positive(D[i, i, :, :])
-        L[i, i, :, :] = matrix_sqrt(D[i, i, :, :])
-
-        # update row i of pseudo-inverse of D
-        D_pinv[i, i, :, :] = matrix_pinv(D[i, i, :, :])
-
-        # update row i of inverse of unitL
-        for j in range(i):
-            for k in range(j, i):
-                unitR[i, j, :, :] -= unitL[i, k, :, :] @ unitR[k, j, :, :]
-
-    return L
-
-
-def algorithm_bipartite_cholesky_modified(
-    C: torch.Tensor,
-    /,
-) -> tuple[
-    torch.Tensor,
-    torch.Tensor,
-]:
-    """
-    Bi-partite Cholesky decomposition with improved numerical stability
-
-    This is equivalent to the algorithm presented in the paper,
-    modified by keeping track of the pseudo-inverse of L instead of the inverse of unitL.
-
-    NOTE: This algorithm has been obtained by simple change of variables,
-    and by observing that
-
-    ```py
-    unitL[i, j] = ... · D[j, j]
-    L[j, j] = √D[j, j]
-    ```
-
-    whence
-
-    ```py
-    unitL[i, j] = unitL[i, j] · L[j, j] · √D[j, j]
-    ```
-
-    must hold.
-    """
-    n1, _, _, _ = C.shape
-    Zero = torch.zeros(C.shape).type_as(C)
-    L = Zero.clone()  # NOTE: L = L · √D
-    R = Zero.clone()  # NOTE: will compute R := L^† = √D^† · hatL^-1
-    D = Zero.clone()
-
-    for i in range(n1):
-        # update row i of L, D
-        D[i, i, :, :] = C[i, i, :, :]
-        for j in range(i):
-            for k in range(i):
-                L[i, j, :, :] += C[i, k, :, :] @ matrix_conj(R[j, k, :, :])
-
-            D[i, i, :, :] = D[i, i, :, :] - L[i, j, :, :] @ matrix_conj(L[i, j, :, :])
-
-        # NOTE: due to numerical imperfections, raise error if not positive within tolerance, otherwise coerce to positive
-        D[i, i, :, :] = assert_matrix_positive(D[i, i, :, :])
-        L[i, i, :, :] = matrix_sqrt(D[i, i, :, :])
-
-        # update row i of R
-        R[i, i, :, :] = matrix_pinv(L[i, i, :, :])
-        for j in range(i):
-            for k in range(j, i):
-                R[i, j, :, :] -= L[i, k, :, :] @ R[k, j, :, :]
-
-            R[i, j, :, :] = R[i, i, :, :] @ R[i, j, :, :]
-
-    return L
-
-
-# ----------------------------------------------------------------
-# TERTIARY METHODS
-# ----------------------------------------------------------------
-
-
-def random_cp(
-    d1: int,
-    d2: int,
-    /,
-    *,
-    precision: torch.dtype = torch.double,
-) -> torch.Tensor:
-    """
-    Produces the Choi-matrix of a 'random' (contractive) CP-map
-    """
-    while True:
-        U = torch.randn(size=(d1, d1, d2, d2), dtype=precision)
-        C = tensor_multiply(U, tensor_conj(U))
-        scale = scipy.linalg.norm(C)
-        if scale > 0:
-            C = C / scale
-            break
-
-    return C
-
-
-def random_cptp_choicholesky(
-    d1: int,
-    d2: int,
-    /,
-    *,
-    precision: torch.dtype = torch.double,
-) -> torch.Tensor:
-    """
-    Uses the resolution lemma (cf. Lemma 2.10 in §2.5)
-    to determine the Choi matrix of a 'random' CPTP-map
-
-    NOTE: The other way to generate Choi matrices of random CPTP-maps
-    is to keep iterating over randomly generated matrices until the 2nd partial trace
-    equals the identity (cf. Table 1 in §2.1).
-    However, this is verify inefficient.
-    """
-    while True:
-        # create random operators
-        U = torch.randn(size=(d1, d1, d2, d2), dtype=precision)
-
-        # use Gram-Schmidt to ensure that rows are tr-orthonormal
-        success = True
-        for i in range(d1):
-            u = U[i, :, :, :]
-
-            # ensure orthogonality
-            for ii in range(i):
-                v = U[ii, :, :, :]
-                scale = torch.einsum(r"jrs,jsr", u, tensor_conj_2(v))
-                u = u - scale * v
-
-            # ensure normality
-            scale = torch.einsum(r"jrs,jsr", u, tensor_conj_2(u))
-            if not (scale > 0):
-                success = False
-                break
-
-            # update row i of bi-partite lower triangular operator
-            u /= math.sqrt(scale)
-            U[i, :, :, :] = u
-
-        if success:
-            break
-
-    C = tensor_multiply(U, tensor_conj(U))
-    return C
-
-
-# ----------------------------------------------------------------
-# AUXILIARY METHODS
-# ----------------------------------------------------------------
-
-
-def tensor_identity(d1: int, d2: int, /) -> torch.Tensor:
-    return torch.einsum(r"ij,kl->ijkl", torch.eye(d1), torch.eye(d2))
-
-
-def tensor_multiply(A: torch.Tensor, B: torch.Tensor, /) -> torch.Tensor:
-    return torch.einsum(r"ijrs,jkst->ikrt", A, B)
-
-
-def tensor_conj(A: torch.Tensor, /) -> torch.Tensor:
-    return A.conj().transpose(0, 1).transpose(2, 3)
-
-
-def tensor_conj_1(A: torch.Tensor, /) -> torch.Tensor:
-    return A.conj().transpose(0, 1)
-
-
-def tensor_conj_2(A: torch.Tensor, /) -> torch.Tensor:
-    k = len(A.shape)
-    return A.conj().transpose(-2, -1)
-
-
-def matrix_conj(A: torch.Tensor, /) -> torch.Tensor:
-    return tensor_conj_1(A)
-
-
-def matrix_multiply(A: torch.Tensor, B: torch.Tensor, /) -> torch.Tensor:
-    return torch.einsum(r"ik,kj->ij", A, B)
-
-
-def assert_matrix_positive(
-    A: torch.Tensor,
-    /,
-    *,
-    tol: float = TOL,
-) -> torch.Tensor:
-    """
-    Checks if self-adjoint matrix is positive within tolerance and forces it to be so
-    """
-    # check if positive
-    values, U = torch.linalg.eigh(A)
-    eig_min = min(values).item()
-
-    if eig_min < -0.5 * tol:
-        raise ValueError(f"operator is not positive within tolerance {tol:.4g}: min(σ(A)) = {eig_min:.4g}")  # fmt: skip
-
-    if eig_min < 0:
-        logging.warning(f"operator has negative eigenvales but is positive within tolerance {tol:.4g}: min(σ(A)) = {eig_min:.4g}")  # fmt: skip
-
-    # force positivity by setting negative eigenvalues to 0
-    values = torch.maximum(values, torch.zeros_like(values))
-    sqrtD = torch.diag(torch.sqrt(values))
-    sqrtA = U @ sqrtD
-    A = sqrtA @ matrix_conj(sqrtA)
-    return A
-
-
-def matrix_sqrt(A: torch.Tensor, /) -> torch.Tensor:
-    """
-    Computes the sqrt of a positive semi-definite matrix
-    """
-    dtype = A.dtype
-    sqrtA = scipy.linalg.sqrtm(A)
-    sqrtA = torch.asarray(sqrtA, dtype=dtype)
-    return sqrtA
-
-
-def matrix_pinv(A: torch.Tensor, /) -> torch.Tensor:
-    """
-    Computes the Moore-Penrose pseudo-inverse of a semi-definite matrix
-    """
-    A_pinv = torch.pinverse(A)
-    return A_pinv
-
-
-def trace_2(A: torch.Tensor, /) -> torch.Tensor:
-    return torch.einsum(r"ijkk->ij", A)
-
-
-def tensor_diff(
-    X: torch.Tensor,
-    Y: torch.Tensor,
-    /,
-    *,
-    rel: bool = True,
-) -> float:
-    """
-    Computes the relative difference between two tensors
-    """
-    eps = torch.linalg.norm(Y - X)
-    if rel:
-        scale = max(torch.linalg.norm(X), torch.linalg.norm(Y)) or 1.0
-        eps = eps / scale
-    return float(eps)
-
-
-# ----------------------------------------------------------------
-# EXECUTION
-# ----------------------------------------------------------------
-
-if __name__ == "__main__":
-    """
-    Sessions settings including optional seeding
-    """
-    sys.tracebacklimit = 0
-    args = sys.argv[1:]
-    if len(args) > 0:
-        seed = args[0]
-        torch.manual_seed(seed)
-
-    """
-    Set up logging
-    """
-
-    logging.basicConfig(
-        format="%(asctime)s $\x1b[92;1m%(name)s\x1b[0m [\x1b[1m%(levelname)s\x1b[0m] %(message)s",
-        datefmt=r"%Y-%m-%d %H:%M:%S",
-        encoding="utf-8",
-    )
-    getLogger(name="root").setLevel(logging.INFO)
-
-    """
-    Main execution
-    """
-    main(num_experiments=NUM_EXPERIMENTS, map_type=MAP_TYPE, mode=MODE)