# /// script
# requires-python = ">=3.12"
# dependencies = ["symbolica==3.0.2", "numpy==2.4.1"]
# ///
"""Validate a subset, then benchmark five native numerical routes on original F13.

uv run --python 3.12 scripts/benchmarks/f13-output.py --terms 256
uv run --python 3.12 scripts/benchmarks/f13-output.py --terms 0 \
    --output scripts/benchmarks/f13-numerical-results.json

Symbolica scalar x64 ASM and SymJIT, Mathematica Compile-to-C with expression
optimization off/on, FORM O4 followed by C++. C/C++ compilers use -O3 and
-march=native, without fast-math. No input polynomial rewriting is performed.
Requires Linux, Symbolica, an activated Wolfram Engine, FORM and GCC/G++.
All runs use fresh processes and one CPU; the RAM sampler uses a separate CPU.
"""
from __future__ import annotations

import argparse
import ast
import ctypes
from collections import Counter
from datetime import datetime, timezone
from fractions import Fraction
import hashlib
from importlib.metadata import version
import json
import math
import multiprocessing
import os
from pathlib import Path
import platform
import random
import re
import shutil
import signal
import statistics
import subprocess
import sys
import tempfile
import time
from urllib.request import urlopen

NAMES = ['alpha','amuq','ammu','xcp1','e1245','xcp4','e3e2','e1234','e2345',
         'e1235','e1345','amel2','e2e1','e5e2','e4e2','e3e1','e4e1','e5e1',
         'ammu2','amuq2','e5e3','e4e3','x5','x6','x1','x3','x4','xcp3','xcp2']
EXPECTED_SHA256 = '5c5384e0d929d25eabbc32f24be2edd2a6a743d673a5f8a9e2aae4b4598365cf'
ROUTES = ['symbolica_asm', 'symbolica_symjit', 'mathematica_unoptimized', 'mathematica_optimized', 'form_cpp']
DENOMINATOR = 8


def reference_data(path: Path, batch_size: int) -> dict:
    """Independent integer parser/evaluator; exact dyadic inputs and results."""
    source = re.sub(r'\s+', '', path.read_text())
    monomials = re.findall(r'[+-]?[^+-]+', source)
    if ''.join(monomials) != source:
        raise ValueError('Input contains unsupported syntax')
    terms = []
    for monomial in monomials:
        factors = monomial.split('*')
        coefficient = int(factors[0])
        powers = []
        for factor in factors[1:]:
            match = re.fullmatch(r'([a-zA-Z][a-zA-Z0-9]*)(?:\^(\d+))?', factor)
            if not match:
                raise ValueError(f'Unsupported monomial factor: {factor}')
            powers.append((NAMES.index(match[1]), int(match[2] or 1)))
        terms.append((coefficient, powers))
    max_degree = max(sum(n for _, n in powers) for _, powers in terms)
    max_power = max(n for _, powers in terms for _, n in powers)
    rng = random.Random(130013)
    # 32 distinct mixed-sign points, replicated to make the timed batch.
    points = [[rng.choice([-1, 1])*rng.randrange(4, 13) for _ in NAMES] for _ in range(32)]
    references, conditions = [], []
    for point in points:
        powers_table = [[v**n for n in range(max_power + 1)] for v in point]
        exact = absolute_sum = 0
        for coefficient, powers in terms:
            value = coefficient
            degree = 0
            for index, exponent in powers:
                value *= powers_table[index][exponent]
                degree += exponent
            value *= DENOMINATOR**(max_degree-degree)
            exact += value
            absolute_sum += abs(value)
        assert exact != 0
        references.append(float(Fraction(exact, DENOMINATOR**max_degree)))
        conditions.append(float(Fraction(absolute_sum, abs(exact))))
    return {'names': NAMES, 'terms': len(terms), 'max_total_degree': max_degree,
            'input_numerators': points, 'input_denominator': DENOMINATOR,
            'expected': references, 'sum_absolute_terms_over_absolute_result': conditions,
            'batch_size': batch_size, 'seed': 130013}


def check_results(values, expected) -> dict:
    assert len(values) == len(expected)
    errors = []
    for actual, wanted in zip(values, expected):
        assert math.isfinite(actual)
        error = abs(actual-wanted)/abs(wanted)
        assert error < 1e-9, (actual, wanted, error)
        errors.append(error)
    return {'verified': True, 'points': len(expected), 'relative_tolerance': 1e-9,
            'max_relative_error': max(errors), 'values': list(map(float, values))}


def process_tree(pid, excluded=0):
    result = []
    if pid == excluded:
        return result
    try:
        status = Path(f'/proc/{pid}/status').read_text()
        rss = int(re.search(r'^VmRSS:\s+(\d+)', status, re.M)[1])*1024
        result.append((pid, rss))
        children = set()
        for task in Path(f'/proc/{pid}/task').iterdir():
            try:
                children.update(map(int, (task/'children').read_text().split()))
            except (FileNotFoundError, ProcessLookupError):
                pass
    except (FileNotFoundError, ProcessLookupError, TypeError):
        return result
    for child in children:
        result.extend(process_tree(child, excluded))
    return result


def wolfram_operations(program: str) -> dict:
    """Count scalar arithmetic in CompilePrint, rejecting unknown operations.

    N-argument addition/multiplication costs N-1; Square costs one multiply;
    positive integer powers use binary exponentiation. Copies/loads are free.
    This is an IR count, before the native compiler's peephole/FMA passes.
    """
    assert 'MainEvaluate' not in program and 'External' not in program
    integer_constants = {name: int(value) for name, value in
                         re.findall(r'^\s*(I\d+) = (-?\d+)\s*$', program, re.M)}
    counts = Counter()
    classified = 0
    for line in program.splitlines():
        numbered = re.match(r'^\d+\s+(.+)$', line)
        if not numbered or numbered[1].strip() == 'Return':
            continue
        match = re.fullmatch(r'([RI]\d+) = (.+)', numbered[1])
        if not match:
            raise ValueError(f'Unclassified Mathematica instruction: {line}')
        target, rhs = match.groups()
        rhs = rhs.strip()
        classified += 1
        if re.fullmatch(r'[RI]\d+|[-+]?\d+(?:\.\d*)?(?:\*\^[+-]?\d+)?', rhs):
            continue
        if re.fullmatch(r'R\d+(?: \+ R\d+)+', rhs):
            counts['add'] += rhs.count('+')
        elif re.fullmatch(r'R\d+(?: \* R\d+)+', rhs):
            counts['mul'] += rhs.count('*')
        elif re.fullmatch(r'Square\[\s*R\d+\s*\]', rhs):
            counts['mul'] += 1
        elif power := re.fullmatch(r'Power\[\s*R\d+,\s*(I\d+)\s*\]', rhs):
            exponent = integer_constants[power[1]]
            assert exponent >= 1
            counts['mul'] += exponent.bit_length()+exponent.bit_count()-2
        elif re.fullmatch(r'- R\d+', rhs):
            counts['neg'] += 1
        elif re.fullmatch(r'R\d+ - R\d+', rhs):
            counts['add'] += 1
        else:
            raise ValueError(f'Unclassified Mathematica instruction: {target} = {rhs}')
    assert classified > 0
    return {'arithmetic': dict(counts), 'arithmetic_operations': sum(counts.values()),
            'classified_assignment_instructions': classified,
            'compile_print_sha256': hashlib.sha256(program.encode()).hexdigest()}


def wait_file(path, process=None, timeout=120):
    deadline = time.monotonic()+timeout
    while not path.exists():
        if process is not None and process.poll() is not None:
            raise RuntimeError(f'Engine exited before {path.name}; see engine.log')
        if time.monotonic() > deadline:
            raise TimeoutError(f'Waiting for {path}')
        time.sleep(.005)


def gate(work, phase):
    (work/f'{phase}.ready').touch()
    wait_file(work/f'{phase}.start', timeout=6000)


def phase_done(work, phase):
    (work/f'{phase}.done').touch()
    wait_file(work/f'{phase}.ack', timeout=6000)


def throughput(evaluate, batch, repeats):
    evaluate(batch)
    samples = []
    for _ in range(repeats):
        start = time.perf_counter()
        evaluate(batch)
        samples.append(time.perf_counter()-start)
    return samples


def symbolica_route(route, source, reference, work, repeats):
    import numpy as np
    from symbolica import E, S
    expression = E(source.read_text())
    assert len(expression) == reference['terms']
    params = [S(name) for name in NAMES]
    points = np.asarray(reference['input_numerators'], dtype=np.float64)/DENOMINATOR
    batch = np.ascontiguousarray(points[np.arange(reference['batch_size']) % len(points)])
    gate(work, 'generation')
    start = time.perf_counter()
    evaluator = expression.evaluator(params, iterations=10, n_cores=1,
                                    jit_compile=False, direct_translation=True)
    optimized = time.perf_counter()
    if route == 'symbolica_asm':
        executable = evaluator.compile('f13', str(work/'generated.cpp'), str(work/'generated.so'),
                                       'real', inline_asm='x64', optimization_level=3,
                                       compiler_path='g++', native=True,
                                       compiler_flags=['-fno-fast-math', '-fno-unsafe-math-optimizations'])
    else:
        evaluator.jit_compile(True, direct_translation=True, optimization_level=3)
        executable = evaluator
    executable.evaluate(points[:4])  # includes lazy JIT, loading, and first evaluation
    finished = time.perf_counter()
    phase_done(work, 'generation')
    instructions = evaluator.get_instructions()
    assert instructions.input_count == len(NAMES)
    counts = Counter()
    for instruction in instructions.instructions:
        assert instruction[0] in ['add', 'mul']
        counts[instruction[0]] += len(instruction[2])-1
    assert not instructions.sub_evaluators
    verification = check_results(executable.evaluate(points).ravel(), reference['expected'])
    gate(work, 'evaluation')
    samples = throughput(executable.evaluate, batch, repeats)
    phase_done(work, 'evaluation')
    if route == 'symbolica_asm':
        assert '__asm__' in (work/'generated.cpp').read_text()
    return {'version': version('symbolica'), 'numpy_version': np.__version__,
            'argument_count': instructions.input_count, 'kernel_callbacks': 0,
            'optimization_seconds': optimized-start, 'native_compilation_seconds': finished-optimized,
            'total_build_seconds': finished-start, 'arithmetic': dict(counts),
            'arithmetic_operations': sum(counts.values()), 'verification': verification,
            'evaluation_batch_seconds': samples,
            'evaluation_seconds_per_point': statistics.median(samples)/len(batch),
            'settings': {'iterations': 10, 'n_cores': 1, 'direct_translation': True,
                         'backend': 'scalar x64 inline ASM, real' if route == 'symbolica_asm' else 'SymJIT',
                         'optimization_level': 3, 'native_architecture': True,
                         'cpp_flags_override': ['-fno-fast-math', '-fno-unsafe-math-optimizations']}}


WOLFRAM_NATIVE = r'''
$HistoryLength=0; $IterationLimit=10000000;
Needs["CompiledFunctionTools`"]; Needs["CCompilerDriver`"]; Needs["CCompilerDriver`GCCCompiler`"];
commands={}; compilerStart=Null;
recordCompiler[command_]:=(If[!NumberQ[compilerStart],compilerStart=SessionTime[]];AppendTo[commands,command]);
SetOptions[CCompilerDriver`GCCCompiler`GCCCompiler,
 "SystemCompileOptions"->{"-fPIC","-O3","-march=native"},
 "ShellCommandFunction"->recordCompiler];
SetSystemOptions["ParallelOptions"->{"MKLThreadNumber"->1,"ParallelThreadNumber"->1}];
reference=Import["REFERENCE","RawJSON"];
vars=Symbol/@reference["names"];
expression=ToExpression[StringReplace[Import["SOURCE","Text"],{"\n"->"","\r"->""}]];
If[Length[expression]!=reference["terms"],Exit[2]];
points=N[reference["input_numerators"]/reference["input_denominator"]];
batch=Transpose[Table[points[[1+Mod[i-1,Length[points]]]],{i,reference["batch_size"]}]];
Export["WORK/generation.ready","ready","Text"];
While[!FileExistsQ["WORK/generation.start"],Pause[.005]];
commands={}; compilerStart=Null; generationStart=SessionTime[];
{build,cf}=AbsoluteTiming[Compile[ARGUMENTS,Evaluate[expression],
 CompilationTarget->"C",RuntimeOptions->"Speed",RuntimeAttributes->{Listable},
 Parallelization->False,CompilationOptions->{"ExpressionOptimization"->OPTIMIZE}]];
Export["WORK/generation.done","done","Text"];
While[!FileExistsQ["WORK/generation.ack"],Pause[.005]];
If[Head[cf]=!=CompiledFunction || FreeQ[cf,_LibraryFunction] ||
 cf[[2]]=!=ConstantArray[_Real,Length[vars]],Exit[3]];
program=CompiledFunctionTools`CompilePrint[cf];
Export["WORK/instructions.txt",program,"Text"];
If[StringContainsQ[program,"MainEvaluate"|"External"],Exit[4]];
values=Apply[cf,Transpose[points]];
Export["WORK/evaluation.ready","ready","Text"];
While[!FileExistsQ["WORK/evaluation.start"],Pause[.005]];
Apply[cf,batch];
samples=Table[First@AbsoluteTiming[Apply[cf,batch];],{REPEATS}];
Export["WORK/evaluation.done","done","Text"];
While[!FileExistsQ["WORK/evaluation.ack"],Pause[.005]];
Export["WORK/native.json",<|"version"->$Version,"argument_count"->Length[cf[[2]]],
 "kernel_callbacks"->0,"native_library_verified"->True,"total_build_seconds"->build,
 "optimization_seconds"->compilerStart-generationStart,
 "native_compilation_seconds"->build-(compilerStart-generationStart),
 "evaluation_batch_seconds"->samples,"evaluation_seconds_per_point"->Median[samples]/reference["batch_size"],
 "verification_values"->values,"compiler_commands"->Flatten[commands],
 "expression_optimization"->OPTIMIZE|>,"RawJSON"];
'''


def mathematica_route(route, source, reference_path, reference, work, repeats):
    code = WOLFRAM_NATIVE
    substitutions = {'ARGUMENTS': '{'+','.join('{'+name+',_Real}' for name in NAMES)+'}',
                     'OPTIMIZE': 'True' if route == 'mathematica_optimized' else 'False',
                     'REFERENCE': str(reference_path), 'SOURCE': str(source),
                     'WORK': str(work), 'REPEATS': str(repeats)}
    for key, value in substitutions.items():
        code = code.replace(key, value)
    script = work/'run.wls'
    script.write_text(code)
    with (work/'engine.log').open('w') as log:
        process = subprocess.Popen(['wolframscript', '-file', str(script)], stdout=log, stderr=log)
        process.wait()
    if process.returncode or not (work/'native.json').exists():
        raise RuntimeError(f'Wolfram exited {process.returncode}: {(work/"engine.log").read_text()}')
    result = json.loads((work/'native.json').read_text())
    assert result['argument_count'] == len(NAMES) and result['native_library_verified']
    assert result['compiler_commands'], 'No actual C compiler command was captured'
    assert any(re.search(r'(?:^|\s)-O3(?:\s|$)', c) for c in result['compiler_commands'])
    assert not any('-O2' in c for c in result['compiler_commands'])
    result.update(wolfram_operations((work/'instructions.txt').read_text()))
    result['verification'] = check_results(result.pop('verification_values'), reference['expected'])
    return result


def cpp_operations(text):
    counts = Counter()
    def visit(node):
        if isinstance(node, ast.BinOp):
            name = {ast.Add: 'add', ast.Sub: 'add', ast.Mult: 'mul', ast.Div: 'div'}.get(type(node.op))
            assert name, ast.dump(node)
            counts[name] += 1
        elif isinstance(node, ast.Call):
            assert isinstance(node.func, ast.Name) and node.func.id == 'pow' and len(node.args) == 2
            exponent = ast.literal_eval(node.args[1]); assert isinstance(exponent, int) and exponent >= 1
            counts['mul'] += exponent.bit_length()+exponent.bit_count()-2
            visit(node.args[0]); return
        elif isinstance(node, ast.UnaryOp) and isinstance(node.op, ast.USub):
            if not isinstance(node.operand, ast.Constant):
                counts['neg'] += 1
        for child in ast.iter_child_nodes(node):
            visit(child)
    for statement in text.split(';'):
        if '=' not in statement:
            continue
        rhs = statement.split('=', 1)[1]
        visit(ast.parse(re.sub(r'\s+', ' ', rhs).strip(), mode='eval'))
    assert counts
    return {'arithmetic': dict(counts), 'arithmetic_operations': sum(counts.values())}


def form_route(source, reference, work, repeats):
    import numpy as np
    # FORM stops itself after parsing/sorting. This keeps input setup outside generation.
    script = f'''Off Statistics;
Symbols {','.join(NAMES)};
Local F =
#include {source}
;
.sort
#system kill -STOP $PPID
ExtraSymbols,array,w;
Format C;
Format O4;
#optimize F
#write <generated.inc> "double w[1+`optimmaxvar_'];"
#write <generated.inc> "%O"
#write <generated.inc> "double result = %e",F
#write "FORM_OPERATIONS:`optimvalue_'"
.end
'''
    (work/'run.frm').write_text(script)
    log = (work/'engine.log').open('w')
    process = subprocess.Popen(['form', '-q', 'run.frm'], cwd=work, stdout=log, stderr=log)
    deadline = time.monotonic()+120
    while True:
        if process.poll() is not None:
            raise RuntimeError((work/'engine.log').read_text())
        if re.search(r'^State:\s+T', Path(f'/proc/{process.pid}/status').read_text(), re.M):
            break
        if time.monotonic() > deadline:
            raise TimeoutError('FORM input preparation')
        time.sleep(.005)
    points = np.asarray(reference['input_numerators'], dtype=np.float64)/DENOMINATOR
    batch = np.ascontiguousarray(points[np.arange(reference['batch_size']) % len(points)])
    gate(work, 'generation')
    start = time.perf_counter()
    process.send_signal(signal.SIGCONT)
    process.wait()
    log.close()
    assert process.returncode == 0, (work/'engine.log').read_text()
    generated = (work/'generated.inc').read_text()
    declarations = '\n'.join(f'const double {name}=p[{index}];' for index,name in enumerate(NAMES))
    cpp = ('#include <cmath>\n#include <cstddef>\n'
           'static double f13_scalar(const double* p) {\n'+declarations+'\n'+generated+
           '\nreturn result;\n}\nextern "C" void f13_batch(const double* p,double* out,std::size_t n) {\n'
           f'for(std::size_t i=0;i<n;++i) out[i]=f13_scalar(p+{len(NAMES)}*i);\n}}\n')
    (work/'generated.cpp').write_text(cpp)
    optimized = time.perf_counter()
    command = ['g++', '-shared', '-fPIC', '-O3', '-march=native', '-std=c++17',
               str(work/'generated.cpp'), '-o', str(work/'generated.so')]
    compiler = subprocess.run(command, capture_output=True, text=True)
    (work/'compiler.log').write_text(compiler.stdout+compiler.stderr)
    assert compiler.returncode == 0, compiler.stderr
    library = ctypes.CDLL(str(work/'generated.so'))
    function = library.f13_batch
    array = np.ctypeslib.ndpointer(dtype=np.float64, flags='C_CONTIGUOUS')
    function.argtypes = [array, array, ctypes.c_size_t]; function.restype = None
    def evaluate(p):
        out = np.empty(len(p), dtype=np.float64)
        function(p, out, len(p))
        return out
    evaluate(points[:4])
    finished = time.perf_counter()
    phase_done(work, 'generation')
    verification = check_results(evaluate(points), reference['expected'])
    gate(work, 'evaluation')
    samples = throughput(evaluate, batch, repeats)
    phase_done(work, 'evaluation')
    return {'version': subprocess.check_output(['form', '-v'], text=True).splitlines()[0],
            'argument_count': len(NAMES), 'kernel_callbacks': 0,
            'optimization_seconds': optimized-start, 'native_compilation_seconds': finished-optimized,
            'total_build_seconds': finished-start, 'compiler_commands': [command],
            'settings': {'form_optimization': 'O4 (default 1000 SA iterations)', 'cpp_flags': command[1:7]},
            'verification': verification, 'evaluation_batch_seconds': samples,
            'evaluation_seconds_per_point': statistics.median(samples)/len(batch),
            **cpp_operations(generated)}


def stop_tree(pid):
    tree = process_tree(pid)
    for child,_ in reversed(tree):
        try: os.kill(child, signal.SIGSTOP)
        except ProcessLookupError: pass
    for child,_ in reversed(tree):
        try: os.kill(child, signal.SIGKILL)
        except ProcessLookupError: pass


def measure_phase(process, work, phase, timeout, min_available_mib):
    wait_file(work/f'{phase}.ready', process)
    baseline = peak = sum(rss for _,rss in process_tree(process.pid))
    samples = 1
    observed = {}
    peak_processes = []
    started = time.perf_counter()
    last_saved = started
    def snapshot(elapsed):
        memory = {'baseline_rss_bytes': baseline, 'peak_tree_rss_bytes': peak,
                  'peak_increase_bytes': max(0, peak-baseline),
                  'sample_interval_seconds': .01, 'sample_count': samples,
                  'observed_compiler_processes': list(observed.values()),
                  'processes_at_sampled_peak': peak_processes,
                  'includes': 'worker and all live engine/compiler/formatter descendants; profiler excluded'}
        progress = {'phase': phase, 'elapsed_seconds': elapsed, 'status': status, 'memory': memory}
        temporary = work/f'{phase}-measurement.json.tmp'
        temporary.write_text(json.dumps(progress, indent=2)+'\n')
        temporary.replace(work/f'{phase}-measurement.json')
        return memory
    (work/f'{phase}.start').touch()
    status = 'completed'
    while not (work/f'{phase}.done').exists():
        if process.poll() is not None:
            status = 'worker_error'
            break
        tree = process_tree(process.pid)
        rss = sum(memory for _,memory in tree)
        if rss > peak:
            peak = rss
            peak_processes = [{'pid': pid, 'rss_bytes': memory} for pid,memory in tree]
        for pid,_ in tree:
            if pid == process.pid:
                continue
            try:
                command = Path(f'/proc/{pid}/cmdline').read_bytes().decode().strip('\0').split('\0')
            except (FileNotFoundError, ProcessLookupError):
                continue
            if command and Path(command[0]).name in {'g++','gcc','cc1','cc1plus','clang-format','clang','clang++','as','ld','collect2'}:
                observed[(pid,tuple(command))] = {'pid': pid, 'command': command}
        samples += 1
        available = int(next(line.split()[1] for line in Path('/proc/meminfo').read_text().splitlines()
                             if line.startswith('MemAvailable:')))*1024
        elapsed = time.perf_counter()-started
        if time.perf_counter()-last_saved >= 1:
            snapshot(elapsed)
            last_saved = time.perf_counter()
        if timeout and elapsed > timeout:
            status = 'generation_timeout'
        if min_available_mib and available < min_available_mib*2**20:
            status = 'memory_limit'
        if status != 'completed':
            stop_tree(process.pid); process.wait()
            break
        time.sleep(.01)
    elapsed = time.perf_counter()-started
    peak = max(peak, sum(rss for _,rss in process_tree(process.pid)))
    if status == 'completed':
        (work/f'{phase}.ack').touch()
    return status, elapsed, snapshot(elapsed)


def route_fingerprint(route):
    """Identify unchanged engine code when carrying validated runs into a resumed matrix."""
    names = ['check_results', 'throughput', 'gate', 'phase_done']
    if route == 'form_cpp':
        names += ['form_route', 'cpp_operations']
    elif route.startswith('mathematica'):
        names += ['mathematica_route', 'wolfram_operations', 'WOLFRAM_NATIVE']
    else:
        names += ['symbolica_route']
    nodes = []
    for node in ast.parse(Path(__file__).read_text()).body:
        if isinstance(node, (ast.FunctionDef, ast.AsyncFunctionDef)) and node.name in names:
            nodes.append(node)
        if isinstance(node, ast.Assign) and any(isinstance(t, ast.Name) and
                t.id in names+['NAMES', 'DENOMINATOR'] for t in node.targets):
            nodes.append(node)
    return hashlib.sha256('\n'.join(ast.dump(n, include_attributes=False) for n in nodes).encode()).hexdigest()


def worker(args):
    os.sched_setaffinity(0, {args.engine_cpu})
    reference = json.loads(args.reference.read_text())
    work = args.work.resolve()
    source = args.input.resolve()
    assert hashlib.sha256(source.read_bytes()).hexdigest() == reference['input_sha256']
    if args.worker.startswith('symbolica'):
        result = symbolica_route(args.worker, source, reference, work, args.repeats)
    elif args.worker.startswith('mathematica'):
        result = mathematica_route(args.worker, source, args.reference.resolve(), reference, work, args.repeats)
    else:
        result = form_route(source, reference, work, args.repeats)
    assert result['argument_count'] == 29 and result['kernel_callbacks'] == 0
    assert result['optimization_seconds'] >= 0 and result['native_compilation_seconds'] >= 0
    assert abs(result['optimization_seconds']+result['native_compilation_seconds']-result['total_build_seconds']) < 1e-6
    result.update(engine=args.worker, status='completed')
    (work/'result.json').write_text(json.dumps(result, indent=2)+'\n')


def main():
    parser = argparse.ArgumentParser(description='Five-route numerical F13 benchmark; subset validation runs first.')
    parser.add_argument('--input', type=Path, default=Path('posts/expression_evaluation/f13.txt'))
    parser.add_argument('--output', type=Path, default=Path('scripts/benchmarks/f13-subset-results.json'))
    parser.add_argument('--terms', type=int, default=256, help='evenly spaced monomials; 0 selects full F13')
    parser.add_argument('--runs', type=int, default=3)
    parser.add_argument('--repeats', type=int, default=7)
    parser.add_argument('--batch-size', type=int, default=4096)
    parser.add_argument('--timeout', type=float, default=5400)
    parser.add_argument('--min-available-mib', type=float, default=4096)
    parser.add_argument('--artifacts-dir', type=Path, default=Path('.quarto/f13-validation'))
    parser.add_argument('--worker', choices=ROUTES)
    parser.add_argument('--reference', type=Path)
    parser.add_argument('--work', type=Path)
    parser.add_argument('--engine-cpu', type=int, default=0)
    parser.add_argument('--validation', type=Path, default=Path('scripts/benchmarks/f13-subset-results.json'))
    parser.add_argument('--reuse-results', type=Path, action='append', default=[],
                        help='carry validated rows from unchanged engine implementations')
    args = parser.parse_args()
    if args.worker:
        worker(args); return
    if args.runs < 1 or args.repeats < 1 or args.batch_size < 32 or args.terms < 0:
        parser.error('invalid sample counts')
    affinity = sorted(os.sched_getaffinity(0))
    assert len(affinity) >= 2, 'Engine and sampler need separate logical CPUs'
    engine_cpu = args.engine_cpu
    assert engine_cpu in affinity
    profiler_cpu = next(cpu for cpu in affinity if cpu != engine_cpu)
    os.sched_setaffinity(0, {profiler_cpu})
    original = args.input.read_bytes()
    assert hashlib.sha256(original).hexdigest() == EXPECTED_SHA256
    monomials = re.findall(r'[+-]?[^+-]+', re.sub(r'\s+', '', original.decode()))
    assert len(monomials) == 105114
    if not args.terms:
        validation = json.loads(args.validation.read_text())
        assert validation['schema_version'] == 3 and 0 < validation['input']['terms'] < 105114
        assert validation['benchmark_source_sha256'] == hashlib.sha256(Path(__file__).read_bytes()).hexdigest()
        assert {r['engine'] for r in validation['results']} == set(ROUTES)
        assert all(r['status'] == 'completed' and r['verification']['verified']
                   and r['argument_count'] == 29 and r['kernel_callbacks'] == 0 for r in validation['results'])
    count = args.terms or len(monomials)
    assert 1 <= count <= len(monomials)
    indices = ([0] if count == 1 else [i*(len(monomials)-1)//(count-1) for i in range(count)])
    selected = '\n'.join(('-' if monomials[i].startswith('-') else '+')+monomials[i].lstrip('+-') for i in indices).lstrip('+')+'\n'
    root = (args.artifacts_dir/f'terms-{count}-{datetime.now(timezone.utc).strftime("%Y%m%dT%H%M%SZ")}').resolve()
    root.mkdir(parents=True, exist_ok=False)
    source = root/'input.txt'; source.write_text(selected)
    reference = reference_data(source, args.batch_size)
    reference['input_sha256'] = hashlib.sha256(source.read_bytes()).hexdigest()
    reference_path = root/'reference.json'; reference_path.write_text(json.dumps(reference))
    cpu = next(line.split(':',1)[1].strip() for line in Path('/proc/cpuinfo').read_text().splitlines() if line.startswith('model name'))
    data = {'schema_version': 3, 'measured_at_utc': datetime.now(timezone.utc).isoformat(),
            'cpu': cpu, 'platform': platform.platform(), 'python': platform.python_version(),
            'engine_cpu': engine_cpu, 'profiler_cpu': profiler_cpu, 'input': reference,
            'original_input_sha256': EXPECTED_SHA256, 'selected_monomial_indices': indices,
            'benchmark_source_sha256': hashlib.sha256(Path(__file__).read_bytes()).hexdigest(),
            'runs_per_engine': args.runs, 'evaluation_repeats': args.repeats,
            'generation_timeout_seconds': args.timeout, 'min_available_mib': args.min_available_mib,
            'compiler_version': subprocess.check_output(['g++','--version'],text=True).splitlines()[0],
            'artifacts_dir': str(root), 'results': [],
            'method': {'memory': '10 ms external samples of benchmark worker and all engine/compiler/formatter descendants, traversing every thread; baseline after loading input; profiler on a separate CPU',
                       'generation': 'input parsing, process startup, imports, instruction inspection and correctness checks excluded; optimization, export, native compilation and library loading included',
                       'cpp': '-O3 -march=native; fast-math disabled; Symbolica scalar real x64 inline assembly, SymJIT level 3; FORM O4 default 1000 iterations',
                       'verification': '32 fixed mixed-sign dyadic points checked against independent exact integer monomial evaluation; 29 real inputs; no Wolfram callbacks',
                       'evaluation': '4096-point batches, one warmup and 7 timed samples; allocation included; scalar real ASM; one CPU',
                       'order': 'fresh processes, route order reversed on alternate repetitions'}}
    args.output.parent.mkdir(parents=True, exist_ok=True)
    for reusable in args.reuse_results:
        previous = json.loads(reusable.read_text())
        assert previous['original_input_sha256'] == EXPECTED_SHA256
        assert previous['input']['input_sha256'] == reference['input_sha256']
        for row in previous['results']:
            assert row['engine'] in ROUTES
            assert row['route_source_sha256'] == route_fingerprint(row['engine'])
            if row['status'] == 'completed':
                assert row['argument_count'] == 29 and row['kernel_callbacks'] == 0
                assert row['verification']['verified'] and row['verification']['points'] == 32
            else:
                assert row['status'] in ('compiler_error', 'generation_timeout', 'memory_limit')
            assert not any(r['engine'] == row['engine'] and r['run'] == row['run'] for r in data['results'])
            data['results'].append(row)
    args.output.write_text(json.dumps(data, indent=2)+'\n')
    failed_routes = set()
    failed_routes.update(r['engine'] for r in data['results'] if r['status'] != 'completed')
    for run in range(args.runs):
        for route in (ROUTES if run%2 == 0 else ROUTES[::-1]):
            if route in failed_routes:
                continue  # a single DNF is retained; do not repeat an expensive failed build
            if any(r['engine'] == route and r['run'] == run+1 for r in data['results']):
                continue
            work = root/f'{run+1}-{route}'; work.mkdir()
            print(f'Run {run+1}/{args.runs}: {route}, {count} terms', flush=True)
            command = [sys.executable, str(Path(__file__).resolve()), '--worker', route,
                       '--input', str(source), '--reference', str(reference_path), '--work', str(work),
                       '--engine-cpu', str(engine_cpu), '--repeats', str(args.repeats)]
            with (work/'worker.log').open('w') as log:
                process = subprocess.Popen(command, stdout=log, stderr=log, start_new_session=True)
                memory = None
                elapsed = 0
                try:
                    status, elapsed, memory = measure_phase(process, work, 'generation', args.timeout, args.min_available_mib)
                    if status == 'completed':
                        eval_status, _, eval_memory = measure_phase(process, work, 'evaluation', 180, args.min_available_mib)
                        assert eval_status == 'completed'
                        process.wait(timeout=180)
                        assert process.returncode == 0, (work/'worker.log').read_text()
                        result = json.loads((work/'result.json').read_text())
                        result['evaluation_memory'] = eval_memory
                    else:
                        result = {'engine': route, 'status': status, 'observed_generation_seconds': elapsed}
                    result.update(run=run+1, generation_memory=memory, artifacts_dir=str(work))
                except Exception as error:
                    stop_tree(process.pid); process.wait()
                    engine_log = (work/'engine.log').read_text() if (work/'engine.log').exists() else ''
                    status = 'compiler_error' if 'internal compiler error' in engine_log else 'worker_error'
                    if memory is None:
                        saved = work/'generation-measurement.json'
                        if saved.exists():
                            checkpoint = json.loads(saved.read_text())
                            memory, elapsed = checkpoint['memory'], checkpoint['elapsed_seconds']
                        else:
                            memory = {'baseline_rss_bytes': None, 'peak_tree_rss_bytes': None,
                                      'peak_increase_bytes': None}
                    result = {'engine': route, 'status': status, 'observed_generation_seconds': elapsed,
                              'error': str(error), 'engine_log': engine_log,
                              'run': run+1, 'generation_memory': memory, 'artifacts_dir': str(work)}
                except BaseException:
                    stop_tree(process.pid); process.wait(); raise
            result.update(engine_cpu=engine_cpu, profiler_cpu=profiler_cpu,
                          benchmark_source_sha256=data['benchmark_source_sha256'],
                          route_source_sha256=route_fingerprint(route))
            data['results'].append(result)
            args.output.write_text(json.dumps(data, indent=2)+'\n')
            peak = memory['peak_tree_rss_bytes']
            print(f'  {result["status"]}: build {result.get("total_build_seconds",elapsed):.4f} s; '
                  f'peak {peak/2**20 if peak is not None else "unavailable"} MiB', flush=True)
            if args.terms and result['status'] != 'completed':
                raise RuntimeError('Subset validation failed; full run is blocked')
            if result['status'] != 'completed':
                failed_routes.add(route)
    print(f'Saved {args.output}', flush=True)


if __name__ == '__main__':
    main()