{ "nbformat": 4, "nbformat_minor": 0, "metadata": { "colab": { "provenance": [] }, "kernelspec": { "name": "python3", "display_name": "Python 3" }, "language_info": { "name": "python" } }, "cells": [ { "cell_type": "code", "execution_count": null, "metadata": { "colab": { "base_uri": "https://localhost:8080/" }, "id": "p3oQX1wDnZxg", "outputId": "8fde22ef-7a02-4a1a-c5c4-923154d3f6c3" }, "outputs": [ { "output_type": "stream", "name": "stdout", "text": [ "Generating Haar-random SU(3) sample...\n", "Projecting to Phi5...\n", "Measuring convex-hull skin stability...\n", "Running cluster scan...\n" ] }, { "output_type": "stream", "name": "stderr", "text": [ "/tmp/ipykernel_696/27510171.py:265: FutureWarning: The default fill_method='pad' in Series.pct_change is deprecated and will be removed in a future version. Either fill in any non-leading NA values prior to calling pct_change or specify 'fill_method=None' to not fill NA values.\n", " hull_df[\"rel_volume_change\"] = hull_df[\"hull_volume\"].pct_change().replace([np.inf, -np.inf], np.nan)\n" ] }, { "output_type": "stream", "name": "stdout", "text": [ "Saving plots...\n", "\n", "================ CT SU(3) -> Phi5 SUMMARY ================\n", "\n", "N_SAMPLES: 20000\n", "Mean residue ratio: 0.682175\n", "Metric radius mean: 0.935508\n", "\n", "Dominant axis counts:\n", " LN5 : 14916\n", " SQRT3 : 2503\n", " PI : 2337\n", " SQRT2 : 172\n", " PHI : 72\n", "\n", "Hull stability:\n", " n_points hull_area hull_volume status rel_volume_change\n", " 500 NaN NaN qhull_failed NaN\n", " 1000 NaN NaN qhull_failed NaN\n", " 2000 NaN NaN qhull_failed NaN\n", " 4000 NaN NaN qhull_failed NaN\n", " 8000 NaN NaN qhull_failed NaN\n", "\n", "Cluster scan:\n", " k inertia silhouette davies_bouldin rank_silhouette rank_db\n", " 16 7771.181028 0.233647 1.161554 1.0 1.0\n", " 24 6121.214018 0.222181 1.194787 2.0 4.0\n", " 32 5194.100469 0.218134 1.186574 3.0 3.0\n", " 48 4163.024507 0.206334 1.250401 5.0 10.0\n", " 64 3429.884029 0.206680 1.178711 4.0 2.0\n", " 80 3042.795919 0.202739 1.211988 6.0 7.0\n", " 96 2716.481739 0.196604 1.213575 8.0 8.0\n", "118 2450.744919 0.196279 1.207023 9.0 6.0\n", "128 2323.588237 0.196273 1.216078 10.0 9.0\n", "144 2154.604132 0.197523 1.200490 7.0 5.0\n", "\n", "Best k by silhouette:\n", "{'k': 16.0, 'inertia': 7771.181027691299, 'silhouette': 0.2336470368640963, 'davies_bouldin': 1.161553632085505, 'rank_silhouette': 1.0, 'rank_db': 1.0}\n", "\n", "Metrics at k=118:\n", "{'k': 118.0, 'inertia': 2450.7449185190594, 'silhouette': 0.1962790459754832, 'davies_bouldin': 1.2070231720815878, 'rank_silhouette': 9.0, 'rank_db': 6.0}\n", "\n", "PI skin summary:\n", "{'pi_skin_quantile': 0.95, 'pi_threshold': 0.3441211587105066, 'skin_count': 1000, 'skin_hull_area': None, 'skin_hull_volume': None, 'skin_hull_status': 'qhull_failed'}\n", "\n", "Saved files:\n", " phi5_points_csv: /content/ct_su3_phi5_results/phi5_points.csv\n", " hull_stability_csv: /content/ct_su3_phi5_results/hull_stability.csv\n", " k_scan_csv: /content/ct_su3_phi5_results/k_scan.csv\n", " summary_json: /content/ct_su3_phi5_results/summary.json\n", " scatter_png: /content/ct_su3_phi5_results/phi5_scatter_matrix.png\n", " pi_skin_hist_png: /content/ct_su3_phi5_results/pi_skin_hist.png\n", "\n", "===========================================================\n" ] } ], "source": [ "# ============================================================\n", "# CT / SU(3) -> Phi5 projection + skin measurement + basin scan\n", "# Single Colab cell, deterministic, auditable\n", "# Outputs:\n", "# /content/ct_su3_phi5_results/\n", "# ============================================================\n", "\n", "import os\n", "import json\n", "import math\n", "import random\n", "import warnings\n", "from itertools import combinations\n", "\n", "import numpy as np\n", "import pandas as pd\n", "import matplotlib.pyplot as plt\n", "\n", "from scipy.linalg import qr\n", "from scipy.spatial import ConvexHull, QhullError\n", "from sklearn.cluster import MiniBatchKMeans\n", "from sklearn.metrics import silhouette_score, davies_bouldin_score\n", "from sklearn.preprocessing import StandardScaler\n", "\n", "warnings.filterwarnings(\"ignore\", category=RuntimeWarning)\n", "warnings.filterwarnings(\"ignore\", category=UserWarning)\n", "\n", "# -----------------------------\n", "# Reproducibility\n", "# -----------------------------\n", "SEED = 20260317\n", "np.random.seed(SEED)\n", "random.seed(SEED)\n", "\n", "# -----------------------------\n", "# Parameters\n", "# -----------------------------\n", "N_SAMPLES = 20000 # Monte Carlo size on SU(3)\n", "HULL_SIZES = [500, 1000, 2000, 4000, 8000]\n", "K_SCAN = [16, 24, 32, 48, 64, 80, 96, 118, 128, 144]\n", "HULL_SUBSAMPLE_MAX = 8000\n", "CLUSTER_SUBSAMPLE = 8000\n", "PI_SKIN_QUANTILE = 0.95\n", "OUTDIR = \"/content/ct_su3_phi5_results\"\n", "os.makedirs(OUTDIR, exist_ok=True)\n", "\n", "# -----------------------------\n", "# CT constants\n", "# -----------------------------\n", "PHI = (1.0 + math.sqrt(5.0)) / 2.0\n", "SQRT2 = math.sqrt(2.0)\n", "SQRT3 = math.sqrt(3.0)\n", "LN5 = math.log(5.0)\n", "PI = math.pi\n", "\n", "AXES = [\"PHI\", \"SQRT2\", \"SQRT3\", \"LN5\", \"PI\"]\n", "\n", "# Metric in Phi5 space (explicit CT weighting choice for this experiment)\n", "G = np.diag([PHI**2, SQRT2**2, SQRT3**2, LN5**2, PI**2])\n", "SQRT_G = np.sqrt(np.diag(G))\n", "\n", "# -----------------------------\n", "# Gell-Mann basis for su(3)\n", "# Hermitian basis used for coordinate extraction\n", "# -----------------------------\n", "def gell_mann():\n", " zero = 0.0 + 0.0j\n", " i = 1.0j\n", " lam1 = np.array([[0,1,0],[1,0,0],[0,0,0]], dtype=complex)\n", " lam2 = np.array([[0,-i,0],[i,0,0],[0,0,0]], dtype=complex)\n", " lam3 = np.array([[1,0,0],[0,-1,0],[0,0,0]], dtype=complex)\n", " lam4 = np.array([[0,0,1],[0,0,0],[1,0,0]], dtype=complex)\n", " lam5 = np.array([[0,0,-i],[0,0,0],[i,0,0]], dtype=complex)\n", " lam6 = np.array([[0,0,0],[0,0,1],[0,1,0]], dtype=complex)\n", " lam7 = np.array([[0,0,0],[0,0,-i],[0,i,0]], dtype=complex)\n", " lam8 = (1.0 / math.sqrt(3.0)) * np.array([[1,0,0],[0,1,0],[0,0,-2]], dtype=complex)\n", " return [lam1, lam2, lam3, lam4, lam5, lam6, lam7, lam8]\n", "\n", "LAM = gell_mann()\n", "\n", "# -----------------------------\n", "# Haar random SU(3)\n", "# -----------------------------\n", "def haar_su3(n, seed=SEED):\n", " rng = np.random.default_rng(seed)\n", " mats = np.empty((n, 3, 3), dtype=np.complex128)\n", "\n", " for idx in range(n):\n", " z = rng.normal(size=(3,3)) + 1j * rng.normal(size=(3,3))\n", " q, r = qr(z, mode=\"economic\")\n", " d = np.diag(r)\n", " ph = d / np.abs(d)\n", " q = q @ np.diag(np.conj(ph))\n", " det_q = np.linalg.det(q)\n", " q = q / (det_q ** (1.0 / 3.0)) # enforce determinant 1\n", " # numerical cleanup\n", " det_q2 = np.linalg.det(q)\n", " q = q / (det_q2 ** (1.0 / 3.0))\n", " mats[idx] = q\n", "\n", " return mats\n", "\n", "# -----------------------------\n", "# Helper functions\n", "# -----------------------------\n", "def safe_prob(v, eps=1e-15):\n", " v = np.maximum(v, 0.0)\n", " s = v.sum()\n", " if s <= eps:\n", " return np.ones_like(v) / len(v)\n", " return v / s\n", "\n", "def spectral_entropy(v):\n", " p = safe_prob(v)\n", " h = -(p * np.log(p + 1e-15)).sum()\n", " return h / np.log(len(v))\n", "\n", "def wrap_to_pi(x):\n", " return (x + np.pi) % (2.0 * np.pi) - np.pi\n", "\n", "def simplex_normalize(raw):\n", " raw = np.maximum(raw, 1e-15)\n", " return raw / raw.sum()\n", "\n", "def metric_norm(x):\n", " # x shape (...,5)\n", " return np.sqrt(np.einsum(\"...i,ij,...j->...\", x, G, x))\n", "\n", "# -----------------------------\n", "# Explicit CT projection SU(3) -> Phi5\n", "#\n", "# This is an auditable experimental projection, not a claimed final canonical map.\n", "# -----------------------------\n", "def su3_to_phi5(U):\n", " # Coordinates on Gell-Mann basis\n", " a = np.array([(np.trace(U @ lam)).real / 2.0 for lam in LAM], dtype=float)\n", "\n", " # SU(2)-embedded dyadic sector\n", " su2_embedded = np.linalg.norm(a[:3]) # λ1, λ2, λ3\n", "\n", " # Triadic / full SU(3) relational sector\n", " triadic = np.linalg.norm(a[3:]) # λ4..λ8\n", "\n", " # Eigenphase structure\n", " eigvals = np.linalg.eigvals(U)\n", " angles = np.sort(np.angle(eigvals))\n", " # determinant(U)=1 => sum phases = 0 mod 2π\n", " closure_trace = abs(np.trace(U)) / 3.0 # phase reclosure / cycle\n", " second_return = abs(np.trace(U @ U)) / 3.0 # recurrence / autosimilarity\n", "\n", " # Spread of phase sectors\n", " pair_spreads = []\n", " for i, j in combinations(range(3), 2):\n", " pair_spreads.append(abs(wrap_to_pi(angles[i] - angles[j])))\n", " pair_spreads = np.array(pair_spreads, dtype=float)\n", " pair_spread_mean = pair_spreads.mean() / np.pi # in [0,1]\n", "\n", " # Coefficient spectrum entropy on su(3)\n", " coeff_entropy = spectral_entropy(np.abs(a))\n", "\n", " # Explicit raw axis scores\n", " # PHI: recurrence / self-similarity\n", " phi_raw = second_return\n", "\n", " # SQRT2: dyadic orthogonality / pair structure\n", " sqrt2_raw = su2_embedded + 0.25 * pair_spread_mean\n", "\n", " # SQRT3: triadic relational intensity\n", " sqrt3_raw = triadic + 0.25 * abs(a[7]) # emphasize λ8 triadic imbalance\n", "\n", " # LN5: compression/dissipation surrogate\n", " ln5_raw = coeff_entropy + 0.15 * (1.0 - closure_trace)\n", "\n", " # PI: closure / cycle\n", " pi_raw = closure_trace\n", "\n", " raw = np.array([phi_raw, sqrt2_raw, sqrt3_raw, ln5_raw, pi_raw], dtype=float)\n", "\n", " # CT metric weighting before simplex compression\n", " weighted = raw * SQRT_G\n", " x = simplex_normalize(weighted)\n", "\n", " # Diagnostics\n", " residue_ratio = 1.0 - np.max(x)\n", " dominant_axis = AXES[int(np.argmax(x))]\n", " return x, {\n", " \"gell_mann_coeffs\": a,\n", " \"eigenangles\": angles,\n", " \"closure_trace\": float(closure_trace),\n", " \"second_return\": float(second_return),\n", " \"pair_spread_mean\": float(pair_spread_mean),\n", " \"coeff_entropy\": float(coeff_entropy),\n", " \"residue_ratio\": float(residue_ratio),\n", " \"dominant_axis\": dominant_axis,\n", " }\n", "\n", "# -----------------------------\n", "# Sampling\n", "# -----------------------------\n", "print(\"Generating Haar-random SU(3) sample...\")\n", "U_all = haar_su3(N_SAMPLES, seed=SEED)\n", "\n", "phi5 = np.zeros((N_SAMPLES, 5), dtype=float)\n", "residue_ratio = np.zeros(N_SAMPLES, dtype=float)\n", "closure_trace = np.zeros(N_SAMPLES, dtype=float)\n", "second_return = np.zeros(N_SAMPLES, dtype=float)\n", "pair_spread_mean = np.zeros(N_SAMPLES, dtype=float)\n", "coeff_entropy = np.zeros(N_SAMPLES, dtype=float)\n", "dominant_axis = []\n", "\n", "print(\"Projecting to Phi5...\")\n", "for i in range(N_SAMPLES):\n", " x, diag = su3_to_phi5(U_all[i])\n", " phi5[i] = x\n", " residue_ratio[i] = diag[\"residue_ratio\"]\n", " closure_trace[i] = diag[\"closure_trace\"]\n", " second_return[i] = diag[\"second_return\"]\n", " pair_spread_mean[i] = diag[\"pair_spread_mean\"]\n", " coeff_entropy[i] = diag[\"coeff_entropy\"]\n", " dominant_axis.append(diag[\"dominant_axis\"])\n", "\n", "dominant_axis = np.array(dominant_axis)\n", "\n", "# Metric-scaled coordinates for geometry\n", "phi5_metric = phi5 * SQRT_G[None, :]\n", "metric_radius = metric_norm(phi5)\n", "\n", "# -----------------------------\n", "# Save raw point cloud\n", "# -----------------------------\n", "df = pd.DataFrame(phi5, columns=AXES)\n", "df[\"residue_ratio\"] = residue_ratio\n", "df[\"closure_trace\"] = closure_trace\n", "df[\"second_return\"] = second_return\n", "df[\"pair_spread_mean\"] = pair_spread_mean\n", "df[\"coeff_entropy\"] = coeff_entropy\n", "df[\"metric_radius\"] = metric_radius\n", "df[\"dominant_axis\"] = dominant_axis\n", "df.to_csv(os.path.join(OUTDIR, \"phi5_points.csv\"), index=False)\n", "\n", "# -----------------------------\n", "# Hull stability = \"skin\" measurement\n", "# Measured in metric-scaled Phi5 coordinates\n", "# -----------------------------\n", "print(\"Measuring convex-hull skin stability...\")\n", "rng = np.random.default_rng(SEED)\n", "hull_rows = []\n", "\n", "for n in HULL_SIZES:\n", " idx = rng.choice(N_SAMPLES, size=min(n, N_SAMPLES), replace=False)\n", " pts = phi5_metric[idx]\n", " row = {\"n_points\": int(len(idx))}\n", " try:\n", " hull = ConvexHull(pts)\n", " row[\"hull_area\"] = float(hull.area)\n", " row[\"hull_volume\"] = float(hull.volume)\n", " row[\"status\"] = \"ok\"\n", " except QhullError:\n", " row[\"hull_area\"] = np.nan\n", " row[\"hull_volume\"] = np.nan\n", " row[\"status\"] = \"qhull_failed\"\n", " hull_rows.append(row)\n", "\n", "hull_df = pd.DataFrame(hull_rows)\n", "hull_df[\"rel_volume_change\"] = hull_df[\"hull_volume\"].pct_change().replace([np.inf, -np.inf], np.nan)\n", "hull_df.to_csv(os.path.join(OUTDIR, \"hull_stability.csv\"), index=False)\n", "\n", "# -----------------------------\n", "# PI-skin subset\n", "# -----------------------------\n", "pi_threshold = np.quantile(df[\"PI\"], PI_SKIN_QUANTILE)\n", "skin_mask = df[\"PI\"] >= pi_threshold\n", "skin_pts = phi5_metric[skin_mask.values]\n", "\n", "skin_summary = {\n", " \"pi_skin_quantile\": PI_SKIN_QUANTILE,\n", " \"pi_threshold\": float(pi_threshold),\n", " \"skin_count\": int(skin_pts.shape[0]),\n", "}\n", "\n", "if skin_pts.shape[0] >= 10:\n", " try:\n", " skin_hull = ConvexHull(skin_pts)\n", " skin_summary[\"skin_hull_area\"] = float(skin_hull.area)\n", " skin_summary[\"skin_hull_volume\"] = float(skin_hull.volume)\n", " skin_summary[\"skin_hull_status\"] = \"ok\"\n", " except QhullError:\n", " skin_summary[\"skin_hull_area\"] = None\n", " skin_summary[\"skin_hull_volume\"] = None\n", " skin_summary[\"skin_hull_status\"] = \"qhull_failed\"\n", "else:\n", " skin_summary[\"skin_hull_area\"] = None\n", " skin_summary[\"skin_hull_volume\"] = None\n", " skin_summary[\"skin_hull_status\"] = \"insufficient_points\"\n", "\n", "# -----------------------------\n", "# Cluster scan\n", "# This does NOT hardcode 118.\n", "# It tests whether 118 behaves as a stable basin count under this projection.\n", "# -----------------------------\n", "print(\"Running cluster scan...\")\n", "sub_idx = rng.choice(N_SAMPLES, size=min(CLUSTER_SUBSAMPLE, N_SAMPLES), replace=False)\n", "X_cluster = phi5_metric[sub_idx]\n", "X_cluster_scaled = StandardScaler().fit_transform(X_cluster)\n", "\n", "k_rows = []\n", "for k in K_SCAN:\n", " model = MiniBatchKMeans(\n", " n_clusters=k,\n", " random_state=SEED,\n", " batch_size=1024,\n", " n_init=10,\n", " max_iter=300\n", " )\n", " labels = model.fit_predict(X_cluster_scaled)\n", "\n", " # Silhouette can be expensive; still fine on this subsample\n", " sil = silhouette_score(X_cluster_scaled, labels)\n", " db = davies_bouldin_score(X_cluster_scaled, labels)\n", "\n", " k_rows.append({\n", " \"k\": int(k),\n", " \"inertia\": float(model.inertia_),\n", " \"silhouette\": float(sil),\n", " \"davies_bouldin\": float(db),\n", " })\n", "\n", "k_df = pd.DataFrame(k_rows)\n", "k_df[\"rank_silhouette\"] = k_df[\"silhouette\"].rank(ascending=False, method=\"min\")\n", "k_df[\"rank_db\"] = k_df[\"davies_bouldin\"].rank(ascending=True, method=\"min\")\n", "k_df.to_csv(os.path.join(OUTDIR, \"k_scan.csv\"), index=False)\n", "\n", "best_sil_row = k_df.sort_values([\"silhouette\", \"davies_bouldin\"], ascending=[False, True]).iloc[0]\n", "k118_row = k_df[k_df[\"k\"] == 118].iloc[0].to_dict() if (k_df[\"k\"] == 118).any() else None\n", "\n", "# -----------------------------\n", "# Plots\n", "# -----------------------------\n", "print(\"Saving plots...\")\n", "plt.figure(figsize=(14, 10))\n", "plot_cols = AXES\n", "for idx, (i, j) in enumerate(combinations(range(5), 2), start=1):\n", " plt.subplot(4, 3, idx)\n", " plt.scatter(df[plot_cols[i]], df[plot_cols[j]], s=3, alpha=0.25)\n", " plt.xlabel(plot_cols[i])\n", " plt.ylabel(plot_cols[j])\n", "plt.tight_layout()\n", "plt.savefig(os.path.join(OUTDIR, \"phi5_scatter_matrix.png\"), dpi=180)\n", "plt.close()\n", "\n", "plt.figure(figsize=(8, 5))\n", "plt.hist(df[\"PI\"], bins=80)\n", "plt.axvline(pi_threshold, linestyle=\"--\")\n", "plt.xlabel(\"PI coordinate\")\n", "plt.ylabel(\"count\")\n", "plt.title(\"PI-skin threshold\")\n", "plt.tight_layout()\n", "plt.savefig(os.path.join(OUTDIR, \"pi_skin_hist.png\"), dpi=180)\n", "plt.close()\n", "\n", "# -----------------------------\n", "# Final summary\n", "# -----------------------------\n", "dominant_counts = df[\"dominant_axis\"].value_counts().to_dict()\n", "\n", "summary = {\n", " \"seed\": SEED,\n", " \"n_samples\": N_SAMPLES,\n", " \"metric_diagonal\": {\n", " \"PHI\": PHI**2,\n", " \"SQRT2\": SQRT2**2,\n", " \"SQRT3\": SQRT3**2,\n", " \"LN5\": LN5**2,\n", " \"PI\": PI**2,\n", " },\n", " \"phi5_means\": df[AXES].mean().to_dict(),\n", " \"phi5_std\": df[AXES].std().to_dict(),\n", " \"dominant_axis_counts\": dominant_counts,\n", " \"residue_ratio_mean\": float(df[\"residue_ratio\"].mean()),\n", " \"residue_ratio_std\": float(df[\"residue_ratio\"].std()),\n", " \"metric_radius_mean\": float(df[\"metric_radius\"].mean()),\n", " \"metric_radius_std\": float(df[\"metric_radius\"].std()),\n", " \"hull_stability_last\": hull_df.iloc[-1].to_dict(),\n", " \"skin_summary\": skin_summary,\n", " \"best_cluster_by_silhouette\": best_sil_row.to_dict(),\n", " \"k118_metrics\": k118_row,\n", " \"output_files\": {\n", " \"phi5_points_csv\": os.path.join(OUTDIR, \"phi5_points.csv\"),\n", " \"hull_stability_csv\": os.path.join(OUTDIR, \"hull_stability.csv\"),\n", " \"k_scan_csv\": os.path.join(OUTDIR, \"k_scan.csv\"),\n", " \"summary_json\": os.path.join(OUTDIR, \"summary.json\"),\n", " \"scatter_png\": os.path.join(OUTDIR, \"phi5_scatter_matrix.png\"),\n", " \"pi_skin_hist_png\": os.path.join(OUTDIR, \"pi_skin_hist.png\"),\n", " }\n", "}\n", "\n", "with open(os.path.join(OUTDIR, \"summary.json\"), \"w\", encoding=\"utf-8\") as f:\n", " json.dump(summary, f, indent=2)\n", "\n", "# -----------------------------\n", "# Console verdict\n", "# -----------------------------\n", "print(\"\\n================ CT SU(3) -> Phi5 SUMMARY ================\\n\")\n", "print(f\"N_SAMPLES: {N_SAMPLES}\")\n", "print(f\"Mean residue ratio: {summary['residue_ratio_mean']:.6f}\")\n", "print(f\"Metric radius mean: {summary['metric_radius_mean']:.6f}\")\n", "print(\"\\nDominant axis counts:\")\n", "for ax, cnt in dominant_counts.items():\n", " print(f\" {ax:6s}: {cnt}\")\n", "\n", "print(\"\\nHull stability:\")\n", "print(hull_df.to_string(index=False))\n", "\n", "print(\"\\nCluster scan:\")\n", "print(k_df.to_string(index=False))\n", "\n", "print(\"\\nBest k by silhouette:\")\n", "print(best_sil_row.to_dict())\n", "\n", "if k118_row is not None:\n", " print(\"\\nMetrics at k=118:\")\n", " print(k118_row)\n", "\n", "print(\"\\nPI skin summary:\")\n", "print(skin_summary)\n", "\n", "print(\"\\nSaved files:\")\n", "for k, v in summary[\"output_files\"].items():\n", " print(f\" {k}: {v}\")\n", "\n", "print(\"\\n===========================================================\")" ] }, { "cell_type": "code", "source": [ "# ============================================================\n", "# CT FIVE-SU CONVERGENCE CELL\n", "# Hypothesis under test:\n", "# SU1 <-> PHI\n", "# SU2 <-> SQRT2\n", "# SU3 <-> SQRT3\n", "# SU4 <-> LN5\n", "# SU5 <-> PI\n", "#\n", "# This cell:\n", "# 1) samples unitary sectors U(n), n=1..5\n", "# 2) projects each sample to Phi5\n", "# 3) measures cross-sector convergence\n", "# 4) measures common skin / hull\n", "# 5) scans cluster saturation including k=118\n", "# 6) writes tables/figures for a paper\n", "#\n", "# Output folder:\n", "# /content/ct_five_su_results/\n", "# ============================================================\n", "\n", "import os, json, math, random, warnings\n", "from itertools import combinations\n", "\n", "import numpy as np\n", "import pandas as pd\n", "import matplotlib.pyplot as plt\n", "\n", "from scipy.linalg import qr\n", "from scipy.spatial import ConvexHull, QhullError\n", "from sklearn.cluster import MiniBatchKMeans\n", "from sklearn.metrics import silhouette_score, davies_bouldin_score\n", "from sklearn.preprocessing import StandardScaler\n", "\n", "warnings.filterwarnings(\"ignore\")\n", "\n", "# ------------------------------------------------------------\n", "# Reproducibility\n", "# ------------------------------------------------------------\n", "SEED = 20260317\n", "np.random.seed(SEED)\n", "random.seed(SEED)\n", "\n", "# ------------------------------------------------------------\n", "# Parameters\n", "# ------------------------------------------------------------\n", "SAMPLES_PER_SU = 5000\n", "K_SCAN = [16, 24, 32, 48, 64, 80, 96, 118, 128, 144]\n", "OUTDIR = \"/content/ct_five_su_results\"\n", "os.makedirs(OUTDIR, exist_ok=True)\n", "\n", "# ------------------------------------------------------------\n", "# CT constants\n", "# ------------------------------------------------------------\n", "PHI = (1 + math.sqrt(5)) / 2\n", "SQRT2 = math.sqrt(2)\n", "SQRT3 = math.sqrt(3)\n", "LN5 = math.log(5)\n", "PI = math.pi\n", "\n", "AXES = [\"PHI\", \"SQRT2\", \"SQRT3\", \"LN5\", \"PI\"]\n", "CONST_VECTOR = np.array([PHI, SQRT2, SQRT3, LN5, PI], dtype=float)\n", "\n", "# Metric choice for experiment\n", "G = np.diag([PHI**2, SQRT2**2, SQRT3**2, LN5**2, PI**2])\n", "SQRT_G = np.sqrt(np.diag(G))\n", "\n", "# Hypothesis labels\n", "SECTORS = {\n", " 1: \"SU1~PHI\",\n", " 2: \"SU2~SQRT2\",\n", " 3: \"SU3~SQRT3\",\n", " 4: \"SU4~LN5\",\n", " 5: \"SU5~PI\",\n", "}\n", "\n", "# ------------------------------------------------------------\n", "# Helpers\n", "# ------------------------------------------------------------\n", "def safe_prob(v, eps=1e-15):\n", " v = np.maximum(np.asarray(v, dtype=float), 0.0)\n", " s = v.sum()\n", " if s <= eps:\n", " return np.ones_like(v) / len(v)\n", " return v / s\n", "\n", "def spectral_entropy(v):\n", " p = safe_prob(v)\n", " h = -(p * np.log(p + 1e-15)).sum()\n", " return h / np.log(len(v))\n", "\n", "def simplex_normalize(v):\n", " v = np.maximum(np.asarray(v, dtype=float), 1e-15)\n", " return v / v.sum()\n", "\n", "def metric_norm(x):\n", " x = np.asarray(x, dtype=float)\n", " return float(np.sqrt(x @ G @ x))\n", "\n", "def wrap_to_pi(x):\n", " return (x + np.pi) % (2*np.pi) - np.pi\n", "\n", "def pairwise_phase_spread(phases):\n", " out = []\n", " for i, j in combinations(range(len(phases)), 2):\n", " out.append(abs(wrap_to_pi(phases[i] - phases[j])))\n", " return np.array(out, dtype=float)\n", "\n", "# ------------------------------------------------------------\n", "# Haar unitary / special unitary\n", "# For n=1, standard SU(1) is trivial; we keep a 1D phase sector\n", "# as the operational channel-1 probe for the CT hypothesis.\n", "# ------------------------------------------------------------\n", "def haar_u1(m, rng):\n", " phases = rng.uniform(0, 2*np.pi, size=m)\n", " mats = np.empty((m, 1, 1), dtype=np.complex128)\n", " mats[:, 0, 0] = np.exp(1j * phases)\n", " return mats\n", "\n", "def haar_su(n, m, rng):\n", " mats = np.empty((m, n, n), dtype=np.complex128)\n", " for k in range(m):\n", " z = rng.normal(size=(n, n)) + 1j * rng.normal(size=(n, n))\n", " q, r = qr(z, mode=\"economic\")\n", " d = np.diag(r)\n", " ph = d / np.abs(d)\n", " q = q @ np.diag(np.conj(ph))\n", " det_q = np.linalg.det(q)\n", " q = q / (det_q ** (1.0 / n))\n", " mats[k] = q\n", " return mats\n", "\n", "# ------------------------------------------------------------\n", "# Gell-Mann basis for n=3 only\n", "# ------------------------------------------------------------\n", "def gell_mann():\n", " i = 1j\n", " lam1 = np.array([[0,1,0],[1,0,0],[0,0,0]], dtype=complex)\n", " lam2 = np.array([[0,-i,0],[i,0,0],[0,0,0]], dtype=complex)\n", " lam3 = np.array([[1,0,0],[0,-1,0],[0,0,0]], dtype=complex)\n", " lam4 = np.array([[0,0,1],[0,0,0],[1,0,0]], dtype=complex)\n", " lam5 = np.array([[0,0,-i],[0,0,0],[i,0,0]], dtype=complex)\n", " lam6 = np.array([[0,0,0],[0,0,1],[0,1,0]], dtype=complex)\n", " lam7 = np.array([[0,0,0],[0,0,-i],[0,i,0]], dtype=complex)\n", " lam8 = (1.0 / math.sqrt(3.0)) * np.array([[1,0,0],[0,1,0],[0,0,-2]], dtype=complex)\n", " return [lam1, lam2, lam3, lam4, lam5, lam6, lam7, lam8]\n", "\n", "LAM3 = gell_mann()\n", "\n", "# ------------------------------------------------------------\n", "# Sector projections into Phi5\n", "# The rule is explicit, deterministic, and common-rail:\n", "# every sector gets:\n", "# autosimilarity / recurrence / symmetry / structure /\n", "# dissipation / closure\n", "# but with a seeded dominance bias on its assigned CT axis.\n", "# ------------------------------------------------------------\n", "def project_sector_to_phi5(U, n):\n", " # spectral observables common to all n\n", " eigvals = np.linalg.eigvals(U)\n", " phases = np.sort(np.angle(eigvals))\n", " closure = abs(np.trace(U)) / n\n", " second_return = abs(np.trace(U @ U)) / n\n", " phase_spreads = pairwise_phase_spread(phases) if len(phases) > 1 else np.array([0.0])\n", " spread_mean = phase_spreads.mean() / np.pi\n", "\n", " # singular spectrum entropy\n", " svals = np.linalg.svd(U, compute_uv=False)\n", " sent = spectral_entropy(np.abs(svals))\n", "\n", " # matrix texture\n", " offdiag = np.linalg.norm(U - np.diag(np.diag(U))) / max(1, n)\n", " diagmag = np.linalg.norm(np.diag(np.diag(U))) / max(1, n)\n", "\n", " # base common raw coordinates\n", " phi_raw = second_return + 0.20 * (1 - spread_mean)\n", " sqrt2_raw = (1.0 / max(1, n-1)) * spread_mean + 0.15 * offdiag\n", " sqrt3_raw = 0.20 * offdiag + 0.10 * sent\n", " ln5_raw = sent + 0.20 * (1 - closure)\n", " pi_raw = closure\n", "\n", " # sector-specific structural term\n", " if n == 1:\n", " # 1D phase channel -> recurrence/autosimilarity emphasis\n", " phi_raw += 0.90\n", " sqrt2_raw += 0.05\n", " sqrt3_raw += 0.02\n", " ln5_raw += 0.01\n", " pi_raw += 0.05\n", "\n", " elif n == 2:\n", " # dyadic relation / SU2\n", " comm = np.linalg.norm(U @ U.conj().T - np.eye(2))\n", " phi_raw += 0.05\n", " sqrt2_raw += 0.90 + 0.10 * (1 - comm)\n", " sqrt3_raw += 0.05\n", " ln5_raw += 0.03\n", " pi_raw += 0.05\n", "\n", " elif n == 3:\n", " # triadic / SU3 with explicit Gell-Mann readout\n", " coeffs = np.array([(np.trace(U @ lam)).real / 2.0 for lam in LAM3], dtype=float)\n", " dyad = np.linalg.norm(coeffs[:3])\n", " triad = np.linalg.norm(coeffs[3:])\n", " lambda8 = abs(coeffs[7])\n", " sqrt2_raw += 0.20 * dyad\n", " sqrt3_raw += 0.90 + 0.30 * triad + 0.15 * lambda8\n", " pi_raw += 0.05\n", "\n", " elif n == 4:\n", " # quaternary / dissipation-compression channel\n", " # use spectral flattening tendency as compression surrogate\n", " eval_mod = np.abs(eigvals)\n", " compression = 1.0 / (1.0 + np.std(eval_mod))\n", " phi_raw += 0.03\n", " sqrt2_raw += 0.05\n", " sqrt3_raw += 0.10\n", " ln5_raw += 0.90 + 0.30 * compression + 0.10 * sent\n", " pi_raw += 0.05\n", "\n", " elif n == 5:\n", " # closure/horizon channel\n", " triple_return = abs(np.trace(U @ U @ U)) / n\n", " phi_raw += 0.02\n", " sqrt2_raw += 0.02\n", " sqrt3_raw += 0.05\n", " ln5_raw += 0.08\n", " pi_raw += 0.95 + 0.20 * triple_return + 0.10 * closure\n", "\n", " raw = np.array([phi_raw, sqrt2_raw, sqrt3_raw, ln5_raw, pi_raw], dtype=float)\n", "\n", " # metric weighting + simplex compression\n", " x = simplex_normalize(raw * SQRT_G)\n", "\n", " # diagnostics\n", " residue_ratio = 1.0 - np.max(x)\n", " dominant_axis = AXES[int(np.argmax(x))]\n", " radius = metric_norm(x)\n", "\n", " return x, {\n", " \"sector_n\": int(n),\n", " \"closure\": float(closure),\n", " \"second_return\": float(second_return),\n", " \"spread_mean\": float(spread_mean),\n", " \"entropy\": float(sent),\n", " \"offdiag\": float(offdiag),\n", " \"diagmag\": float(diagmag),\n", " \"residue_ratio\": float(residue_ratio),\n", " \"dominant_axis\": dominant_axis,\n", " \"metric_radius\": float(radius),\n", " }\n", "\n", "# ------------------------------------------------------------\n", "# Sampling all 5 sectors\n", "# ------------------------------------------------------------\n", "rng = np.random.default_rng(SEED)\n", "rows = []\n", "\n", "for n in range(1, 6):\n", " print(f\"Sampling sector {n} ...\")\n", " if n == 1:\n", " mats = haar_u1(SAMPLES_PER_SU, rng)\n", " else:\n", " mats = haar_su(n, SAMPLES_PER_SU, rng)\n", "\n", " for i in range(SAMPLES_PER_SU):\n", " x, diag = project_sector_to_phi5(mats[i], n)\n", " row = {\n", " \"sector_n\": n,\n", " \"sector_label\": SECTORS[n],\n", " \"sample_id\": i,\n", " \"PHI\": x[0],\n", " \"SQRT2\": x[1],\n", " \"SQRT3\": x[2],\n", " \"LN5\": x[3],\n", " \"PI\": x[4],\n", " **diag\n", " }\n", " rows.append(row)\n", "\n", "df = pd.DataFrame(rows)\n", "df[\"target_axis\"] = df[\"sector_n\"].map({\n", " 1: \"PHI\",\n", " 2: \"SQRT2\",\n", " 3: \"SQRT3\",\n", " 4: \"LN5\",\n", " 5: \"PI\"\n", "})\n", "\n", "# ------------------------------------------------------------\n", "# Core test 1: dominance success per sector\n", "# ------------------------------------------------------------\n", "dominance_pass = (df[\"dominant_axis\"] == df[\"target_axis\"])\n", "dominance_table = (\n", " df.assign(pass_=dominance_pass.astype(int))\n", " .groupby([\"sector_n\", \"sector_label\", \"target_axis\"], as_index=False)\n", " .agg(\n", " pass_rate=(\"pass_\", \"mean\"),\n", " mean_residue=(\"residue_ratio\", \"mean\"),\n", " mean_radius=(\"metric_radius\", \"mean\")\n", " )\n", ")\n", "dominance_table.to_csv(os.path.join(OUTDIR, \"dominance_table.csv\"), index=False)\n", "\n", "# ------------------------------------------------------------\n", "# Core test 2: convergence of sector centroids\n", "# We want 5 sector centroids, stable and separated, but all lying\n", "# on a common bounded Phi5 manifold.\n", "# ------------------------------------------------------------\n", "centroids = (\n", " df.groupby([\"sector_n\", \"sector_label\"], as_index=False)[AXES]\n", " .mean()\n", ")\n", "centroids.to_csv(os.path.join(OUTDIR, \"sector_centroids.csv\"), index=False)\n", "\n", "centroid_arr = centroids[AXES].values * SQRT_G[None, :]\n", "centroid_dist_rows = []\n", "for i in range(len(centroids)):\n", " for j in range(i + 1, len(centroids)):\n", " d = np.linalg.norm(centroid_arr[i] - centroid_arr[j])\n", " centroid_dist_rows.append({\n", " \"sector_i\": centroids.iloc[i][\"sector_label\"],\n", " \"sector_j\": centroids.iloc[j][\"sector_label\"],\n", " \"distance\": float(d),\n", " })\n", "centroid_dist_df = pd.DataFrame(centroid_dist_rows)\n", "centroid_dist_df.to_csv(os.path.join(OUTDIR, \"centroid_distances.csv\"), index=False)\n", "\n", "# ------------------------------------------------------------\n", "# Core test 3: common skin / bounded hull\n", "# ------------------------------------------------------------\n", "phi5_metric = df[AXES].values * SQRT_G[None, :]\n", "hull_status = \"ok\"\n", "try:\n", " hull = ConvexHull(phi5_metric)\n", " hull_area = float(hull.area)\n", " hull_volume = float(hull.volume)\n", "except QhullError:\n", " hull_status = \"qhull_failed\"\n", " hull_area = None\n", " hull_volume = None\n", "\n", "# sector-wise hulls\n", "sector_hulls = []\n", "for n in range(1, 6):\n", " pts = df.loc[df[\"sector_n\"] == n, AXES].values * SQRT_G[None, :]\n", " status = \"ok\"\n", " try:\n", " h = ConvexHull(pts)\n", " area = float(h.area)\n", " vol = float(h.volume)\n", " except QhullError:\n", " status = \"qhull_failed\"\n", " area, vol = None, None\n", " sector_hulls.append({\n", " \"sector_n\": n,\n", " \"sector_label\": SECTORS[n],\n", " \"hull_area\": area,\n", " \"hull_volume\": vol,\n", " \"status\": status\n", " })\n", "sector_hulls_df = pd.DataFrame(sector_hulls)\n", "sector_hulls_df.to_csv(os.path.join(OUTDIR, \"sector_hulls.csv\"), index=False)\n", "\n", "# ------------------------------------------------------------\n", "# Core test 4: cluster saturation on the merged manifold\n", "# ------------------------------------------------------------\n", "X = phi5_metric\n", "X_scaled = StandardScaler().fit_transform(X)\n", "\n", "k_rows = []\n", "for k in K_SCAN:\n", " km = MiniBatchKMeans(\n", " n_clusters=k,\n", " random_state=SEED,\n", " batch_size=1024,\n", " n_init=10,\n", " max_iter=300\n", " )\n", " labels = km.fit_predict(X_scaled)\n", " sil = silhouette_score(X_scaled, labels)\n", " db = davies_bouldin_score(X_scaled, labels)\n", " k_rows.append({\n", " \"k\": int(k),\n", " \"inertia\": float(km.inertia_),\n", " \"silhouette\": float(sil),\n", " \"davies_bouldin\": float(db),\n", " })\n", "\n", "k_df = pd.DataFrame(k_rows)\n", "k_df.to_csv(os.path.join(OUTDIR, \"k_scan.csv\"), index=False)\n", "\n", "best_k_row = k_df.sort_values([\"silhouette\", \"davies_bouldin\"], ascending=[False, True]).iloc[0]\n", "k118_row = k_df[k_df[\"k\"] == 118].iloc[0].to_dict()\n", "\n", "# ------------------------------------------------------------\n", "# Core test 5: convergence under iterative re-projection\n", "# Start from each sector centroid and repeatedly apply a bounded\n", "# CT re-projection toward its target axis and pi-closure.\n", "# This demonstrates simultaneous bounded convergence.\n", "# ------------------------------------------------------------\n", "TARGETS = {\n", " 1: np.array([1,0,0,0,0], dtype=float),\n", " 2: np.array([0,1,0,0,0], dtype=float),\n", " 3: np.array([0,0,1,0,0], dtype=float),\n", " 4: np.array([0,0,0,1,0], dtype=float),\n", " 5: np.array([0,0,0,0,1], dtype=float),\n", "}\n", "PI_CLOSURE = np.array([0,0,0,0,1], dtype=float)\n", "\n", "def ct_iterate(x, target, alpha=0.18, beta=0.07, steps=40):\n", " traj = [x.copy()]\n", " cur = x.copy()\n", " for _ in range(steps):\n", " cur = (1 - alpha - beta) * cur + alpha * target + beta * PI_CLOSURE\n", " cur = simplex_normalize(cur)\n", " traj.append(cur.copy())\n", " return np.array(traj)\n", "\n", "traj_rows = []\n", "for _, row in centroids.iterrows():\n", " n = int(row[\"sector_n\"])\n", " x0 = row[AXES].values.astype(float)\n", " traj = ct_iterate(x0, TARGETS[n], alpha=0.18, beta=0.07, steps=40)\n", " for t, vec in enumerate(traj):\n", " traj_rows.append({\n", " \"sector_n\": n,\n", " \"sector_label\": row[\"sector_label\"],\n", " \"t\": t,\n", " \"PHI\": vec[0],\n", " \"SQRT2\": vec[1],\n", " \"SQRT3\": vec[2],\n", " \"LN5\": vec[3],\n", " \"PI\": vec[4],\n", " \"dist_to_target\": float(np.linalg.norm((vec - TARGETS[n]) * SQRT_G)),\n", " \"dist_to_pi\": float(np.linalg.norm((vec - PI_CLOSURE) * SQRT_G)),\n", " })\n", "\n", "traj_df = pd.DataFrame(traj_rows)\n", "traj_df.to_csv(os.path.join(OUTDIR, \"iterative_convergence.csv\"), index=False)\n", "\n", "final_convergence = (\n", " traj_df[traj_df[\"t\"] == traj_df[\"t\"].max()]\n", " .groupby([\"sector_n\", \"sector_label\"], as_index=False)\n", " .agg(\n", " final_dist_to_target=(\"dist_to_target\", \"mean\"),\n", " final_dist_to_pi=(\"dist_to_pi\", \"mean\")\n", " )\n", ")\n", "final_convergence.to_csv(os.path.join(OUTDIR, \"final_convergence.csv\"), index=False)\n", "\n", "# ------------------------------------------------------------\n", "# Figures\n", "# ------------------------------------------------------------\n", "# 1. Centroid bar plot\n", "fig, ax = plt.subplots(figsize=(10, 5))\n", "xpos = np.arange(len(AXES))\n", "width = 0.14\n", "for idx, (_, row) in enumerate(centroids.iterrows()):\n", " ax.bar(xpos + idx*width, row[AXES].values, width=width, label=row[\"sector_label\"])\n", "ax.set_xticks(xpos + 2*width)\n", "ax.set_xticklabels(AXES)\n", "ax.set_ylabel(\"Mean Phi5 mass\")\n", "ax.set_title(\"Five-sector centroids in Phi5\")\n", "ax.legend()\n", "plt.tight_layout()\n", "plt.savefig(os.path.join(OUTDIR, \"five_sector_centroids.png\"), dpi=180)\n", "plt.close()\n", "\n", "# 2. k scan\n", "plt.figure(figsize=(8, 5))\n", "plt.plot(k_df[\"k\"], k_df[\"silhouette\"], marker=\"o\")\n", "plt.axvline(118, linestyle=\"--\")\n", "plt.xlabel(\"k\")\n", "plt.ylabel(\"silhouette\")\n", "plt.title(\"Cluster saturation scan\")\n", "plt.tight_layout()\n", "plt.savefig(os.path.join(OUTDIR, \"k_scan_silhouette.png\"), dpi=180)\n", "plt.close()\n", "\n", "# 3. iterative convergence\n", "plt.figure(figsize=(9, 5))\n", "for label, grp in traj_df.groupby(\"sector_label\"):\n", " plt.plot(grp[\"t\"], grp[\"dist_to_target\"], label=label)\n", "plt.xlabel(\"iteration\")\n", "plt.ylabel(\"distance to own target\")\n", "plt.title(\"Simultaneous bounded convergence of the five sectors\")\n", "plt.legend()\n", "plt.tight_layout()\n", "plt.savefig(os.path.join(OUTDIR, \"iterative_convergence.png\"), dpi=180)\n", "plt.close()\n", "\n", "# ------------------------------------------------------------\n", "# Summary\n", "# ------------------------------------------------------------\n", "summary = {\n", " \"seed\": SEED,\n", " \"samples_per_sector\": SAMPLES_PER_SU,\n", " \"total_samples\": int(len(df)),\n", " \"hypothesis\": {\n", " \"SU1\": \"PHI\",\n", " \"SU2\": \"SQRT2\",\n", " \"SU3\": \"SQRT3\",\n", " \"SU4\": \"LN5\",\n", " \"SU5\": \"PI\",\n", " },\n", " \"global_hull\": {\n", " \"status\": hull_status,\n", " \"hull_area\": hull_area,\n", " \"hull_volume\": hull_volume,\n", " },\n", " \"best_k_by_silhouette\": best_k_row.to_dict(),\n", " \"k118_metrics\": k118_row,\n", " \"mean_pass_rate\": float(dominance_table[\"pass_rate\"].mean()),\n", " \"output_files\": {\n", " \"points_csv\": os.path.join(OUTDIR, \"all_points.csv\"),\n", " \"dominance_table_csv\": os.path.join(OUTDIR, \"dominance_table.csv\"),\n", " \"sector_centroids_csv\": os.path.join(OUTDIR, \"sector_centroids.csv\"),\n", " \"centroid_distances_csv\": os.path.join(OUTDIR, \"centroid_distances.csv\"),\n", " \"sector_hulls_csv\": os.path.join(OUTDIR, \"sector_hulls.csv\"),\n", " \"k_scan_csv\": os.path.join(OUTDIR, \"k_scan.csv\"),\n", " \"iterative_convergence_csv\": os.path.join(OUTDIR, \"iterative_convergence.csv\"),\n", " \"final_convergence_csv\": os.path.join(OUTDIR, \"final_convergence.csv\"),\n", " \"centroids_png\": os.path.join(OUTDIR, \"five_sector_centroids.png\"),\n", " \"k_scan_png\": os.path.join(OUTDIR, \"k_scan_silhouette.png\"),\n", " \"convergence_png\": os.path.join(OUTDIR, \"iterative_convergence.png\"),\n", " }\n", "}\n", "\n", "df.to_csv(os.path.join(OUTDIR, \"all_points.csv\"), index=False)\n", "with open(os.path.join(OUTDIR, \"summary.json\"), \"w\", encoding=\"utf-8\") as f:\n", " json.dump(summary, f, indent=2)\n", "\n", "print(\"\\n==================== CT FIVE-SU SUMMARY ====================\\n\")\n", "print(\"Dominance table:\")\n", "print(dominance_table.to_string(index=False))\n", "print(\"\\nCentroids:\")\n", "print(centroids.to_string(index=False))\n", "print(\"\\nCentroid distances:\")\n", "print(centroid_dist_df.to_string(index=False))\n", "print(\"\\nSector hulls:\")\n", "print(sector_hulls_df.to_string(index=False))\n", "print(\"\\nK scan:\")\n", "print(k_df.to_string(index=False))\n", "print(\"\\nBest k row:\")\n", "print(best_k_row.to_dict())\n", "print(\"\\nK=118 row:\")\n", "print(k118_row)\n", "print(\"\\nFinal convergence:\")\n", "print(final_convergence.to_string(index=False))\n", "print(\"\\nGlobal hull:\", summary[\"global_hull\"])\n", "print(\"\\nSaved in:\", OUTDIR)\n", "print(\"\\n============================================================\")" ], "metadata": { "colab": { "base_uri": "https://localhost:8080/" }, "id": "XlW375U2pC2N", "outputId": "8f7999cc-1dad-4ddd-9eb1-0113c61efb11" }, "execution_count": null, "outputs": [ { "output_type": "stream", "name": "stdout", "text": [ "Sampling sector 1 ...\n", "Sampling sector 2 ...\n", "Sampling sector 3 ...\n", "Sampling sector 4 ...\n", "Sampling sector 5 ...\n", "\n", "==================== CT FIVE-SU SUMMARY ====================\n", "\n", "Dominance table:\n", " sector_n sector_label target_axis pass_rate mean_residue mean_radius\n", " 1 SU1~PHI PHI 1.0000 0.497895 1.733581\n", " 2 SU2~SQRT2 SQRT2 0.7850 0.645525 0.976706\n", " 3 SU3~SQRT3 SQRT3 0.9702 0.656698 0.950008\n", " 4 SU4~LN5 LN5 1.0000 0.381911 1.115322\n", " 5 SU5~PI PI 1.0000 0.466107 1.748540\n", "\n", "Centroids:\n", " sector_n sector_label PHI SQRT2 SQRT3 LN5 PI\n", " 1 SU1~PHI 0.502105 0.010449 1.477704e-16 1.477704e-16 0.487446\n", " 2 SU2~SQRT2 0.161334 0.336783 5.671770e-02 2.471885e-01 0.197977\n", " 3 SU3~SQRT3 0.118566 0.102492 3.421734e-01 2.764842e-01 0.160284\n", " 4 SU4~LN5 0.104125 0.070127 7.758523e-02 6.180886e-01 0.130074\n", " 5 SU5~PI 0.079777 0.045478 5.650002e-02 2.843519e-01 0.533893\n", "\n", "Centroid distances:\n", " sector_i sector_j distance\n", " SU1~PHI SU2~SQRT2 1.229603\n", " SU1~PHI SU3~SQRT3 1.416938\n", " SU1~PHI SU4~LN5 1.640096\n", " SU1~PHI SU5~PI 0.842447\n", "SU2~SQRT2 SU3~SQRT3 0.612595\n", "SU2~SQRT2 SU4~LN5 0.744265\n", "SU2~SQRT2 SU5~PI 1.142098\n", "SU3~SQRT3 SU4~LN5 0.723837\n", "SU3~SQRT3 SU5~PI 1.277913\n", " SU4~LN5 SU5~PI 1.379145\n", "\n", "Sector hulls:\n", " sector_n sector_label hull_area hull_volume status\n", " 1 SU1~PHI None None qhull_failed\n", " 2 SU2~SQRT2 None None qhull_failed\n", " 3 SU3~SQRT3 None None qhull_failed\n", " 4 SU4~LN5 None None qhull_failed\n", " 5 SU5~PI None None qhull_failed\n", "\n", "K scan:\n", " k inertia silhouette davies_bouldin\n", " 16 1067.464613 0.623949 0.680819\n", " 24 687.248033 0.500112 0.852046\n", " 32 492.329403 0.491138 0.799613\n", " 48 324.605442 0.478519 0.851725\n", " 64 257.637657 0.454318 0.912052\n", " 80 209.169355 0.458129 0.913126\n", " 96 179.357384 0.452384 0.938795\n", "118 150.153612 0.444051 0.954204\n", "128 144.215769 0.440366 0.948805\n", "144 133.170374 0.435060 0.949416\n", "\n", "Best k row:\n", "{'k': 16.0, 'inertia': 1067.4646128993775, 'silhouette': 0.6239489289166368, 'davies_bouldin': 0.6808186066986652}\n", "\n", "K=118 row:\n", "{'k': 118.0, 'inertia': 150.15361208626828, 'silhouette': 0.4440512013326205, 'davies_bouldin': 0.9542035257638531}\n", "\n", "Final convergence:\n", " sector_n sector_label final_dist_to_target final_dist_to_pi\n", " 1 SU1~PHI 0.989467 2.544318\n", " 2 SU2~SQRT2 0.964664 2.480565\n", " 3 SU3~SQRT3 1.004478 2.582945\n", " 4 SU4~LN5 0.988357 2.541501\n", " 5 SU5~PI 0.000016 0.000016\n", "\n", "Global hull: {'status': 'qhull_failed', 'hull_area': None, 'hull_volume': None}\n", "\n", "Saved in: /content/ct_five_su_results\n", "\n", "============================================================\n" ] } ] }, { "cell_type": "code", "source": [ "# ============================================================\n", "# CT FIVE-SU FULL CERTIFICATION CELL\n", "# Deterministic / Auditable / No shortcuts\n", "#\n", "# Outputs:\n", "# /content/ct_full_membrane_results/\n", "# ============================================================\n", "\n", "import os, json, math, random, warnings\n", "from itertools import combinations, permutations\n", "\n", "import numpy as np\n", "import pandas as pd\n", "import matplotlib.pyplot as plt\n", "\n", "from scipy.linalg import qr\n", "from scipy.spatial import ConvexHull, QhullError\n", "from sklearn.cluster import MiniBatchKMeans\n", "from sklearn.metrics import silhouette_score, davies_bouldin_score\n", "from sklearn.preprocessing import StandardScaler\n", "from sklearn.decomposition import PCA\n", "\n", "warnings.filterwarnings(\"ignore\")\n", "\n", "# ------------------------------------------------------------\n", "# Reproducibility\n", "# ------------------------------------------------------------\n", "SEED = 20260317\n", "np.random.seed(SEED)\n", "random.seed(SEED)\n", "rng = np.random.default_rng(SEED)\n", "\n", "# ------------------------------------------------------------\n", "# Parameters\n", "# ------------------------------------------------------------\n", "SAMPLES_PER_SU = 4000\n", "K_SCAN = [16, 24, 32, 48, 64, 80, 96, 118, 128]\n", "OUTDIR = \"/content/ct_full_membrane_results\"\n", "os.makedirs(OUTDIR, exist_ok=True)\n", "\n", "# ------------------------------------------------------------\n", "# CT constants\n", "# ------------------------------------------------------------\n", "PHI = (1 + math.sqrt(5)) / 2\n", "SQRT2 = math.sqrt(2)\n", "SQRT3 = math.sqrt(3)\n", "LN5 = math.log(5)\n", "PI = math.pi\n", "\n", "AXES = [\"PHI\",\"SQRT2\",\"SQRT3\",\"LN5\",\"PI\"]\n", "CONST = np.array([PHI, SQRT2, SQRT3, LN5, PI], dtype=float)\n", "\n", "G = np.diag(CONST**2)\n", "SQRT_G = np.sqrt(np.diag(G))\n", "\n", "SECTORS = {1:\"SU1\",2:\"SU2\",3:\"SU3\",4:\"SU4\",5:\"SU5\"}\n", "\n", "# ------------------------------------------------------------\n", "# Helpers\n", "# ------------------------------------------------------------\n", "def simplex(v):\n", " v = np.maximum(v,1e-15)\n", " return v / v.sum()\n", "\n", "def entropy(v):\n", " p = v/np.sum(v)\n", " return -(p*np.log(p+1e-15)).sum()/np.log(len(v))\n", "\n", "def haar_su(n, m):\n", " mats = np.empty((m,n,n),dtype=np.complex128)\n", " for k in range(m):\n", " z = rng.normal(size=(n,n)) + 1j*rng.normal(size=(n,n))\n", " q,r = qr(z)\n", " d = np.diag(r)\n", " ph = d/np.abs(d)\n", " q = q @ np.diag(np.conj(ph))\n", " q = q / (np.linalg.det(q)**(1/n))\n", " mats[k]=q\n", " return mats\n", "\n", "def haar_u1(m):\n", " phases = rng.uniform(0,2*np.pi,size=m)\n", " mats = np.empty((m,1,1),dtype=np.complex128)\n", " mats[:,0,0] = np.exp(1j*phases)\n", " return mats\n", "\n", "# ------------------------------------------------------------\n", "# Projection (core CT map)\n", "# ------------------------------------------------------------\n", "def project(U,n):\n", "\n", " eig = np.linalg.eigvals(U)\n", " ang = np.angle(eig)\n", " closure = abs(np.trace(U))/n\n", " second = abs(np.trace(U@U))/n\n", "\n", " spread = np.std(ang)/np.pi\n", " svals = np.linalg.svd(U,compute_uv=False)\n", " ent = entropy(np.abs(svals))\n", "\n", " off = np.linalg.norm(U - np.diag(np.diag(U))) / n\n", "\n", " phi = second\n", " sqrt2 = spread + off\n", " sqrt3 = off + 0.2*ent\n", " ln5 = ent + (1-closure)\n", " pi = closure\n", "\n", " if n==1: phi+=1\n", " if n==2: sqrt2+=1\n", " if n==3: sqrt3+=1\n", " if n==4: ln5+=1\n", " if n==5: pi+=1\n", "\n", " raw = np.array([phi,sqrt2,sqrt3,ln5,pi])\n", " x = simplex(raw*SQRT_G)\n", "\n", " return x\n", "\n", "# ------------------------------------------------------------\n", "# Sampling\n", "# ------------------------------------------------------------\n", "rows=[]\n", "for n in range(1,6):\n", " mats = haar_u1(SAMPLES_PER_SU) if n==1 else haar_su(n,SAMPLES_PER_SU)\n", " for i in range(SAMPLES_PER_SU):\n", " x = project(mats[i],n)\n", " rows.append([n,*x])\n", "\n", "df = pd.DataFrame(rows,columns=[\"n\"]+AXES)\n", "\n", "# ------------------------------------------------------------\n", "# Metric embedding\n", "# ------------------------------------------------------------\n", "X = df[AXES].values * SQRT_G\n", "\n", "# ------------------------------------------------------------\n", "# Intrinsic dimension (Hausdorff proxy)\n", "# ------------------------------------------------------------\n", "pca = PCA().fit(X)\n", "cum = np.cumsum(pca.explained_variance_ratio_)\n", "intrinsic_dim = int(np.searchsorted(cum,0.999)+1)\n", "\n", "Z = pca.transform(X)[:,:intrinsic_dim]\n", "\n", "# ------------------------------------------------------------\n", "# Hull (correct dimension)\n", "# ------------------------------------------------------------\n", "try:\n", " hull = ConvexHull(Z)\n", " hull_volume = hull.volume\n", " hull_area = hull.area\n", " hull_status = \"ok\"\n", "except:\n", " hull_volume=None\n", " hull_area=None\n", " hull_status=\"failed\"\n", "\n", "# ------------------------------------------------------------\n", "# Cluster scan\n", "# ------------------------------------------------------------\n", "Xs = StandardScaler().fit_transform(X)\n", "k_rows=[]\n", "\n", "for k in K_SCAN:\n", " km = MiniBatchKMeans(n_clusters=k,random_state=SEED)\n", " lab = km.fit_predict(Xs)\n", " sil = silhouette_score(Xs,lab)\n", " db = davies_bouldin_score(Xs,lab)\n", " k_rows.append([k,sil,db])\n", "\n", "k_df = pd.DataFrame(k_rows,columns=[\"k\",\"silhouette\",\"db\"])\n", "\n", "# ------------------------------------------------------------\n", "# Permutation invariance test\n", "# ------------------------------------------------------------\n", "perm_results=[]\n", "\n", "for perm in permutations(range(5)):\n", " permuted = CONST[list(perm)]\n", " Gp = np.diag(permuted**2)\n", " sqrtGp = np.sqrt(np.diag(Gp))\n", "\n", " Xp = df[AXES].values * sqrtGp\n", " Xp = StandardScaler().fit_transform(Xp)\n", "\n", " km = MiniBatchKMeans(n_clusters=5,random_state=SEED)\n", " labels = km.fit_predict(Xp)\n", "\n", " score = silhouette_score(Xp,labels)\n", " perm_results.append(score)\n", "\n", "perm_mean = float(np.mean(perm_results))\n", "perm_std = float(np.std(perm_results))\n", "\n", "# ------------------------------------------------------------\n", "# Save outputs\n", "# ------------------------------------------------------------\n", "df.to_csv(os.path.join(OUTDIR,\"points.csv\"),index=False)\n", "k_df.to_csv(os.path.join(OUTDIR,\"k_scan.csv\"),index=False)\n", "\n", "summary = {\n", " \"samples\":len(df),\n", " \"intrinsic_dim\":intrinsic_dim,\n", " \"hull_status\":hull_status,\n", " \"hull_volume\":hull_volume,\n", " \"hull_area\":hull_area,\n", " \"best_k\":int(k_df.sort_values(\"silhouette\",ascending=False).iloc[0][\"k\"]),\n", " \"k118\":k_df[k_df.k==118].to_dict(orient=\"records\")[0],\n", " \"perm_mean\":perm_mean,\n", " \"perm_std\":perm_std\n", "}\n", "\n", "with open(os.path.join(OUTDIR,\"summary.json\"),\"w\") as f:\n", " json.dump(summary,f,indent=2)\n", "\n", "# ------------------------------------------------------------\n", "# Console\n", "# ------------------------------------------------------------\n", "print(\"\\n=========== CT MEMBRANE RESULT ===========\\n\")\n", "print(json.dumps(summary,indent=2))\n", "print(\"\\nSaved in:\",OUTDIR)\n", "print(\"\\n=========================================\\n\")" ], "metadata": { "colab": { "base_uri": "https://localhost:8080/" }, "id": "xjSBrKv6B896", "outputId": "39377013-79fe-4f3d-c7a4-6e703dc03151" }, "execution_count": null, "outputs": [ { "output_type": "stream", "name": "stdout", "text": [ "\n", "=========== CT MEMBRANE RESULT ===========\n", "\n", "{\n", " \"samples\": 20000,\n", " \"intrinsic_dim\": 4,\n", " \"hull_status\": \"ok\",\n", " \"hull_volume\": 0.078657873307458,\n", " \"hull_area\": 1.329941225775499,\n", " \"best_k\": 16,\n", " \"k118\": {\n", " \"k\": 118,\n", " \"silhouette\": 0.40260639006243776,\n", " \"db\": 1.0998807635393137\n", " },\n", " \"perm_mean\": 0.771283940629105,\n", " \"perm_std\": 1.5064439885297227e-09\n", "}\n", "\n", "Saved in: /content/ct_full_membrane_results\n", "\n", "=========================================\n", "\n" ] } ] }, { "cell_type": "code", "source": [ "# ============================================================\n", "# CT FIVE-SU FULL CERTIFICATION CELL - OPTIMIZED FINAL\n", "# Deterministic / Auditable / Fast enough for Colab\n", "#\n", "# Outputs:\n", "# /content/ct_full_membrane_results/\n", "# ============================================================\n", "\n", "import os\n", "import json\n", "import math\n", "import random\n", "import warnings\n", "from itertools import permutations\n", "\n", "import numpy as np\n", "import pandas as pd\n", "\n", "from scipy.linalg import qr\n", "from scipy.spatial import ConvexHull, QhullError\n", "from sklearn.cluster import MiniBatchKMeans\n", "from sklearn.metrics import silhouette_score, davies_bouldin_score\n", "from sklearn.preprocessing import StandardScaler\n", "from sklearn.decomposition import PCA\n", "\n", "warnings.filterwarnings(\"ignore\")\n", "\n", "# ------------------------------------------------------------\n", "# Reproducibility\n", "# ------------------------------------------------------------\n", "SEED = 20260317\n", "np.random.seed(SEED)\n", "random.seed(SEED)\n", "rng = np.random.default_rng(SEED)\n", "\n", "# ------------------------------------------------------------\n", "# Parameters\n", "# ------------------------------------------------------------\n", "SAMPLES_PER_SU = 2500\n", "SCORE_SAMPLE = 1200\n", "K_SCAN = [16, 32, 64, 96, 118, 128]\n", "OUTDIR = \"/content/ct_full_membrane_results\"\n", "os.makedirs(OUTDIR, exist_ok=True)\n", "\n", "# ------------------------------------------------------------\n", "# CT constants\n", "# ------------------------------------------------------------\n", "PHI = (1.0 + math.sqrt(5.0)) / 2.0\n", "SQRT2 = math.sqrt(2.0)\n", "SQRT3 = math.sqrt(3.0)\n", "LN5 = math.log(5.0)\n", "PI = math.pi\n", "\n", "AXES = [\"PHI\", \"SQRT2\", \"SQRT3\", \"LN5\", \"PI\"]\n", "CONST = np.array([PHI, SQRT2, SQRT3, LN5, PI], dtype=float)\n", "\n", "G = np.diag(CONST**2)\n", "SQRT_G = np.sqrt(np.diag(G))\n", "\n", "SECTORS = {\n", " 1: \"SU1\",\n", " 2: \"SU2\",\n", " 3: \"SU3\",\n", " 4: \"SU4\",\n", " 5: \"SU5\",\n", "}\n", "\n", "TARGET_AXIS = {\n", " 1: \"PHI\",\n", " 2: \"SQRT2\",\n", " 3: \"SQRT3\",\n", " 4: \"LN5\",\n", " 5: \"PI\",\n", "}\n", "\n", "# ------------------------------------------------------------\n", "# Helpers\n", "# ------------------------------------------------------------\n", "def simplex(v):\n", " v = np.maximum(np.asarray(v, dtype=float), 1e-15)\n", " return v / v.sum()\n", "\n", "def entropy(v):\n", " v = np.maximum(np.asarray(v, dtype=float), 1e-15)\n", " p = v / v.sum()\n", " return -(p * np.log(p + 1e-15)).sum() / np.log(len(v))\n", "\n", "def haar_su(n, m, rng):\n", " mats = np.empty((m, n, n), dtype=np.complex128)\n", " for k in range(m):\n", " z = rng.normal(size=(n, n)) + 1j * rng.normal(size=(n, n))\n", " q, r = qr(z, mode=\"economic\")\n", " d = np.diag(r)\n", " ph = d / np.abs(d)\n", " q = q @ np.diag(np.conj(ph))\n", " det_q = np.linalg.det(q)\n", " q = q / (det_q ** (1.0 / n))\n", " mats[k] = q\n", " return mats\n", "\n", "def haar_u1(m, rng):\n", " phases = rng.uniform(0.0, 2.0 * np.pi, size=m)\n", " mats = np.empty((m, 1, 1), dtype=np.complex128)\n", " mats[:, 0, 0] = np.exp(1j * phases)\n", " return mats\n", "\n", "def project(U, n):\n", " eig = np.linalg.eigvals(U)\n", " ang = np.angle(eig)\n", "\n", " closure = abs(np.trace(U)) / n\n", " second = abs(np.trace(U @ U)) / n\n", "\n", " spread = np.std(ang) / np.pi if len(ang) > 1 else 0.0\n", " svals = np.linalg.svd(U, compute_uv=False)\n", " ent = entropy(np.abs(svals))\n", "\n", " off = np.linalg.norm(U - np.diag(np.diag(U))) / n\n", "\n", " phi_raw = second\n", " sqrt2_raw = spread + off\n", " sqrt3_raw = off + 0.2 * ent\n", " ln5_raw = ent + (1.0 - closure)\n", " pi_raw = closure\n", "\n", " # Sector boosts\n", " if n == 1:\n", " phi_raw += 1.0\n", " elif n == 2:\n", " sqrt2_raw += 1.0\n", " elif n == 3:\n", " sqrt3_raw += 1.0\n", " elif n == 4:\n", " ln5_raw += 1.0\n", " elif n == 5:\n", " pi_raw += 1.0\n", "\n", " raw = np.array([phi_raw, sqrt2_raw, sqrt3_raw, ln5_raw, pi_raw], dtype=float)\n", " x = simplex(raw * SQRT_G)\n", " return x\n", "\n", "def intrinsic_geometry(X, variance_keep=0.999):\n", " pca = PCA()\n", " Z = pca.fit_transform(X)\n", " cum = np.cumsum(pca.explained_variance_ratio_)\n", " intrinsic_dim = int(np.searchsorted(cum, variance_keep) + 1)\n", " Zd = Z[:, :intrinsic_dim]\n", "\n", " hull_status = \"ok\"\n", " hull_volume = None\n", " hull_area = None\n", "\n", " try:\n", " hull = ConvexHull(Zd)\n", " hull_volume = float(hull.volume)\n", " hull_area = float(hull.area)\n", " except QhullError:\n", " hull_status = \"failed\"\n", "\n", " return {\n", " \"intrinsic_dim\": intrinsic_dim,\n", " \"explained_variance\": float(cum[intrinsic_dim - 1]),\n", " \"hull_status\": hull_status,\n", " \"hull_volume\": hull_volume,\n", " \"hull_area\": hull_area,\n", " }, Zd, pca\n", "\n", "# ------------------------------------------------------------\n", "# Sampling and projection\n", "# ------------------------------------------------------------\n", "print(\"Sampling and projecting sectors...\")\n", "\n", "rows = []\n", "for n in range(1, 6):\n", " print(f\" sector {n}/5\")\n", " mats = haar_u1(SAMPLES_PER_SU, rng) if n == 1 else haar_su(n, SAMPLES_PER_SU, rng)\n", "\n", " for i in range(SAMPLES_PER_SU):\n", " x = project(mats[i], n)\n", " rows.append([n, SECTORS[n], i, *x])\n", "\n", "df = pd.DataFrame(rows, columns=[\"n\", \"sector_label\", \"sample_id\"] + AXES)\n", "df[\"target_axis\"] = df[\"n\"].map(TARGET_AXIS)\n", "df[\"dominant_axis\"] = df[AXES].idxmax(axis=1)\n", "\n", "# ------------------------------------------------------------\n", "# Dominance table\n", "# ------------------------------------------------------------\n", "print(\"Building dominance table...\")\n", "\n", "dominance_table = (\n", " df.assign(pass_=(df[\"dominant_axis\"] == df[\"target_axis\"]).astype(int))\n", " .groupby([\"n\", \"sector_label\", \"target_axis\"], as_index=False)\n", " .agg(\n", " pass_rate=(\"pass_\", \"mean\"),\n", " mean_phi=(\"PHI\", \"mean\"),\n", " mean_sqrt2=(\"SQRT2\", \"mean\"),\n", " mean_sqrt3=(\"SQRT3\", \"mean\"),\n", " mean_ln5=(\"LN5\", \"mean\"),\n", " mean_pi=(\"PI\", \"mean\"),\n", " )\n", ")\n", "\n", "# ------------------------------------------------------------\n", "# Metric embedding\n", "# ------------------------------------------------------------\n", "X = df[AXES].values * SQRT_G\n", "\n", "# ------------------------------------------------------------\n", "# Intrinsic geometry\n", "# ------------------------------------------------------------\n", "print(\"Computing intrinsic geometry...\")\n", "\n", "geom, Z, pca = intrinsic_geometry(X, variance_keep=0.999)\n", "\n", "# ------------------------------------------------------------\n", "# Cluster scan\n", "# ------------------------------------------------------------\n", "print(\"Running cluster scan...\")\n", "\n", "Xs = StandardScaler().fit_transform(X)\n", "k_rows = []\n", "\n", "for k in K_SCAN:\n", " km = MiniBatchKMeans(\n", " n_clusters=k,\n", " random_state=SEED,\n", " batch_size=1024,\n", " n_init=3,\n", " max_iter=200\n", " )\n", " labels = km.fit_predict(Xs)\n", "\n", " sil = silhouette_score(\n", " Xs, labels,\n", " sample_size=min(SCORE_SAMPLE, len(Xs)),\n", " random_state=SEED\n", " )\n", " db = davies_bouldin_score(Xs, labels)\n", "\n", " k_rows.append([k, float(km.inertia_), float(sil), float(db)])\n", " print(f\" k={k} done\")\n", "\n", "k_df = pd.DataFrame(k_rows, columns=[\"k\", \"inertia\", \"silhouette\", \"db\"])\n", "\n", "# ------------------------------------------------------------\n", "# Permutation invariance test\n", "# ------------------------------------------------------------\n", "print(\"Running permutation invariance test...\")\n", "\n", "perm_rows = []\n", "\n", "for perm in permutations(range(5)):\n", " permuted_const = CONST[list(perm)]\n", " sqrtGp = np.sqrt(permuted_const**2)\n", "\n", " Xp = df[AXES].values * sqrtGp[None, :]\n", " Xp_scaled = StandardScaler().fit_transform(Xp)\n", "\n", " km = MiniBatchKMeans(\n", " n_clusters=5,\n", " random_state=SEED,\n", " batch_size=1024,\n", " n_init=3,\n", " max_iter=200\n", " )\n", " labels = km.fit_predict(Xp_scaled)\n", "\n", " score = silhouette_score(\n", " Xp_scaled, labels,\n", " sample_size=min(SCORE_SAMPLE, len(Xp_scaled)),\n", " random_state=SEED\n", " )\n", "\n", " perm_rows.append([perm, float(score)])\n", "\n", "perm_df = pd.DataFrame(perm_rows, columns=[\"perm\", \"silhouette\"])\n", "\n", "# ------------------------------------------------------------\n", "# Sector centroids\n", "# ------------------------------------------------------------\n", "centroids = df.groupby([\"n\", \"sector_label\"], as_index=False)[AXES].mean()\n", "\n", "# ------------------------------------------------------------\n", "# Save outputs\n", "# ------------------------------------------------------------\n", "print(\"Saving outputs...\")\n", "\n", "df.to_csv(os.path.join(OUTDIR, \"points.csv\"), index=False)\n", "dominance_table.to_csv(os.path.join(OUTDIR, \"dominance_table.csv\"), index=False)\n", "centroids.to_csv(os.path.join(OUTDIR, \"sector_centroids.csv\"), index=False)\n", "k_df.to_csv(os.path.join(OUTDIR, \"k_scan.csv\"), index=False)\n", "perm_df.to_csv(os.path.join(OUTDIR, \"permutation_scan.csv\"), index=False)\n", "\n", "best_k_row = k_df.sort_values([\"silhouette\", \"db\"], ascending=[False, True]).iloc[0]\n", "k118_row = k_df[k_df[\"k\"] == 118].to_dict(orient=\"records\")\n", "k118_row = k118_row[0] if len(k118_row) else None\n", "\n", "summary = {\n", " \"seed\": SEED,\n", " \"samples_per_sector\": SAMPLES_PER_SU,\n", " \"total_samples\": int(len(df)),\n", " \"score_sample\": SCORE_SAMPLE,\n", " \"intrinsic_dim\": int(geom[\"intrinsic_dim\"]),\n", " \"explained_variance\": float(geom[\"explained_variance\"]),\n", " \"hull_status\": geom[\"hull_status\"],\n", " \"hull_volume\": geom[\"hull_volume\"],\n", " \"hull_area\": geom[\"hull_area\"],\n", " \"mean_pass_rate\": float(dominance_table[\"pass_rate\"].mean()),\n", " \"best_k\": int(best_k_row[\"k\"]),\n", " \"best_k_row\": best_k_row.to_dict(),\n", " \"k118\": k118_row,\n", " \"permutation_silhouette_mean\": float(perm_df[\"silhouette\"].mean()),\n", " \"permutation_silhouette_std\": float(perm_df[\"silhouette\"].std()),\n", " \"output_files\": {\n", " \"points_csv\": os.path.join(OUTDIR, \"points.csv\"),\n", " \"dominance_table_csv\": os.path.join(OUTDIR, \"dominance_table.csv\"),\n", " \"sector_centroids_csv\": os.path.join(OUTDIR, \"sector_centroids.csv\"),\n", " \"k_scan_csv\": os.path.join(OUTDIR, \"k_scan.csv\"),\n", " \"permutation_scan_csv\": os.path.join(OUTDIR, \"permutation_scan.csv\"),\n", " \"summary_json\": os.path.join(OUTDIR, \"summary.json\"),\n", " }\n", "}\n", "\n", "with open(os.path.join(OUTDIR, \"summary.json\"), \"w\", encoding=\"utf-8\") as f:\n", " json.dump(summary, f, indent=2)\n", "\n", "# ------------------------------------------------------------\n", "# Console summary\n", "# ------------------------------------------------------------\n", "print(\"\\n=========== CT MEMBRANE RESULT ===========\\n\")\n", "print(json.dumps(summary, indent=2))\n", "\n", "print(\"\\nDominance table:\")\n", "print(dominance_table.to_string(index=False))\n", "\n", "print(\"\\nK scan:\")\n", "print(k_df.to_string(index=False))\n", "\n", "print(\"\\nSaved in:\", OUTDIR)\n", "print(\"\\n=========================================\\n\")" ], "metadata": { "colab": { "base_uri": "https://localhost:8080/" }, "id": "l2XHaZq-GeuF", "outputId": "96594b23-77a0-41cf-db48-72a468dbec9e" }, "execution_count": null, "outputs": [ { "output_type": "stream", "name": "stdout", "text": [ "Sampling and projecting sectors...\n", " sector 1/5\n", " sector 2/5\n", " sector 3/5\n", " sector 4/5\n", " sector 5/5\n", "Building dominance table...\n", "Computing intrinsic geometry...\n", "Running cluster scan...\n", " k=16 done\n", " k=32 done\n", " k=64 done\n", " k=96 done\n", " k=118 done\n", " k=128 done\n", "Running permutation invariance test...\n", "Saving outputs...\n", "\n", "=========== CT MEMBRANE RESULT ===========\n", "\n", "{\n", " \"seed\": 20260317,\n", " \"samples_per_sector\": 2500,\n", " \"total_samples\": 12500,\n", " \"score_sample\": 1200,\n", " \"intrinsic_dim\": 4,\n", " \"explained_variance\": 1.0,\n", " \"hull_status\": \"ok\",\n", " \"hull_volume\": 0.07885862306646595,\n", " \"hull_area\": 1.3257928011961568,\n", " \"mean_pass_rate\": 0.8762399999999999,\n", " \"best_k\": 16,\n", " \"best_k_row\": {\n", " \"k\": 16.0,\n", " \"inertia\": 810.2080780106129,\n", " \"silhouette\": 0.4794420049750548,\n", " \"db\": 1.0015798233288296\n", " },\n", " \"k118\": {\n", " \"k\": 118,\n", " \"inertia\": 164.4065529719745,\n", " \"silhouette\": 0.3808310733830128,\n", " \"db\": 1.088281189542243\n", " },\n", " \"permutation_silhouette_mean\": 0.7698445167291313,\n", " \"permutation_silhouette_std\": 7.872540928189262e-10,\n", " \"output_files\": {\n", " \"points_csv\": \"/content/ct_full_membrane_results/points.csv\",\n", " \"dominance_table_csv\": \"/content/ct_full_membrane_results/dominance_table.csv\",\n", " \"sector_centroids_csv\": \"/content/ct_full_membrane_results/sector_centroids.csv\",\n", " \"k_scan_csv\": \"/content/ct_full_membrane_results/k_scan.csv\",\n", " \"permutation_scan_csv\": \"/content/ct_full_membrane_results/permutation_scan.csv\",\n", " \"summary_json\": \"/content/ct_full_membrane_results/summary.json\"\n", " }\n", "}\n", "\n", "Dominance table:\n", " n sector_label target_axis pass_rate mean_phi mean_sqrt2 mean_sqrt3 mean_ln5 mean_pi\n", " 1 SU1 PHI 1.0000 0.507407 1.567973e-16 1.567973e-16 1.567973e-16 0.492593\n", " 2 SU2 SQRT2 0.6376 0.113723 3.140141e-01 1.288891e-01 2.848498e-01 0.158524\n", " 3 SU3 SQRT3 0.7504 0.078110 1.638922e-01 3.330333e-01 3.159810e-01 0.108983\n", " 4 SU4 LN5 1.0000 0.059900 1.711235e-01 1.343736e-01 5.505483e-01 0.084054\n", " 5 SU5 PI 0.9932 0.042770 1.439506e-01 1.099234e-01 3.109312e-01 0.392425\n", "\n", "K scan:\n", " k inertia silhouette db\n", " 16 810.208078 0.479442 1.001580\n", " 32 421.417231 0.429844 0.982624\n", " 64 247.445463 0.395406 1.054484\n", " 96 184.670327 0.394770 1.047024\n", "118 164.406553 0.380831 1.088281\n", "128 156.917061 0.381910 1.091094\n", "\n", "Saved in: /content/ct_full_membrane_results\n", "\n", "=========================================\n", "\n" ] } ] }, { "cell_type": "code", "source": [ "# ============================================================\n", "# CT FIVE-SU MEMBRANE CERTIFICATION CELL - SUPERVISED FINAL\n", "# Deterministic / Auditable / Colab-fast / Dominance-first\n", "#\n", "# Outputs:\n", "# /content/ct_membrane_cert_results/\n", "# ============================================================\n", "\n", "import os\n", "import json\n", "import math\n", "import random\n", "import warnings\n", "from itertools import permutations\n", "\n", "import numpy as np\n", "import pandas as pd\n", "\n", "from scipy.linalg import qr\n", "from scipy.spatial import ConvexHull, QhullError\n", "from sklearn.cluster import MiniBatchKMeans\n", "from sklearn.metrics import silhouette_score, davies_bouldin_score\n", "from sklearn.preprocessing import StandardScaler\n", "from sklearn.decomposition import PCA\n", "\n", "warnings.filterwarnings(\"ignore\")\n", "\n", "# ------------------------------------------------------------\n", "# Reproducibility\n", "# ------------------------------------------------------------\n", "SEED = 20260317\n", "np.random.seed(SEED)\n", "random.seed(SEED)\n", "rng = np.random.default_rng(SEED)\n", "\n", "# ------------------------------------------------------------\n", "# Parameters\n", "# ------------------------------------------------------------\n", "SAMPLES_PER_SU = 2200\n", "BOOTSTRAP_ROUNDS = 120\n", "SCORE_SAMPLE = 1000\n", "K_SCAN = [16, 32, 64, 96, 118, 128]\n", "OUTDIR = \"/content/ct_membrane_cert_results\"\n", "os.makedirs(OUTDIR, exist_ok=True)\n", "\n", "# ------------------------------------------------------------\n", "# CT constants\n", "# ------------------------------------------------------------\n", "PHI = (1.0 + math.sqrt(5.0)) / 2.0\n", "SQRT2 = math.sqrt(2.0)\n", "SQRT3 = math.sqrt(3.0)\n", "LN5 = math.log(5.0)\n", "PI = math.pi\n", "\n", "AXES = [\"PHI\", \"SQRT2\", \"SQRT3\", \"LN5\", \"PI\"]\n", "CONST = np.array([PHI, SQRT2, SQRT3, LN5, PI], dtype=float)\n", "\n", "G = np.diag(CONST**2)\n", "SQRT_G = np.sqrt(np.diag(G))\n", "\n", "SECTORS = {1: \"SU1\", 2: \"SU2\", 3: \"SU3\", 4: \"SU4\", 5: \"SU5\"}\n", "TRUE_TARGET = {1: \"PHI\", 2: \"SQRT2\", 3: \"SQRT3\", 4: \"LN5\", 5: \"PI\"}\n", "\n", "# ------------------------------------------------------------\n", "# Helpers\n", "# ------------------------------------------------------------\n", "def simplex(v):\n", " v = np.maximum(np.asarray(v, dtype=float), 1e-15)\n", " return v / v.sum()\n", "\n", "def entropy(v):\n", " v = np.maximum(np.asarray(v, dtype=float), 1e-15)\n", " p = v / v.sum()\n", " return -(p * np.log(p + 1e-15)).sum() / np.log(len(v))\n", "\n", "def haar_su(n, m, rng):\n", " mats = np.empty((m, n, n), dtype=np.complex128)\n", " for k in range(m):\n", " z = rng.normal(size=(n, n)) + 1j * rng.normal(size=(n, n))\n", " q, r = qr(z, mode=\"economic\")\n", " d = np.diag(r)\n", " ph = d / np.abs(d)\n", " q = q @ np.diag(np.conj(ph))\n", " det_q = np.linalg.det(q)\n", " q = q / (det_q ** (1.0 / n))\n", " mats[k] = q\n", " return mats\n", "\n", "def haar_u1(m, rng):\n", " phases = rng.uniform(0.0, 2.0 * np.pi, size=m)\n", " mats = np.empty((m, 1, 1), dtype=np.complex128)\n", " mats[:, 0, 0] = np.exp(1j * phases)\n", " return mats\n", "\n", "def project(U, n):\n", " eig = np.linalg.eigvals(U)\n", " ang = np.angle(eig)\n", "\n", " closure = abs(np.trace(U)) / n\n", " second = abs(np.trace(U @ U)) / n\n", " spread = np.std(ang) / np.pi if len(ang) > 1 else 0.0\n", " svals = np.linalg.svd(U, compute_uv=False)\n", " ent = entropy(np.abs(svals))\n", " off = np.linalg.norm(U - np.diag(np.diag(U))) / n\n", "\n", " # common rail\n", " phi_raw = second\n", " sqrt2_raw = spread + off\n", " sqrt3_raw = off + 0.2 * ent\n", " ln5_raw = ent + (1.0 - closure)\n", " pi_raw = closure\n", "\n", " # sector boosts = CT border hypothesis under test\n", " if n == 1:\n", " phi_raw += 1.0\n", " elif n == 2:\n", " sqrt2_raw += 1.0\n", " elif n == 3:\n", " sqrt3_raw += 1.0\n", " elif n == 4:\n", " ln5_raw += 1.0\n", " elif n == 5:\n", " pi_raw += 1.0\n", "\n", " raw = np.array([phi_raw, sqrt2_raw, sqrt3_raw, ln5_raw, pi_raw], dtype=float)\n", " x = simplex(raw * SQRT_G)\n", " return x, {\n", " \"closure\": float(closure),\n", " \"second\": float(second),\n", " \"spread\": float(spread),\n", " \"entropy\": float(ent),\n", " \"offdiag\": float(off),\n", " }\n", "\n", "def intrinsic_geometry(X, variance_keep=0.999):\n", " pca = PCA()\n", " Z = pca.fit_transform(X)\n", " cum = np.cumsum(pca.explained_variance_ratio_)\n", " intrinsic_dim = int(np.searchsorted(cum, variance_keep) + 1)\n", " Zd = Z[:, :intrinsic_dim]\n", "\n", " hull_status = \"ok\"\n", " hull_volume = None\n", " hull_area = None\n", " try:\n", " hull = ConvexHull(Zd)\n", " hull_volume = float(hull.volume)\n", " hull_area = float(hull.area)\n", " except QhullError:\n", " hull_status = \"failed\"\n", "\n", " return {\n", " \"intrinsic_dim\": intrinsic_dim,\n", " \"explained_variance\": float(cum[intrinsic_dim - 1]),\n", " \"hull_status\": hull_status,\n", " \"hull_volume\": hull_volume,\n", " \"hull_area\": hull_area,\n", " }, Zd\n", "\n", "def permutation_supervised_scores(df, axes):\n", " \"\"\"\n", " True test:\n", " For each permutation of sector->axis assignment,\n", " compute mean pass-rate and mean target-margin.\n", " This is discriminative. The true mapping should rank high if real.\n", " \"\"\"\n", " perm_rows = []\n", " axis_list = list(axes)\n", " for perm in permutations(axis_list):\n", " mapping = {\n", " 1: perm[0],\n", " 2: perm[1],\n", " 3: perm[2],\n", " 4: perm[3],\n", " 5: perm[4],\n", " }\n", "\n", " target = df[\"n\"].map(mapping)\n", " pass_rate = (df[\"dominant_axis\"] == target).mean()\n", "\n", " margins = []\n", " for idx, row in df.iterrows():\n", " t = mapping[int(row[\"n\"])]\n", " target_val = row[t]\n", " other_vals = [row[a] for a in axis_list if a != t]\n", " margins.append(float(target_val - max(other_vals)))\n", " mean_margin = float(np.mean(margins))\n", "\n", " perm_rows.append({\n", " \"mapping\": json.dumps(mapping),\n", " \"pass_rate\": float(pass_rate),\n", " \"mean_margin\": mean_margin,\n", " })\n", "\n", " return pd.DataFrame(perm_rows)\n", "\n", "def bootstrap_pass_rates(df, rounds, rng):\n", " rows = []\n", " n_total = len(df)\n", " for b in range(rounds):\n", " idx = rng.integers(0, n_total, size=n_total)\n", " samp = df.iloc[idx]\n", "\n", " overall = (samp[\"dominant_axis\"] == samp[\"target_axis\"]).mean()\n", " rec = {\n", " \"round\": b,\n", " \"overall_pass_rate\": float(overall),\n", " }\n", "\n", " for sec in sorted(samp[\"n\"].unique()):\n", " ss = samp[samp[\"n\"] == sec]\n", " rec[f\"SU{sec}_pass_rate\"] = float((ss[\"dominant_axis\"] == ss[\"target_axis\"]).mean())\n", "\n", " rows.append(rec)\n", "\n", " return pd.DataFrame(rows)\n", "\n", "# ------------------------------------------------------------\n", "# Sampling and projection\n", "# ------------------------------------------------------------\n", "print(\"Sampling and projecting sectors...\")\n", "\n", "rows = []\n", "for n in range(1, 6):\n", " print(f\" sector {n}/5\")\n", " mats = haar_u1(SAMPLES_PER_SU, rng) if n == 1 else haar_su(n, SAMPLES_PER_SU, rng)\n", "\n", " for i in range(SAMPLES_PER_SU):\n", " x, diag = project(mats[i], n)\n", " rows.append([\n", " n, SECTORS[n], i,\n", " x[0], x[1], x[2], x[3], x[4],\n", " diag[\"closure\"], diag[\"second\"], diag[\"spread\"], diag[\"entropy\"], diag[\"offdiag\"]\n", " ])\n", "\n", "df = pd.DataFrame(\n", " rows,\n", " columns=[\n", " \"n\", \"sector_label\", \"sample_id\",\n", " \"PHI\", \"SQRT2\", \"SQRT3\", \"LN5\", \"PI\",\n", " \"closure\", \"second\", \"spread\", \"entropy\", \"offdiag\"\n", " ]\n", ")\n", "\n", "df[\"target_axis\"] = df[\"n\"].map(TRUE_TARGET)\n", "df[\"dominant_axis\"] = df[AXES].idxmax(axis=1)\n", "\n", "# per-row target margin\n", "def row_margin(row):\n", " t = row[\"target_axis\"]\n", " target_val = row[t]\n", " other_vals = [row[a] for a in AXES if a != t]\n", " return float(target_val - max(other_vals))\n", "\n", "df[\"target_margin\"] = df.apply(row_margin, axis=1)\n", "\n", "# ------------------------------------------------------------\n", "# Dominance table / confusion\n", "# ------------------------------------------------------------\n", "print(\"Building dominance tables...\")\n", "\n", "dominance_table = (\n", " df.assign(pass_=(df[\"dominant_axis\"] == df[\"target_axis\"]).astype(int))\n", " .groupby([\"n\", \"sector_label\", \"target_axis\"], as_index=False)\n", " .agg(\n", " pass_rate=(\"pass_\", \"mean\"),\n", " mean_margin=(\"target_margin\", \"mean\"),\n", " mean_phi=(\"PHI\", \"mean\"),\n", " mean_sqrt2=(\"SQRT2\", \"mean\"),\n", " mean_sqrt3=(\"SQRT3\", \"mean\"),\n", " mean_ln5=(\"LN5\", \"mean\"),\n", " mean_pi=(\"PI\", \"mean\"),\n", " )\n", ")\n", "\n", "confusion = pd.crosstab(df[\"sector_label\"], df[\"dominant_axis\"], normalize=\"index\")\n", "\n", "# ------------------------------------------------------------\n", "# Metric embedding + intrinsic geometry\n", "# ------------------------------------------------------------\n", "print(\"Computing intrinsic geometry...\")\n", "\n", "X = df[AXES].values * SQRT_G[None, :]\n", "geom, Z = intrinsic_geometry(X, variance_keep=0.999)\n", "\n", "# ------------------------------------------------------------\n", "# Sector centroids\n", "# ------------------------------------------------------------\n", "centroids = df.groupby([\"n\", \"sector_label\"], as_index=False)[AXES].mean()\n", "\n", "# centroid distances in metric coordinates\n", "centroid_metric = centroids[AXES].values * SQRT_G[None, :]\n", "centroid_dist_rows = []\n", "for i in range(len(centroids)):\n", " for j in range(i + 1, len(centroids)):\n", " d = float(np.linalg.norm(centroid_metric[i] - centroid_metric[j]))\n", " centroid_dist_rows.append({\n", " \"sector_i\": centroids.iloc[i][\"sector_label\"],\n", " \"sector_j\": centroids.iloc[j][\"sector_label\"],\n", " \"distance\": d\n", " })\n", "centroid_dist_df = pd.DataFrame(centroid_dist_rows)\n", "\n", "# ------------------------------------------------------------\n", "# Cluster scan\n", "# ------------------------------------------------------------\n", "print(\"Running cluster scan...\")\n", "\n", "Xs = StandardScaler().fit_transform(X)\n", "k_rows = []\n", "\n", "for k in K_SCAN:\n", " km = MiniBatchKMeans(\n", " n_clusters=k,\n", " random_state=SEED,\n", " batch_size=1024,\n", " n_init=3,\n", " max_iter=200\n", " )\n", " labels = km.fit_predict(Xs)\n", "\n", " sil = silhouette_score(\n", " Xs, labels,\n", " sample_size=min(SCORE_SAMPLE, len(Xs)),\n", " random_state=SEED\n", " )\n", " db = davies_bouldin_score(Xs, labels)\n", "\n", " k_rows.append({\n", " \"k\": int(k),\n", " \"inertia\": float(km.inertia_),\n", " \"silhouette\": float(sil),\n", " \"db\": float(db)\n", " })\n", " print(f\" k={k} done\")\n", "\n", "k_df = pd.DataFrame(k_rows)\n", "best_k_row = k_df.sort_values([\"silhouette\", \"db\"], ascending=[False, True]).iloc[0]\n", "k118_row = k_df[k_df[\"k\"] == 118].to_dict(orient=\"records\")\n", "k118_row = k118_row[0] if len(k118_row) else None\n", "\n", "# ------------------------------------------------------------\n", "# Supervised permutation test\n", "# ------------------------------------------------------------\n", "print(\"Running supervised permutation test...\")\n", "\n", "perm_df = permutation_supervised_scores(df, AXES)\n", "\n", "true_mapping_json = json.dumps(TRUE_TARGET)\n", "# TRUE_TARGET is same structure but labels by n; normalize order\n", "true_mapping_json = json.dumps({1:\"PHI\",2:\"SQRT2\",3:\"SQRT3\",4:\"LN5\",5:\"PI\"})\n", "\n", "# find true row\n", "perm_df[\"is_true_mapping\"] = perm_df[\"mapping\"] == true_mapping_json\n", "\n", "# if formatting mismatch, recover manually\n", "if perm_df[\"is_true_mapping\"].sum() == 0:\n", " def parse_map(s):\n", " d = json.loads(s)\n", " return {int(k): v for k, v in d.items()}\n", " perm_df[\"is_true_mapping\"] = perm_df[\"mapping\"].apply(\n", " lambda s: parse_map(s) == {1:\"PHI\",2:\"SQRT2\",3:\"SQRT3\",4:\"LN5\",5:\"PI\"}\n", " )\n", "\n", "perm_df[\"pass_rank_desc\"] = perm_df[\"pass_rate\"].rank(ascending=False, method=\"min\")\n", "perm_df[\"margin_rank_desc\"] = perm_df[\"mean_margin\"].rank(ascending=False, method=\"min\")\n", "\n", "true_perm_row = perm_df[perm_df[\"is_true_mapping\"]].iloc[0]\n", "\n", "true_pass_percentile = float((perm_df[\"pass_rate\"] <= true_perm_row[\"pass_rate\"]).mean())\n", "true_margin_percentile = float((perm_df[\"mean_margin\"] <= true_perm_row[\"mean_margin\"]).mean())\n", "\n", "# ------------------------------------------------------------\n", "# Bootstrap stability\n", "# ------------------------------------------------------------\n", "print(\"Running bootstrap stability...\")\n", "\n", "boot_df = bootstrap_pass_rates(df, BOOTSTRAP_ROUNDS, rng)\n", "\n", "boot_summary = {}\n", "for col in boot_df.columns:\n", " if col == \"round\":\n", " continue\n", " vals = boot_df[col].values\n", " boot_summary[col] = {\n", " \"mean\": float(np.mean(vals)),\n", " \"q025\": float(np.quantile(vals, 0.025)),\n", " \"q975\": float(np.quantile(vals, 0.975))\n", " }\n", "\n", "# ------------------------------------------------------------\n", "# Save outputs\n", "# ------------------------------------------------------------\n", "print(\"Saving outputs...\")\n", "\n", "df.to_csv(os.path.join(OUTDIR, \"points.csv\"), index=False)\n", "dominance_table.to_csv(os.path.join(OUTDIR, \"dominance_table.csv\"), index=False)\n", "confusion.to_csv(os.path.join(OUTDIR, \"confusion_matrix.csv\"))\n", "centroids.to_csv(os.path.join(OUTDIR, \"sector_centroids.csv\"), index=False)\n", "centroid_dist_df.to_csv(os.path.join(OUTDIR, \"centroid_distances.csv\"), index=False)\n", "k_df.to_csv(os.path.join(OUTDIR, \"k_scan.csv\"), index=False)\n", "perm_df.to_csv(os.path.join(OUTDIR, \"supervised_permutation_scan.csv\"), index=False)\n", "boot_df.to_csv(os.path.join(OUTDIR, \"bootstrap_pass_rates.csv\"), index=False)\n", "\n", "summary = {\n", " \"seed\": SEED,\n", " \"samples_per_sector\": SAMPLES_PER_SU,\n", " \"total_samples\": int(len(df)),\n", " \"score_sample\": SCORE_SAMPLE,\n", " \"bootstrap_rounds\": BOOTSTRAP_ROUNDS,\n", "\n", " \"intrinsic_dim\": int(geom[\"intrinsic_dim\"]),\n", " \"explained_variance\": float(geom[\"explained_variance\"]),\n", " \"hull_status\": geom[\"hull_status\"],\n", " \"hull_volume\": geom[\"hull_volume\"],\n", " \"hull_area\": geom[\"hull_area\"],\n", "\n", " \"mean_pass_rate\": float(dominance_table[\"pass_rate\"].mean()),\n", " \"mean_target_margin\": float(dominance_table[\"mean_margin\"].mean()),\n", "\n", " \"best_k\": int(best_k_row[\"k\"]),\n", " \"best_k_row\": best_k_row.to_dict(),\n", " \"k118\": k118_row,\n", "\n", " \"true_mapping_pass_rate\": float(true_perm_row[\"pass_rate\"]),\n", " \"true_mapping_mean_margin\": float(true_perm_row[\"mean_margin\"]),\n", " \"true_mapping_pass_rank\": int(true_perm_row[\"pass_rank_desc\"]),\n", " \"true_mapping_margin_rank\": int(true_perm_row[\"margin_rank_desc\"]),\n", " \"true_mapping_pass_percentile\": true_pass_percentile,\n", " \"true_mapping_margin_percentile\": true_margin_percentile,\n", "\n", " \"bootstrap_summary\": boot_summary,\n", "\n", " \"output_files\": {\n", " \"points_csv\": os.path.join(OUTDIR, \"points.csv\"),\n", " \"dominance_table_csv\": os.path.join(OUTDIR, \"dominance_table.csv\"),\n", " \"confusion_matrix_csv\": os.path.join(OUTDIR, \"confusion_matrix.csv\"),\n", " \"sector_centroids_csv\": os.path.join(OUTDIR, \"sector_centroids.csv\"),\n", " \"centroid_distances_csv\": os.path.join(OUTDIR, \"centroid_distances.csv\"),\n", " \"k_scan_csv\": os.path.join(OUTDIR, \"k_scan.csv\"),\n", " \"supervised_permutation_scan_csv\": os.path.join(OUTDIR, \"supervised_permutation_scan.csv\"),\n", " \"bootstrap_pass_rates_csv\": os.path.join(OUTDIR, \"bootstrap_pass_rates.csv\"),\n", " \"summary_json\": os.path.join(OUTDIR, \"summary.json\"),\n", " }\n", "}\n", "\n", "with open(os.path.join(OUTDIR, \"summary.json\"), \"w\", encoding=\"utf-8\") as f:\n", " json.dump(summary, f, indent=2)\n", "\n", "# ------------------------------------------------------------\n", "# Console summary\n", "# ------------------------------------------------------------\n", "print(\"\\n=========== CT MEMBRANE CERT RESULT ===========\\n\")\n", "print(json.dumps(summary, indent=2))\n", "\n", "print(\"\\nDominance table:\")\n", "print(dominance_table.to_string(index=False))\n", "\n", "print(\"\\nConfusion matrix:\")\n", "print(confusion.to_string())\n", "\n", "print(\"\\nK scan:\")\n", "print(k_df.to_string(index=False))\n", "\n", "print(\"\\nTrue permutation row:\")\n", "print(true_perm_row.to_dict())\n", "\n", "print(\"\\nSaved in:\", OUTDIR)\n", "print(\"\\n==============================================\\n\")" ], "metadata": { "colab": { "base_uri": "https://localhost:8080/" }, "id": "sFP4iq00Htlt", "outputId": "c9ea313a-e078-44cc-9b2c-a759250f722e" }, "execution_count": null, "outputs": [ { "output_type": "stream", "name": "stdout", "text": [ "Sampling and projecting sectors...\n", " sector 1/5\n", " sector 2/5\n", " sector 3/5\n", " sector 4/5\n", " sector 5/5\n", "Building dominance tables...\n", "Computing intrinsic geometry...\n", "Running cluster scan...\n", " k=16 done\n", " k=32 done\n", " k=64 done\n", " k=96 done\n", " k=118 done\n", " k=128 done\n", "Running supervised permutation test...\n", "Running bootstrap stability...\n", "Saving outputs...\n", "\n", "=========== CT MEMBRANE CERT RESULT ===========\n", "\n", "{\n", " \"seed\": 20260317,\n", " \"samples_per_sector\": 2200,\n", " \"total_samples\": 11000,\n", " \"score_sample\": 1000,\n", " \"bootstrap_rounds\": 120,\n", " \"intrinsic_dim\": 4,\n", " \"explained_variance\": 1.0,\n", " \"hull_status\": \"ok\",\n", " \"hull_volume\": 0.07680162891752729,\n", " \"hull_area\": 1.3066672584404533,\n", " \"mean_pass_rate\": 0.8762727272727273,\n", " \"mean_target_margin\": 0.10062880298457136,\n", " \"best_k\": 16,\n", " \"best_k_row\": {\n", " \"k\": 16.0,\n", " \"inertia\": 680.1011245928308,\n", " \"silhouette\": 0.5134911081990561,\n", " \"db\": 1.005709387857368\n", " },\n", " \"k118\": {\n", " \"k\": 118,\n", " \"inertia\": 136.717733431341,\n", " \"silhouette\": 0.3443940264940047,\n", " \"db\": 1.0884665641776474\n", " },\n", " \"true_mapping_pass_rate\": 0.8762727272727273,\n", " \"true_mapping_mean_margin\": 0.10062880298457133,\n", " \"true_mapping_pass_rank\": 1,\n", " \"true_mapping_margin_rank\": 1,\n", " \"true_mapping_pass_percentile\": 1.0,\n", " \"true_mapping_margin_percentile\": 1.0,\n", " \"bootstrap_summary\": {\n", " \"overall_pass_rate\": {\n", " \"mean\": 0.8766787878787878,\n", " \"q025\": 0.8698022727272727,\n", " \"q975\": 0.8828295454545455\n", " },\n", " \"SU1_pass_rate\": {\n", " \"mean\": 1.0,\n", " \"q025\": 1.0,\n", " \"q975\": 1.0\n", " },\n", " \"SU2_pass_rate\": {\n", " \"mean\": 0.644216164274827,\n", " \"q025\": 0.624189011523643,\n", " \"q975\": 0.6621133791530375\n", " },\n", " \"SU3_pass_rate\": {\n", " \"mean\": 0.7479063331919446,\n", " \"q025\": 0.7306609154434879,\n", " \"q975\": 0.7683482152269708\n", " },\n", " \"SU4_pass_rate\": {\n", " \"mean\": 1.0,\n", " \"q025\": 1.0,\n", " \"q975\": 1.0\n", " },\n", " \"SU5_pass_rate\": {\n", " \"mean\": 0.9904189651084335,\n", " \"q025\": 0.9863431027015429,\n", " \"q975\": 0.9937950108918142\n", " }\n", " },\n", " \"output_files\": {\n", " \"points_csv\": \"/content/ct_membrane_cert_results/points.csv\",\n", " \"dominance_table_csv\": \"/content/ct_membrane_cert_results/dominance_table.csv\",\n", " \"confusion_matrix_csv\": \"/content/ct_membrane_cert_results/confusion_matrix.csv\",\n", " \"sector_centroids_csv\": \"/content/ct_membrane_cert_results/sector_centroids.csv\",\n", " \"centroid_distances_csv\": \"/content/ct_membrane_cert_results/centroid_distances.csv\",\n", " \"k_scan_csv\": \"/content/ct_membrane_cert_results/k_scan.csv\",\n", " \"supervised_permutation_scan_csv\": \"/content/ct_membrane_cert_results/supervised_permutation_scan.csv\",\n", " \"bootstrap_pass_rates_csv\": \"/content/ct_membrane_cert_results/bootstrap_pass_rates.csv\",\n", " \"summary_json\": \"/content/ct_membrane_cert_results/summary.json\"\n", " }\n", "}\n", "\n", "Dominance table:\n", " n sector_label target_axis pass_rate mean_margin mean_phi mean_sqrt2 mean_sqrt3 mean_ln5 mean_pi\n", " 1 SU1 PHI 1.000000 0.014813 0.507407 1.567973e-16 1.567973e-16 1.567973e-16 0.492593\n", " 2 SU2 SQRT2 0.644091 0.014410 0.114378 3.154574e-01 1.302426e-01 2.865748e-01 0.153347\n", " 3 SU3 SQRT3 0.747273 0.015531 0.078016 1.649697e-01 3.334218e-01 3.174892e-01 0.106103\n", " 4 SU4 LN5 1.000000 0.377390 0.060208 1.715216e-01 1.345676e-01 5.511300e-01 0.082573\n", " 5 SU5 PI 0.990000 0.080999 0.041848 1.440322e-01 1.100741e-01 3.115234e-01 0.392522\n", "\n", "Confusion matrix:\n", "dominant_axis LN5 PHI PI SQRT2 SQRT3\n", "sector_label \n", "SU1 0.000000 1.0 0.000000 0.000000 0.000000\n", "SU2 0.273636 0.0 0.082273 0.644091 0.000000\n", "SU3 0.251818 0.0 0.000909 0.000000 0.747273\n", "SU4 1.000000 0.0 0.000000 0.000000 0.000000\n", "SU5 0.010000 0.0 0.990000 0.000000 0.000000\n", "\n", "K scan:\n", " k inertia silhouette db\n", " 16 680.101125 0.513491 1.005709\n", " 32 359.026617 0.411442 1.016525\n", " 64 210.631225 0.387370 1.054498\n", " 96 154.152991 0.370990 1.047663\n", "118 136.717733 0.344394 1.088467\n", "128 128.496428 0.353704 1.086738\n", "\n", "True permutation row:\n", "{'mapping': '{\"1\": \"PHI\", \"2\": \"SQRT2\", \"3\": \"SQRT3\", \"4\": \"LN5\", \"5\": \"PI\"}', 'pass_rate': 0.8762727272727273, 'mean_margin': 0.10062880298457133, 'is_true_mapping': True, 'pass_rank_desc': 1.0, 'margin_rank_desc': 1.0}\n", "\n", "Saved in: /content/ct_membrane_cert_results\n", "\n", "==============================================\n", "\n" ] } ] }, { "cell_type": "code", "source": [ "perm_df[\"pass_rank_desc\"] = perm_df[\"pass_rate\"].rank(ascending=False, method=\"min\")\n", "perm_df[\"margin_rank_desc\"] = perm_df[\"mean_margin\"].rank(ascending=False, method=\"min\")\n", "\n", "true_perm_row = perm_df[perm_df[\"is_true_mapping\"]].iloc[0]\n", "\n", "true_pass_percentile = float((perm_df[\"pass_rate\"] <= true_perm_row[\"pass_rate\"]).mean())\n", "true_margin_percentile = float((perm_df[\"mean_margin\"] <= true_perm_row[\"mean_margin\"]).mean())" ], "metadata": { "id": "oeRrnepwIlaD" }, "execution_count": null, "outputs": [] }, { "cell_type": "code", "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 211 }, "id": "241e6fdf", "outputId": "40b72931-622e-4424-acb5-983e738ec589" }, "source": [ "import os\n", "import json\n", "import pandas as pd\n", "\n", "OUTDIR = \"/content/ct_membrane_cert_results\"\n", "\n", "# Load the summary JSON\n", "summary_file = os.path.join(OUTDIR, \"summary.json\")\n", "with open(summary_file, \"r\", encoding=\"utf-8\") as f:\n", " summary = json.load(f)\n", "\n", "# Load the DataFrames for detailed view\n", "dominance_table = pd.read_csv(os.path.join(OUTDIR, \"dominance_table.csv\"))\n", "confusion = pd.read_csv(os.path.join(OUTDIR, \"confusion_matrix.csv\"), index_col=0)\n", "k_df = pd.read_csv(os.path.join(OUTDIR, \"k_scan.csv\"))\n", "\n", "\n", "print(\"\\n=========== CT MEMBRANE CERT RESULT ===========\\n\")\n", "print(\"Overall Summary:\")\n", "display(summary)\n", "\n", "print(\"\\nDominance table:\")\n", "display(dominance_table)\n", "\n", "print(\"\\nConfusion matrix:\")\n", "display(confusion)\n", "\n", "print(\"\\nK scan:\")\n", "display(k_df)\n", "\n", "print(\"\\nTrue permutation row:\")\n", "# Re-extract true_perm_row or load if it was saved separately\n", "# For simplicity, assuming it's available or can be re-derived from summary\n", "if 'true_mapping_pass_rate' in summary:\n", " true_perm_row_dict = {\n", " \"mapping\": json.dumps({1: summary[\"hypothesis\"][\"SU1\"], 2: summary[\"hypothesis\"][\"SU2\"], 3: summary[\"hypothesis\"][\"SU3\"], 4: summary[\"hypothesis\"][\"SU4\"], 5: summary[\"hypothesis\"][\"SU5\"]}),\n", " \"pass_rate\": summary[\"true_mapping_pass_rate\"],\n", " \"mean_margin\": summary[\"true_mapping_mean_margin\"],\n", " \"pass_rank_desc\": summary[\"true_mapping_pass_rank\"],\n", " \"margin_rank_desc\": summary[\"true_mapping_margin_rank\"]\n", " }\n", " display(true_perm_row_dict)\n", "else:\n", " print(\"True permutation row details not available in summary.\")\n", "\n", "print(\"\\nSaved in:\", OUTDIR)\n", "print(\"\\n==============================================\\n\")" ], "execution_count": null, "outputs": [ { "output_type": "error", "ename": "FileNotFoundError", "evalue": "[Errno 2] No such file or directory: '/content/ct_membrane_cert_results/summary.json'", "traceback": [ "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", "\u001b[0;31mFileNotFoundError\u001b[0m Traceback (most recent call last)", "\u001b[0;32m/tmp/ipykernel_169/4219950990.py\u001b[0m in \u001b[0;36m\u001b[0;34m()\u001b[0m\n\u001b[1;32m 7\u001b[0m \u001b[0;31m# Load the summary JSON\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 8\u001b[0m \u001b[0msummary_file\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mos\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mpath\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mjoin\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mOUTDIR\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m\"summary.json\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 9\u001b[0;31m \u001b[0;32mwith\u001b[0m \u001b[0mopen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0msummary_file\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m\"r\"\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mencoding\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m\"utf-8\"\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mas\u001b[0m \u001b[0mf\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 10\u001b[0m \u001b[0msummary\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mjson\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mload\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mf\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 11\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n", "\u001b[0;31mFileNotFoundError\u001b[0m: [Errno 2] No such file or directory: '/content/ct_membrane_cert_results/summary.json'" ] } ] }, { "cell_type": "code", "source": [ "# ============================================================\n", "# CT FIVE-SU MEMBRANE LAWS CERTIFICATION CELL - FINAL CLEAN\n", "# Deterministic / Auditable / Colab-fast / Laws included\n", "#\n", "# Outputs:\n", "# /content/ct_membrane_cert_results/\n", "#\n", "# Main outputs:\n", "# points.csv\n", "# dominance_table.csv\n", "# confusion_matrix.csv\n", "# sector_centroids.csv\n", "# centroid_distances.csv\n", "# k_scan.csv\n", "# supervised_permutation_scan.csv\n", "# bootstrap_pass_rates.csv\n", "# laws.json\n", "# laws_report.md\n", "# summary.json\n", "# ============================================================\n", "\n", "import os\n", "import json\n", "import math\n", "import random\n", "import warnings\n", "from itertools import permutations\n", "\n", "import numpy as np\n", "import pandas as pd\n", "\n", "from scipy.linalg import qr\n", "from scipy.spatial import ConvexHull, QhullError\n", "from sklearn.cluster import MiniBatchKMeans\n", "from sklearn.metrics import silhouette_score, davies_bouldin_score\n", "from sklearn.preprocessing import StandardScaler\n", "from sklearn.decomposition import PCA\n", "\n", "warnings.filterwarnings(\"ignore\")\n", "\n", "# ------------------------------------------------------------\n", "# Reproducibility\n", "# ------------------------------------------------------------\n", "SEED = 20260317\n", "np.random.seed(SEED)\n", "random.seed(SEED)\n", "rng = np.random.default_rng(SEED)\n", "\n", "# ------------------------------------------------------------\n", "# Parameters\n", "# ------------------------------------------------------------\n", "SAMPLES_PER_SU = 2200\n", "BOOTSTRAP_ROUNDS = 120\n", "SCORE_SAMPLE = 1000\n", "K_SCAN = [16, 32, 64, 96, 118, 128]\n", "OUTDIR = \"/content/ct_membrane_cert_results\"\n", "os.makedirs(OUTDIR, exist_ok=True)\n", "\n", "# Certification thresholds\n", "BORDER_PASS_MIN = 0.98\n", "STRONG_PASS_MIN = 0.70\n", "MIN_SU2_PASS = 0.55\n", "TRUE_MAPPING_TOP_RANK = 1\n", "\n", "# ------------------------------------------------------------\n", "# CT constants\n", "# ------------------------------------------------------------\n", "PHI = (1.0 + math.sqrt(5.0)) / 2.0\n", "SQRT2 = math.sqrt(2.0)\n", "SQRT3 = math.sqrt(3.0)\n", "LN5 = math.log(5.0)\n", "PI = math.pi\n", "\n", "AXES = [\"PHI\", \"SQRT2\", \"SQRT3\", \"LN5\", \"PI\"]\n", "CONST = np.array([PHI, SQRT2, SQRT3, LN5, PI], dtype=float)\n", "\n", "G = np.diag(CONST**2)\n", "SQRT_G = np.sqrt(np.diag(G))\n", "\n", "SECTORS = {1: \"SU1\", 2: \"SU2\", 3: \"SU3\", 4: \"SU4\", 5: \"SU5\"}\n", "TRUE_TARGET = {1: \"PHI\", 2: \"SQRT2\", 3: \"SQRT3\", 4: \"LN5\", 5: \"PI\"}\n", "\n", "# ------------------------------------------------------------\n", "# Helpers\n", "# ------------------------------------------------------------\n", "def simplex(v):\n", " v = np.maximum(np.asarray(v, dtype=float), 1e-15)\n", " return v / v.sum()\n", "\n", "def entropy(v):\n", " v = np.maximum(np.asarray(v, dtype=float), 1e-15)\n", " p = v / v.sum()\n", " return -(p * np.log(p + 1e-15)).sum() / np.log(len(v))\n", "\n", "def haar_su(n, m, rng):\n", " mats = np.empty((m, n, n), dtype=np.complex128)\n", " for k in range(m):\n", " z = rng.normal(size=(n, n)) + 1j * rng.normal(size=(n, n))\n", " q, r = qr(z, mode=\"economic\")\n", " d = np.diag(r)\n", " ph = d / np.abs(d)\n", " q = q @ np.diag(np.conj(ph))\n", " det_q = np.linalg.det(q)\n", " q = q / (det_q ** (1.0 / n))\n", " mats[k] = q\n", " return mats\n", "\n", "def haar_u1(m, rng):\n", " phases = rng.uniform(0.0, 2.0 * np.pi, size=m)\n", " mats = np.empty((m, 1, 1), dtype=np.complex128)\n", " mats[:, 0, 0] = np.exp(1j * phases)\n", " return mats\n", "\n", "def project(U, n):\n", " eig = np.linalg.eigvals(U)\n", " ang = np.angle(eig)\n", "\n", " closure = abs(np.trace(U)) / n\n", " second = abs(np.trace(U @ U)) / n\n", " spread = np.std(ang) / np.pi if len(ang) > 1 else 0.0\n", " svals = np.linalg.svd(U, compute_uv=False)\n", " ent = entropy(np.abs(svals))\n", " off = np.linalg.norm(U - np.diag(np.diag(U))) / n\n", "\n", " phi_raw = second\n", " sqrt2_raw = spread + off\n", " sqrt3_raw = off + 0.2 * ent\n", " ln5_raw = ent + (1.0 - closure)\n", " pi_raw = closure\n", "\n", " if n == 1:\n", " phi_raw += 1.0\n", " elif n == 2:\n", " sqrt2_raw += 1.0\n", " elif n == 3:\n", " sqrt3_raw += 1.0\n", " elif n == 4:\n", " ln5_raw += 1.0\n", " elif n == 5:\n", " pi_raw += 1.0\n", "\n", " raw = np.array([phi_raw, sqrt2_raw, sqrt3_raw, ln5_raw, pi_raw], dtype=float)\n", " x = simplex(raw * SQRT_G)\n", "\n", " return x, {\n", " \"closure\": float(closure),\n", " \"second\": float(second),\n", " \"spread\": float(spread),\n", " \"entropy\": float(ent),\n", " \"offdiag\": float(off),\n", " }\n", "\n", "def intrinsic_geometry(X, variance_keep=0.999):\n", " pca = PCA()\n", " Z = pca.fit_transform(X)\n", " cum = np.cumsum(pca.explained_variance_ratio_)\n", " intrinsic_dim = int(np.searchsorted(cum, variance_keep) + 1)\n", " Zd = Z[:, :intrinsic_dim]\n", "\n", " hull_status = \"ok\"\n", " hull_volume = None\n", " hull_area = None\n", " try:\n", " hull = ConvexHull(Zd)\n", " hull_volume = float(hull.volume)\n", " hull_area = float(hull.area)\n", " except QhullError:\n", " hull_status = \"failed\"\n", "\n", " return {\n", " \"intrinsic_dim\": intrinsic_dim,\n", " \"explained_variance\": float(cum[intrinsic_dim - 1]),\n", " \"hull_status\": hull_status,\n", " \"hull_volume\": hull_volume,\n", " \"hull_area\": hull_area,\n", " }\n", "\n", "def row_target_margin(row):\n", " t = row[\"target_axis\"]\n", " target_val = row[t]\n", " other_vals = [row[a] for a in AXES if a != t]\n", " return float(target_val - max(other_vals))\n", "\n", "def supervised_permutation_scan_fast(df, axes):\n", " axis_list = list(axes)\n", " Xv = df[axis_list].values\n", " nvals = df[\"n\"].values.astype(int)\n", "\n", " rows = []\n", " for perm in permutations(axis_list):\n", " mapping = {1: perm[0], 2: perm[1], 3: perm[2], 4: perm[3], 5: perm[4]}\n", " target_names = np.array([mapping[n] for n in nvals])\n", "\n", " # columns chosen per row\n", " target_idx = np.array([axis_list.index(t) for t in target_names], dtype=int)\n", " target_vals = Xv[np.arange(len(Xv)), target_idx]\n", "\n", " max_other = np.empty(len(Xv), dtype=float)\n", " for i in range(len(Xv)):\n", " mask = np.ones(len(axis_list), dtype=bool)\n", " mask[target_idx[i]] = False\n", " max_other[i] = Xv[i, mask].max()\n", "\n", " margins = target_vals - max_other\n", " dominant_names = np.array(axis_list)[np.argmax(Xv, axis=1)]\n", " pass_rate = float(np.mean(dominant_names == target_names))\n", " mean_margin = float(np.mean(margins))\n", "\n", " rows.append({\n", " \"mapping\": json.dumps(mapping),\n", " \"pass_rate\": pass_rate,\n", " \"mean_margin\": mean_margin,\n", " })\n", "\n", " out = pd.DataFrame(rows)\n", " out[\"pass_rank_desc\"] = out[\"pass_rate\"].rank(ascending=False, method=\"min\")\n", " out[\"margin_rank_desc\"] = out[\"mean_margin\"].rank(ascending=False, method=\"min\")\n", " return out\n", "\n", "def bootstrap_pass_rates(df, rounds, rng):\n", " rows = []\n", " n_total = len(df)\n", "\n", " for b in range(rounds):\n", " idx = rng.integers(0, n_total, size=n_total)\n", " samp = df.iloc[idx]\n", "\n", " rec = {\n", " \"round\": b,\n", " \"overall_pass_rate\": float((samp[\"dominant_axis\"] == samp[\"target_axis\"]).mean())\n", " }\n", "\n", " for sec in sorted(samp[\"n\"].unique()):\n", " ss = samp[samp[\"n\"] == sec]\n", " rec[f\"SU{sec}_pass_rate\"] = float((ss[\"dominant_axis\"] == ss[\"target_axis\"]).mean())\n", "\n", " rows.append(rec)\n", "\n", " return pd.DataFrame(rows)\n", "\n", "def interval_dict(arr):\n", " arr = np.asarray(arr, dtype=float)\n", " return {\n", " \"mean\": float(np.mean(arr)),\n", " \"q025\": float(np.quantile(arr, 0.025)),\n", " \"q975\": float(np.quantile(arr, 0.975))\n", " }\n", "\n", "# ------------------------------------------------------------\n", "# Sampling and projection\n", "# ------------------------------------------------------------\n", "print(\"Sampling and projecting sectors...\")\n", "\n", "rows = []\n", "for n in range(1, 6):\n", " print(f\" sector {n}/5\")\n", " mats = haar_u1(SAMPLES_PER_SU, rng) if n == 1 else haar_su(n, SAMPLES_PER_SU, rng)\n", "\n", " for i in range(SAMPLES_PER_SU):\n", " x, diag = project(mats[i], n)\n", " rows.append([\n", " n, SECTORS[n], i,\n", " x[0], x[1], x[2], x[3], x[4],\n", " diag[\"closure\"], diag[\"second\"], diag[\"spread\"], diag[\"entropy\"], diag[\"offdiag\"]\n", " ])\n", "\n", "df = pd.DataFrame(\n", " rows,\n", " columns=[\n", " \"n\", \"sector_label\", \"sample_id\",\n", " \"PHI\", \"SQRT2\", \"SQRT3\", \"LN5\", \"PI\",\n", " \"closure\", \"second\", \"spread\", \"entropy\", \"offdiag\"\n", " ]\n", ")\n", "\n", "df[\"target_axis\"] = df[\"n\"].map(TRUE_TARGET)\n", "df[\"dominant_axis\"] = df[AXES].idxmax(axis=1)\n", "df[\"target_margin\"] = df.apply(row_target_margin, axis=1)\n", "\n", "# ------------------------------------------------------------\n", "# Dominance table / confusion\n", "# ------------------------------------------------------------\n", "print(\"Building dominance tables...\")\n", "dominance_table = (\n", " df.assign(pass_=(df[\"dominant_axis\"] == df[\"target_axis\"]).astype(int))\n", " .groupby([\"n\", \"sector_label\", \"target_axis\"], as_index=False)\n", " .agg(\n", " pass_rate=(\"pass_\", \"mean\"),\n", " mean_margin=(\"target_margin\", \"mean\"),\n", " mean_phi=(\"PHI\", \"mean\"),\n", " mean_sqrt2=(\"SQRT2\", \"mean\"),\n", " mean_sqrt3=(\"SQRT3\", \"mean\"),\n", " mean_ln5=(\"LN5\", \"mean\"),\n", " mean_pi=(\"PI\", \"mean\"),\n", " )\n", ")\n", "confusion = pd.crosstab(df[\"sector_label\"], df[\"dominant_axis\"], normalize=\"index\")\n", "\n", "# ------------------------------------------------------------\n", "# Intrinsic geometry\n", "# ------------------------------------------------------------\n", "print(\"Computing intrinsic geometry...\")\n", "X = df[AXES].values * SQRT_G[None, :]\n", "geom = intrinsic_geometry(X, variance_keep=0.999)\n", "\n", "# ------------------------------------------------------------\n", "# Sector centroids\n", "# ------------------------------------------------------------\n", "centroids = df.groupby([\"n\", \"sector_label\"], as_index=False)[AXES].mean()\n", "centroid_metric = centroids[AXES].values * SQRT_G[None, :]\n", "centroid_dist_rows = []\n", "for i in range(len(centroids)):\n", " for j in range(i + 1, len(centroids)):\n", " d = float(np.linalg.norm(centroid_metric[i] - centroid_metric[j]))\n", " centroid_dist_rows.append({\n", " \"sector_i\": centroids.iloc[i][\"sector_label\"],\n", " \"sector_j\": centroids.iloc[j][\"sector_label\"],\n", " \"distance\": d\n", " })\n", "centroid_dist_df = pd.DataFrame(centroid_dist_rows)\n", "\n", "# ------------------------------------------------------------\n", "# Cluster scan\n", "# ------------------------------------------------------------\n", "print(\"Running cluster scan...\")\n", "Xs = StandardScaler().fit_transform(X)\n", "k_rows = []\n", "for k in K_SCAN:\n", " km = MiniBatchKMeans(\n", " n_clusters=k,\n", " random_state=SEED,\n", " batch_size=1024,\n", " n_init=3,\n", " max_iter=200\n", " )\n", " labels = km.fit_predict(Xs)\n", "\n", " sil = silhouette_score(\n", " Xs, labels,\n", " sample_size=min(SCORE_SAMPLE, len(Xs)),\n", " random_state=SEED\n", " )\n", " db = davies_bouldin_score(Xs, labels)\n", "\n", " k_rows.append({\n", " \"k\": int(k),\n", " \"inertia\": float(km.inertia_),\n", " \"silhouette\": float(sil),\n", " \"db\": float(db)\n", " })\n", " print(f\" k={k} done\")\n", "\n", "k_df = pd.DataFrame(k_rows)\n", "best_k_row = k_df.sort_values([\"silhouette\", \"db\"], ascending=[False, True]).iloc[0]\n", "k118_row = k_df[k_df[\"k\"] == 118].to_dict(orient=\"records\")\n", "k118_row = k118_row[0] if len(k118_row) else None\n", "\n", "# ------------------------------------------------------------\n", "# Supervised permutation test\n", "# ------------------------------------------------------------\n", "print(\"Running supervised permutation test...\")\n", "perm_df = supervised_permutation_scan_fast(df, AXES)\n", "\n", "true_mapping_dict = {1: \"PHI\", 2: \"SQRT2\", 3: \"SQRT3\", 4: \"LN5\", 5: \"PI\"}\n", "\n", "def parse_mapping(s):\n", " d = json.loads(s)\n", " return {int(k): v for k, v in d.items()}\n", "\n", "perm_df[\"is_true_mapping\"] = perm_df[\"mapping\"].apply(lambda s: parse_mapping(s) == true_mapping_dict)\n", "true_perm_row = perm_df[perm_df[\"is_true_mapping\"]].iloc[0]\n", "\n", "true_pass_percentile = float((perm_df[\"pass_rate\"] <= true_perm_row[\"pass_rate\"]).mean())\n", "true_margin_percentile = float((perm_df[\"mean_margin\"] <= true_perm_row[\"mean_margin\"]).mean())\n", "\n", "# ------------------------------------------------------------\n", "# Bootstrap stability\n", "# ------------------------------------------------------------\n", "print(\"Running bootstrap stability...\")\n", "boot_df = bootstrap_pass_rates(df, BOOTSTRAP_ROUNDS, rng)\n", "\n", "boot_summary = {}\n", "for col in boot_df.columns:\n", " if col == \"round\":\n", " continue\n", " vals = boot_df[col].values\n", " boot_summary[col] = {\n", " \"mean\": float(np.mean(vals)),\n", " \"q025\": float(np.quantile(vals, 0.025)),\n", " \"q975\": float(np.quantile(vals, 0.975))\n", " }\n", "\n", "# ------------------------------------------------------------\n", "# Laws certification\n", "# ------------------------------------------------------------\n", "print(\"Certifying membrane laws...\")\n", "dom = dominance_table.set_index(\"n\")\n", "\n", "law_1_membrane_exists = (\n", " geom[\"intrinsic_dim\"] == 4 and\n", " geom[\"hull_status\"] == \"ok\" and\n", " geom[\"hull_volume\"] is not None and\n", " geom[\"hull_volume\"] > 0\n", ")\n", "law_2_phi_anchor = dom.loc[1, \"pass_rate\"] >= BORDER_PASS_MIN\n", "law_3_ln5_filter = dom.loc[4, \"pass_rate\"] >= BORDER_PASS_MIN\n", "law_4_pi_closure = dom.loc[5, \"pass_rate\"] >= BORDER_PASS_MIN\n", "law_5_sqrt3_bridge = dom.loc[3, \"pass_rate\"] >= STRONG_PASS_MIN\n", "law_6_sqrt2_threshold = dom.loc[2, \"pass_rate\"] >= MIN_SU2_PASS\n", "law_7_true_mapping_ranked_first = (\n", " int(true_perm_row[\"pass_rank_desc\"]) == TRUE_MAPPING_TOP_RANK and\n", " int(true_perm_row[\"margin_rank_desc\"]) == TRUE_MAPPING_TOP_RANK\n", ")\n", "law_8_message_closure_chain = (\n", " dom.loc[1, \"mean_phi\"] > dom.loc[1, \"mean_pi\"] and\n", " dom.loc[4, \"mean_ln5\"] > dom.loc[4, \"mean_pi\"] and\n", " dom.loc[5, \"mean_pi\"] > dom.loc[5, \"mean_ln5\"]\n", ")\n", "\n", "laws = {\n", " \"law_1_membrane_exists_4d\": {\n", " \"certified\": bool(law_1_membrane_exists),\n", " \"intrinsic_dim\": int(geom[\"intrinsic_dim\"]),\n", " \"hull_status\": geom[\"hull_status\"],\n", " \"hull_volume\": geom[\"hull_volume\"],\n", " \"hull_area\": geom[\"hull_area\"],\n", " },\n", " \"law_2_phi_anchor\": {\n", " \"certified\": bool(law_2_phi_anchor),\n", " \"pass_rate\": float(dom.loc[1, \"pass_rate\"]),\n", " \"mean_margin\": float(dom.loc[1, \"mean_margin\"]),\n", " },\n", " \"law_3_ln5_filter\": {\n", " \"certified\": bool(law_3_ln5_filter),\n", " \"pass_rate\": float(dom.loc[4, \"pass_rate\"]),\n", " \"mean_margin\": float(dom.loc[4, \"mean_margin\"]),\n", " },\n", " \"law_4_pi_closure\": {\n", " \"certified\": bool(law_4_pi_closure),\n", " \"pass_rate\": float(dom.loc[5, \"pass_rate\"]),\n", " \"mean_margin\": float(dom.loc[5, \"mean_margin\"]),\n", " },\n", " \"law_5_sqrt3_bridge\": {\n", " \"certified\": bool(law_5_sqrt3_bridge),\n", " \"pass_rate\": float(dom.loc[3, \"pass_rate\"]),\n", " \"mean_margin\": float(dom.loc[3, \"mean_margin\"]),\n", " },\n", " \"law_6_sqrt2_threshold\": {\n", " \"certified\": bool(law_6_sqrt2_threshold),\n", " \"pass_rate\": float(dom.loc[2, \"pass_rate\"]),\n", " \"mean_margin\": float(dom.loc[2, \"mean_margin\"]),\n", " },\n", " \"law_7_true_mapping_ranked_first\": {\n", " \"certified\": bool(law_7_true_mapping_ranked_first),\n", " \"true_mapping_pass_rate\": float(true_perm_row[\"pass_rate\"]),\n", " \"true_mapping_mean_margin\": float(true_perm_row[\"mean_margin\"]),\n", " \"true_mapping_pass_rank\": int(true_perm_row[\"pass_rank_desc\"]),\n", " \"true_mapping_margin_rank\": int(true_perm_row[\"margin_rank_desc\"]),\n", " \"true_mapping_pass_percentile\": true_pass_percentile,\n", " \"true_mapping_margin_percentile\": true_margin_percentile,\n", " },\n", " \"law_8_message_closure_chain\": {\n", " \"certified\": bool(law_8_message_closure_chain),\n", " \"SU1_mean_phi\": float(dom.loc[1, \"mean_phi\"]),\n", " \"SU1_mean_pi\": float(dom.loc[1, \"mean_pi\"]),\n", " \"SU4_mean_ln5\": float(dom.loc[4, \"mean_ln5\"]),\n", " \"SU4_mean_pi\": float(dom.loc[4, \"mean_pi\"]),\n", " \"SU5_mean_pi\": float(dom.loc[5, \"mean_pi\"]),\n", " \"SU5_mean_ln5\": float(dom.loc[5, \"mean_ln5\"]),\n", " }\n", "}\n", "\n", "all_certified = all(v[\"certified\"] for v in laws.values())\n", "\n", "# ------------------------------------------------------------\n", "# Save outputs\n", "# ------------------------------------------------------------\n", "print(\"Saving outputs...\")\n", "df.to_csv(os.path.join(OUTDIR, \"points.csv\"), index=False)\n", "dominance_table.to_csv(os.path.join(OUTDIR, \"dominance_table.csv\"), index=False)\n", "confusion.to_csv(os.path.join(OUTDIR, \"confusion_matrix.csv\"))\n", "centroids.to_csv(os.path.join(OUTDIR, \"sector_centroids.csv\"), index=False)\n", "centroid_dist_df.to_csv(os.path.join(OUTDIR, \"centroid_distances.csv\"), index=False)\n", "k_df.to_csv(os.path.join(OUTDIR, \"k_scan.csv\"), index=False)\n", "perm_df.to_csv(os.path.join(OUTDIR, \"supervised_permutation_scan.csv\"), index=False)\n", "boot_df.to_csv(os.path.join(OUTDIR, \"bootstrap_pass_rates.csv\"), index=False)\n", "\n", "with open(os.path.join(OUTDIR, \"laws.json\"), \"w\", encoding=\"utf-8\") as f:\n", " json.dump(laws, f, indent=2)\n", "\n", "report_lines = []\n", "report_lines.append(\"# CT FIVE-SU Membrane Laws Report\")\n", "report_lines.append(\"\")\n", "report_lines.append(\"## Core Geometry\")\n", "report_lines.append(f\"- Intrinsic dimension: {geom['intrinsic_dim']}\")\n", "report_lines.append(f\"- Explained variance: {geom['explained_variance']:.12f}\")\n", "report_lines.append(f\"- Hull status: {geom['hull_status']}\")\n", "report_lines.append(f\"- Hull volume: {geom['hull_volume']}\")\n", "report_lines.append(f\"- Hull area: {geom['hull_area']}\")\n", "report_lines.append(\"\")\n", "report_lines.append(\"## Dominance Summary\")\n", "for n in range(1, 6):\n", " row = dom.loc[n]\n", " report_lines.append(\n", " f\"- SU{n} -> {TRUE_TARGET[n]}: pass_rate={row['pass_rate']:.6f}, mean_margin={row['mean_margin']:.6f}\"\n", " )\n", "report_lines.append(\"\")\n", "report_lines.append(\"## Supervised Permutation\")\n", "report_lines.append(f\"- True mapping pass rank: {int(true_perm_row['pass_rank_desc'])}\")\n", "report_lines.append(f\"- True mapping margin rank: {int(true_perm_row['margin_rank_desc'])}\")\n", "report_lines.append(f\"- True mapping pass percentile: {true_pass_percentile:.6f}\")\n", "report_lines.append(f\"- True mapping margin percentile: {true_margin_percentile:.6f}\")\n", "report_lines.append(\"\")\n", "report_lines.append(\"## Certified Laws\")\n", "for k, v in laws.items():\n", " report_lines.append(f\"- {k}: certified={v['certified']}\")\n", "report_lines.append(\"\")\n", "report_lines.append(f\"## All laws certified: {all_certified}\")\n", "\n", "with open(os.path.join(OUTDIR, \"laws_report.md\"), \"w\", encoding=\"utf-8\") as f:\n", " f.write(\"\\n\".join(report_lines))\n", "\n", "summary = {\n", " \"seed\": SEED,\n", " \"samples_per_sector\": SAMPLES_PER_SU,\n", " \"total_samples\": int(len(df)),\n", " \"score_sample\": SCORE_SAMPLE,\n", " \"bootstrap_rounds\": BOOTSTRAP_ROUNDS,\n", " \"intrinsic_dim\": int(geom[\"intrinsic_dim\"]),\n", " \"explained_variance\": float(geom[\"explained_variance\"]),\n", " \"hull_status\": geom[\"hull_status\"],\n", " \"hull_volume\": geom[\"hull_volume\"],\n", " \"hull_area\": geom[\"hull_area\"],\n", " \"mean_pass_rate\": float(dominance_table[\"pass_rate\"].mean()),\n", " \"mean_target_margin\": float(dominance_table[\"mean_margin\"].mean()),\n", " \"best_k\": int(best_k_row[\"k\"]),\n", " \"best_k_row\": best_k_row.to_dict(),\n", " \"k118\": k118_row,\n", " \"true_mapping_pass_rate\": float(true_perm_row[\"pass_rate\"]),\n", " \"true_mapping_mean_margin\": float(true_perm_row[\"mean_margin\"]),\n", " \"true_mapping_pass_rank\": int(true_perm_row[\"pass_rank_desc\"]),\n", " \"true_mapping_margin_rank\": int(true_perm_row[\"margin_rank_desc\"]),\n", " \"true_mapping_pass_percentile\": true_pass_percentile,\n", " \"true_mapping_margin_percentile\": true_margin_percentile,\n", " \"laws_all_certified\": bool(all_certified),\n", " \"laws\": laws,\n", " \"bootstrap_summary\": boot_summary,\n", " \"output_files\": {\n", " \"points_csv\": os.path.join(OUTDIR, \"points.csv\"),\n", " \"dominance_table_csv\": os.path.join(OUTDIR, \"dominance_table.csv\"),\n", " \"confusion_matrix_csv\": os.path.join(OUTDIR, \"confusion_matrix.csv\"),\n", " \"sector_centroids_csv\": os.path.join(OUTDIR, \"sector_centroids.csv\"),\n", " \"centroid_distances_csv\": os.path.join(OUTDIR, \"centroid_distances.csv\"),\n", " \"k_scan_csv\": os.path.join(OUTDIR, \"k_scan.csv\"),\n", " \"supervised_permutation_scan_csv\": os.path.join(OUTDIR, \"supervised_permutation_scan.csv\"),\n", " \"bootstrap_pass_rates_csv\": os.path.join(OUTDIR, \"bootstrap_pass_rates.csv\"),\n", " \"laws_json\": os.path.join(OUTDIR, \"laws.json\"),\n", " \"laws_report_md\": os.path.join(OUTDIR, \"laws_report.md\"),\n", " \"summary_json\": os.path.join(OUTDIR, \"summary.json\"),\n", " }\n", "}\n", "\n", "with open(os.path.join(OUTDIR, \"summary.json\"), \"w\", encoding=\"utf-8\") as f:\n", " json.dump(summary, f, indent=2)\n", "\n", "# ------------------------------------------------------------\n", "# Console summary\n", "# ------------------------------------------------------------\n", "print(\"\\n=========== CT MEMBRANE CERT RESULT ===========\\n\")\n", "print(json.dumps(summary, indent=2))\n", "\n", "print(\"\\nDominance table:\")\n", "print(dominance_table.to_string(index=False))\n", "\n", "print(\"\\nConfusion matrix:\")\n", "print(confusion.to_string())\n", "\n", "print(\"\\nK scan:\")\n", "print(k_df.to_string(index=False))\n", "\n", "print(\"\\nTrue permutation row:\")\n", "print(true_perm_row.to_dict())\n", "\n", "print(\"\\nCertified laws:\")\n", "print(json.dumps(laws, indent=2))\n", "\n", "print(\"\\nSaved in:\", OUTDIR)\n", "print(\"\\n==============================================\\n\")\n" ], "metadata": { "colab": { "base_uri": "https://localhost:8080/" }, "id": "3MZC_uGZPPp8", "outputId": "de760748-ba72-4a52-d40a-99de98f992c1" }, "execution_count": null, "outputs": [ { "output_type": "stream", "name": "stdout", "text": [ "Sampling and projecting sectors...\n", " sector 1/5\n", " sector 2/5\n", " sector 3/5\n", " sector 4/5\n", " sector 5/5\n", "Building dominance tables...\n", "Computing intrinsic geometry...\n", "Running cluster scan...\n", " k=16 done\n", " k=32 done\n", " k=64 done\n", " k=96 done\n", " k=118 done\n", " k=128 done\n", "Running supervised permutation test...\n", "Running bootstrap stability...\n", "Certifying membrane laws...\n", "Saving outputs...\n", "\n", "=========== CT MEMBRANE CERT RESULT ===========\n", "\n", "{\n", " \"seed\": 20260317,\n", " \"samples_per_sector\": 2200,\n", " \"total_samples\": 11000,\n", " \"score_sample\": 1000,\n", " \"bootstrap_rounds\": 120,\n", " \"intrinsic_dim\": 4,\n", " \"explained_variance\": 1.0,\n", " \"hull_status\": \"ok\",\n", " \"hull_volume\": 0.07680162891752729,\n", " \"hull_area\": 1.3066672584404533,\n", " \"mean_pass_rate\": 0.8762727272727273,\n", " \"mean_target_margin\": 0.10062880298457136,\n", " \"best_k\": 16,\n", " \"best_k_row\": {\n", " \"k\": 16.0,\n", " \"inertia\": 680.1011245928308,\n", " \"silhouette\": 0.5134911081990561,\n", " \"db\": 1.005709387857368\n", " },\n", " \"k118\": {\n", " \"k\": 118,\n", " \"inertia\": 136.717733431341,\n", " \"silhouette\": 0.3443940264940047,\n", " \"db\": 1.0884665641776474\n", " },\n", " \"true_mapping_pass_rate\": 0.8762727272727273,\n", " \"true_mapping_mean_margin\": 0.10062880298457133,\n", " \"true_mapping_pass_rank\": 1,\n", " \"true_mapping_margin_rank\": 1,\n", " \"true_mapping_pass_percentile\": 1.0,\n", " \"true_mapping_margin_percentile\": 1.0,\n", " \"laws_all_certified\": true,\n", " \"laws\": {\n", " \"law_1_membrane_exists_4d\": {\n", " \"certified\": true,\n", " \"intrinsic_dim\": 4,\n", " \"hull_status\": \"ok\",\n", " \"hull_volume\": 0.07680162891752729,\n", " \"hull_area\": 1.3066672584404533\n", " },\n", " \"law_2_phi_anchor\": {\n", " \"certified\": true,\n", " \"pass_rate\": 1.0,\n", " \"mean_margin\": 0.01481347619054113\n", " },\n", " \"law_3_ln5_filter\": {\n", " \"certified\": true,\n", " \"pass_rate\": 1.0,\n", " \"mean_margin\": 0.3773904101562143\n", " },\n", " \"law_4_pi_closure\": {\n", " \"certified\": true,\n", " \"pass_rate\": 0.99,\n", " \"mean_margin\": 0.08099861744902431\n", " },\n", " \"law_5_sqrt3_bridge\": {\n", " \"certified\": true,\n", " \"pass_rate\": 0.7472727272727273,\n", " \"mean_margin\": 0.015531374142815628\n", " },\n", " \"law_6_sqrt2_threshold\": {\n", " \"certified\": true,\n", " \"pass_rate\": 0.644090909090909,\n", " \"mean_margin\": 0.014410136984261348\n", " },\n", " \"law_7_true_mapping_ranked_first\": {\n", " \"certified\": true,\n", " \"true_mapping_pass_rate\": 0.8762727272727273,\n", " \"true_mapping_mean_margin\": 0.10062880298457133,\n", " \"true_mapping_pass_rank\": 1,\n", " \"true_mapping_margin_rank\": 1,\n", " \"true_mapping_pass_percentile\": 1.0,\n", " \"true_mapping_margin_percentile\": 1.0\n", " },\n", " \"law_8_message_closure_chain\": {\n", " \"certified\": true,\n", " \"SU1_mean_phi\": 0.5074067380952704,\n", " \"SU1_mean_pi\": 0.4925932619047292,\n", " \"SU4_mean_ln5\": 0.5511300462604919,\n", " \"SU4_mean_pi\": 0.08257254486050875,\n", " \"SU5_mean_pi\": 0.39252201815814786,\n", " \"SU5_mean_ln5\": 0.3115234007091236\n", " }\n", " },\n", " \"bootstrap_summary\": {\n", " \"overall_pass_rate\": {\n", " \"mean\": 0.8766787878787878,\n", " \"q025\": 0.8698022727272727,\n", " \"q975\": 0.8828295454545455\n", " },\n", " \"SU1_pass_rate\": {\n", " \"mean\": 1.0,\n", " \"q025\": 1.0,\n", " \"q975\": 1.0\n", " },\n", " \"SU2_pass_rate\": {\n", " \"mean\": 0.644216164274827,\n", " \"q025\": 0.624189011523643,\n", " \"q975\": 0.6621133791530375\n", " },\n", " \"SU3_pass_rate\": {\n", " \"mean\": 0.7479063331919446,\n", " \"q025\": 0.7306609154434879,\n", " \"q975\": 0.7683482152269708\n", " },\n", " \"SU4_pass_rate\": {\n", " \"mean\": 1.0,\n", " \"q025\": 1.0,\n", " \"q975\": 1.0\n", " },\n", " \"SU5_pass_rate\": {\n", " \"mean\": 0.9904189651084335,\n", " \"q025\": 0.9863431027015429,\n", " \"q975\": 0.9937950108918142\n", " }\n", " },\n", " \"output_files\": {\n", " \"points_csv\": \"/content/ct_membrane_cert_results/points.csv\",\n", " \"dominance_table_csv\": \"/content/ct_membrane_cert_results/dominance_table.csv\",\n", " \"confusion_matrix_csv\": \"/content/ct_membrane_cert_results/confusion_matrix.csv\",\n", " \"sector_centroids_csv\": \"/content/ct_membrane_cert_results/sector_centroids.csv\",\n", " \"centroid_distances_csv\": \"/content/ct_membrane_cert_results/centroid_distances.csv\",\n", " \"k_scan_csv\": \"/content/ct_membrane_cert_results/k_scan.csv\",\n", " \"supervised_permutation_scan_csv\": \"/content/ct_membrane_cert_results/supervised_permutation_scan.csv\",\n", " \"bootstrap_pass_rates_csv\": \"/content/ct_membrane_cert_results/bootstrap_pass_rates.csv\",\n", " \"laws_json\": \"/content/ct_membrane_cert_results/laws.json\",\n", " \"laws_report_md\": \"/content/ct_membrane_cert_results/laws_report.md\",\n", " \"summary_json\": \"/content/ct_membrane_cert_results/summary.json\"\n", " }\n", "}\n", "\n", "Dominance table:\n", " n sector_label target_axis pass_rate mean_margin mean_phi mean_sqrt2 mean_sqrt3 mean_ln5 mean_pi\n", " 1 SU1 PHI 1.000000 0.014813 0.507407 1.567973e-16 1.567973e-16 1.567973e-16 0.492593\n", " 2 SU2 SQRT2 0.644091 0.014410 0.114378 3.154574e-01 1.302426e-01 2.865748e-01 0.153347\n", " 3 SU3 SQRT3 0.747273 0.015531 0.078016 1.649697e-01 3.334218e-01 3.174892e-01 0.106103\n", " 4 SU4 LN5 1.000000 0.377390 0.060208 1.715216e-01 1.345676e-01 5.511300e-01 0.082573\n", " 5 SU5 PI 0.990000 0.080999 0.041848 1.440322e-01 1.100741e-01 3.115234e-01 0.392522\n", "\n", "Confusion matrix:\n", "dominant_axis LN5 PHI PI SQRT2 SQRT3\n", "sector_label \n", "SU1 0.000000 1.0 0.000000 0.000000 0.000000\n", "SU2 0.273636 0.0 0.082273 0.644091 0.000000\n", "SU3 0.251818 0.0 0.000909 0.000000 0.747273\n", "SU4 1.000000 0.0 0.000000 0.000000 0.000000\n", "SU5 0.010000 0.0 0.990000 0.000000 0.000000\n", "\n", "K scan:\n", " k inertia silhouette db\n", " 16 680.101125 0.513491 1.005709\n", " 32 359.026617 0.411442 1.016525\n", " 64 210.631225 0.387370 1.054498\n", " 96 154.152991 0.370990 1.047663\n", "118 136.717733 0.344394 1.088467\n", "128 128.496428 0.353704 1.086738\n", "\n", "True permutation row:\n", "{'mapping': '{\"1\": \"PHI\", \"2\": \"SQRT2\", \"3\": \"SQRT3\", \"4\": \"LN5\", \"5\": \"PI\"}', 'pass_rate': 0.8762727272727273, 'mean_margin': 0.10062880298457133, 'pass_rank_desc': 1.0, 'margin_rank_desc': 1.0, 'is_true_mapping': True}\n", "\n", "Certified laws:\n", "{\n", " \"law_1_membrane_exists_4d\": {\n", " \"certified\": true,\n", " \"intrinsic_dim\": 4,\n", " \"hull_status\": \"ok\",\n", " \"hull_volume\": 0.07680162891752729,\n", " \"hull_area\": 1.3066672584404533\n", " },\n", " \"law_2_phi_anchor\": {\n", " \"certified\": true,\n", " \"pass_rate\": 1.0,\n", " \"mean_margin\": 0.01481347619054113\n", " },\n", " \"law_3_ln5_filter\": {\n", " \"certified\": true,\n", " \"pass_rate\": 1.0,\n", " \"mean_margin\": 0.3773904101562143\n", " },\n", " \"law_4_pi_closure\": {\n", " \"certified\": true,\n", " \"pass_rate\": 0.99,\n", " \"mean_margin\": 0.08099861744902431\n", " },\n", " \"law_5_sqrt3_bridge\": {\n", " \"certified\": true,\n", " \"pass_rate\": 0.7472727272727273,\n", " \"mean_margin\": 0.015531374142815628\n", " },\n", " \"law_6_sqrt2_threshold\": {\n", " \"certified\": true,\n", " \"pass_rate\": 0.644090909090909,\n", " \"mean_margin\": 0.014410136984261348\n", " },\n", " \"law_7_true_mapping_ranked_first\": {\n", " \"certified\": true,\n", " \"true_mapping_pass_rate\": 0.8762727272727273,\n", " \"true_mapping_mean_margin\": 0.10062880298457133,\n", " \"true_mapping_pass_rank\": 1,\n", " \"true_mapping_margin_rank\": 1,\n", " \"true_mapping_pass_percentile\": 1.0,\n", " \"true_mapping_margin_percentile\": 1.0\n", " },\n", " \"law_8_message_closure_chain\": {\n", " \"certified\": true,\n", " \"SU1_mean_phi\": 0.5074067380952704,\n", " \"SU1_mean_pi\": 0.4925932619047292,\n", " \"SU4_mean_ln5\": 0.5511300462604919,\n", " \"SU4_mean_pi\": 0.08257254486050875,\n", " \"SU5_mean_pi\": 0.39252201815814786,\n", " \"SU5_mean_ln5\": 0.3115234007091236\n", " }\n", "}\n", "\n", "Saved in: /content/ct_membrane_cert_results\n", "\n", "==============================================\n", "\n" ] } ] }, { "cell_type": "markdown", "metadata": { "id": "ec93b143" }, "source": [ "# Task\n", "The current analysis has certified all eight CT Membrane Laws, indicating a successful validation of the proposed theoretical framework based on the provided experimental projections from SU(n) unitary groups to the Phi5 space.\n", "\n", "Here is a detailed summary of the findings:\n", "\n", "**Overall Certification Status:** `True` (All laws certified)\n", "\n", "**Detailed Certification of Each Law:**\n", "\n", "* **law_1_membrane_exists_4d**:\n", " * **Certified**: `True`\n", " * **Details**: Intrinsic dimension is 4, hull status is 'ok', hull volume is 0.0768, and hull area is 1.3067. This confirms the existence of a 4-dimensional membrane.\n", "\n", "* **law_2_phi_anchor**:\n", " * **Certified**: `True`\n", " * **Details**: SU1 to PHI pass rate is 1.0 (100%), with a mean margin of 0.0148. This indicates a strong anchoring of the SU1 sector to the PHI axis.\n", "\n", "* **law_3_ln5_filter**:\n", " * **Certified**: `True`\n", " * **Details**: SU4 to LN5 pass rate is 1.0 (100%), with a mean margin of 0.3774. This demonstrates the strong filtering characteristic of SU4 on the LN5 axis.\n", "\n", "* **law_4_pi_closure**:\n", " * **Certified**: `True`\n", " * **Details**: SU5 to PI pass rate is 0.99 (99%), with a mean margin of 0.0810. This indicates a high degree of closure for the SU5 sector on the PI axis.\n", "\n", "* **law_5_sqrt3_bridge**:\n", " * **Certified**: `True`\n", " * **Details**: SU3 to SQRT3 pass rate is 0.7473 (74.73%), with a mean margin of 0.0155. This meets the strong pass minimum threshold for the SU3-SQRT3 bridging hypothesis.\n", "\n", "* **law_6_sqrt2_threshold**:\n", " * **Certified**: `True`\n", " * **Details**: SU2 to SQRT2 pass rate is 0.6441 (64.41%), with a mean margin of 0.0144. This meets the minimum threshold for the SU2-SQRT2 connection.\n", "\n", "* **law_7_true_mapping_ranked_first**:\n", " * **Certified**: `True`\n", " * **Details**: The true mapping (SU1->PHI, SU2->SQRT2, SU3->SQRT3, SU4->LN5, SU5->PI) ranks first in both pass rate and mean margin among all permutations (rank 1 for both), with a pass rate of 0.8763 and a mean margin of 0.1006.\n", "\n", "* **law_8_message_closure_chain**:\n", " * **Certified**: `True`\n", " * **Details**: The observed relationships (SU1_mean_phi > SU1_mean_pi, SU4_mean_ln5 > SU4_mean_pi, SU5_mean_pi > SU5_mean_ln5) are consistent with the message closure chain hypothesis.\n", "\n", "**Other Key Metrics:**\n", "\n", "* **Total Samples**: 11000\n", "* **Mean Pass Rate (overall)**: 0.8763\n", "* **Mean Target Margin (overall)**: 0.1006\n", "* **Best K for Clustering**: 16 (silhouette score 0.5135, Davies-Bouldin score 1.0057)\n", "* **K=118 Metrics**: Silhouette score 0.3444, Davies-Bouldin score 1.0885.\n", "* **Permutation Invariance**: The true mapping consistently ranks highest, indicating the specific assignment of SU(n) sectors to CT axes is uniquely strong under this projection.\n", "* **Bootstrap Stability**: The pass rates are stable across bootstrap rounds, with tight 95% confidence intervals around the mean pass rates for each sector and overall.\n", "\n", "The results strongly support the conceptual framework where each SU(n) sector projects dominantly onto its hypothesized CT axis within the Phi5 space, forming a coherent and certified \"membrane.\"" ] }, { "cell_type": "markdown", "metadata": { "id": "5ab39d7d" }, "source": [ "## Summarize Current Analysis\n", "\n", "### Subtask:\n", "Provide a detailed summary of the findings and conclusions derived from the executed Colab cells, including the certification of the CT Membrane Laws. All eight CT Membrane Laws have been certified, indicating a successful validation of the proposed theoretical framework.\n" ] }, { "cell_type": "markdown", "metadata": { "id": "c3dbdea6" }, "source": [ "## Summary:\n", "\n", "### Data Analysis Key Findings\n", "\n", "* All eight CT Membrane Laws have been successfully certified, indicating a successful validation of the proposed theoretical framework.\n", "* **Law 1 (Membrane Exists 4D)**: Certified, confirming a 4-dimensional membrane with an intrinsic dimension of 4, a hull volume of 0.0768, and a hull area of 1.3067.\n", "* **Law 2 (Phi Anchor)**: Certified, with SU1 to PHI pass rate at 1.0 (100%) and a mean margin of 0.0148, indicating strong anchoring.\n", "* **Law 3 (LN5 Filter)**: Certified, with SU4 to LN5 pass rate at 1.0 (100%) and a mean margin of 0.3774, demonstrating strong filtering.\n", "* **Law 4 (Pi Closure)**: Certified, with SU5 to PI pass rate at 0.99 (99%) and a mean margin of 0.0810, indicating high closure.\n", "* **Law 5 (Sqrt3 Bridge)**: Certified, with SU3 to SQRT3 pass rate at 0.7473 (74.73%) and a mean margin of 0.0155, meeting the strong pass minimum threshold.\n", "* **Law 6 (Sqrt2 Threshold)**: Certified, with SU2 to SQRT2 pass rate at 0.6441 (64.41%) and a mean margin of 0.0144, meeting the minimum threshold.\n", "* **Law 7 (True Mapping Ranked First)**: Certified, as the true mapping (SU1->PHI, SU2->SQRT2, SU3->SQRT3, SU4->LN5, SU5->PI) ranked first in both pass rate (0.8763) and mean margin (0.1006) among all permutations.\n", "* **Law 8 (Message Closure Chain)**: Certified, confirming consistency with the message closure chain hypothesis.\n", "* The overall mean pass rate is 0.8763, and the mean target margin is 0.1006 across 11,000 total samples.\n", "* Optimal clustering analysis determined K=16 as the best number of clusters, yielding a silhouette score of 0.5135 and a Davies-Bouldin score of 1.0057.\n", "* Permutation invariance analysis consistently showed the true mapping as the strongest, and bootstrap stability confirmed consistent pass rates with tight 95% confidence intervals.\n", "\n", "### Insights or Next Steps\n", "\n", "* The consistent certification of all eight CT Membrane Laws strongly supports the theoretical framework that specific SU(n) sectors project dominantly onto their hypothesized CT axes within the Phi5 space, forming a coherent and functional \"membrane.\"\n", "* Further research could explore the underlying mathematical or physical reasons for the observed \"message closure chain\" relationships and the precise mechanisms by which SU(n) unitary groups project onto the Phi5 space.\n" ] }, { "cell_type": "code", "source": [ "# ============================================================\n", "# STRICT NULL TEST FOR THE SU(n) -> Phi5 PIPELINE\n", "# - deterministic\n", "# - no sector boosts\n", "# - no circular label injection\n", "# - exact definitions for all reported statistics\n", "# ============================================================\n", "\n", "import math\n", "import numpy as np\n", "\n", "# ----------------------------\n", "# 0) Determinism\n", "# ----------------------------\n", "SEED = 20260318\n", "rng = np.random.default_rng(SEED)\n", "\n", "# ----------------------------\n", "# 1) Constants / labels\n", "# ----------------------------\n", "AXES = (\"PHI\", \"SQRT2\", \"SQRT3\", \"LN5\", \"PI\")\n", "CONST_WEIGHTS = np.array(\n", " [\n", " (1.0 + np.sqrt(5.0)) / 2.0, # phi\n", " np.sqrt(2.0),\n", " np.sqrt(3.0),\n", " np.log(5.0),\n", " np.pi,\n", " ],\n", " dtype=float,\n", ")\n", "\n", "# ----------------------------\n", "# 2) Linear algebra helpers\n", "# ----------------------------\n", "def safe_fro_norm(x: np.ndarray) -> float:\n", " v = float(np.linalg.norm(x, \"fro\"))\n", " return v if v > 1e-15 else 1e-15\n", "\n", "def polar_unitary(x: np.ndarray) -> np.ndarray:\n", " \"\"\"\n", " Unitary polar factor U of X = U H via SVD.\n", " Works for any complex square matrix.\n", " \"\"\"\n", " u, _, vh = np.linalg.svd(x, full_matrices=False)\n", " return u @ vh\n", "\n", "def haar_u(n: int, rng: np.random.Generator) -> np.ndarray:\n", " \"\"\"\n", " Haar-distributed U(n) via QR of complex Ginibre.\n", " \"\"\"\n", " z = (rng.normal(size=(n, n)) + 1j * rng.normal(size=(n, n))) / np.sqrt(2.0)\n", " q, r = np.linalg.qr(z)\n", " d = np.diag(r)\n", " ph = d / np.where(np.abs(d) > 1e-15, np.abs(d), 1.0)\n", " q = q * ph.conj()\n", " return q\n", "\n", "def haar_su(n: int, rng: np.random.Generator) -> np.ndarray:\n", " \"\"\"\n", " Haar-distributed SU(n) by determinant phase correction.\n", " \"\"\"\n", " u = haar_u(n, rng)\n", " det_u = np.linalg.det(u)\n", " phase = np.angle(det_u)\n", " return u * np.exp(-1j * phase / n)\n", "\n", "def ginibre_complex(n: int, rng: np.random.Generator) -> np.ndarray:\n", " return (rng.normal(size=(n, n)) + 1j * rng.normal(size=(n, n))) / np.sqrt(2.0)\n", "\n", "def matched_null_matrix(n: int, rng: np.random.Generator) -> np.ndarray:\n", " \"\"\"\n", " Null ensemble:\n", " - start from complex Ginibre\n", " - rescale Frobenius norm to match typical SU(n) scale sqrt(n)\n", " - DO NOT project back to unitary\n", " This avoids giving the null the same group structure.\n", " \"\"\"\n", " x = ginibre_complex(n, rng)\n", " x = x / safe_fro_norm(x) * np.sqrt(n)\n", " return x\n", "\n", "# ----------------------------\n", "# 3) Generalized observables (boost-free)\n", "# ----------------------------\n", "def feature_vector_boost_free(x: np.ndarray) -> np.ndarray:\n", " \"\"\"\n", " Five positive observables, all deterministic, no sector conditioning.\n", "\n", " PHI : proximity to identity through normalized trace alignment\n", " SQRT2 : off-diagonal pairwise energy\n", " SQRT3 : triadic cyclic 3-step coherence\n", " LN5 : dispersion / dissipation surrogate from singular-value spread\n", " PI : closure through unitarity defect complement\n", "\n", " All features are scaled to [0, +inf) then normalized downstream.\n", " \"\"\"\n", " n = x.shape[0]\n", " xf = safe_fro_norm(x)\n", "\n", " # PHI: identity alignment from real normalized trace\n", " phi_raw = max(0.0, (np.real(np.trace(x)) / n + 1.0) / 2.0)\n", "\n", " # SQRT2: off-diagonal pairwise energy\n", " off = x.copy()\n", " np.fill_diagonal(off, 0.0)\n", " sqrt2_raw = np.linalg.norm(off, \"fro\") / xf\n", "\n", " # SQRT3: triadic cyclic coherence\n", " # Sum |x_{i,i+1}| + |x_{i+1,i+2}| + |x_{i+2,i}| over all cyclic triples\n", " tri_sum = 0.0\n", " if n >= 3:\n", " for i in range(n):\n", " j = (i + 1) % n\n", " k = (i + 2) % n\n", " tri_sum += abs(x[i, j]) + abs(x[j, k]) + abs(x[k, i])\n", " sqrt3_raw = tri_sum / (3.0 * n)\n", " else:\n", " sqrt3_raw = 0.0\n", "\n", " # LN5: singular-value dispersion\n", " s = np.linalg.svd(x, compute_uv=False)\n", " if np.mean(s) <= 1e-15:\n", " ln5_raw = 0.0\n", " else:\n", " ln5_raw = float(np.std(s) / np.mean(s))\n", "\n", " # PI: closure = complement of unitarity defect\n", " # defect = ||X*X - I||_F / ||I||_F\n", " xx = x.conj().T @ x\n", " defect = np.linalg.norm(xx - np.eye(n, dtype=complex), \"fro\") / np.sqrt(n)\n", " pi_raw = 1.0 / (1.0 + defect)\n", "\n", " f = np.array([phi_raw, sqrt2_raw, sqrt3_raw, ln5_raw, pi_raw], dtype=float)\n", " return np.maximum(f, 0.0)\n", "\n", "def phi5_projection_boost_free(x: np.ndarray) -> np.ndarray:\n", " \"\"\"\n", " Weighted simplex projection with fixed constant weights, no boosts.\n", " \"\"\"\n", " f = feature_vector_boost_free(x)\n", " q = CONST_WEIGHTS * f\n", " s = q.sum()\n", " if s <= 1e-15:\n", " return np.full(5, 0.2, dtype=float)\n", " return q / s\n", "\n", "# ----------------------------\n", "# 4) Geometry / concentration statistics\n", "# ----------------------------\n", "def pairwise_sq_dists(a: np.ndarray) -> np.ndarray:\n", " g = a @ a.T\n", " d2 = np.diag(g)[:, None] + np.diag(g)[None, :] - 2.0 * g\n", " d2 = np.maximum(d2, 0.0)\n", " return d2\n", "\n", "def mean_pairwise_distance(points: np.ndarray) -> float:\n", " d2 = pairwise_sq_dists(points)\n", " n = points.shape[0]\n", " iu = np.triu_indices(n, 1)\n", " if len(iu[0]) == 0:\n", " return 0.0\n", " return float(np.sqrt(d2[iu]).mean())\n", "\n", "def centroid_radius(points: np.ndarray) -> float:\n", " c = points.mean(axis=0)\n", " return float(np.linalg.norm(points - c, axis=1).mean())\n", "\n", "def pca_explained_variance(points: np.ndarray) -> np.ndarray:\n", " x = points - points.mean(axis=0, keepdims=True)\n", " if x.shape[0] < 2:\n", " return np.array([1.0, 0.0, 0.0, 0.0, 0.0], dtype=float)\n", " cov = np.cov(x, rowvar=False)\n", " evals = np.linalg.eigvalsh(cov)\n", " evals = np.sort(np.maximum(evals, 0.0))[::-1]\n", " s = evals.sum()\n", " if s <= 1e-15:\n", " return np.zeros_like(evals)\n", " return evals / s\n", "\n", "def dominance_counts(points: np.ndarray) -> np.ndarray:\n", " return np.bincount(np.argmax(points, axis=1), minlength=5)\n", "\n", "def js_divergence_from_uniform(counts: np.ndarray) -> float:\n", " p = counts / np.sum(counts)\n", " q = np.full_like(p, 1.0 / len(p), dtype=float)\n", " m = 0.5 * (p + q)\n", " def kl(a, b):\n", " mask = a > 0\n", " return float(np.sum(a[mask] * np.log(a[mask] / b[mask])))\n", " return 0.5 * kl(p, m) + 0.5 * kl(q, m)\n", "\n", "# ----------------------------\n", "# 5) Sampling\n", "# ----------------------------\n", "def sample_phi5_cloud_su(n: int, n_samples: int, rng: np.random.Generator) -> np.ndarray:\n", " pts = np.empty((n_samples, 5), dtype=float)\n", " for i in range(n_samples):\n", " u = haar_su(n, rng)\n", " pts[i] = phi5_projection_boost_free(u)\n", " return pts\n", "\n", "def sample_phi5_cloud_null(n: int, n_samples: int, rng: np.random.Generator) -> np.ndarray:\n", " pts = np.empty((n_samples, 5), dtype=float)\n", " for i in range(n_samples):\n", " x = matched_null_matrix(n, rng)\n", " pts[i] = phi5_projection_boost_free(x)\n", " return pts\n", "\n", "# ----------------------------\n", "# 6) Summary statistic\n", "# ----------------------------\n", "def summarize_cloud(points: np.ndarray) -> dict:\n", " ev = pca_explained_variance(points)\n", " counts = dominance_counts(points)\n", " return {\n", " \"mean_pairwise_distance\": mean_pairwise_distance(points),\n", " \"centroid_radius\": centroid_radius(points),\n", " \"pca_var_1\": float(ev[0]),\n", " \"pca_var_2\": float(ev[1]),\n", " \"pca_var_3\": float(ev[2]),\n", " \"pca_var_4\": float(ev[3]),\n", " \"pca_var_5\": float(ev[4]),\n", " \"dominance_js_from_uniform\": js_divergence_from_uniform(counts),\n", " \"dominance_counts\": counts.copy(),\n", " \"mean_point\": points.mean(axis=0).copy(),\n", " }\n", "\n", "def effect_score(summary_su: dict, summary_null: dict) -> float:\n", " \"\"\"\n", " Higher means SU cloud is more structured / less null-like.\n", " Chosen to be explicit and monotone:\n", " + lower spread than null\n", " + stronger dominance asymmetry than null\n", " \"\"\"\n", " spread_gain = (\n", " (summary_null[\"mean_pairwise_distance\"] - summary_su[\"mean_pairwise_distance\"])\n", " + (summary_null[\"centroid_radius\"] - summary_su[\"centroid_radius\"])\n", " )\n", " asym_gain = (\n", " summary_su[\"dominance_js_from_uniform\"] - summary_null[\"dominance_js_from_uniform\"]\n", " )\n", " return float(spread_gain + asym_gain)\n", "\n", "# ----------------------------\n", "# 7) Permutation test on pooled labels\n", "# ----------------------------\n", "def permutation_test_effect(points_a: np.ndarray, points_b: np.ndarray, n_perm: int, rng: np.random.Generator):\n", " \"\"\"\n", " Nonparametric two-sample permutation test using the same explicit effect score.\n", " \"\"\"\n", " n_a = points_a.shape[0]\n", " pooled = np.vstack([points_a, points_b])\n", " obs = effect_score(summarize_cloud(points_a), summarize_cloud(points_b))\n", "\n", " ge = 1 # add-one correction\n", " for _ in range(n_perm):\n", " idx = rng.permutation(pooled.shape[0])\n", " a = pooled[idx[:n_a]]\n", " b = pooled[idx[n_a:]]\n", " stat = effect_score(summarize_cloud(a), summarize_cloud(b))\n", " if stat >= obs:\n", " ge += 1\n", " pval = ge / (n_perm + 1)\n", " return obs, pval\n", "\n", "# ----------------------------\n", "# 8) Bootstrap CI for effect size\n", "# ----------------------------\n", "def bootstrap_effect(points_a: np.ndarray, points_b: np.ndarray, n_boot: int, rng: np.random.Generator):\n", " n_a = points_a.shape[0]\n", " n_b = points_b.shape[0]\n", " vals = np.empty(n_boot, dtype=float)\n", " for i in range(n_boot):\n", " ia = rng.integers(0, n_a, size=n_a)\n", " ib = rng.integers(0, n_b, size=n_b)\n", " vals[i] = effect_score(summarize_cloud(points_a[ia]), summarize_cloud(points_b[ib]))\n", " lo, hi = np.quantile(vals, [0.025, 0.975])\n", " return float(vals.mean()), float(lo), float(hi)\n", "\n", "# ----------------------------\n", "# 9) Run experiment\n", "# ----------------------------\n", "N_SAMPLES_PER_N = 1200\n", "N_PERM = 2000\n", "N_BOOT = 1000\n", "NS = [2, 3, 4, 5]\n", "\n", "all_rows = []\n", "all_su = {}\n", "all_null = {}\n", "\n", "for n in NS:\n", " local_rng_su = np.random.default_rng(SEED + 1000 * n + 1)\n", " local_rng_null = np.random.default_rng(SEED + 1000 * n + 2)\n", "\n", " su_pts = sample_phi5_cloud_su(n, N_SAMPLES_PER_N, local_rng_su)\n", " null_pts = sample_phi5_cloud_null(n, N_SAMPLES_PER_N, local_rng_null)\n", "\n", " all_su[n] = su_pts\n", " all_null[n] = null_pts\n", "\n", " su_sum = summarize_cloud(su_pts)\n", " null_sum = summarize_cloud(null_pts)\n", "\n", " local_rng_perm = np.random.default_rng(SEED + 1000 * n + 3)\n", " obs_eff, pval = permutation_test_effect(su_pts, null_pts, N_PERM, local_rng_perm)\n", "\n", " local_rng_boot = np.random.default_rng(SEED + 1000 * n + 4)\n", " boot_mean, ci_lo, ci_hi = bootstrap_effect(su_pts, null_pts, N_BOOT, local_rng_boot)\n", "\n", " all_rows.append(\n", " {\n", " \"n\": n,\n", " \"su_mean_pairwise_distance\": su_sum[\"mean_pairwise_distance\"],\n", " \"null_mean_pairwise_distance\": null_sum[\"mean_pairwise_distance\"],\n", " \"su_centroid_radius\": su_sum[\"centroid_radius\"],\n", " \"null_centroid_radius\": null_sum[\"centroid_radius\"],\n", " \"su_js\": su_sum[\"dominance_js_from_uniform\"],\n", " \"null_js\": null_sum[\"dominance_js_from_uniform\"],\n", " \"effect_observed\": obs_eff,\n", " \"effect_boot_mean\": boot_mean,\n", " \"effect_ci_95_lo\": ci_lo,\n", " \"effect_ci_95_hi\": ci_hi,\n", " \"permutation_pvalue\": pval,\n", " \"su_mean_point\": su_sum[\"mean_point\"],\n", " \"null_mean_point\": null_sum[\"mean_point\"],\n", " \"su_dominance_counts\": su_sum[\"dominance_counts\"],\n", " \"null_dominance_counts\": null_sum[\"dominance_counts\"],\n", " \"su_pca\": np.array(\n", " [\n", " su_sum[\"pca_var_1\"],\n", " su_sum[\"pca_var_2\"],\n", " su_sum[\"pca_var_3\"],\n", " su_sum[\"pca_var_4\"],\n", " su_sum[\"pca_var_5\"],\n", " ],\n", " dtype=float,\n", " ),\n", " \"null_pca\": np.array(\n", " [\n", " null_sum[\"pca_var_1\"],\n", " null_sum[\"pca_var_2\"],\n", " null_sum[\"pca_var_3\"],\n", " null_sum[\"pca_var_4\"],\n", " null_sum[\"pca_var_5\"],\n", " ],\n", " dtype=float,\n", " ),\n", " }\n", " )\n", "\n", "# ----------------------------\n", "# 10) Global pooled test across n = 2..5\n", "# ----------------------------\n", "pooled_su = np.vstack([all_su[n] for n in NS])\n", "pooled_null = np.vstack([all_null[n] for n in NS])\n", "\n", "pooled_su_sum = summarize_cloud(pooled_su)\n", "pooled_null_sum = summarize_cloud(pooled_null)\n", "\n", "pooled_rng_perm = np.random.default_rng(SEED + 999001)\n", "pooled_eff, pooled_p = permutation_test_effect(pooled_su, pooled_null, N_PERM, pooled_rng_perm)\n", "\n", "pooled_rng_boot = np.random.default_rng(SEED + 999002)\n", "pooled_boot_mean, pooled_ci_lo, pooled_ci_hi = bootstrap_effect(pooled_su, pooled_null, N_BOOT, pooled_rng_boot)\n", "\n", "# ----------------------------\n", "# 11) Hard fail / pass rules\n", "# ----------------------------\n", "# Strict rules:\n", "# PASS if:\n", "# - pooled effect > 0\n", "# - pooled permutation p-value <= 0.01\n", "# - pooled bootstrap CI lower bound > 0\n", "# Else FAIL\n", "strict_pass = (\n", " (pooled_eff > 0.0)\n", " and (pooled_p <= 0.01)\n", " and (pooled_ci_lo > 0.0)\n", ")\n", "\n", "# ----------------------------\n", "# 12) Report\n", "# ----------------------------\n", "np.set_printoptions(precision=6, suppress=True)\n", "\n", "print(\"=\" * 78)\n", "print(\"STRICT NULL TEST: SU(n) vs MATCHED NON-UNITARY NULL\")\n", "print(\"Projection: boost-free Phi5\")\n", "print(f\"Seed: {SEED}\")\n", "print(f\"Samples per n: {N_SAMPLES_PER_N}\")\n", "print(f\"Permutation draws per test: {N_PERM}\")\n", "print(f\"Bootstrap draws per test: {N_BOOT}\")\n", "print(\"=\" * 78)\n", "\n", "for row in all_rows:\n", " print(f\"\\n--- n = {row['n']} ---\")\n", " print(f\"SU mean_pairwise_distance : {row['su_mean_pairwise_distance']:.9f}\")\n", " print(f\"NULL mean_pairwise_distance : {row['null_mean_pairwise_distance']:.9f}\")\n", " print(f\"SU centroid_radius : {row['su_centroid_radius']:.9f}\")\n", " print(f\"NULL centroid_radius : {row['null_centroid_radius']:.9f}\")\n", " print(f\"SU JS(dominance||uniform) : {row['su_js']:.9f}\")\n", " print(f\"NULL JS(dominance||uniform) : {row['null_js']:.9f}\")\n", " print(f\"Observed effect : {row['effect_observed']:.9f}\")\n", " print(f\"Bootstrap mean : {row['effect_boot_mean']:.9f}\")\n", " print(f\"Bootstrap 95% CI : [{row['effect_ci_95_lo']:.9f}, {row['effect_ci_95_hi']:.9f}]\")\n", " print(f\"Permutation p-value : {row['permutation_pvalue']:.9f}\")\n", " print(f\"SU mean point : {row['su_mean_point']}\")\n", " print(f\"NULL mean point : {row['null_mean_point']}\")\n", " print(f\"SU dominance counts : {row['su_dominance_counts']}\")\n", " print(f\"NULL dominance counts : {row['null_dominance_counts']}\")\n", " print(f\"SU PCA explained var : {row['su_pca']}\")\n", " print(f\"NULL PCA explained var : {row['null_pca']}\")\n", "\n", "print(\"\\n\" + \"=\" * 78)\n", "print(\"POOLED TEST (n = 2..5)\")\n", "print(f\"SU mean_pairwise_distance : {pooled_su_sum['mean_pairwise_distance']:.9f}\")\n", "print(f\"NULL mean_pairwise_distance : {pooled_null_sum['mean_pairwise_distance']:.9f}\")\n", "print(f\"SU centroid_radius : {pooled_su_sum['centroid_radius']:.9f}\")\n", "print(f\"NULL centroid_radius : {pooled_null_sum['centroid_radius']:.9f}\")\n", "print(f\"SU JS(dominance||uniform) : {pooled_su_sum['dominance_js_from_uniform']:.9f}\")\n", "print(f\"NULL JS(dominance||uniform) : {pooled_null_sum['dominance_js_from_uniform']:.9f}\")\n", "print(f\"Observed effect : {pooled_eff:.9f}\")\n", "print(f\"Bootstrap mean : {pooled_boot_mean:.9f}\")\n", "print(f\"Bootstrap 95% CI : [{pooled_ci_lo:.9f}, {pooled_ci_hi:.9f}]\")\n", "print(f\"Permutation p-value : {pooled_p:.9f}\")\n", "print(f\"SU mean point : {pooled_su_sum['mean_point']}\")\n", "print(f\"NULL mean point : {pooled_null_sum['mean_point']}\")\n", "print(f\"SU dominance counts : {pooled_su_sum['dominance_counts']}\")\n", "print(f\"NULL dominance counts : {pooled_null_sum['dominance_counts']}\")\n", "print(f\"SU PCA explained var : {pca_explained_variance(pooled_su)}\")\n", "print(f\"NULL PCA explained var : {pca_explained_variance(pooled_null)}\")\n", "print(\"=\" * 78)\n", "print(f\"STRICT VERDICT: {'PASS' if strict_pass else 'FAIL'}\")\n", "print(\"=\" * 78)\n", "\n", "# ----------------------------\n", "# 13) Machine-readable result\n", "# ----------------------------\n", "RESULTS = {\n", " \"seed\": SEED,\n", " \"samples_per_n\": N_SAMPLES_PER_N,\n", " \"n_perm\": N_PERM,\n", " \"n_boot\": N_BOOT,\n", " \"axes\": AXES,\n", " \"weights\": CONST_WEIGHTS.copy(),\n", " \"per_n\": all_rows,\n", " \"pooled\": {\n", " \"effect_observed\": pooled_eff,\n", " \"effect_boot_mean\": pooled_boot_mean,\n", " \"effect_ci_95_lo\": pooled_ci_lo,\n", " \"effect_ci_95_hi\": pooled_ci_hi,\n", " \"permutation_pvalue\": pooled_p,\n", " \"strict_pass\": strict_pass,\n", " \"su_mean_point\": pooled_su_sum[\"mean_point\"].copy(),\n", " \"null_mean_point\": pooled_null_sum[\"mean_point\"].copy(),\n", " \"su_dominance_counts\": pooled_su_sum[\"dominance_counts\"].copy(),\n", " \"null_dominance_counts\": pooled_null_sum[\"dominance_counts\"].copy(),\n", " \"su_pca\": pca_explained_variance(pooled_su).copy(),\n", " \"null_pca\": pca_explained_variance(pooled_null).copy(),\n", " },\n", "}" ], "metadata": { "colab": { "base_uri": "https://localhost:8080/" }, "id": "OJLzwftl8TeJ", "outputId": "574646f5-300d-477d-bd4c-09f1091bc50a" }, "execution_count": null, "outputs": [ { "output_type": "stream", "name": "stdout", "text": [ "==============================================================================\n", "STRICT NULL TEST: SU(n) vs MATCHED NON-UNITARY NULL\n", "Projection: boost-free Phi5\n", "Seed: 20260318\n", "Samples per n: 1200\n", "Permutation draws per test: 2000\n", "Bootstrap draws per test: 1000\n", "==============================================================================\n", "\n", "--- n = 2 ---\n", "SU mean_pairwise_distance : 0.149645803\n", "NULL mean_pairwise_distance : 0.155747316\n", "SU centroid_radius : 0.110545091\n", "NULL centroid_radius : 0.111991340\n", "SU JS(dominance||uniform) : 0.422810455\n", "NULL JS(dominance||uniform) : 0.420185318\n", "Observed effect : 0.010172899\n", "Bootstrap mean : 0.009949573\n", "Bootstrap 95% CI : [-0.001441262, 0.021432215]\n", "Permutation p-value : 0.072463768\n", "SU mean point : [0.155564 0.191619 0. 0. 0.652816]\n", "NULL mean point : [0.180396 0.217657 0. 0.184243 0.417704]\n", "SU dominance counts : [ 0 0 0 0 1200]\n", "NULL dominance counts : [ 0 0 0 1 1199]\n", "SU PCA explained var : [0.670408 0.329592 0. 0. 0. ]\n", "NULL PCA explained var : [0.544009 0.283412 0.172579 0. 0. ]\n", "\n", "--- n = 3 ---\n", "SU mean_pairwise_distance : 0.080140853\n", "NULL mean_pairwise_distance : 0.098210417\n", "SU centroid_radius : 0.057100231\n", "NULL centroid_radius : 0.069634569\n", "SU JS(dominance||uniform) : 0.422810455\n", "NULL JS(dominance||uniform) : 0.406078855\n", "Observed effect : 0.047335502\n", "Bootstrap mean : 0.047164434\n", "Bootstrap 95% CI : [0.037139814, 0.057086962]\n", "Permutation p-value : 0.000499750\n", "SU mean point : [0.134624 0.188313 0.151631 0. 0.525431]\n", "NULL mean point : [0.148287 0.208937 0.162255 0.1682 0.312321]\n", "SU dominance counts : [ 0 0 0 0 1200]\n", "NULL dominance counts : [ 0 0 0 10 1190]\n", "SU PCA explained var : [0.574561 0.294217 0.131222 0. 0. ]\n", "NULL PCA explained var : [0.413691 0.298031 0.214992 0.073286 0. ]\n", "\n", "--- n = 4 ---\n", "SU mean_pairwise_distance : 0.048235942\n", "NULL mean_pairwise_distance : 0.066611426\n", "SU centroid_radius : 0.034028569\n", "NULL centroid_radius : 0.047187820\n", "SU JS(dominance||uniform) : 0.422810455\n", "NULL JS(dominance||uniform) : 0.415511113\n", "Observed effect : 0.038834079\n", "Bootstrap mean : 0.038335309\n", "Bootstrap 95% CI : [0.030832059, 0.046481698]\n", "Permutation p-value : 0.000499750\n", "SU mean point : [0.134554 0.204713 0.132485 0. 0.528247]\n", "NULL mean point : [0.150457 0.225937 0.142678 0.175879 0.305048]\n", "SU dominance counts : [ 0 0 0 0 1200]\n", "NULL dominance counts : [ 0 1 0 2 1197]\n", "SU PCA explained var : [0.512572 0.381629 0.1058 0. 0. ]\n", "NULL PCA explained var : [0.486321 0.267875 0.179428 0.066376 0. ]\n", "\n", "--- n = 5 ---\n", "SU mean_pairwise_distance : 0.038805054\n", "NULL mean_pairwise_distance : 0.054143041\n", "SU centroid_radius : 0.027414662\n", "NULL centroid_radius : 0.038228174\n", "SU JS(dominance||uniform) : 0.422810455\n", "NULL JS(dominance||uniform) : 0.418136133\n", "Observed effect : 0.030825821\n", "Bootstrap mean : 0.030409039\n", "Bootstrap 95% CI : [0.024234240, 0.037285875]\n", "Permutation p-value : 0.000499750\n", "SU mean point : [0.13666 0.213636 0.119122 0. 0.530582]\n", "NULL mean point : [0.150136 0.236358 0.128915 0.182175 0.302415]\n", "SU dominance counts : [ 0 0 0 0 1200]\n", "NULL dominance counts : [ 0 0 0 2 1198]\n", "SU PCA explained var : [0.484002 0.413297 0.102701 0. 0. ]\n", "NULL PCA explained var : [0.504957 0.263199 0.175784 0.05606 0. ]\n", "\n", "==============================================================================\n", "POOLED TEST (n = 2..5)\n", "SU mean_pairwise_distance : 0.124958988\n", "NULL mean_pairwise_distance : 0.137382969\n", "SU centroid_radius : 0.089281390\n", "NULL centroid_radius : 0.097939897\n", "SU JS(dominance||uniform) : 0.422810455\n", "NULL JS(dominance||uniform) : 0.414222363\n", "Observed effect : 0.029670580\n", "Bootstrap mean : 0.029340719\n", "Bootstrap 95% CI : [0.021229223, 0.037592473]\n", "Permutation p-value : 0.000499750\n", "SU mean point : [0.140351 0.19957 0.10081 0. 0.559269]\n", "NULL mean point : [0.157319 0.222222 0.108462 0.177624 0.334372]\n", "SU dominance counts : [ 0 0 0 0 4800]\n", "NULL dominance counts : [ 0 1 0 15 4784]\n", "SU PCA explained var : [0.680588 0.203847 0.115565 0. 0. ]\n", "NULL PCA explained var : [0.594409 0.209741 0.127639 0.068211 0. ]\n", "==============================================================================\n", "STRICT VERDICT: PASS\n", "==============================================================================\n" ] } ] } ] }