"""Standalone event likelihood; no Lab code, private inputs or runtime imports."""
from pathlib import Path
import ctypes
import hashlib
import json
import platform
import subprocess
import numpy as np
from scipy.interpolate import CubicHermiteSpline
from scipy.optimize import minimize_scalar


def load_observations(folder, manifest):
    tables = []; station_names = []; band_gammas = []
    for station, item in enumerate(manifest['files']):
        content = (folder / item['name']).read_bytes()
        assert len(content) == item['bytes'] and hashlib.sha256(content).hexdigest() == item['sha256']
        values = np.loadtxt(folder / item['name'])
        assert values.shape == (item['rows'], 3) and np.all(np.isfinite(values)) and np.all(values[:, 2] > 0)
        tables.append(np.column_stack((values, np.full(len(values), station))))
        station_names.append(item['name'])
        band_gammas.append(.626 if item['band'] == 'R' else .538)
    table = np.concatenate(tables)
    return {'time': table[:, 0], 'magnitude': table[:, 1], 'sigma': table[:, 2],
            'station': table[:, 3].astype(int), 'station_names': station_names,
            'gamma': np.array(band_gammas)[table[:, 3].astype(int)]}


def compile_reference(folder):
    reference = folder / 'reference'
    library = folder / ('event-reference.dylib' if platform.system() == 'Darwin' else 'event-reference.so')
    command = ['c++', '-O3', '-std=c++11', '-I' + str(reference),
               '-dynamiclib' if platform.system() == 'Darwin' else '-shared', '-fPIC',
               str(folder / 'event_reference.cpp'), str(reference / 'VBBinaryLensingLibrary.cpp'), '-o', str(library)]
    subprocess.run(command, check=True)
    return library, command


class EventReference:
    def __init__(self, library):
        self.library = ctypes.CDLL(str(library.resolve()))
        pointer = np.ctypeslib.ndpointer(dtype=np.float64, flags='C_CONTIGUOUS')
        self.library.ac23_magnifications.argtypes = [ctypes.c_int, ctypes.c_double, ctypes.c_double,
                                                     ctypes.c_double, pointer, pointer, pointer,
                                                     ctypes.c_double, pointer]
        self.library.ac23_magnifications.restype = ctypes.c_int

    def binary(self, source, separation, ratio, radius, gamma, tolerance=1e-5):
        source = np.asarray(source)
        x = np.ascontiguousarray(source.real); y = np.ascontiguousarray(source.imag)
        gamma = np.ascontiguousarray(gamma, dtype=np.float64)
        result = np.empty(len(source), dtype=np.float64)
        error = self.library.ac23_magnifications(len(source), separation, ratio, radius, x, y, gamma, tolerance, result)
        if error:
            raise ArithmeticError('Reference calculation failed at sample ' + str(error - 1))
        return result


def projected_earth(folder, times, reference_time=3582.731, stride=1, time_shift_days=0):
    text = (folder / 'earth-horizons.txt').read_text()
    rows = [line.split(',') for line in text.split('$$SOE')[1].split('$$EOE')[0].strip().splitlines()]
    ephemeris = np.array([[float(row[0]), *map(float, row[2:8])] for row in rows])[::stride]
    spline = CubicHermiteSpline(ephemeris[:, 0] - 2450000, ephemeris[:, 1:4], ephemeris[:, 4:7])
    shifted = np.asarray(times) + time_shift_days
    if np.min(shifted) < ephemeris[0, 0] - 2450000 or np.max(shifted) > ephemeris[-1, 0] - 2450000:
        raise ValueError('No Earth-ephemeris extrapolation')
    shifted_reference = reference_time + time_shift_days
    displacement = spline(shifted) - spline(shifted_reference) - (shifted - shifted_reference)[:, None] * spline(shifted_reference, 1)
    # Fixed J2000 source position from the discovery identification, explicit axes.
    right_ascension = np.deg2rad(15 * (17 + 54 / 60 + 19.2 / 3600))
    declination = np.deg2rad(-(30 + 22 / 60 + 38 / 3600))
    east = np.array([-np.sin(right_ascension), np.cos(right_ascension), 0])
    north = np.array([-np.sin(declination) * np.cos(right_ascension), -np.sin(declination) * np.sin(right_ascension), np.cos(declination)])
    return displacement @ north, displacement @ east


def source_trajectory(times, parameters, earth=None):
    t0, log_time, impact, separation, log_ratio, angle, log_radius = parameters[:7]
    tau = (np.asarray(times) - t0) / np.exp(log_time)
    beta = np.full_like(tau, impact)
    if earth is not None and len(parameters) >= 9:
        north, east = earth
        tau = tau + parameters[7] * north + parameters[8] * east
        beta = beta + parameters[7] * east - parameters[8] * north
    return np.exp(1j * angle) * (tau + 1j * beta)


def profile_magnitudes(observations, magnification):
    prediction = np.empty_like(magnification); stations = []
    for index, name in enumerate(observations['station_names']):
        selected = observations['station'] == index
        magnitude = observations['magnitude'][selected]
        weights = 1 / observations['sigma'][selected] ** 2
        extra = magnification[selected] - 1
        def evaluate(fraction):
            offset = 2.5 * np.log10(1 + fraction * extra)
            baseline = np.sum(weights * (magnitude + offset)) / np.sum(weights)
            residual = (baseline - offset - magnitude) * np.sqrt(weights)
            return float(np.dot(residual, residual)), float(baseline), residual
        interior = minimize_scalar(lambda value: evaluate(value)[0], bounds=(0, 1), method='bounded', options={'xatol': 1e-10})
        choices = [(evaluate(value), value) for value in (0., float(interior.x), 1.)]
        (chi_square, baseline, residual), fraction = min(choices, key=lambda row: row[0][0])
        prediction[selected] = observations['magnitude'][selected] + residual * observations['sigma'][selected]
        stations.append({'name': name, 'baseline_magnitude': baseline, 'source_baseline_fraction': fraction,
                         'chi_square': chi_square, 'samples': int(np.sum(selected)), 'blend_boundary': fraction == 0 or fraction == 1})
    residual = (prediction - observations['magnitude']) / observations['sigma']
    return residual, stations, prediction


def finite_point_lens(source, radius, gamma, radial_points=16, angular_points=32):
    nodes, weights = np.polynomial.legendre.leggauss(radial_points)
    mu = (nodes + 1) / 2
    radial = radius * np.sqrt(1 - mu ** 2)
    angles = (np.arange(angular_points) + .5) * 2 * np.pi / angular_points
    separations = abs(np.asarray(source)[:, None, None] + radial[None, :, None] * np.exp(1j * angles)[None, None, :])
    point = (separations ** 2 + 2) / (separations * np.sqrt(separations ** 2 + 4))
    intensity = 1 - np.asarray(gamma)[:, None] + 1.5 * np.asarray(gamma)[:, None] * mu[None, :]
    return np.sum(np.mean(point, axis=2) * (weights * mu)[None, :] * intensity, axis=1)