#!/usr/bin/env python3
"""
PROVENANCE: SURROGATE

Computation 130 -- the A6 bias residual c_u is a cell-volume EQUIDISTRIBUTION defect of
the emergence cloud; matched-scaling (the A1 count) does NOT force it (open_research §1 A6)
=================================================================================
Backs the cycle-6 (shape-regularity crux) result, which CORRECTED the cycle-5 reading:
the deterministic constant
    c_u = (1/|D_n|) sum_a |grad u(rho_n(a))|^2  /  ((1/|K|) int_K |grad u|^2)
is a UNIFORM-WEIGHT NODAL / cell-volume equidistribution defect of the point cloud --
NOT the Delaunay finite-element interpolation error that Brenner-Scott bounds by mesh
shape-regularity. c_u -> 1 iff the emergence nodes volume-equidistribute (each node
represents an equal Voronoi volume); shape-regularity is neither necessary nor sufficient.

The A6 matched-scaling assumption A1 is a DENSITY / COUNTING condition (|D_n|/V_M ~ d_0^{-4}),
fixing only the AVERAGE spacing -- not a geometric regularity condition on the realised
cloud. This script shows that A1 (a fixed total count) does NOT force c_u -> 1: a cloud with
the SAME count but a non-uniform local density gives a biased c_u that does NOT decay, while
only a volume-equidistributed cloud gives c_u -> 1. Hence the bias half is conditional on a
deterministic equidistribution lemma, the sibling of the open A2 weak-discrepancy, and not on
the point count. (PST-specific aggravation: the signed-tension metric manufactures vanishing-separation pairs
on the null cone, so A1 does not even force a separation lower bound.)

CHECKS
------
  (1) VOLUME-EQUIDISTRIBUTED cloud (deterministic jittered grid): c_u -> 1, defect -> 0 as
      the count grows -- equidistribution gives the closure.
  (2) NON-EQUIDISTRIBUTED cloud (uniform on [0,1]^2 minus a central disk void in a low-f
      region): c_u converges to a BIASED value != 1 that does NOT decay under refinement --
      the asymptotic density is correct everywhere outside the void, only equidistribution
      fails. The bias equals the deterministic void defect <f>_{K\\void}/<f>_K, f = |grad u|^2.
      This is exactly the matched-count near-null void/sliver family of the crux's L1/L3.
  (3) The defect grows monotonically with the void radius (the non-equidistribution amplitude),
      matching the analytic void bias -- isolating cell-volume equidistribution (not the count,
      not shape-regularity) as the binding property.

A "pass": (1) defect decays to ~0 for the uniform cloud; (2) the void cloud's defect is flat
(refinement-independent) and matches the analytic void bias; (3) defect grows with the void radius.
"""

from __future__ import annotations
import numpy as np

CONT = np.pi ** 2 / 2.0   # (1/|K|) int_{[0,1]^2} |grad sin(pi x) sin(pi y)|^2 = pi^2/2


def gradf(P):
    x, y = P[:, 0], P[:, 1]
    return np.pi ** 2 * (np.cos(np.pi * x) ** 2 * np.sin(np.pi * y) ** 2
                         + np.sin(np.pi * x) ** 2 * np.cos(np.pi * y) ** 2)


def c_u(P):
    """uniform-weight nodal quadrature defect c_u = mean_a f(rho_a) / continuum-average(f)."""
    return gradf(P).mean() / CONT


def uniform_cloud(N, rng):
    """deterministic volume-equidistributed cloud: jittered ~square grid on [0,1]^2."""
    m = int(round(np.sqrt(N)))
    xs = (np.arange(m) + 0.5) / m
    X, Y = np.meshgrid(xs, xs)
    P = np.column_stack([X.ravel(), Y.ravel()])
    P += rng.uniform(-0.4, 0.4, P.shape) / m
    return np.clip(P, 1e-6, 1 - 1e-6)


def void_cloud(N, r, rng):
    """jittered ~square grid on [0,1]^2 with a central disk of radius r emptied -- a
    non-equidistributed (count-deficient-then-rebalanced) cloud; the void sits in a low-f
    region, exactly the matched-count near-null void/sliver family of the crux's L1/L3."""
    m = int(round(np.sqrt(N)))
    xs = (np.arange(m) + 0.5) / m
    X, Y = np.meshgrid(xs, xs)
    P = np.column_stack([X.ravel(), Y.ravel()])
    P += rng.uniform(-0.4, 0.4, P.shape) / m
    d = (P[:, 0] - 0.5) ** 2 + (P[:, 1] - 0.5) ** 2
    return np.clip(P[d > r * r], 1e-6, 1 - 1e-6)


def analytic_void_bias(r):
    """c_u limit for uniform nodes on [0,1]^2 minus a central disk: <f>_{K\\void} / <f>_K."""
    g = np.linspace(0, 1, 3000)
    X, Y = np.meshgrid(g, g)
    f = gradf(np.column_stack([X.ravel(), Y.ravel()])).reshape(X.shape)
    void = (X - 0.5) ** 2 + (Y - 0.5) ** 2 <= r * r
    return f[~void].mean() / f.mean()


def main():
    print("=" * 100)
    print("  Computation 130 -- A6 c_u is a cell-volume equidistribution defect; the A1 count does not force it")
    print("=" * 100)
    print()
    rng = np.random.default_rng(20260628)

    print("  (1) VOLUME-EQUIDISTRIBUTED cloud (jittered grid): c_u -> 1 as the count grows")
    print(f"      {'N':>7} {'|c_u-1|':>10}")
    d_uni = []
    for N in (400, 1600, 6400, 25600):
        d = abs(c_u(uniform_cloud(N, rng)) - 1.0)
        d_uni.append(d)
        print(f"      {N:>7} {d:>10.5f}")
    p1 = d_uni[-1] < 0.01 and d_uni[-1] < d_uni[0]
    print(f"      -> defect decays toward 0 (equidistribution closes it): {p1}")
    print()

    print("  (2) NON-EQUIDISTRIBUTED (central void r=0.2, low-f region emptied): c_u biased, FLAT in N")
    print(f"      {'N':>7} {'pts':>7} {'|c_u-1|':>10} {'analytic bias':>14}")
    bias = abs(analytic_void_bias(0.2) - 1.0)
    d_void = []
    for N in (1600, 6400, 25600, 102400):
        P = void_cloud(N, 0.2, rng)
        d = abs(c_u(P) - 1.0)
        d_void.append(d)
        print(f"      {N:>7} {len(P):>7} {d:>10.5f} {bias:>14.5f}")
    p2 = (d_void[-1] > 0.03) and (abs(d_void[-1] - bias) < 0.01)
    print(f"      -> defect FLAT in N (refinement does NOT close it) and matches the analytic void bias {bias:.4f}: {p2}")
    print(f"         the cloud has the right asymptotic density everywhere outside the void; only equidistribution fails.")
    print()

    print("  (3) defect grows with the void radius (cell-volume equidistribution is the binding property)")
    print(f"      {'r':>6} {'|c_u-1| (N=25600)':>18} {'analytic':>10}")
    ds = []
    for r in (0.0, 0.1, 0.15, 0.2, 0.25):
        d = abs(c_u(void_cloud(25600, r, rng)) - 1.0)
        ds.append(d)
        print(f"      {r:>6.2f} {d:>18.5f} {abs(analytic_void_bias(r)-1):>10.5f}")
    p3 = all(ds[i] <= ds[i + 1] + 0.004 for i in range(len(ds) - 1)) and ds[-1] > ds[0] + 0.02
    print(f"      -> monotone in the void radius (non-equidistribution amplitude): {p3}")
    print()

    print("=" * 100)
    print("  ASSESSMENT")
    print("=" * 100)
    print(f"  (1) equidistributed cloud: c_u -> 1                              : {p1}")
    print(f"  (2) matched COUNT but non-equidistributed: c_u biased, flat in N : {p2}")
    print(f"  (3) defect set by density non-uniformity (not count/shape)       : {p3}")
    ok = p1 and p2 and p3
    print()
    msg = ("CONFIRMS the cycle-6 correction -- c_u is a cell-volume equidistribution defect; the matched-scaling "
           "A1 count does NOT force c_u -> 1, so the A6 bias half is conditional on a deterministic "
           "equidistribution lemma (the sibling of the open A2 weak-discrepancy), not on the point count") if ok else "MISMATCH -- revise"
    print("  RESULT: " + msg)


if __name__ == "__main__":
    main()