diff --git a/docs/examples_dilations.ipynb b/docs/examples_dilations.ipynb
new file mode 100644
index 0000000..1a0df70
--- /dev/null
+++ b/docs/examples_dilations.ipynb
@@ -0,0 +1,229 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# Examples - non-dilatable $d$-parameter $C_{0}$-semigroups #\n",
+ "\n",
+ "This notebook provides supplementary material to the research paper <https://doi.org/10.48550/arXiv.2210.02353>.\n",
+ "\n",
+ "In §5.3 a construction is given which yields for any $d \\geq 2$\n",
+ "a $d$-parameter $C_{0}$-semigroup, $T$ satisfying:\n",
+ "\n",
+ "- $T$ is contractive (equivalently its marginals $T_{i}$ are for each $i$);\n",
+ "- the generator $A_{i}$ of $T_{i}$ has strictly negative spectral bound for each $i$;\n",
+ "\n",
+ "and such that $T$ is __not__ **completely dissipative**.\n",
+ "Thus by the classification Theorem (Thm 1.1), $T$ does not have a regular unitary dilation.\n",
+ "\n",
+ "This Notebook demonstrates this general result empirically."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 12,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import os;\n",
+ "import sys;\n",
+ "\n",
+ "# NOTE: need this to force jupyter to reload imports:\n",
+ "for key in list(sys.modules.keys()):\n",
+ " if key.startswith('src.'):\n",
+ " del sys.modules[key];\n",
+ "\n",
+ "os.chdir(os.path.dirname(_dh[0]));\n",
+ "sys.path.insert(0, os.getcwd());\n",
+ "\n",
+ "from src.examples_dilations import *;"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 13,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# User input:\n",
+ "N = 4; # dimension of the Hilbert space.\n",
+ "d = 4;\n",
+ "\n",
+ "# If you ensure that the failure of S_{T,K} >= 0 only occurs for K = {1,2,...,d}\n",
+ "# + you want S_TK > 0 (strictly) for all K ≠ {1,2,...,d}:\n",
+ "alpha = 1/math.sqrt(d - 0.5);\n",
+ "\n",
+ "# If you ensure that the failure of S_{T,K} >= 0 only occurs for K = {1,2,...,d}\n",
+ "# alpha = 1/math.sqrt(d - 1);\n",
+ "\n",
+ "# Otherwise:\n",
+ "# alpha = 1;"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 14,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "\n",
+ "The marginal semigroups T_i and their generators A_i:\n",
+ "\n",
+ " i spec bound of A_i A_i dissipative (<==> T_i contractive)?\n",
+ "--- ------------------- -----------------------------------------\n",
+ " 1 -0.465478 True\n",
+ " 2 -0.465478 True\n",
+ " 3 -0.465478 True\n",
+ " 4 -0.465478 True\n"
+ ]
+ }
+ ],
+ "source": [
+ "# create the generators `A_i` of the marginal semigroups `T_i`:`\n",
+ "A = [\n",
+ " generate_semigroup_generator(\n",
+ " shape = [N, N],\n",
+ " rational = True,\n",
+ " base = 100,\n",
+ " alpha = alpha,\n",
+ " )\n",
+ " for _ in range(d)\n",
+ "];\n",
+ "\n",
+ "data = [];\n",
+ "for i, A_i in enumerate(A):\n",
+ " omega_Re = spec_bounds((1/2)*(A_i + A_i.T.conj()));\n",
+ " omega = spec_bounds(A_i);\n",
+ " data.append((i+1, omega_Re, True if (omega_Re <= 0) else False));\n",
+ "\n",
+ "repr = tabulate(\n",
+ " tabular_data = data,\n",
+ " headers = ['i', 'spec bound of A_i', 'A_i dissipative (<==> T_i contractive)?'],\n",
+ " showindex = False,\n",
+ " floatfmt = '.6f',\n",
+ " colalign = ['center', 'center', 'center'],\n",
+ " tablefmt = 'simple',\n",
+ ");\n",
+ "print(f'\\nThe marginal semigroups T_i and their generators A_i:\\n\\n{repr}');"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 15,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "\n",
+ "Dissipation operators:\n",
+ "\n",
+ " K min σ(S_{T,K}) S_{T,K} >= 0?\n",
+ "------------ ---------------- ---------------\n",
+ " [] 1.000000 True\n",
+ " [3] 0.465478 True\n",
+ " [2] 0.465478 True\n",
+ " [1] 0.465478 True\n",
+ " [0] 0.465478 True\n",
+ " [2, 3] 0.212675 True\n",
+ " [1, 3] 0.183738 True\n",
+ " [1, 2] 0.217453 True\n",
+ " [0, 3] 0.200206 True\n",
+ " [0, 2] 0.301332 True\n",
+ " [0, 1] 0.215681 True\n",
+ " [1, 2, 3] 0.058427 True\n",
+ " [0, 2, 3] 0.075037 True\n",
+ " [0, 1, 3] 0.056030 True\n",
+ " [0, 1, 2] 0.077350 True\n",
+ "[0, 1, 2, 3] -0.082403 False\n"
+ ]
+ }
+ ],
+ "source": [
+ "# compute the dissipation operators `S_TK` for each `K ⊆ {1,2,...,d}``:\n",
+ "S, beta_T = dissipation_operators(shape=[N, N], A=A);\n",
+ "\n",
+ "data = [];\n",
+ "for K, S_TK, b in sorted(S, key=lambda x: len(x[0])):\n",
+ " data.append((K, b, True if b >= -MACHINE_EPS else False))\n",
+ "\n",
+ "repr = tabulate(\n",
+ " tabular_data = data,\n",
+ " headers = ['K', 'min σ(S_{T,K})', 'S_{T,K} >= 0?'],\n",
+ " showindex = False,\n",
+ " floatfmt = '.6f',\n",
+ " colalign = ['center', 'center', 'center'],\n",
+ " tablefmt = 'simple',\n",
+ ");\n",
+ "print(f'\\nDissipation operators:\\n\\n{repr}');"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 16,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "\n",
+ "β_T = min_K min σ(S_{T,K}) = -0.082403\n",
+ "⟹ T is not compeletely dissipative.\n",
+ "⟹ (by Thm 1.1) T does not have a regular unitary dilation.\n"
+ ]
+ }
+ ],
+ "source": [
+ "# Display summary:\n",
+ "print('')\n",
+ "print(f'β_T = min_K min σ(S_{{T,K}}) = {beta_T:.6f}');\n",
+ "if beta_T == 0:\n",
+ " print('⟹ T is compeletely dissipative.');\n",
+ " print('⟹ (by Thm 1.1) T has a regular unitary dilation.');\n",
+ "elif beta_T > 0:\n",
+ " print('⟹ T is compeletely super dissipative.');\n",
+ " print('⟹ (by Thm 1.1) T has a regular unitary dilation.');\n",
+ "else:\n",
+ " print('⟹ T is not compeletely dissipative.');\n",
+ " print('⟹ (by Thm 1.1) T does not have a regular unitary dilation.');"
+ ]
+ }
+ ],
+ "metadata": {
+ "jupytext": {
+ "cell_metadata_filter": "-all"
+ },
+ "kernelspec": {
+ "display_name": "Python 3.10.6 64-bit",
+ "language": "python",
+ "name": "python3"
+ },
+ "language_info": {
+ "codemirror_mode": {
+ "name": "ipython",
+ "version": 3
+ },
+ "file_extension": ".py",
+ "mimetype": "text/x-python",
+ "name": "python",
+ "nbconvert_exporter": "python",
+ "pygments_lexer": "ipython3",
+ "version": "3.10.6"
+ },
+ "widgets": {
+ "application/vnd.jupyter.widget-state+json": {
+ "state": {},
+ "version_major": 2,
+ "version_minor": 0
+ }
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}
diff --git a/src/__init__.py b/src/__init__.py
new file mode 100644
index 0000000..e69de29
diff --git a/src/examples_dilations.py b/src/examples_dilations.py
new file mode 100644
index 0000000..7db2272
--- /dev/null
+++ b/src/examples_dilations.py
@@ -0,0 +1,288 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+# AUTHOR: Raj Dahya
+# CREATED: 05.09.2022
+# EDITED: 08.10.2022
+# TYPE: Source code
+# DEPARTMENT: Fakult\"at for Mathematik und Informatik
+# INSTITUTE: Universit\"at Leipzig
+# NOTE: Collapse this block.
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+# IMPORTS
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+import re;
+import math
+from tkinter import E;
+import numpy as np;
+from typing import Generator;
+import itertools;
+from tabulate import tabulate;
+
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+# CONSTANTS + SETTINGS
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+np.set_printoptions(precision=3);
+MACHINE_EPS = 0.5e-12;
+
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+# SECONDARY METHODS - Generate C0-Semigroup Generator
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+def generate_semigroup_generator(
+ shape: list[int, int],
+ rational: bool = False,
+ base: int = 2,
+ alpha: float = 1,
+) -> np.ndarray:
+ '''
+ Creates bounded generators for d commuting C0-semigroups as in the paper (cf. §5.3).
+
+ @inputs
+ - `shape` - <int, int> Desired shape of output.
+ - `rational` - <bool> Whether or not entries are to be rational.
+ - `base` - <int> If `rational = True`, fixes the denominator of the rational numbers.
+ - `alpha` - <float> Additional parameter to scale the D_i operators.
+
+ NOTE: in the paper `α = 1` was chosen. However one can choose any value in `(1/√d, \infty)`.
+ By choosing any value `α ≥ 1/√(d-1)`, by the computations in Proposition 5.3
+ one can force that the S_TK operators only fail to be positive when |K| > d-1.
+
+ @returns
+ Randomly created bounded generators for C0-semigroups as in the paper (cf. §5.3).
+ '''
+ sh = tuple(shape)
+ m, n = sh;
+ assert m == n, 'must be square';
+ assert m % 2 == 0, 'dimensions must be divisible by 2';
+ N = int(m/2);
+ I = generate_op_identity(shape=[2*N, 2*N]);
+ O = generate_op_zero(shape=[N, N]);
+ U = generate_op_unitary(
+ shape = [N, N],
+ rational = rational,
+ base = base,
+ );
+ D = alpha * U;
+ A = -I + np.vstack([
+ np.hstack([ O, -2*D ]),
+ np.hstack([ O, O ]),
+ ]);
+ return A;
+
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+# SECONDARY METHODS - Compute spectral bounds
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+def spec_bounds(A: np.ndarray) -> float:
+ '''
+ Computes the spectral bounds of an operator:
+ ```tex
+ \max_{i} \max \{ Re \lambda : \lambda \in \sigma(A_{i}) \}
+ ```
+ '''
+ sigma_A = np.linalg.eigvals(A);
+ return max(np.real(sigma_A));
+
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+# SECONDARY METHODS - Generate Dissipation Operator
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+def dissipation_operators(
+ shape: list[int, int],
+ A: list[np.ndarray],
+) -> tuple[
+ list[tuple[
+ list[int],
+ np.ndarray,
+ float,
+ ]],
+ float,
+]:
+ '''
+ Computes the dissipation operators and their minimum spectral values.
+ '''
+ d = len(A);
+ indexes = list(range(d));
+ data = [];
+ beta = np.inf;
+ for K in power_set(indexes):
+ S_TK, eig_K = dissipation_operator(shape=shape, A=A, K=K);
+ data.append((K, S_TK, eig_K));
+ beta = min(eig_K, beta);
+ return data, beta;
+
+def dissipation_operator(
+ shape: list[int, int],
+ A: list[np.ndarray],
+ K: list[int],
+) -> tuple[np.ndarray, float]:
+ '''
+ Computes the dissipation operator `S_{T,K}` and its minimum spectral values.
+ '''
+ sh = tuple(shape)
+ m, n = sh;
+ assert m == n, 'must be square';
+ N = m;
+ O = generate_op_zero(shape=[N, N]);
+ I = generate_op_identity(shape=[N, N]);
+ S = O;
+ for C, C_ in partitions_two(K):
+ A_C = I;
+ A_C_ = I;
+ for i in C:
+ A_C = A_C @ -A[i];
+ for i in C_:
+ A_C_ = A_C_ @ -A[i];
+ S += A_C.conj().T @ A_C_;
+ S = (1/2**len(K)) * S;
+ eig = min(np.real(np.linalg.eigvals(S)));
+ return S, eig;
+
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+# MISCELLANEOUS METHODS - SETS
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+def power_set(K: list[int]) -> Generator[list[int], None, None]:
+ '''
+ Iterates through all subsets `C` of a set `K`.
+ '''
+ n = len(K);
+ for ch in itertools.product(*([[0,1]]*n)):
+ yield [ K[i] for i in range(n) if ch[i] == 1 ];
+
+def partitions_two(K: list[int]) -> Generator[tuple[list[int],list[int]], None, None]:
+ '''
+ Iterates through all partitions `[C1, C2]` of a set `K` of size 2.
+ '''
+ n = len(K);
+ for ch in itertools.product(*([[0,1]]*n)):
+ yield (
+ [ K[i] for i in range(n) if ch[i] == 1 ],
+ [ K[i] for i in range(n) if ch[i] == 0 ],
+ );
+
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+# MISCELLANEOUS METHODS - BASIC OPERATORS
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+def generate_op_zero(shape: list[int, int]) -> np.ndarray:
+ '''
+ Returns the zero operator with the right shape and type.
+ '''
+ sh = tuple(shape)
+ O = np.zeros(sh).astype(complex);
+ return O;
+
+def generate_op_identity(shape: list[int, int]) -> np.ndarray:
+ '''
+ Returns the identity operator with the right shape and type.
+ '''
+ sh = tuple(shape)
+ m, n = sh;
+ assert m == n, 'must be square';
+ I = np.eye(m).astype(complex);
+ return I;
+
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+# MISCELLANEOUS METHODS - GENERATE RANDOM MATRICES
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+def generate_op_selfadjoint(
+ shape: list[int, int],
+ bound: float = 1.,
+) -> np.ndarray:
+ '''
+ Generates a self-adjoint operator.
+
+ @inputs
+ - `shape` - <int, int> Desired shape of the output.
+ - `bound` - <float> The desired norm of output operator.
+
+ @returns
+ A randomly generated self-adjoint operator as a matrix.
+ '''
+ sh = tuple(shape)
+ m, n = sh;
+ assert m == n, 'must be square';
+ A = np.random.random(sh) + 1j * np.random.random(sh);
+ S = (A + A.conj().T)/2.;
+ M = np.linalg.norm(S);
+ if M > 0:
+ S = (bound/M) * S;
+ return S;
+
+def generate_onb(shape: list[int, int]) -> np.ndarray:
+ '''
+ Generates an ONB.
+
+ @inputs
+ - `shape` - <int, int> Desired shape of matrices.
+
+ @returns
+ An ONB for an n-dim complex Hilbert space as a matrix.
+ '''
+ sh = tuple(shape)
+ m, n = sh;
+ assert m == n, 'must be square';
+ N = m;
+ O = generate_op_zero(shape=[N, N]);
+ P = O;
+ for k in range(N):
+ while True:
+ v = np.random.random((N,)) + 1j * np.random.random((N,));
+ for i in range(k):
+ # NOTE: np.inner does not work.
+ coeff = np.sum([ v[j] * P[j, i].conj() for j in range(N)]);
+ v -= coeff * P[:, i];
+ norm_v = np.linalg.norm(v);
+ if norm_v > 0.:
+ P[:, k] = v / norm_v;
+ break;
+ return P;
+
+def generate_op_unitary(
+ shape: list[int, int],
+ rational: bool = False,
+ base: int = 2,
+) -> np.ndarray:
+ '''
+ Generates a unitary operator.
+
+ @inputs
+ - `shape` - <int, int> Desired shape of output.
+ - `rational` - <bool> Whether or not entries are to be rational.
+ - `base` - <int> If `rational = True`, fixes the denominator of the rational numbers.
+
+ @returns
+ A randomly generated unitary operator as a matrix.
+ '''
+ sh = tuple(shape)
+ m, n = sh;
+ assert m == n, 'must be square';
+ N = m;
+ P = generate_onb(shape=shape);
+ if rational:
+ theta = 2 * np.pi * np.random.randint(low=0, high=base, size=N) / base;
+ else:
+ theta = 2 * np.pi * np.random.random(size=N);
+ D = np.diag(np.exp(1j * theta));
+ U = P @ D @ P.conj().T;
+ return U;
+
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+# MISCELLANEOUS METHODS - Get user input
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+def ask_confirmation(message: str):
+ while True:
+ answer = input(f'{message} (y/n): ');
+ if re.match(pattern='^y(es)?|1$', string=answer, flags=re.IGNORECASE):
+ return True;
+ elif re.match(pattern='^n(o)?|0$', string=answer, flags=re.IGNORECASE):
+ return False;