#!/usr/bin/env python3
"""Independent empirical-PSF injections test the local Gaussian registration.

Images are downloaded/read locally; output contains positions and numeric
residual summaries only. Injecting onto a frozen real background is a
conditional algorithm test, not a new independent astronomical exposure.
"""
import argparse
import json
import math
from pathlib import Path
import numpy as np
from gpi_data import load_selected
from gpi_reduction import registered_science, reduce_registered, sky_measurements
from image_tools import insert_template, remeasure_spot


def validate(manifest, cache, output):
    dataset = load_selected(manifest, cache)
    trials = []
    for exposure in range(3):
        for channel in (0, 4, 9):
            image = dataset["science"][exposure, channel]
            star_center = np.mean(
                [spot["measured"] for spot in dataset["spots"][exposure][channel]],
                axis=0,
            )
            first = dataset["spots"][exposure][channel][0]
            reference_angle = math.atan2(
                first["measured"][1] - star_center[1],
                first["measured"][0] - star_center[0],
            )
            radius = float(np.linalg.norm(np.asarray(first["measured"]) - star_center))
            for location in range(4):
                angle = reference_angle + math.pi / 4 + location * math.pi / 2
                x, y = (
                    star_center[0] + radius * math.cos(angle) + 0.3,
                    star_center[1] + radius * math.sin(angle) - 0.2,
                )
                for spot in (0, 3):
                    for scale in (0.3, 1):
                        template = dataset["individual_templates"][
                            exposure, channel, spot
                        ]
                        amplitude = (
                            np.mean(
                                [
                                    item["aperture_adu"]
                                    for item in dataset["spots"][exposure][channel]
                                ]
                            )
                            * scale
                        )
                        added = insert_template(image.shape, template, x, y) * amplitude
                        trial = {
                            "exposure": exposure,
                            "channel_index": channel,
                            "original_channel": dataset["manifest"]["channels"][
                                channel
                            ],
                            "template_spot": spot,
                            "relative_spot_flux": scale,
                            "true_position": [x, y],
                        }
                        try:
                            fit = remeasure_spot(image + added, x, y)
                            trial.update(
                                status="accepted",
                                measured_position=fit,
                                bias_pixels=math.hypot(fit[0] - x, fit[1] - y),
                            )
                        except ValueError as error:
                            trial.update(status="failed-local-model", reason=str(error))
                        trials.append(trial)
    strict = {}
    data, _ = registered_science(dataset, strict_quality=True)
    try:
        images = reduce_registered(data, dataset)
        value = sky_measurements(images, 430, 211.6)
        strict.update(status="available", aperture_measurements_adu=value.tolist())
    except ValueError as error:
        strict.update(status="not-identifiable-with-strict-mask", reason=str(error))
    valid = [row["bias_pixels"] for row in trials if row["status"] == "accepted"]
    result = {
        "trials": trials,
        "accepted": len(valid),
        "failed": len(trials) - len(valid),
        "accepted_bias_median_pixels": float(np.median(valid)),
        "accepted_bias_95pct_pixels": float(np.quantile(valid, 0.95)),
        "accepted_bias_max_pixels": max(valid),
        "strict_any_nonzero_DQ_sensitivity": strict,
        "interpretation": "A local one-pixel Gaussian-plus-plane search can fail or be biased on a real structured halo. Failed boundary fits are retained. These numbers are conditional on injected empirical spot templates and frozen observed backgrounds, not astrometric error bars for the real planet.",
    }
    output.write_text(json.dumps(result, indent=2) + "\n")
    print(
        "COMPLETE registration injections",
        result["accepted"],
        "accepted",
        result["failed"],
        "failed; 95pct bias",
        result["accepted_bias_95pct_pixels"],
        flush=True,
    )


if __name__ == "__main__":
    parser = argparse.ArgumentParser(description=__doc__)
    parser.add_argument(
        "--manifest", type=Path, default=Path("data-source-manifest.json")
    )
    parser.add_argument("--cache", type=Path, default=Path("gpi-cache"))
    parser.add_argument("--output", type=Path, default=Path("gpi-registration.json"))
    args = parser.parse_args()
    validate(args.manifest, args.cache, args.output)