"""Reproduce a frozen Perseus catalogue slice with Python's standard library."""
import argparse
from decimal import Decimal, getcontext
import hashlib
import json
from pathlib import Path
import sys
getcontext().prec = 50
D = Decimal


def calculate(data):
    if data['schemaVersion'] != 1 or data['source']['arxiv'] != 'astro-ph/0605323v2':
        raise ValueError('SOURCE_VERSION')
    if data['units']['pV'] != '1e58 erg' or data['units']['age'] != '1e7 yr' or data['units']['luminosity'] != '1e42 erg/s':
        raise ValueError('UNIT_MISMATCH')
    if data['aperture']['kind'] != 'detector-chip-edge-truncated':
        raise ValueError('APERTURE_FOOTNOTE_MISSING')
    factor = D(data['tutorial']['enthalpyFactor'])
    if factor != 4:
        raise ValueError('RELATIVISTIC_ENTHALPY_FACTOR')
    selected = [row for row in data['cavities'] if row['id'] in data['tutorial']['selectedCavityIds']]
    if len(selected) != 2 or any(row['generation'] != 'inner' for row in selected):
        raise ValueError('DO_NOT_MIX_CAVITY_GENERATIONS')
    for row in selected:
        if D(row['pV']) <= 0 or D(row['pV']) - D(row['pVMinus']) <= 0:
            raise ValueError('NONPOSITIVE_CAVITY_PV')
        if any(D(row[key]) <= 0 for key in ['soundAge', 'refillAge', 'buoyancyAge']):
            raise ValueError('NONPOSITIVE_AGE')
    aperture = data['aperture']
    loss = D(aperture['Lx']) - D(aperture['Lcool'])
    loss_low = D(aperture['Lx']) - D(aperture['LxMinus']) - D(aperture['Lcool']) - D(aperture['LcoolPlus'])
    loss_high = D(aperture['Lx']) + D(aperture['LxPlus']) - D(aperture['Lcool']) + D(aperture['LcoolMinus'])
    if loss_low <= 0:
        raise ValueError('NONPOSITIVE_HEATING_REQUIREMENT')
    energy_unit_j = D('1e58') * D(data['units']['ergJoules'])
    age_unit_s = D('1e7') * D(data['units']['yearSeconds'])
    luminosity_unit_w = D('1e42') * D(data['units']['ergJoules'])
    energy = factor * sum(D(row['pV']) for row in selected) * energy_unit_j
    energy_low = factor * sum(D(row['pV']) - D(row['pVMinus']) for row in selected) * energy_unit_j
    energy_high = factor * sum(D(row['pV']) + D(row['pVPlus']) for row in selected) * energy_unit_j
    rows = []
    for age in ['buoyancyAge', 'soundAge', 'refillAge']:
        individual = [factor * D(row['pV']) * energy_unit_j / (D(row[age]) * age_unit_s) for row in selected]
        power = sum(individual)
        power_low = sum(factor * (D(row['pV']) - D(row['pVMinus'])) * energy_unit_j / (D(row[age]) * age_unit_s) for row in selected)
        power_high = sum(factor * (D(row['pV']) + D(row['pVPlus'])) * energy_unit_j / (D(row[age]) * age_unit_s) for row in selected)
        rows.append({'ageColumn': age, 'individualPowerW': list(map(str, individual)),
            'pairConditionalPowerW': str(power), 'powerOverRequiredHeating': str(power / (loss * luminosity_unit_w)),
            'fixedAgeEndpointSensitivityRatio': [str(power_low / (loss_high * luminosity_unit_w)), str(power_high / (loss_low * luminosity_unit_w))],
            'effectiveEnergyWeightedAgeMyr': str(energy / power / D(data['units']['yearSeconds']) / D('1e6'))})
    published_power_w = D(data['publishedAggregateReferenceOnly']['cavityPower']) * luminosity_unit_w
    # This mismatch is deliberately retained: different selection/time averaging is not a rounding error.
    return {'status': 'PASS_FROZEN_INNER_PAIR_REPLAY_NOT_GLOBAL_BALANCE', 'sourceVersion': data['source']['arxiv'],
        'selectedIds': [row['id'] for row in selected], 'innerPairEnthalpyJ': str(energy),
        'pVEndpointEnthalpySensitivityJ': [str(energy_low), str(energy_high)],
        'requiredHeatingW': str(loss * luminosity_unit_w),
        'luminosityEndpointHeatingSensitivityW': [str(loss_low * luminosity_unit_w), str(loss_high * luminosity_unit_w)],
        'sameEnthalpyRadiativeReplacementTimeMyr': str(energy / (loss * luminosity_unit_w) / D(data['units']['yearSeconds']) / D('1e6')),
        'ageComparisons': rows, 'apertureRadiusKpc': aperture['radiusKpc'], 'apertureIsChipEdgeTruncated': True,
        'publishedAggregatePowerWReferenceOnly': str(published_power_w),
        'publishedAggregatePowerOverRequiredHeating': str(published_power_w / (loss * luminosity_unit_w)),
        'tutorialToPublishedAggregatePowerRatio': str(D(rows[0]['pairConditionalPowerW']) / published_power_w),
        'publishedAggregateEstimatorReproduced': False,
        'uncertaintyConvention': data['tutorial']['uncertaintyConvention'], 'boundary': data['tutorial']['scientificBoundary']}


def main():
    parser = argparse.ArgumentParser()
    parser.add_argument('--input', type=Path, default=Path(__file__).with_name('perseus-input.json'))
    parser.add_argument('--output', type=Path)
    args = parser.parse_args()
    raw = args.input.read_bytes()
    result = calculate(json.loads(raw))
    result['inputSHA256'] = hashlib.sha256(raw).hexdigest()
    result['scriptSHA256'] = hashlib.sha256(Path(__file__).read_bytes()).hexdigest()
    result['actualPython'] = sys.version
    output = json.dumps(result, indent=2, allow_nan=False) + '\n'
    if args.output:
        with args.output.open('x') as stream:
            stream.write(output)
    print(output, end='')


if __name__ == '__main__':
    main()