#!/usr/bin/env python3
"""Analyze three exact upstream GPI exposures locally; no observed pixels are redistributed.

NumPy only. This is a sparse demonstration, not a paper-spectrum replication.
The source pixels and locally derived image figures stay in the reader's folder.
"""
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,
    real_comparison,
    real_response,
)
from imaging_reduction import fit_spectrum
from gpi_profile import position_profile


def json_default(value):
    if isinstance(value, np.ndarray):
        return value.tolist()
    if isinstance(value, np.generic):
        return value.item()
    raise TypeError(str(type(value)))


def main(manifest, cache, output, algorithm="adi", local_images=False):
    dataset = load_selected(manifest, cache)
    data, centers = registered_science(dataset)
    images = reduce_registered(data, dataset, algorithm)
    # Published historical astrometry defines a disclosed LOCAL region of
    # interest; this is not a blind discovery or a prior-independent detection.
    separation, pa = 430.0, 211.6
    observed = sky_measurements(images, separation, pa)
    mean, covariance, comparison_angles = real_comparison(images, separation, pa)
    response = real_response(data, dataset, separation, pa, separation, pa, algorithm)
    fit = fit_spectrum(observed - mean, response, covariance)
    calibration_fraction = 0.4 * math.log(10) * 0.15
    calibration_covariance = (
        np.outer(fit["contrast"], fit["contrast"]) * calibration_fraction**2
    )
    variants = []
    for variant in range(4):
        matrix = real_response(
            data, dataset, separation, pa, separation, pa, algorithm, variant=variant
        )
        try:
            variants.append(
                fit_spectrum(observed - mean, matrix, covariance)["contrast"]
            )
        except ValueError:
            variants.append(None)
    try:
        profile_status = position_profile(data, dataset, images, separation, pa, algorithm)
    except ValueError as error:
        profile_status = {"status": "not-identifiable", "reason": str(error), "points": []}
    # This observes the actual maximum in a disclosed annulus using CD-derived
    # north-up/east-left coordinates; publication PA is an independent comparison.
    white = np.mean(
        images / dataset["star_aperture_adu"].T.mean(axis=1)[:, None, None], axis=0
    )
    rows, columns = np.indices(white.shape)
    radius = np.hypot(columns - 55, rows - 55)
    masked = np.where((radius > 25) & (radius < 36), white, -np.inf)
    row, column = np.unravel_index(np.nanargmax(masked), masked.shape)
    east, north = -(column - 55) * 14.14, (55 - row) * 14.14
    peak_pa = math.degrees(math.atan2(east, north)) % 360
    result = {
        "scope": "Three actual reduced GPI exposures and ten selected native channels; historical-ROI conditional extraction",
        "algorithm": algorithm,
        "channels_zero_based": dataset["manifest"]["channels"],
        "wavelengths_micron": dataset["manifest"]["files"][0]["observed"][
            "selected_wavelengths_micron"
        ],
        "native_pixel_units": "ADU per coadd",
        "contrast": fit["contrast"],
        "conditional_noise_covariance": fit["covariance"],
        "common_spot_calibration_covariance": calibration_covariance,
        "spot_to_star_delta_mag": 9.36,
        "spot_calibration_sigma_mag": 0.15,
        "comparison_count": len(comparison_angles),
        "comparison_angles_degrees": comparison_angles,
        "response": response,
        "response_singular_values": np.linalg.svd(response, compute_uv=False),
        "individual_spot_template_sensitivity": variants,
        "spots": dataset["spots"],
        "remeasured_centers_native": centers,
        "position_spectrum_profile": profile_status,
        "orientation_check": {
            "annular_peak_pixel": [int(column), int(row)],
            "header_cd_peak_pa_degrees": peak_pa,
            "header_cd_peak_separation_mas": math.hypot(east, north),
            "historical_published_pa_degrees": 211.6,
            "historical_published_sigma_pa_degrees": 1.3,
            "interpretation": "Image peak comparison checks orientation, not an independent high-precision astrometric calibration.",
        },
        "limits": [
            "No exact nineteen-exposure membership or published 2014 reduction reproduced; supplied files were reduced in 2016.",
            "No ERR extension or calibrated gain: no invented Poisson electron error for reduced ADU.",
            "A constant historical J-band spot ratio is adopted across selected reduced channels; no precision chromatic calibration or uncertain bare-table off-axis throughput correction is applied.",
            "Observed-pixel and derived-image redistribution remains unresolved; this JSON stores numerical facts, not images.",
            "Local Gaussian-plus-plane registration and finite satellite templates are approximations; individual-template sensitivity is retained.",
            "The common 0.15-mag spot calibration is historical and not an independent per-channel noise term.",
            "AVPARANG/CD conventions are preserved. Original repeated SATS and PSFCENT headers are never overwritten.",
            "Within-exposure field smear and upstream lenslet spectral covariance are not fully modeled.",
            "An expanded and refined local position profile with fixed empirical covariance is not a posterior, a calibrated interval, or a global false-alarm probability.",
        ],
    }
    output.mkdir(parents=True, exist_ok=True)
    (output / ("gpi-" + algorithm + "-results.json")).write_text(
        json.dumps(result, default=json_default, indent=2) + "\n"
    )
    if local_images:
        np.savez_compressed(
            output / ("LOCAL-ONLY-gpi-" + algorithm + "-images.npz"),
            raw_registered=data,
            residuals=images,
        )
    print(
        "COMPLETE sparse real extraction",
        algorithm,
        "contrast",
        fit["contrast"],
        "orientation PA",
        peak_pa,
        "profile",
        profile_status["status"],
        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("local-results"))
    parser.add_argument("--algorithm", choices=["adi", "asdi"], default="adi")
    parser.add_argument(
        "--local-images",
        action="store_true",
        help="Write local derived pixels; these are not website assets",
    )
    args = parser.parse_args()
    main(args.manifest, args.cache, args.output, args.algorithm, args.local_images)