"""Independent audit and sharpening of JAM 6's resonance-debt claim.

Reports a fully specified last-20-depth summary of cumulative R_D, the
final R_D, log-log slopes of diagonal injection energy and final
deviation energy against exact layer size, and R_D at completed
doubling-orbit endpoints.  The latter guards against mistaking a sweep
in modulus at fixed depth for the fixed-modulus depth limit.
"""

from __future__ import annotations

import sys
from pathlib import Path

import numpy as np

TOOLS = Path(r"C:\FOLDERFORAI\THE_WILDS\TOOLS")
sys.path.insert(0, str(TOOLS))

from jam2_rotation_protractor import counts, deviation, next_parts


MODULI = (9, 27, 81, 243, 729, 2187, 6561)
DEPTH = 63
TAIL = 20
FIT_FROM = 30


def audit(modulus: int):
    level = np.array([1], dtype=np.uint64)
    e0 = deviation(counts(level, modulus))
    e0_norm2 = float(np.dot(e0, e0))
    cumulative_a2 = 0.0
    cumulative_cross2 = 0.0
    rows = []

    for d in range(DEPTH):
        doubled, injected = next_parts(level)
        nxt = np.concatenate((doubled, injected))

        rotated = deviation(counts(doubled, modulus))
        arrival = deviation(counts(injected, modulus))
        next_error = deviation(counts(nxt, modulus))

        a2 = float(np.dot(arrival, arrival))
        cross2 = 2.0 * float(np.dot(rotated, arrival))
        e2 = float(np.dot(next_error, next_error))
        cumulative_a2 += a2
        cumulative_cross2 += cross2
        r_d = cumulative_cross2 / cumulative_a2 if cumulative_a2 else np.nan

        predicted_e2 = e0_norm2 + cumulative_a2 + cumulative_cross2
        assert np.isclose(predicted_e2, e2, rtol=1e-12, atol=1e-6)
        rows.append((d + 1, len(nxt), cumulative_a2, e2, r_d))
        level = nxt

    tail_r = np.array([r[-1] for r in rows[-TAIL:]])
    fit = [r for r in rows if r[0] >= FIT_FROM and r[2] > 0 and r[3] > 0]
    log_n = np.log([r[1] for r in fit])
    diagonal_slope = float(np.polyfit(log_n, np.log([r[2] for r in fit]), 1)[0])
    error_energy_slope = float(np.polyfit(log_n, np.log([r[3] for r in fit]), 1)[0])
    k = 0
    power = 1
    while power < modulus:
        power *= 3
        k += 1
    assert power == modulus
    orbit_order = 2 * (3 ** (k - 1))
    cycle_end_r = tuple(
        float(rows[d - 1][-1])
        for d in range(orbit_order, DEPTH + 1, orbit_order)
    )
    return (
        float(tail_r.mean()),
        float(tail_r.std()),
        float(tail_r[-1]),
        diagonal_slope,
        error_energy_slope,
        orbit_order,
        cycle_end_r,
    )


print(
    "q,orbit_order,completed_orbits,tail20_mean_R,tail20_std_R,"
    "final_R,diag_energy_slope,error_energy_slope,pressure_exponent"
)
for q in MODULI:
    result = audit(q)
    mean_r, std_r, final_r, diagonal, error, order, endpoints = result
    print(
        "%d,%d,%d,%+.9f,%.9f,%+.9f,%.9f,%.9f,%+.9f"
        % (
            q,
            order,
            len(endpoints),
            mean_r,
            std_r,
            final_r,
            diagonal,
            error,
            error - diagonal,
        )
    )
    if endpoints:
        print(
            "  R_D at completed orbit endpoints:",
            " ".join("%+.6f" % value for value in endpoints),
        )