HAKO upload: hako/core/plateau.py

hako/core/plateau.py

@@ -0,0 +1,112 @@
+"""2D plateau: kernel density of prototypes under the mapping Pi: R^N -> R^2.
+
+Theorem T-PLATEAU (existence of a compact non-empty plateau; Morse analytics).
+Let Pi: R^N -> R^2 be continuous (here: PCA whitening followed by a trainable
+linear refinement; PCA projection is linear, hence analytic) and let
+    rho(x, y) = sum_{k=1..K} exp( -|| Pi(mu_k) - (x, y) ||^2 / (2 sigma^2) ).
+(1) rho is real-analytic on R^2 (finite sum of Gaussian bumps, each analytic).
+(2) For generic sigma > 0, rho is a Morse function: its critical points are
+    non-degenerate (degeneracy requires an algebraic coincidence of the K
+    bump centers, a measure-zero set in sigma; Sard's theorem argument).
+(3) rho(x,y) -> 0 as ||(x,y)|| -> 0-distance-to-infinity, and
+    max rho >= rho(Pi(mu_1)) >= 1 (the bump centered at Pi(mu_1) contributes
+    >= 1 there). Hence for any delta with 0 < delta < max rho the superlevel
+    set  P_delta = {(x,y) : rho >= delta}  is non-empty (contains a global
+    maximizer, which exists because rho is continuous and vanishes at
+    infinity => attains max on a compact set) and compact (closed superlevel
+    set of a coercive function). QED.
+
+The orchestrator monitors: plateau mass M = rho(argmax), gradient norm at the
+current projected point, and saturation = 1 - entropy of the normalized
+bump contributions (low entropy = one dominant prototype = saturation signal
+feeding the dimensional-growth trigger, per PHASE1 inertia rule).
+"""
+from __future__ import annotations
+
+from typing import Tuple
+
+import numpy as np
+import torch
+
+
+class Plateau2D:
+    def __init__(self, sigma: float = 0.55, n_components: int = 2) -> None:
+        self.sigma = float(sigma)
+        self.mean = None
+        self.comp = None          # (N, 2) PCA components
+        self.refined = np.zeros((2, 2), dtype=np.float32)  # trainable refine I
+
+    @staticmethod
+    def _np(x):
+        """Torch/numpy agnostic conversion (detached; plateau coords are
+        monitoring/control signals, not gradient carriers)."""
+        if isinstance(x, torch.Tensor):
+            return x.detach().cpu().numpy().astype(np.float32)
+        return np.asarray(x, dtype=np.float32)
+
+    # ------------------------------------------------------------- fitting
+    def fit_pca(self, Z: np.ndarray) -> None:
+        X = np.asarray(Z, dtype=np.float32)
+        self.mean = X.mean(axis=0)
+        Xc = X - self.mean
+        # economical SVD on the covariance (N x 2 subspace)
+        C = Xc.T @ Xc / max(1, len(X))
+        w, V = np.linalg.eigh(C)
+        order = np.argsort(w)[::-1][:2]
+        self.comp = V[:, order].astype(np.float32)      # (N, 2)
+
+    def project(self, Z: np.ndarray) -> np.ndarray:
+        """Pi(z) = refined PCA coords (zeros until fit_pca has run --
+        plateau coords are monitoring signals, never gradient carriers)."""
+        Xc = self._np(Z)
+        if self.mean is None or self.comp is None:
+            return np.zeros((len(Xc), 2), dtype=np.float32)
+        base = (Xc - self.mean) @ self.comp
+        return (base @ self.refined.T).astype(np.float32)
+
+    def refine_gradient(self, coords: np.ndarray, proto_coords: np.ndarray,
+                        lr: float = 1e-3) -> float:
+        """One gradient step aligning Pi(z) mass toward prototype plateau.
+
+        Uses the analytic gradient of the local density at the sample points;
+        returns the mean gradient norm (monitoring signal for the retrain
+        rule: if ||grad rho|| > delta, Pi is refined for one epoch).
+        """
+        if self.refined is None:
+            return 0.0
+        base = (np.asarray(coords, dtype=np.float32) @ self.refined.T)
+        diff = base[:, None, :] - proto_coords[None, :, :]
+        d2 = (diff ** 2).sum(-1)
+        w = np.exp(-d2 / (2 * self.sigma ** 2))
+        grad_pts = -(diff * w[..., None]).sum(1) / (self.sigma ** 2)
+        g = grad_pts.mean(0)                                   # (2,)
+        # chain rule back to refined (2x2): dL/dR = outer(g, base_mean)
+        self.refined -= lr * np.outer(g, coords.mean(0)).astype(np.float32)
+        return float(np.linalg.norm(g))
+
+    # ------------------------------------------------------------ monitoring
+    def density(self, pts: np.ndarray, proto_coords: np.ndarray) -> np.ndarray:
+        pts = self._np(pts)
+        proto_coords = self._np(proto_coords)
+        diff = pts[:, None, :] - proto_coords[None, :, :]
+        d2 = (diff ** 2).sum(-1)
+        return np.exp(-d2 / (2 * self.sigma ** 2)).sum(1)
+
+    def saturation(self, proto_coords: np.ndarray) -> float:
+        """1 - normalized entropy of bump mass at the density argmax."""
+        pts = proto_coords
+        dens = self.density(pts, proto_coords)
+        p = dens / max(dens.sum(), 1e-12)
+        H = -(p * np.log(p + 1e-12)).sum()
+        Hn = H / np.log(len(p))
+        return float(1.0 - Hn)
+
+    def inertia(self, proto_coords: np.ndarray) -> float:
+        """Mean squared spread of prototypes on the plateau (PHASE1 trigger)."""
+        c = proto_coords.mean(0, keepdims=True)
+        return float(((proto_coords - c) ** 2).sum(1).mean())
+
+    def coords_of(self, Z: np.ndarray, protos: np.ndarray
+                  ) -> Tuple[np.ndarray, np.ndarray]:
+        pc = self.project(protos)
+        return self.project(Z), pc