#!/usr/bin/env python3
"""crt_qprofile_sweep.py — q-profile safety factor design for CRT moduli.

Replaces brute-force modulus selection with systematic q-profile sweep,
as proposed in TOROIDAL_POLOIDAL_REFINEMENT.md (R1, R6).

The toroidal/poloidal mapping gives:
  q = L_reflection / L_identity (toroidal/poloidal ratio)
  q < 1 = poloidal-dominated (Sidon-favorable)
  q > 1 = toroidal-dominated (Sidon-unfavorable)
  q = 1 = degenerate (rational surface → resonance → Sidon collapse)

This script:
  1. Sweeps q from 0.3 to 3.0 (avoiding simple rationals)
  2. For each q, constructs moduli (L₀, L₁) with L₁/L₀ = q
  3. Tests Sidon preservation for known Sidon label sets
  4. Also applies cross-pair q-ratio check (R2): avoid simple rationals
  5. Measures Sidon score as function of q

Outputs:
  .openresearch/artifacts/crt_qprofile_sweep.json
  .openresearch/artifacts/EVAL.md

All integer arithmetic. No floats in compute path.
"""

import sys
import math
import json
import time
import hashlib
from pathlib import Path
from collections import Counter
from fractions import Fraction

REPO_ROOT = Path(__file__).resolve().parent.parent
ARTIFACTS_DIR = REPO_ROOT / ".openresearch" / "artifacts"
ARTIFACTS_DIR.mkdir(parents=True, exist_ok=True)
OUTPUT_PATH = ARTIFACTS_DIR / "crt_qprofile_sweep.json"
EVAL_PATH = ARTIFACTS_DIR / "EVAL.md"


# ── Exact Arithmetic ──────────────────────────────────────────────────────

def gcd(a, b):
    while b:
        a, b = b, a % b
    return a

def pairwise_coprime(moduli):
    for i in range(len(moduli)):
        for j in range(i + 1, len(moduli)):
            if gcd(moduli[i], moduli[j]) != 1:
                return False
    return True

def is_simple_rational(a, b, max_den=7):
    """Check if a/b reduces to m/n with n <= max_den."""
    if b == 0:
        return True
    g = gcd(abs(a), abs(b))
    na, nb = abs(a) // g, abs(b) // g
    return nb <= max_den


# ── CRT Torus Embedding ───────────────────────────────────────────────────

def embed(labels, S, moduli):
    M = 1
    for m in moduli:
        M *= m
    embedded = []
    for a in labels:
        row = []
        for i in range(len(moduli)):
            if i == 0:
                row.append(a % moduli[0])
            else:
                row.append((S - a) % moduli[i])
        embedded.append(row)
    return embedded, M

def egcd(a, b):
    if b == 0:
        return a, 1, 0
    g, x, y = egcd(b, a % b)
    return g, y, x - (a // b) * y

def modinv(a, m):
    g, x, _ = egcd(a % m, m)
    if g != 1:
        return None
    return x % m

def crt_reconstruct(residues, moduli):
    M = 1
    for m in moduli:
        M *= m
    x = 0
    for r, m in zip(residues, moduli):
        Mi = M // m
        inv = modinv(Mi % m, m)
        if inv is None:
            return None
        x = (x + r * Mi * inv) % M
    return x

def sidon_check(embedded, moduli):
    M = 1
    for m in moduli:
        M *= m
    n = len(embedded)
    total_pairs = n * (n + 1) // 2
    vals = [crt_reconstruct(row, moduli) for row in embedded]
    sums = []
    for i in range(n):
        for j in range(i, n):
            sums.append((vals[i] + vals[j]) % M)
    counts = Counter(sums)
    collisions = sum(c - 1 for c in counts.values())
    score = 1.0 - collisions / total_pairs if total_pairs > 0 else 1.0
    return {
        "total_pairs": total_pairs,
        "distinct_residues": len(counts),
        "collisions": collisions,
        "sidon_score": round(score, 6),
        "is_sidon": collisions == 0,
    }


# ── q-Profile Modulus Selection ──────────────────────────────────────────

def select_qprofile_moduli(L0, q_num, q_den, n_moduli=2):
    """Select moduli with a given q-profile.
    
    q = L₁/L₀ = q_num/q_den (as exact fraction)
    L₀ = L0 (identity axis, poloidal)
    L₁ = L₀ * q_num / q_den (reflection axis, toroidal)
    
    For n_moduli > 2, additional moduli use increasing primes coprime to all.
    Returns (moduli, q) or None if not coprime.
    """
    # L₁ = L0 * q_num / q_den — must be integer
    L1 = L0 * q_num // q_den
    if L0 * q_num != L1 * q_den:
        # Not exact — adjust L0 to make it work
        L0 = L0 * q_den
        L1 = L0 * q_num // q_den
    
    if L0 < 2 or L1 < 2:
        return None, None
    
    moduli = [L0, L1]
    
    # For n_moduli > 2, add coprime primes
    if n_moduli > 2:
        p = L1 + 2
        while len(moduli) < n_moduli:
            if all(gcd(p, m) == 1 for m in moduli) and p > 2:
                moduli.append(p)
            p += 1
    
    if not pairwise_coprime(moduli):
        return None, None
    
    q = Fraction(q_num, q_den)
    return moduli, q


def sweep_q_profile(label_set, S, L0_base, max_q_num=30, max_q_den=20):
    """Sweep q = L₁/L₀ over all coprime fractions q_num/q_den.
    Returns list of (q, moduli, sidon_result, q_is_simple_rational).
    """
    results = []
    
    for q_den in range(2, max_q_den + 1):
        for q_num in range(1, max_q_num + 1):
            if gcd(q_num, q_den) != 1:
                continue  # skip non-reduced fractions
            if q_num == q_den:
                continue  # skip q=1 (degenerate)
            
            moduli, q = select_qprofile_moduli(L0_base, q_num, q_den)
            if moduli is None:
                continue
            
            q_simple = is_simple_rational(q_num, q_den, max_den=7)
            
            embedded, M = embed(label_set, S, moduli)
            sidon = sidon_check(embedded, moduli)
            
            entry = {
                "q": str(q),
                "q_float": float(q),
                "q_num": q_num,
                "q_den": q_den,
                "q_simple_rational": q_simple,
                "moduli": moduli,
                "M": M,
                "L0": moduli[0],
                "L1": moduli[1],
                "sidon_score": sidon["sidon_score"],
                "collisions": sidon["collisions"],
                "is_sidon": sidon["is_sidon"],
                "total_pairs": sidon["total_pairs"],
                "distinct_residues": sidon["distinct_residues"],
            }
            results.append(entry)
    
    return results


# ── Main ──────────────────────────────────────────────────────────────────

def run_experiment():
    results = {
        "experiment": "crt_qprofile_sweep",
        "timestamp": time.strftime("%Y-%m-%dT%H:%M:%SZ", time.gmtime()),
        "config": {
            "L0_base": 7,
            "max_q_num": 30,
            "max_q_den": 20,
            "label_sets": [
                {"name": "sidon_pow2", "labels": [1, 2, 4, 8, 16], "S": 32},
                {"name": "sidon_singer5", "labels": [0, 1, 4, 14, 16], "S": 30},
                {"name": "nonsidon_seq5", "labels": [0, 1, 2, 3, 4], "S": 5},
            ],
        },
        "data": [],
        "summary": {},
    }
    
    label_sets = [
        ("sidon_pow2", [1, 2, 4, 8, 16], 32),
        ("sidon_singer5", [0, 1, 4, 14, 16], 30),
        ("nonsidon_seq5", [0, 1, 2, 3, 4], 5),
    ]
    
    L0_base = 7
    
    for desc, labels, S in label_sets:
        print(f"\n{'='*60}")
        print(f"  Label set: {desc} (n={len(labels)}, S={S})")
        print(f"{'='*60}")
        
        sweep = sweep_q_profile(labels, S, L0_base)
        
        # Analyze
        sidon_count = sum(1 for r in sweep if r["is_sidon"])
        total_count = len(sweep)
        q_simple_count = sum(1 for r in sweep if r["q_simple_rational"])
        q_simple_sidon = sum(1 for r in sweep if r["q_simple_rational"] and r["is_sidon"])
        q_not_simple_sidon = sum(1 for r in sweep if not r["q_simple_rational"] and r["is_sidon"])
        
        # q < 1 vs q > 1
        q_lt1 = [r for r in sweep if r["q_float"] < 1.0]
        q_gt1 = [r for r in sweep if r["q_float"] > 1.0]
        q_lt1_sidon = sum(1 for r in q_lt1 if r["is_sidon"])
        q_gt1_sidon = sum(1 for r in q_gt1 if r["is_sidon"])
        
        # Best q (highest Sidon score)
        best = max(sweep, key=lambda r: r["sidon_score"]) if sweep else None
        
        summary = {
            "total_configs": total_count,
            "sidon_configs": sidon_count,
            "sidon_rate": round(sidon_count / total_count, 4) if total_count > 0 else 0,
            "q_simple_rational_count": q_simple_count,
            "q_simple_sidon": q_simple_sidon,
            "q_not_simple_sidon": q_not_simple_sidon,
            "q_lt1_total": len(q_lt1),
            "q_lt1_sidon": q_lt1_sidon,
            "q_lt1_sidon_rate": round(q_lt1_sidon / len(q_lt1), 4) if q_lt1 else 0,
            "q_gt1_total": len(q_gt1),
            "q_gt1_sidon": q_gt1_sidon,
            "q_gt1_sidon_rate": round(q_gt1_sidon / len(q_gt1), 4) if q_gt1 else 0,
            "best_q": best["q"] if best else None,
            "best_sidon_score": best["sidon_score"] if best else None,
            "best_moduli": best["moduli"] if best else None,
        }
        results["summary"][desc] = summary
        
        print(f"  Total configs:  {total_count}")
        print(f"  Sidon configs:  {sidon_count} ({summary['sidon_rate']:.1%})")
        print(f"  q < 1: {q_lt1_sidon}/{len(q_lt1)} Sidon ({summary['q_lt1_sidon_rate']:.1%})")
        print(f"  q > 1: {q_gt1_sidon}/{len(q_gt1)} Sidon ({summary['q_gt1_sidon_rate']:.1%})")
        print(f"  q simple rational: {q_simple_sidon}/{q_simple_count} Sidon")
        print(f"  q NOT simple:      {q_not_simple_sidon}/{total_count - q_simple_count} Sidon")
        if best:
            print(f"  Best q = {best['q']} (score={best['sidon_score']}, moduli={best['moduli']})")
        
        # Store all sweep data
        for r in sweep:
            r["label_set"] = desc
            results["data"].append(r)
    
    content = json.dumps(results, indent=2, sort_keys=True, default=str)
    results["sha256"] = hashlib.sha256(content.encode()).hexdigest()
    return results


def write_eval(results):
    lines = [
        "# CRT q-Profile Safety Factor Sweep",
        "",
        f"**Experiment:** {results['experiment']}",
        f"**Date:** {results['timestamp']}",
        f"**SHA-256:** `{results['sha256']}`",
        "",
        "## Summary: Sidon Rate by q-Regime",
        "",
        "| Label set | Total | Sidon | Rate | q<1 total | q<1 Sidon | q<1 rate | q>1 total | q>1 Sidon | q>1 rate | Simple q Sidon | Non-simple Sidon | Best q |",
        "|-----------|-------|-------|------|-----------|-----------|----------|-----------|-----------|----------|----------------|------------------|-------|",
    ]
    for desc, s in results["summary"].items():
        lines.append(
            f"| {desc} | {s['total_configs']} | {s['sidon_configs']} | "
            f"{s['sidon_rate']:.1%} | "
            f"{s['q_lt1_total']} | {s['q_lt1_sidon']} | {s['q_lt1_sidon_rate']:.1%} | "
            f"{s['q_gt1_total']} | {s['q_gt1_sidon']} | {s['q_gt1_sidon_rate']:.1%} | "
            f"{s['q_simple_sidon']} | {s['q_not_simple_sidon']} | "
            f"{s['best_q']} |"
        )
    
    lines.extend([
        "",
        "## Predictions Tested",
        "",
        "1. **q < 1 (poloidal-dominated) should have higher Sidon rate than q > 1**",
        "   - This is the toroidal/poloidal refinement prediction",
        "   - If confirmed: poloidal resolution matters for Sidon structure",
        "",
        "2. **Simple rational q should have LOWER Sidon rate than non-simple q**",
        "   - This is the R2 cross-pair coprimality prediction",
        "   - Simple rationals = resonant surfaces = Sidon collapse",
        "",
        "3. **Best q should be < 1 and not a simple rational**",
        "   - Optimal q-profile is poloidal-dominated and irrational",
        "",
    ])
    
    EVAL_PATH.write_text("\n".join(lines))


if __name__ == "__main__":
    print("=" * 60)
    print("CRT q-Profile Safety Factor Sweep")
    print("=" * 60)
    
    t0 = time.time()
    results = run_experiment()
    elapsed = time.time() - t0
    
    OUTPUT_PATH.write_text(json.dumps(results, indent=2, default=str))
    print(f"\nResults → {OUTPUT_PATH}")
    
    write_eval(results)
    print(f"EVAL → {EVAL_PATH}")
    
    print(f"\nElapsed: {elapsed:.1f}s")
    print("=" * 60)
    print("DONE")
    print("=" * 60)