#!/usr/bin/env python3
"""
EMPIRICAL VALIDATION SUITE
===========================
The Davis Duality of Approximation and Obstruction
Davis Geometric — Bee Rosa Davis — April 2026

Tests INTRINSIC curvature (function complexity) not extrinsic embedding.

Key design principle: curvature = local nonlinearity of the TARGET FUNCTION.
The metric tensor comes from the function's Hessian, not from the surface shape.
ML sees flat coordinates (x,y,...) and tries to learn f(x,y,...).
Curvature measures how hard f is to learn locally.

NON-CIRCULARITY GUARANTEES:
  - Curvature computed from ANALYTIC FORMULAS on known functions
  - ML models trained on COORDINATES ONLY (never see curvature)
  - Curvature-error correlations computed POST HOC
  - Target functions are INDEPENDENT of curvature computation
  - Multiple architectures prevent overfitting to one model class
"""

import numpy as np
from scipy.stats import spearmanr, pearsonr
from sklearn.ensemble import GradientBoostingRegressor, RandomForestRegressor
from sklearn.neighbors import KNeighborsRegressor
from sklearn.neural_network import MLPRegressor
from sklearn.model_selection import cross_val_predict
from sklearn.metrics import r2_score, mean_squared_error
from sklearn.preprocessing import StandardScaler
import warnings
warnings.filterwarnings('ignore')
np.random.seed(42)

PASS = 0
FAIL = 0
RESULTS = {}

def report(test_num, name, passed, detail=""):
    global PASS, FAIL
    status = "✓ PASS" if passed else "✗ FAIL"
    if passed: PASS += 1
    else: FAIL += 1
    RESULTS[test_num] = {'name': name, 'passed': passed, 'detail': detail}
    print(f"\n{'='*70}")
    print(f"TEST {test_num}: {name}")
    print(f"  {status}")
    if detail:
        for line in detail.split('\n'):
            print(f"  {line}")
    print(f"{'='*70}")

def function_curvature(f_vals, x, y, h=0.05):
    """Estimate ||Hess f|| at each point via finite differences.
    This is the INTRINSIC curvature: how nonlinear f is locally.
    """
    n = len(x)
    hess_norm = np.zeros(n)
    for i in range(n):
        # Find nearby points for local Hessian estimation
        dx = x - x[i]
        dy = y - y[i]
        dist = np.sqrt(dx**2 + dy**2)
        mask = (dist > 0) & (dist < h * 3)
        if mask.sum() < 6:
            # Use wider neighborhood
            mask = (dist > 0) & (dist < h * 6)
        if mask.sum() < 6:
            hess_norm[i] = 0
            continue
        
        # Local quadratic fit: f ≈ a + bx + cy + dx² + exy + fy²
        local_dx = dx[mask]
        local_dy = dy[mask]
        local_df = f_vals[mask] - f_vals[i]
        
        A = np.column_stack([
            local_dx, local_dy,
            0.5*local_dx**2, local_dx*local_dy, 0.5*local_dy**2
        ])
        
        try:
            coeffs, _, _, _ = np.linalg.lstsq(A, local_df, rcond=None)
            # Hessian entries: [[d, e], [e, f]] from coefficients [b,c,d,e,f]
            H = np.array([[coeffs[2], coeffs[3]], [coeffs[3], coeffs[4]]])
            hess_norm[i] = np.linalg.norm(H, 'fro')
        except Exception:
            hess_norm[i] = 0
    
    return hess_norm

print("#"*70)
print("# EMPIRICAL VALIDATION SUITE")
print("# The Davis Duality of Approximation and Obstruction")
print("# Testing: do ML errors concentrate where the TARGET is locally complex?")
print("#"*70)

# ============================================================
# TEST 1: Function with spatially varying Hessian
# ============================================================
print("\n--- TEST 1: Errors track function Hessian norm ---")

n = 2000
x = np.random.uniform(-4, 4, n)
y = np.random.uniform(-4, 4, n)

# Function with HIGH curvature near origin, LOW curvature far out
# f = exp(-r²/σ²) * sin(ωr) — rapidly oscillating near center
r = np.sqrt(x**2 + y**2)
f_true = np.exp(-r**2/8) * np.sin(3*r) * np.cos(2*np.arctan2(y, x))
f_noisy = f_true + np.random.normal(0, 0.02, n)

# INTRINSIC curvature: ||Hess f|| computed from the function itself
K_intrinsic = function_curvature(f_true, x, y, h=0.3)

# ML: sees (x,y) only, tries to learn f
X = np.column_stack([x, y])
model = GradientBoostingRegressor(n_estimators=300, max_depth=5, random_state=42)
pred = cross_val_predict(model, X, f_noisy, cv=5)
errors = np.abs(f_noisy - pred)

# Remove points where curvature estimation failed
valid = K_intrinsic > 0
rho, p = spearmanr(K_intrinsic[valid], errors[valid])

# Quintile analysis
edges = np.percentile(K_intrinsic[valid], np.linspace(0, 100, 6))
bins = np.digitize(K_intrinsic[valid], edges[1:-1])
q_err = [np.mean(errors[valid][bins==q]) for q in range(5)]
ratio = q_err[-1] / q_err[0] if q_err[0] > 0 else 0

detail = f"ρ(||Hess f||, |error|) = {rho:.4f}, p = {p:.4e}\n"
detail += f"Q5/Q1 error ratio: {ratio:.2f}x\n"
detail += f"Quintile errors: {[f'{e:.4f}' for e in q_err]}"
report(1, "Errors track function Hessian norm", rho > 0.05 and p < 0.05, detail)

# ============================================================
# TEST 2: Smooth vs oscillatory regions — same feature space
# ============================================================
print("\n--- TEST 2: Smooth vs oscillatory regions ---")

n = 2000
x = np.random.uniform(-6, 6, n)
y = np.random.uniform(-6, 6, n)

# Left half: smooth. Right half: oscillatory.
# This creates a STEP FUNCTION in curvature at x=0.
f_smooth = 0.5 * np.sin(0.5*x) * np.cos(0.5*y)
f_oscillatory = 0.5 * np.sin(3*x) * np.cos(3*y)
f_true = np.where(x < 0, f_smooth, f_oscillatory)
f_noisy = f_true + np.random.normal(0, 0.02, n)

# Intrinsic curvature: smooth region has low Hessian, oscillatory has high
# Analytic: for sin(ωx)cos(ωy), ||Hess|| ~ ω²
K_analytic = np.where(x < 0, 0.25, 9.0)  # ω²: 0.5² vs 3²

X = np.column_stack([x, y])
model = GradientBoostingRegressor(n_estimators=300, max_depth=5, random_state=42)
pred = cross_val_predict(model, X, f_noisy, cv=5)
errors = np.abs(f_noisy - pred)

err_smooth = np.mean(errors[x < 0])
err_oscillatory = np.mean(errors[x >= 0])
ratio = err_oscillatory / err_smooth

rho, p = spearmanr(K_analytic, errors)
detail = f"Mean error in smooth region: {err_smooth:.4f}\n"
detail += f"Mean error in oscillatory region: {err_oscillatory:.4f}\n"
detail += f"Ratio (oscillatory/smooth): {ratio:.2f}x\n"
detail += f"ρ(K_analytic, |error|) = {rho:.4f}, p = {p:.4e}"
report(2, "Oscillatory regions have higher error", ratio > 1.3 and p < 0.01, detail)

# ============================================================
# TEST 3: Frequency gradient — continuously varying curvature
# ============================================================
print("\n--- TEST 3: Frequency gradient (continuous curvature variation) ---")

n = 2000
x = np.random.uniform(0, 10, n)
y = np.random.uniform(-3, 3, n)

# Frequency increases linearly with x: ω(x) = 0.5 + 0.3*x
# ||Hess f|| ~ ω(x)² — curvature increases quadratically
omega = 0.5 + 0.3 * x
f_true = np.sin(omega * x) * np.cos(omega * y) * np.exp(-y**2/8)
f_noisy = f_true + np.random.normal(0, 0.02, n)
K_analytic = omega**2  # ||Hess f|| ~ ω²

X = np.column_stack([x, y])
model = GradientBoostingRegressor(n_estimators=300, max_depth=5, random_state=42)
pred = cross_val_predict(model, X, f_noisy, cv=5)
errors = np.abs(f_noisy - pred)

rho, p = spearmanr(K_analytic, errors)
edges = np.percentile(K_analytic, np.linspace(0, 100, 6))
bins = np.digitize(K_analytic, edges[1:-1])
q_err = [np.mean(errors[bins==q]) for q in range(5)]
ratio = q_err[-1] / q_err[0] if q_err[0] > 0 else 0

detail = f"ρ(ω², |error|) = {rho:.4f}, p = {p:.4e}\n"
detail += f"Q5/Q1 error ratio: {ratio:.2f}x\n"
detail += f"Quintile errors: {[f'{e:.4f}' for e in q_err]}"
report(3, "Errors increase with local frequency (curvature)", rho > 0.1 and p < 0.01, detail)

# ============================================================
# TEST 4: Architecture independence on intrinsic curvature
# ============================================================
print("\n--- TEST 4: Architecture independence (intrinsic curvature) ---")

# Reuse the frequency gradient data from Test 3
X = np.column_stack([x, y])
arch_results = {}
for name, model in [
    ("GBT", GradientBoostingRegressor(n_estimators=200, max_depth=4, random_state=42)),
    ("RF", RandomForestRegressor(n_estimators=200, max_depth=6, random_state=42)),
    ("KNN", KNeighborsRegressor(n_neighbors=5)),
    ("MLP", MLPRegressor(hidden_layer_sizes=(128,64,32), max_iter=2000, random_state=42)),
]:
    pred = cross_val_predict(model, X, f_noisy, cv=5)
    errs = np.abs(f_noisy - pred)
    rho_a, p_a = spearmanr(K_analytic, errs)
    arch_results[name] = (rho_a, p_a)

detail = f"{'Arch':>6} {'ρ(ω²,|error|)':>14} {'p-value':>12}\n"
pos_count = 0
for nm, (rr, pp) in arch_results.items():
    sig = "***" if pp<.001 else ("**" if pp<.01 else ("*" if pp<.05 else "n.s."))
    detail += f"{nm:>6} {rr:>14.4f} {pp:>10.2e} {sig}\n"
    if rr > 0 and pp < 0.1:
        pos_count += 1

report(4, "Architecture independence (intrinsic curvature)", pos_count >= 3, detail)

# ============================================================
# TEST 5: Scaling with capacity at fixed curvature
# ============================================================
print("\n--- TEST 5: Error decreases with capacity, modulated by curvature ---")

n = 1500
x = np.random.uniform(-4, 4, n)
y = np.random.uniform(-4, 4, n)
omega_local = 0.5 + 0.4 * np.abs(x)
f_true = np.sin(omega_local * x) * np.cos(omega_local * y)
f_noisy = f_true + np.random.normal(0, 0.02, n)
K_local = omega_local**2
X = np.column_stack([x, y])

# Split into low-K and high-K regions
median_K = np.median(K_local)
low_K = K_local < median_K
high_K = K_local >= median_K

capacities = [50, 100, 200, 400, 800]
ratios_at_capacity = []

detail = f"{'Capacity':>10} {'err_lowK':>10} {'err_highK':>10} {'ratio':>8}\n"
for n_est in capacities:
    model = GradientBoostingRegressor(n_estimators=n_est, max_depth=4, random_state=42)
    pred = cross_val_predict(model, X, f_noisy, cv=5)
    errs = np.abs(f_noisy - pred)
    
    err_low = np.mean(errs[low_K])
    err_high = np.mean(errs[high_K])
    ratio_val = err_high / err_low if err_low > 0 else 0
    ratios_at_capacity.append(ratio_val)
    detail += f"{n_est:>10} {err_low:>10.4f} {err_high:>10.4f} {ratio_val:>8.2f}x\n"

# High-K region should consistently have higher error
avg_ratio = np.mean(ratios_at_capacity)
detail += f"Average high-K/low-K error ratio: {avg_ratio:.2f}x"
report(5, "High-curvature regions consistently worse across capacities", avg_ratio > 1.2, detail)

# ============================================================
# TEST 6: Wrong-target discrimination (intrinsic version)
# ============================================================
print("\n--- TEST 6: Wrong-target discrimination (intrinsic) ---")

n = 1000
x = np.random.uniform(-4, 4, n)
y = np.random.uniform(-4, 4, n)

# True target: a function with meaningful Hessian structure
f_true = np.sin(2*x) * np.exp(-y**2/4)
K_true = function_curvature(f_true, x, y, h=0.3)

# Wrong targets
wrong1 = np.random.normal(0, 0.5, n)
wrong2 = np.sin(2*np.pi*np.arange(n)/50)
wrong3 = np.cumsum(np.random.normal(0, 0.1, n))
wrong3 = (wrong3 - wrong3.mean()) / wrong3.std() * np.std(f_true)

valid = K_true > 0
r_true = abs(pearsonr(K_true[valid], np.abs(f_true[valid]))[0])
r_wrongs = [abs(pearsonr(K_true[valid], np.abs(w[valid]))[0]) for w in [wrong1, wrong2, wrong3]]
best_wrong = max(r_wrongs)
disc = r_true / (best_wrong + 1e-8)

detail = f"|r| with true target: {r_true:.4f}\n"
detail += f"|r| with best wrong target: {best_wrong:.4f}\n"
detail += f"Discrimination ratio: {disc:.2f}x"
report(6, "Curvature discriminates true vs wrong targets", disc > 1.5, detail)

# ============================================================
# TEST 7: Neural network (MLP) specifically — the theorem's actual subject
# ============================================================
print("\n--- TEST 7: Neural network errors track function curvature ---")

n = 2000
x = np.random.uniform(0, 10, n)
y = np.random.uniform(-3, 3, n)
omega = 0.5 + 0.3 * x
f_true = np.sin(omega * x) * np.cos(omega * y) * np.exp(-y**2/8)
f_noisy = f_true + np.random.normal(0, 0.02, n)
K_analytic = omega**2

X = StandardScaler().fit_transform(np.column_stack([x, y]))

model = MLPRegressor(hidden_layer_sizes=(256, 128, 64), max_iter=3000,
                     learning_rate_init=0.001, random_state=42)
pred = cross_val_predict(model, X, f_noisy, cv=5)
errors = np.abs(f_noisy - pred)

rho, p = spearmanr(K_analytic, errors)
edges = np.percentile(K_analytic, np.linspace(0, 100, 6))
bins = np.digitize(K_analytic, edges[1:-1])
q_err = [np.mean(errors[bins==q]) for q in range(5)]
ratio = q_err[-1] / q_err[0] if q_err[0] > 0 else 0

detail = f"ρ(ω², |error|) = {rho:.4f}, p = {p:.4e}\n"
detail += f"Q5/Q1 error ratio: {ratio:.2f}x\n"
detail += f"Quintile errors: {[f'{e:.4f}' for e in q_err]}"
report(7, "Neural network errors track function curvature", rho > 0.1 and p < 0.01, detail)

# ============================================================
# TEST 8: Regime transition — increase frequency until ML breaks
# ============================================================
print("\n--- TEST 8: Regime transition at curvature threshold ---")

n = 800
x = np.random.uniform(-3, 3, n)
y = np.random.uniform(-3, 3, n)

frequencies = [0.5, 1.0, 2.0, 3.0, 5.0, 8.0, 12.0, 16.0]
r2_by_freq = []

for omega in frequencies:
    f_true = np.sin(omega * x) * np.cos(omega * y)
    f_noisy = f_true + np.random.normal(0, 0.02, n)
    X = np.column_stack([x, y])
    model = GradientBoostingRegressor(n_estimators=200, max_depth=5, random_state=42)
    pred = cross_val_predict(model, X, f_noisy, cv=5)
    r2 = r2_score(f_noisy, pred)
    r2_by_freq.append(r2)

rho, p = spearmanr(frequencies, r2_by_freq)
detail = f"ρ(frequency, R²) = {rho:.4f}, p = {p:.4e}\n"
detail += f"{'ω':>6} {'κ~ω²':>8} {'R²':>8}\n"
for om, r2 in zip(frequencies, r2_by_freq):
    detail += f"{om:>6.1f} {om**2:>8.1f} {r2:>8.4f}\n"

report(8, "Performance degrades with increasing frequency", rho < -0.7 and p < 0.05, detail)

# ============================================================
# TEST 9: Coordinate invariance — same function, rotated coords
# ============================================================
print("\n--- TEST 9: Error is coordinate-invariant ---")

n = 1500
x = np.random.uniform(-3, 3, n)
y = np.random.uniform(-3, 3, n)
f_true = np.sin(2*x) * np.cos(3*y) + 0.5*x*y
f_noisy = f_true + np.random.normal(0, 0.03, n)

# Coordinate system 1: (x, y)
X1 = np.column_stack([x, y])
model = GradientBoostingRegressor(n_estimators=200, max_depth=4, random_state=42)
pred1 = cross_val_predict(model, X1, f_noisy, cv=5)
err1 = np.mean(np.abs(f_noisy - pred1))

# Coordinate system 2: rotated by 37 degrees
theta = 37 * np.pi / 180
x2 = x * np.cos(theta) + y * np.sin(theta)
y2 = -x * np.sin(theta) + y * np.cos(theta)
X2 = np.column_stack([x2, y2])
model2 = GradientBoostingRegressor(n_estimators=200, max_depth=4, random_state=42)
pred2 = cross_val_predict(model2, X2, f_noisy, cv=5)
err2 = np.mean(np.abs(f_noisy - pred2))

# Coordinate system 3: nonlinear (polar-ish)
x3 = np.sqrt(x**2 + y**2)
y3 = np.arctan2(y, x)
X3 = np.column_stack([x3, y3])
model3 = GradientBoostingRegressor(n_estimators=200, max_depth=4, random_state=42)
pred3 = cross_val_predict(model3, X3, f_noisy, cv=5)
err3 = np.mean(np.abs(f_noisy - pred3))

# Error should change with coordinates (ML is NOT coordinate-invariant)
# BUT the curvature-error correlation should be preserved
detail = f"Error in (x,y) coordinates: {err1:.4f}\n"
detail += f"Error in rotated coordinates: {err2:.4f}\n"
detail += f"Error in polar coordinates: {err3:.4f}\n"
detail += f"Ratio rotated/original: {err2/err1:.3f}\n"
detail += f"Ratio polar/original: {err3/err1:.3f}"

# Rotation should give similar error (isometry). Polar might differ.
rotated_close = abs(err2/err1 - 1.0) < 0.15
report(9, "Rotation invariance (isometry preserves error)", rotated_close, detail)

# ============================================================
# TEST 10: Multi-frequency decomposition
# ============================================================
print("\n--- TEST 10: Error scales with local spectral complexity ---")

n = 2000
x = np.random.uniform(-5, 5, n)
y = np.random.uniform(-5, 5, n)

# Compose function with 3 frequency bands, spatially separated
f_low = 0.5 * np.sin(0.5*x) * np.cos(0.5*y)   # ω=0.5, κ~0.25
f_mid = 0.5 * np.sin(2*x) * np.cos(2*y)        # ω=2, κ~4
f_high = 0.5 * np.sin(6*x) * np.cos(6*y)       # ω=6, κ~36

# Weight by spatial position: low-freq on left, high-freq on right
w_low = np.clip(1 - (x + 5)/5, 0, 1)
w_mid = np.clip(1 - np.abs(x)/3, 0, 1)
w_high = np.clip(1 - (5 - x)/5, 0, 1)
weights_sum = w_low + w_mid + w_high + 1e-8
w_low /= weights_sum; w_mid /= weights_sum; w_high /= weights_sum

f_true = w_low * f_low + w_mid * f_mid + w_high * f_high
f_noisy = f_true + np.random.normal(0, 0.02, n)

# Local spectral complexity ~ weighted frequency²
K_spectral = w_low * 0.25 + w_mid * 4.0 + w_high * 36.0

X = np.column_stack([x, y])
model = GradientBoostingRegressor(n_estimators=300, max_depth=5, random_state=42)
pred = cross_val_predict(model, X, f_noisy, cv=5)
errors = np.abs(f_noisy - pred)

rho, p = spearmanr(K_spectral, errors)
edges = np.percentile(K_spectral, np.linspace(0, 100, 6))
bins = np.digitize(K_spectral, edges[1:-1])
q_err = [np.mean(errors[bins==q]) for q in range(5)]
ratio = q_err[-1] / q_err[0] if q_err[0] > 0 else 0

detail = f"ρ(spectral complexity, |error|) = {rho:.4f}, p = {p:.4e}\n"
detail += f"Q5/Q1 error ratio: {ratio:.2f}x\n"
detail += f"Quintile errors: {[f'{e:.4f}' for e in q_err]}"
report(10, "Error scales with local spectral complexity", rho > 0.1 and p < 0.01, detail)

# ============================================================
# TEST 11: The curvature prefactor — high-K regions need more parameters
# ============================================================
print("\n--- TEST 11: High-K regions need more parameters to reach same error ---")

n = 1200
x = np.random.uniform(-4, 4, n)
y = np.random.uniform(-4, 4, n)

# Two regions: low-freq left, high-freq right
f_left = 0.5 * np.sin(0.7*x) * np.cos(0.7*y)
f_right = 0.5 * np.sin(4*x) * np.cos(4*y)
f_true = np.where(x < 0, f_left, f_right)
f_noisy = f_true + np.random.normal(0, 0.02, n)

left_mask = x < 0
right_mask = x >= 0

capacities = [30, 60, 120, 240, 480]
detail_lines = f"{'N_est':>6} {'err_lowK':>10} {'err_highK':>10} {'ratio':>8}\n"
ratios = []

for n_est in capacities:
    X = np.column_stack([x, y])
    model = GradientBoostingRegressor(n_estimators=n_est, max_depth=4, random_state=42)
    pred = cross_val_predict(model, X, f_noisy, cv=5)
    errs = np.abs(f_noisy - pred)
    
    el = np.mean(errs[left_mask])
    er = np.mean(errs[right_mask])
    ratio_val = er / el if el > 0 else 0
    ratios.append(ratio_val)
    detail_lines += f"{n_est:>6} {el:>10.4f} {er:>10.4f} {ratio_val:>8.2f}x\n"

avg_ratio = np.mean(ratios)
detail_lines += f"Average ratio: {avg_ratio:.2f}x"
report(11, "High-K regions consistently need more capacity", avg_ratio > 1.3, detail_lines)

# ============================================================
# FINAL SCORECARD
# ============================================================
print("\n" + "#"*70)
print("# FINAL SCORECARD")
print("#"*70)
print(f"\n  PASSED: {PASS}/11")
print(f"  FAILED: {FAIL}/11")
print(f"  PASS RATE: {100*PASS/11:.0f}%")

print(f"\n{'Test':>4} {'Result':>8} {'Name'}")
print("-"*70)
for t in sorted(RESULTS.keys()):
    r = RESULTS[t]
    status = "  ✓" if r['passed'] else "  ✗"
    print(f"  {t:>2}   {status:>6}   {r['name']}")

print(f"""
{'='*70}
NON-CIRCULARITY GUARANTEES:
  Curvature = function Hessian norm (INTRINSIC), not surface shape
  ML sees flat coordinates only — curvature is HIDDEN in the function
  Local frequency (omega) serves as analytically tractable curvature proxy
  Tests include actual neural networks (MLP), not just tree models
  Coordinate invariance tested explicitly
{'='*70}
Davis Geometric — C = tau/K
The Davis Duality of Approximation and Obstruction
{'='*70}
""")