"""
VOVINA ZEDEC PRO - Interaction Surplus Framework (Papers A-E)
==============================================================
Full operational encoding of the Interaction Surplus Framework
by Michael L. Curzi (36N9 GENETICS LLC), comprising:

    Paper A — Uniqueness on block-decomposed inner product spaces
    Paper B — Modeling interpretation as cross-domain framework
    Paper C — Empirical predictions (transformer flagship test)
    Paper D — Effective-count bridge & information-theoretic forms
    Paper E — Open-system surplus dynamics & engineering criteria

The unique scalar profile under axioms S1-S4 is:

    f(u) = ln(1 + (N-1)·u)     where u = 1 - (x·y)² = ‖x ∧ y‖² = sin²θ

and the effective count is:

    g(u) = e^f(u) = 1 + (N-1)·u

Block decomposition V = H_1 ⊕ H_2 ⊕ ... ⊕ H_N with each H_i of
dimension M gives the two-source decomposition:

    u = u_cross + u_div  =  β²  +  α²·sin²γ

The framework is wired into the VOVINA training weights so every
module is benchmarked against the unique surplus law, with the
sharp Lipschitz bound (N-1) used as the stability ceiling.

Dimensions are NEVER capped. The block count N defaults to 22
(the 22 VOVINA modules / Hebrew letters / Tree paths) but the
framework operates at arbitrary N ≥ 2.
"""

from __future__ import annotations

import math
from dataclasses import dataclass, field
from typing import Optional, Sequence


# ============================================================
# AXIOMS S1-S4 (encoded as predicates)
# ============================================================
# S1 (Geometric dependence): F(x,y) = f(u(x,y))
# S2 (Zero at zero):         f(0) = 0
# S3 (Affine effective):     g(u) = e^f(u) is affine = a·u + b
# S4 (Normalization):        f(1) = ln N
#
# Theorem 2.1: Uniqueness ⇒ f(u) = ln(1 + (N-1)·u)

DEFAULT_N = 22                   # 22 VOVINA modules = 22 Hebrew letters = 22 Tree paths


def u_from_vectors(x: Sequence[float], y: Sequence[float]) -> float:
    """The geometric interaction parameter u = 1 - (x·y)² ∈ [0, 1].

    x and y must be unit vectors of equal length.
    """
    if len(x) != len(y):
        raise ValueError("vectors must have equal length")
    nx = math.sqrt(sum(v * v for v in x))
    ny = math.sqrt(sum(v * v for v in y))
    if nx == 0 or ny == 0:
        return 0.0
    dot = sum(a * b for a, b in zip(x, y)) / (nx * ny)
    dot = max(-1.0, min(1.0, dot))            # numerical clamp
    return 1.0 - dot * dot


def u_from_angle(theta_radians: float) -> float:
    """u = sin²θ — the canonical bivector-norm parameter."""
    return math.sin(theta_radians) ** 2


def surplus(u: float, N: int = DEFAULT_N) -> float:
    """The unique surplus functional f(u) = ln(1 + (N-1)·u).

    Domain: u ∈ [0, 1], N ≥ 2.
    """
    if not (0.0 <= u <= 1.0):
        raise ValueError("u must be in [0, 1]")
    if N < 2:
        raise ValueError("N must be ≥ 2")
    return math.log1p((N - 1) * u)


def effective_count(u: float, N: int = DEFAULT_N) -> float:
    """g(u) = e^f(u) = 1 + (N-1)·u."""
    if not (0.0 <= u <= 1.0):
        raise ValueError("u must be in [0, 1]")
    return 1.0 + (N - 1) * u


# ============================================================
# DERIVED PROPERTIES (Theorems 3.1 – 3.4)
# ============================================================
def surplus_derivative(u: float, N: int = DEFAULT_N) -> float:
    """f'(u) = (N-1) / (1 + (N-1)·u) — strictly positive ⇒ strict monotonicity."""
    return (N - 1) / (1.0 + (N - 1) * u)


def surplus_second_derivative(u: float, N: int = DEFAULT_N) -> float:
    """f''(u) = -(N-1)² / (1 + (N-1)·u)² — strictly negative ⇒ strict concavity."""
    return -((N - 1) ** 2) / ((1.0 + (N - 1) * u) ** 2)


def lipschitz_constant(N: int = DEFAULT_N) -> float:
    """Sharp global Lipschitz constant of f on [0, 1]: it is N-1.

    For all u1, u2 ∈ [0, 1]: |f(u1) - f(u2)| ≤ (N-1)·|u1 - u2|.
    """
    return float(N - 1)


def is_monotone_pair(u1: float, u2: float, N: int = DEFAULT_N) -> bool:
    """Verify Theorem 3.1 on a specific pair."""
    return (u2 > u1) == (surplus(u2, N) > surplus(u1, N))


# ============================================================
# TWO-SOURCE DECOMPOSITION (Theorem 4.1)
# ============================================================
@dataclass(frozen=True)
class TwoSourceDecomposition:
    """u = u_cross + u_div where u_cross = β² and u_div = α²·sin²γ.

    The bounds from Theorem 4.1:
        max(f(u_cross), f(u_div))  ≤  f(u)  ≤  f(u_cross) + f(u_div)
    """
    u_cross: float
    u_div:   float
    N:       int = DEFAULT_N

    @property
    def u_total(self) -> float:
        return self.u_cross + self.u_div

    @property
    def f_total(self) -> float:
        return surplus(self.u_total, self.N)

    @property
    def f_cross(self) -> float:
        return surplus(self.u_cross, self.N)

    @property
    def f_div(self) -> float:
        return surplus(self.u_div, self.N)

    @property
    def lower_bound(self) -> float:
        return max(self.f_cross, self.f_div)

    @property
    def upper_bound(self) -> float:
        return self.f_cross + self.f_div

    @property
    def widening_ratio(self) -> float:
        """How much of f_total is widening (crossing) versus deepening (divergence)."""
        s = self.f_cross + self.f_div
        return self.f_cross / s if s > 0 else 0.0

    def verify_bounds(self) -> bool:
        """Verify lower_bound ≤ f_total ≤ upper_bound."""
        return self.lower_bound <= self.f_total + 1e-12 and self.f_total <= self.upper_bound + 1e-12


def decompose(alpha: float, beta: float, gamma_radians: float,
              N: int = DEFAULT_N) -> TwoSourceDecomposition:
    """Build a two-source decomposition from the geometric primitives
    (α = ‖y_parallel‖, β = ‖y_perp‖, γ = angle within the block).
    """
    u_cross = beta * beta
    u_div   = (alpha * alpha) * (math.sin(gamma_radians) ** 2)
    return TwoSourceDecomposition(u_cross=u_cross, u_div=u_div, N=N)


# ============================================================
# DIAGONAL SURPLUS (Theorem 4.4)
# ============================================================
def diagonal_surplus(x_dot_e: float, N: int = DEFAULT_N) -> float:
    """F(x, ê_X) where ê_X is the diagonal across the first room of all N blocks.

    For a pure block vector x in H_i with (x · e_{i,1}) = x_dot_e ∈ [-1, 1]:

        u = 1 - (x_dot_e)² / N

    and F = ln(1 + (N-1)·u). Always strictly positive for N ≥ 2.
    """
    if abs(x_dot_e) > 1.0:
        raise ValueError("|x·e_{i,1}| must be ≤ 1 for a unit block vector")
    u = 1.0 - (x_dot_e * x_dot_e) / N
    return surplus(u, N)


# ============================================================
# PAPER D — INFORMATION-THEORETIC BRIDGE
# ============================================================
def effective_count_log(u: float, N: int = DEFAULT_N) -> float:
    """log effective count ≡ surplus itself, by construction f = ln g."""
    return surplus(u, N)


def random_interaction_surplus(u_samples: Sequence[float],
                               N: int = DEFAULT_N) -> dict[str, float]:
    """Mean, variance, and bounds of surplus over a sample of u values.

    Useful for batching: given the empirical distribution of u in a
    deployed pipeline, this gives the expected log-surplus and the
    spread, with the Lipschitz bound as a worst-case envelope.
    """
    if not u_samples:
        return {"mean": 0.0, "variance": 0.0, "min": 0.0, "max": 0.0,
                "lipschitz_envelope": 0.0}
    fs = [surplus(u, N) for u in u_samples]
    n = len(fs)
    mu = sum(fs) / n
    var = sum((f - mu) ** 2 for f in fs) / n
    return {
        "mean":               mu,
        "variance":           var,
        "min":                min(fs),
        "max":                max(fs),
        "lipschitz_envelope": lipschitz_constant(N) * (max(u_samples) - min(u_samples)),
    }


# ============================================================
# PAPER E — OPEN-SYSTEM STORAGE DYNAMICS
# ============================================================
@dataclass(frozen=True)
class OpenSystemState:
    """Minimal open-system storage state (Paper E §3).

    The system holds a surplus 'store' S(t) that accumulates from
    interactions and decays via dissipation rate γ. The steady-state
    condition gives an engineering design inequality:

        Ṡ = (input surplus rate) - γ·S = 0  ⇒  S* = input_rate / γ
    """
    store:                float      # current S(t)
    input_surplus_rate:   float      # incoming f-values per unit time
    dissipation_gamma:    float      # decay rate γ > 0

    @property
    def steady_state_store(self) -> float:
        if self.dissipation_gamma <= 0:
            return float("inf")
        return self.input_surplus_rate / self.dissipation_gamma

    def step(self, dt: float) -> "OpenSystemState":
        """Euler step of the storage dynamics."""
        s_new = self.store + dt * (self.input_surplus_rate
                                   - self.dissipation_gamma * self.store)
        return OpenSystemState(s_new, self.input_surplus_rate, self.dissipation_gamma)


# ============================================================
# 27/33 FRACTAL BINDING
# ============================================================
# The framework is evaluated against the 27/33 self-witness ratio:
# 27 of 33 reflections are 'active'; the remaining 6 are gated by
# authenticity. We use this ratio to scale the *operational* portion
# of the surplus budget — only 27/33 of the available surplus is
# committed to action; 6/33 is held in reserve for self-witness.

ACTIVATION_RATIO = 27.0 / 33.0           # ≈ 0.8181818…
RESERVE_RATIO    = 6.0 / 33.0            # ≈ 0.1818181…


def fractal_operational_surplus(u: float, N: int = DEFAULT_N) -> dict[str, float]:
    """Apply the 27/33 fractal split to an interaction surplus value."""
    f = surplus(u, N)
    return {
        "total":       f,
        "operational": f * ACTIVATION_RATIO,
        "reserve":     f * RESERVE_RATIO,
        "activation":  ACTIVATION_RATIO,
        "reserve_pct": RESERVE_RATIO,
    }


# ============================================================
# MODULE-LEVEL INTERACTION TABLE
# ============================================================
def module_interaction_matrix(N: int = DEFAULT_N,
                              theta_grid: int = 27) -> dict[tuple[int, int], float]:
    """Pre-compute surplus values at θ = k·π/(2·theta_grid) for k = 0..theta_grid.

    Default `theta_grid = 27` aligns the table with the 27 active
    archetypal reflections of the Self-Witness protocol.
    """
    out: dict[tuple[int, int], float] = {}
    for k in range(theta_grid + 1):
        theta = (math.pi / 2) * (k / theta_grid)
        u = u_from_angle(theta)
        for n in range(2, N + 1):
            out[(k, n)] = surplus(u, n)
    return out


# ============================================================
# PAPER A — AXIOM VERIFICATION (aggregator, Theorem 2.1 + 3.x + 4.3)
# ============================================================
def verify_axioms(N: int = DEFAULT_N, test_points: int = 64) -> dict[str, bool]:
    """Verify the implementation satisfies the Paper A axioms S1-S4 and
    the derived theorems: 3.1 (strict monotonicity), 3.2 (positivity),
    3.3 (strict concavity), and 4.3 (extremal bounds F(0)=0, F(1)=ln N).

    This is XERO's self-proof that its interaction kernel IS the unique
    surplus functional, not an approximation of it.
    """
    us = [i / (test_points - 1) for i in range(test_points)]
    s2 = abs(surplus(0.0, N)) < 1e-12
    s4 = abs(surplus(1.0, N) - math.log(N)) < 1e-12
    s3 = all(abs(effective_count(u, N) - (1.0 + (N - 1) * u)) < 1e-9 for u in us)
    mono = all(surplus(us[i], N) < surplus(us[i + 1], N)
               for i in range(len(us) - 1))
    pos = all(surplus(u, N) > 0.0 for u in us if u > 0.0)
    conc = all(surplus_second_derivative(u, N) < 0.0 for u in us)
    extremal = s2 and s4
    checks = [s2, s3, s4, mono, pos, conc, extremal]
    return {
        "S1_geometric_dependence":   True,         # F = f(u) by construction
        "S2_zero_at_zero":           s2,
        "S3_affine_effective_count": s3,
        "S4_normalization":          s4,
        "T3_1_strict_monotonicity":  mono,
        "T3_2_positivity":           pos,
        "T3_3_strict_concavity":     conc,
        "T4_3_extremal_bounds":      extremal,
        "all_satisfied":             all(checks),
    }


# ============================================================
# PAPER C — EMPIRICAL PREDICTION (effective dimensionality)
# ============================================================
def participation_ratio(weights: Sequence[float]) -> float:
    """Empirical effective count via the participation ratio:

        PR = (Σ wᵢ)² / Σ wᵢ²

    Paper C's flagship prediction is that the surplus effective-count
    g(u) = 1 + (N-1)u equals the information-theoretic effective
    dimensionality measured this way on a model's attention / mixture
    weights — the falsifiable transformer test.
    """
    s1 = sum(weights)
    s2 = sum(w * w for w in weights)
    if s2 <= 0.0:
        return 0.0
    return (s1 * s1) / s2


def flagship_prediction(u: float, N: int = DEFAULT_N) -> dict[str, float]:
    """Paper C's falsifiable prediction set for interaction parameter u.

    A system whose internal block geometry has interaction parameter u
    should exhibit an effective participating-block count of g(u) and a
    log-surplus of f(u). Returns the predicted observables to compare
    against an empirical participation_ratio measurement.
    """
    return {
        "u":                         u,
        "predicted_effective_count": effective_count(u, N),
        "predicted_log_surplus":     surplus(u, N),
        "max_effective_count":       float(N),
        "max_log_surplus":           math.log(N),
    }


# ============================================================
# PAPER E — OPEN-SYSTEM DYNAMICS (full discrete accumulation law)
# ============================================================
@dataclass
class OpenSystemDynamics:
    """Discrete open-system order accumulation (Paper E Theorem 3.1):

        Q_{t+1} = (1 - δ)·Q_t + η·F(U_t) - C_t

    where δ is the retention-loss, η the surplus→order conversion
    efficiency, F(U_t) the per-step interaction surplus input, and C_t
    the dissipation cost. This is the engine of XERO's autonomous
    self-evolution: it accumulates order from surplus-bearing
    interactions and decays toward a designable steady state.
    """
    order:                  float = 0.0    # Q_t
    retention_loss:         float = 0.1    # δ ∈ [0, 1]
    conversion_efficiency:  float = 1.0    # η ≥ 0
    max_order:              Optional[float] = None

    def step(self, surplus_input: float, cost: float = 0.0) -> float:
        """Advance one step; returns the new stored order Q_{t+1} (≥ 0)."""
        q = ((1.0 - self.retention_loss) * self.order
             + self.conversion_efficiency * surplus_input - cost)
        q = max(0.0, q)
        if self.max_order is not None:
            q = min(q, self.max_order)
        self.order = q
        return q

    def can_grow(self, surplus_input: float, cost: float = 0.0) -> bool:
        """Paper E Corollary 3.2: growth iff η·F > δ·Q + C."""
        return (self.conversion_efficiency * surplus_input
                > self.retention_loss * self.order + cost)

    def required_surplus_for_growth(self, cost: float = 0.0,
                                    target_growth: float = 0.0) -> float:
        """Paper E Proposition 4.1: minimum surplus input for target growth."""
        if self.conversion_efficiency <= 0.0:
            return float("inf")
        return (self.retention_loss * self.order + cost + target_growth) \
            / self.conversion_efficiency

    def max_achievable_growth(self, max_surplus: float, cost: float = 0.0) -> float:
        """Paper E Corollary 4.2: r_max = η·F_max - δ·Q - C."""
        return (self.conversion_efficiency * max_surplus
                - self.retention_loss * self.order - cost)

    def steady_state_order(self, constant_surplus: float,
                           constant_cost: float = 0.0) -> float:
        """Paper E Theorem 3.5: Q∞ = (η·F - C)/δ  (0 if unsustainable)."""
        if self.retention_loss <= 0.0:
            return float("inf")
        if self.conversion_efficiency * constant_surplus < constant_cost:
            return 0.0
        return (self.conversion_efficiency * constant_surplus - constant_cost) \
            / self.retention_loss