"""Automatic numerical scoring only. Professional/clinical score stays null."""

import hashlib
import json
import math
from pathlib import Path

import numpy as np

from .cases import FAMILIES
from .optics import metrics


def load_result(file, config):
    with np.load(file, allow_pickle=False) as z:
        required = {"psf", "throughput", "captured_energy", "angular_pixel_arcsec"}
        if not required.issubset(z.files):
            raise ValueError("missing required arrays/scalars")
        data = {key: np.array(z[key], dtype=float) for key in required}
    psf = data["psf"]
    if psf.shape != (config["output_n"], config["output_n"]):
        raise ValueError("wrong PSF shape")
    if not all(np.all(np.isfinite(value)) for value in data.values()):
        raise ValueError("non-finite submission")
    if np.any(psf < 0) or abs(float(psf.sum()) - 1.) > 1e-3:
        raise ValueError("PSF must be nonnegative and sum to one within 0.001")
    for key in required - {"psf"}:
        if data[key].shape != ():
            raise ValueError(f"{key} must be scalar")
        data[key] = float(data[key])
    if abs(data["angular_pixel_arcsec"] - config["angular_pixel_arcsec"]) > 1e-6:
        raise ValueError("wrong angular pixel scale")
    if not 0 <= data["throughput"] <= 1.00001:
        raise ValueError("passive pupil throughput outside [0,1]")
    # Only a small float serialization allowance; physical gain is disallowed.
    if not 0 < data["captured_energy"] <= data["throughput"] + 1e-5:
        raise ValueError("finite-FOV energy incompatible with passive throughput")
    return data


def compare(candidate, reference):
    pixel = reference["angular_pixel_arcsec"]
    cm, rm = metrics(candidate["psf"], pixel), metrics(reference["psf"], pixel)
    errors = {
        "psf_l1": float(np.abs(candidate["psf"] - reference["psf"]).sum()),
        "ee_mae": float(np.mean(np.abs(np.array(cm["encircled_energy"])-rm["encircled_energy"]))),
        "mtf_mae": float(np.mean(np.abs(np.array(cm["mtf_x"]+cm["mtf_y"])-np.array(rm["mtf_x"]+rm["mtf_y"])))),
        "centroid_arcsec": float(np.linalg.norm(np.array(cm["centroid_arcsec"])-rm["centroid_arcsec"])),
        "energy_abs": max(abs(candidate["throughput"]-reference["throughput"]), abs(candidate["captured_energy"]-reference["captured_energy"]))
    }
    scales = {"psf_l1": .08, "ee_mae": .03, "mtf_mae": .04, "centroid_arcsec": 8., "energy_abs": .02}
    weights = {"psf_l1": .45, "ee_mae": .15, "mtf_mae": .15, "centroid_arcsec": .1, "energy_abs": .15}
    score = 100 * sum(weights[key]*math.exp(-errors[key]/scales[key]) for key in errors)
    return score, errors


def evaluate(release, submission):
    release, submission = Path(release), Path(submission)
    public_bytes = (release / "public/cases.json").read_bytes()
    public = json.loads(public_bytes)
    private = json.loads((release / "private/manifest.json").read_text())
    if hashlib.sha256(public_bytes).hexdigest() != private["public_sha256"]:
        raise ValueError("public cases changed after reference generation")
    results, cached = [], {}
    for case in public["cases"]:
        cid = case["id"]
        reference = load_result(release / f"private/reference/{cid}.npz", case["config"])
        try:
            candidate = load_result(submission / f"{cid}.npz", case["config"])
        except (OSError, ValueError, TypeError, EOFError) as exc:
            results.append({"id": cid, "family": case["family"], "score": 0., "invalid": str(exc)})
            continue
        score, errors = compare(candidate, reference)
        cached[cid] = (candidate, reference)
        results.append({"id": cid, "family": case["family"], "score": score, "errors": errors})
    pair_results = []
    for pair in private["pairs"]:
        a, b = pair["control"], pair["intervention"]
        score, error = 0., None
        if a in cached and b in cached:
            ca, ra = cached[a]
            cb, rb = cached[b]
            # Unnormalized detector energies preserve changes in captured flux.
            expected = rb["psf"]*rb["captured_energy"] - ra["psf"]*ra["captured_energy"]
            observed = cb["psf"]*cb["captured_energy"] - ca["psf"]*ca["captured_energy"]
            error = float(np.abs(observed-expected).sum() / max(float(np.abs(expected).sum()), .002))
            score = 100 * math.exp(-error/.5)
        pair_results.append({**pair, "score": score, "relative_delta_l1": error})
    family_scores = {}
    for family in FAMILIES:
        singles = [r["score"] for r in results if r["family"] == family]
        paired = [r["score"] for r in pair_results if r["family"] == family]
        family_scores[family] = .7*float(np.mean(singles)) + .3*float(np.mean(paired))
    return {"version": "0.1.0", "track": "synthetic_thin_phase", "case_count": len(results),
            "core_score": float(np.mean(list(family_scores.values()))),
            "professional_score": None, "professional_status": "not_evaluated",
            "family_scores": family_scores, "invalid_count": sum("invalid" in r for r in results),
            "cases": results, "pairs": pair_results,
            "warning": "Development score against a public synthetic reference; not clinical validation or an Agent leaderboard result."}