#!/usr/bin/env python3 """Independent complete noisy acquisitions, re-estimated covariance and position profile. No chi-square tuning: raw trials, boundary minima and failures are retained. A synthetic instrument with fixed prescribed OPD is the population here, not real atmospheric speckles or a claim of survey false-alarm calibration. """ import argparse import json import math from pathlib import Path import numpy as np from imaging_reduction import ( reduce_cube, background_subtracted_weights, spatial_weights, measurements, comparison_covariance, finite_response, fit_spectrum, ) from synthetic_optics import source_templates, stellar_expectation from noise_covariance import unit_source_covariances, bootstrap_source_covariance def serializable(value): if isinstance(value, np.ndarray): return value.tolist() if isinstance(value, np.generic): return value.item() raise TypeError(type(value).__name__) def validate(materials, output, draws=64, algorithm="adi", position_half_width=0.002): metadata = json.loads((materials / "synthetic-input.json").read_text()) arrays = np.load(materials / "synthetic-observation.npz", allow_pickle=False) settings, angles = metadata["settings"], np.asarray(metadata["angles"]) wavelengths = metadata["camera"]["wavelengthsMicrons"] templates = arrays["source_templates"] truth = settings["contrast"] * np.asarray([1, 0.9, 0.62, 0.48, 0.76]) angle = math.radians(settings["positionAngleDegrees"]) x, y = settings["separation"] * math.cos(angle), settings["separation"] * math.sin( angle ) weights = spatial_weights(80, x, y) calibration_weights = background_subtracted_weights(80, settings, True) expected = ( stellar_expectation(settings, angles, wavelengths) + 3e7 * truth[:, None, None, None] * templates ) unit_cov = ( unit_source_covariances(templates, calibration_weights, angles, wavelengths, algorithm) if algorithm != "klip" else None ) offsets = ( [ -position_half_width, -position_half_width / 2, 0, position_half_width / 2, position_half_width, ] if algorithm != "klip" else [-position_half_width, 0, position_half_width] ) grid = [] for dx in offsets: for dy in offsets: trial_templates = ( templates if dx == dy == 0 else source_templates(settings, angles, wavelengths, x + dx, y + dy) ) response = finite_response( arrays["data"], trial_templates, angles, wavelengths, calibration_weights, algorithm ) grid.append( {"dx": dx, "dy": dy, "response": response, "templates": trial_templates} ) center_index = next( index for index, point in enumerate(grid) if point["dx"] == point["dy"] == 0 ) trials, raw_vectors = [], [] for trial in range(draws): rng = np.random.default_rng(22022000 + trial) data = rng.poisson(expected) + rng.normal(0, 2, expected.shape) images = reduce_cube(data, angles, wavelengths, algorithm) mean, comparison, count = comparison_covariance(images, settings, spatial=True) observed = measurements(images, weights) - mean raw_vectors.append(observed) response = ( grid[center_index]["response"] if algorithm != "klip" else finite_response( data, templates, angles, wavelengths, calibration_weights, algorithm ) ) preliminary = fit_spectrum(observed, response, comparison) if unit_cov is not None: photon_cov = np.einsum( "c,cij->ij", np.maximum(preliminary["contrast"], 0), unit_cov ) else: photon_cov = bootstrap_source_covariance( data, templates, preliminary["contrast"], calibration_weights, angles, wavelengths, np.random.default_rng(22023000 + trial), ) covariance = comparison + photon_cov points = [] for point in grid: this_response = ( point["response"] if algorithm != "klip" else finite_response( data, point["templates"], angles, wavelengths, calibration_weights, algorithm ) ) fit = fit_spectrum(observed, this_response, covariance) points.append(fit) best_index = min( range(len(points)), key=lambda index: points[index]["negative2loglikelihood"], ) center, best = points[center_index], points[best_index] diagonal_errors = np.sqrt(np.diag(center["covariance"])) row = { "seed": 22022000 + trial, "contrast": center["contrast"], "standard_errors": diagonal_errors, "covered_1sigma": np.abs(center["contrast"] - truth) <= diagonal_errors, "covered_2sigma": np.abs(center["contrast"] - truth) <= 2 * diagonal_errors, "center_chi_square": center["chi_square"], "profile_delta_truth": center["negative2loglikelihood"] - best["negative2loglikelihood"], "best_offset": [grid[best_index]["dx"], grid[best_index]["dy"]], "boundary_minimum": abs(grid[best_index]["dx"]) == max(offsets) or abs(grid[best_index]["dy"]) == max(offsets), "comparison_rank": int(np.linalg.matrix_rank(comparison)), "comparison_condition": float(np.linalg.cond(comparison)), "comparison_count": count, "full_covariance_condition": center["condition"], "photon_covariance_representatives": [ photon_cov[0, 0], photon_cov[0, 1], photon_cov[0, 3], photon_cov[6, 7], ], "covariance_representatives": [ covariance[0, 0], covariance[0, 1], covariance[0, 3], covariance[6, 7], ], } if algorithm == "klip" and trial < 4: fixed = bootstrap_source_covariance( data, templates, preliminary["contrast"], calibration_weights, angles, wavelengths, np.random.default_rng(22023000 + trial), fixed_basis=True, ) row["fixed_vs_reestimated_basis_diagonal_ratio"] = np.diag(fixed) / np.diag( photon_cov ) trials.append(row) print(algorithm, "complete acquisition", trial + 1, "/", draws, flush=True) estimates = np.asarray([row["contrast"] for row in trials]) result = { "algorithm": algorithm, "draws": draws, "truth": truth, "profile_grid_offsets": offsets, "coverage_1sigma_counts": np.sum( [row["covered_1sigma"] for row in trials], axis=0 ), "coverage_2sigma_counts": np.sum( [row["covered_2sigma"] for row in trials], axis=0 ), "estimate_mean": estimates.mean(axis=0), "estimate_empirical_covariance": np.cov(estimates, rowvar=False), "measurement_empirical_covariance": np.cov(raw_vectors, rowvar=False), "profile_delta_quantiles": np.quantile( [row["profile_delta_truth"] for row in trials], [0.5, 0.68, 0.95] ), "boundary_minima": sum(row["boundary_minimum"] for row in trials), "trials": trials, "interpretation": "Conditional Monte Carlo with new Poisson/read noise, unchanged prescribed OPD and fresh comparison covariance in every acquisition. Grid boundary cases and poor coverage are retained. Finite sample covariance, PSF mismatch, source-position grid and KLIP basis changes prevent an exact full-pipeline error claim.", } output.write_text(json.dumps(result, default=serializable, indent=2) + "\n") print( "COMPLETE", algorithm, result["coverage_1sigma_counts"], result["coverage_2sigma_counts"], "boundary", result["boundary_minima"], flush=True, ) if __name__ == "__main__": parser = argparse.ArgumentParser(description=__doc__) parser.add_argument("--materials", type=Path, default=Path(".")) parser.add_argument("--output", type=Path, default=Path("noise-coverage.json")) parser.add_argument("--draws", type=int, default=64) parser.add_argument("--position-half-width", type=float, default=0.002) parser.add_argument("--algorithm", choices=["adi", "asdi", "klip"], default="adi") args = parser.parse_args() validate( args.materials, args.output, args.draws, args.algorithm, args.position_half_width, )