# SilverSight NUVMAP — Julia Cross-Validation Test
#
# Runs the same test data as tests/test_nuvmap_equivalence.py and
# tests/test_nuvmap_r_ref.R, then asserts the Julia Q16_16 output
# matches Python/R within Q16_16 quantization tolerance.

include("/home/allaun/SilverSight/julia/CoreFormalism/Q16_16.jl")
include("/home/allaun/SilverSight/julia/nuvmap/projection_engine.jl")

using .Q16_16
using .NUVMAP
using Printf

# ── Test data ───────────────────────────────────────────────────────

test_data = [
    Dict{Symbol, Any}(
        :equation_id => 1, :amvr => 0.8, :avmr => 0.75, :cr => 0.05, :cs => "achiral_stable"
    ),
    Dict{Symbol, Any}(
        :equation_id => 2, :amvr => 0.3, :avmr => 0.4, :cr => 0.2, :cs => "left_handed_mass_bias"
    ),
    Dict{Symbol, Any}(
        :equation_id => 3, :amvr => 0.1, :avmr => 0.15, :cr => 0.6, :cs => "chiral_scarred"
    ),
    Dict{Symbol, Any}(
        :equation_id => 4, :amvr => 0.5, :avmr => 0.5, :cr => 0.1, :cs => "right_handed_vector_bias"
    ),
    Dict{Symbol, Any}(
        :equation_id => 5, :amvr => 0.9, :avmr => 0.85, :cr => 0.02, :cs => "achiral_stable"
    ),
]

function convert_to_q16(data)
    return [Dict{Symbol, Any}(
        :equation_id => d[:equation_id],
        # Use round() to match Python f2q (banker's rounding) exactly.
        # Julia of_float uses floor (matching Lean ofFloat), but
        # the test data was generated by Python f2q's round().
        :amvr_q16 => Q16_16.of_raw_int(Int32(round(d[:amvr] * 65536.0))),
        :avmr_q16 => Q16_16.of_raw_int(Int32(round(d[:avmr] * 65536.0))),
        :chiral_residual_q16 => Q16_16.of_raw_int(Int32(round(d[:cr] * 65536.0))),
        :chiral_state => d[:cs],
    ) for d in data]
end

# ── Known Python Q16_16 outputs (from test run) ─────────────────────

# From tests/test_nuvmap_equivalence.py output:
#   Cell 0: q_i=73,  E_i=7.748764, R_i=0.114288
#   Cell 1: q_i=0,   E_i=0.466660, R_i=0.599991
#   Cell 2: q_i=0,   E_i=0.033340, R_i=1.285721
#   Cell 3: q_i=2,   E_i=0.933289, R_i=0.428574
#   Cell 4: q_i=10780, E_i=99.902435, R_i=0.009995
#   total_qubits=10855, B=21.816895

expected_q_i    = [73, 0, 0, 2, 10780]
expected_E_i    = [7.748764, 0.466660, 0.033340, 0.933289, 99.902435]
expected_R_i    = [0.114288, 0.599991, 1.285721, 0.428574, 0.009995]
expected_qubits = 10855
expected_B     = 21.816895

# ── Run Julia projector ─────────────────────────────────────────────

println("="^60)
println("NUVMAP Julia — Cross-Validation Test")
println("="^60)

data_q16 = convert_to_q16(test_data)

engine = NUVMAPProjectionEngine(
    total_qubit_budget = 0,
    chi_max_q16 = Q16_16.half,
    R_max_q16 = Q16_16.half,
    # Integer division on raw Q16_16 value: 65536 ÷ 10 = 6553 = 0.1
    landauer_threshold_q16 = Q16_16.of_raw_int(Q16_16.to_int(Q16_16.one) ÷ 10)
)

surface = NUVMAP.project(engine, data_q16)

println("\nResults:")
println("  Cells: $(length(surface.cells))")
println("  Total qubits: $(surface.total_qubits)  (expected $expected_qubits)")
println("  Bekenstein: $(Q16_16.to_float(surface.bekenstein_bound))  (expected $expected_B)")

# ── Compare ─────────────────────────────────────────────────────────

println("\nCell-by-cell comparison:")
println("  cell  q_i       E_i           R_i")
println("  ----  -------   -----------   -----------")
all_ok = true

for i in 1:length(surface.cells)
    c = surface.cells[i]

    q_ok = c.q_i == expected_q_i[i]
    e_diff = abs(Q16_16.to_float(c.E_i) - expected_E_i[i])
    e_ok = e_diff < 0.02
    r_diff = abs(Q16_16.to_float(c.R_i) - expected_R_i[i])
    r_ok = r_diff < 0.0001

    status = (q_ok && e_ok && r_ok) ? "✓" : "✗"
    println("  [$(@sprintf("%d", i-1))] q=$(c.q_i)$(@sprintf("%+d", c.q_i - expected_q_i[i]))  E=$(Q16_16.to_float(c.E_i))  R=$(Q16_16.to_float(c.R_i))  $status")

    if !(q_ok && e_ok && r_ok)
        println("        expected: q=$(@sprintf("%d", expected_q_i[i])) E=$(expected_E_i[i]) R=$(expected_R_i[i])")
        all_ok = false
    end
end

# Surface-level checks
qt_diff = abs(surface.total_qubits - expected_qubits)
if qt_diff > 1
    println("\n✗ total_qubits: $(surface.total_qubits) vs expected $expected_qubits")
    all_ok = false
end

b_diff = abs(Q16_16.to_float(surface.bekenstein_bound) - expected_B)
if b_diff > 0.02
    println("✗ bekenstein_bound: $(Q16_16.to_float(surface.bekenstein_bound)) vs expected $expected_B")
    all_ok = false
end

println("\n" * "="^60)
if all_ok
    println("JULIA VALIDATION: PASS — matches Python/R Q16_16 output ✓")
    exit(0)
else
    println("JULIA VALIDATION: FAIL — divergence detected ✗")
    exit(1)
end