"""Reintegrate exported native current-day spectra on an independent polar chart.

These are instantaneous Gaussian mixtures. The GL9 saved exposure observations and
inferred period/age are separate products; no period truth enters this check.
"""
from pathlib import Path
import argparse
import json
import numpy as np
from rotation_reference import disk_quadrature


def polar_line(row, order):
    x, y, z, quiet_weight = disk_quadrature(*order)
    settings = row['surfaceSettings']; time = row['timeDays']
    inclination = np.deg2rad(settings['spinInclinationDegrees'])
    latitude = np.deg2rad(settings['spotLatitudeDegrees'])
    shear = settings.get('differentialRotationShear', 0)
    phase = np.deg2rad(settings['spotLongitudeDegrees']) + 2*np.pi*time/settings['rotationDays']*(1-shear*np.sin(latitude)**2)
    direction = [np.cos(latitude)*np.sin(phase), np.sin(latitude)*np.sin(inclination)-np.cos(latitude)*np.cos(phase)*np.cos(inclination), -np.sin(latitude)*np.cos(inclination)-np.cos(latitude)*np.cos(phase)*np.sin(inclination)]
    amplitude = np.exp(-.5*((time-settings.get('spotEvolutionPeakDay', 0))/settings['spotEvolutionSigmaDays'])**2) if 'spotEvolutionSigmaDays' in settings else 1
    region = amplitude*np.exp(-(1-x*direction[0]-y*direction[1]-z*direction[2])/np.deg2rad(settings['spotWidthDegrees'])**2)
    active_weight = quiet_weight*(1-settings['spotContrast']*region)
    quiet_velocity = 2*np.pi*695700000/(settings['rotationDays']*86400)*np.sin(inclination)*x*(1-shear*(y*np.sin(inclination)-z*np.cos(inclination))**2)
    active_velocity = quiet_velocity+settings['convectiveSuppressionMps']*region
    indices = np.unique(np.linspace(0, len(row['line']['samples'])-1, 61).astype(int))
    velocities = np.array([row['line']['samples'][index]['velocityMps'] for index in indices])
    active = np.array([.6*np.sum(active_weight*np.exp(-.5*((velocity-active_velocity)/2500)**2))/np.sum(active_weight) for velocity in velocities])
    quiet = np.array([.6*np.sum(quiet_weight*np.exp(-.5*((velocity-quiet_velocity)/2500)**2))/np.sum(quiet_weight) for velocity in velocities])
    centroid = np.sum(active_weight*active_velocity)/np.sum(active_weight)
    variance = np.sum(active_weight*(active_velocity-centroid)**2)/np.sum(active_weight)
    moments = np.array([np.sum(active_weight)/np.sum(quiet_weight), centroid, variance])
    return indices, active, quiet, moments


def verify(native):
    assert native['schema'] == 'ac27-instantaneous-line-native-v1'
    assert native['savedExposureBoundary']['unchanged'] and native['savedExposureBoundary']['gaussNodes'] == 9
    maximum = np.zeros(6); results = []
    # For linear velocity deposition, |G''| <= 1/sigma^2 gives .6*h^2/(8*sigma^2).
    deposition_bound = .6*native['velocityStepMps']**2/(8*native['localSigmaMps']**2)
    for row in native['cases']:
        indices, active, quiet, moments = polar_line(row, (192, 960))
        _, coarse_active, coarse_quiet, coarse_moments = polar_line(row, (96, 480))
        native_active = np.array([row['line']['samples'][index]['deficit'] for index in indices])
        native_quiet = np.array([row['quiet']['samples'][index]['deficit'] for index in indices])
        point_error = max(np.max(np.abs(active-native_active)), np.max(np.abs(quiet-native_quiet)))
        quadrature_change = max(np.max(np.abs(active-coarse_active)), np.max(np.abs(quiet-coarse_quiet)))
        flux_error, centroid_error, variance_error = np.abs(moments-[row['surfaceFlux'], row['sameSurfaceCentroidMps'], row['sameSurfaceVarianceMpsSquared']])
        variance_relative = variance_error/max(1, moments[2])
        # Positive linear deposition adds at most h^2/4 to the velocity variance.
        variance_excess = row['discreteVarianceMpsSquared']-(2500**2+row['sameSurfaceVarianceMpsSquared'])
        line_centroid_error = abs(row['line']['centroidMps']-row['sameSurfaceCentroidMps'])
        assert point_error < deposition_bound+2e-6, (row['name'], point_error, deposition_bound)
        assert quadrature_change < 1e-6, (row['name'], quadrature_change)
        assert flux_error < 1e-6 and centroid_error < .001 and variance_relative < 1e-6, (row['name'], flux_error, centroid_error, variance_relative)
        assert -1e-6 <= variance_excess <= native['velocityStepMps']**2/4+1e-5, (row['name'], variance_excess)
        assert line_centroid_error < 1e-8, (row['name'], line_centroid_error)
        assert row['cellCount'] == row['surfaceOrder']**2
        for index, difference in enumerate(row['difference']):
            assert abs(difference-(row['line']['samples'][index]['deficit']-row['quiet']['samples'][index]['deficit'])) < 1e-15
        errors = np.array([point_error, quadrature_change, flux_error, centroid_error, variance_relative, line_centroid_error])
        maximum = np.maximum(maximum, errors)
        results.append({'name': row['name'], 'day': row['timeDays'], 'surface_order': row['surfaceOrder'],
                        'maximum_profile_absolute_error': float(point_error), 'independent_order_refinement': float(quadrature_change),
                        'positive_deposition_variance_excess_mps_squared': float(variance_excess),
                        'coarse_to_fine_centroid_change_mps': float(abs(moments[1]-coarse_moments[1])),
                        'difference_range': float(np.ptp(row['difference']))})
    by_name = {row['name']: row for row in native['cases']}
    def profile_distance(a, b):
        left, right = by_name[a]['line']['samples'], by_name[b]['line']['samples']
        assert [row['velocityMps'] for row in left] == [row['velocityMps'] for row in right]
        return max(abs(x['deficit']-y['deficit']) for x, y in zip(left, right))
    coupling = {'latitude': profile_distance('latitude-zero', 'latitude-north'),
                'shear': profile_distance('shear-zero', 'shear-strong'),
                'lifetime': profile_distance('lifetime-short', 'lifetime-long')}
    assert min(coupling.values()) > 1e-4, coupling
    assert np.ptp(by_name['default']['difference']) > .001
    assert np.ptp(by_name['pole']['difference']) < 1e-10
    assert np.ptp(by_name['faded']['difference']) < 1e-10
    # A pole-on convective perturbation need not have the quiet normalized shape.
    assert np.ptp(by_name['pole-convective']['difference']) > 1e-5
    return {'status': 'PASS actual exported native current samples against independent polar Gaussian mixtures and moments',
            'cases': results, 'maximum_errors_profile_refinement_flux_centroid_relativevariance_linecentroid': maximum.tolist(),
            'continuum_absolute_deposition_bound': deposition_bound, 'actual_profile_coupling': coupling,
            'reference': '61 actual native velocity bins per case; independent positive polar quadrature 96x480 to192x960; direct exp Gaussian mixture, no velocity deposition/convolution',
            'discretization': 'Native profile uses80m/s linear velocity deposition: pointwise interpolation bound and positive variance excess are checked separately from surface quadrature.',
            'boundary': 'Own instantaneous synthetic records only. Finite floating-point tests are not a formal continuum certificate, observed spectrum reduction, exposure validation or age-accuracy proof.'}


def main():
    parser = argparse.ArgumentParser()
    parser.add_argument('--native', type=Path, default=Path('instantaneous-line-injections.json'))
    parser.add_argument('--output', type=Path, default=Path('results/instantaneous-line-verification.json'))
    arguments = parser.parse_args()
    report = verify(json.loads(arguments.native.read_text()))
    arguments.output.parent.mkdir(parents=True, exist_ok=True)
    arguments.output.write_text(json.dumps(report, indent=2)+'\n')
    print(json.dumps({'status': report['status'], 'cases': len(report['cases']), 'maximum_errors': report['maximum_errors_profile_refinement_flux_centroid_relativevariance_linecentroid'], 'coupling': report['actual_profile_coupling']}))


if __name__ == '__main__': main()