wyj
0a46a7095b
Add error-controlled DP5(4) integration and trusted first-crossing localization, including non-monotonic energy thresholds and representable-time stepping. Preserve adaptive state and independent step/time retry grants across RayPool, refinement and movie scheduling. Expose numerical controls, record actual persistent-sample costs, and add v3 lens-map provenance with legacy v2 RK4 import. Use DP54 by default and select an 8M Schwarzschild maximum step from bounded scans and a two-run 4K comparison. Retain the conservative minimum-step guard and document critical-ray and backend capability limits. Archive self-contained benchmark inputs and raw output; keep fixed RK4 HDR references explicit. Validation: make -B -j4 BUILD_TYPE=Debug test passed; explicit RK4 HDR references have zero differences. Bounded convergence checks, benchmark reproduction, Release build and focused reviews passed. No numerical-relativity backend is added.
78 lines
3.4 KiB
Python
78 lines
3.4 KiB
Python
#!/usr/bin/env python3
|
|
"""Pure GRLENS v3 map load/compare helpers.
|
|
|
|
No renderer, no environment probing, no data generation: this module is only
|
|
imported by the benchmark scripts under
|
|
benchmarks/adaptive_step_bounds_2026-10-05/. Keeping it here means the
|
|
benchmark never imports an untracked module from local/.
|
|
"""
|
|
import collections
|
|
import hashlib
|
|
import math
|
|
import struct
|
|
import zlib
|
|
|
|
|
|
def load(path):
|
|
"""Parse one GRLENS v3 map; returns (info dict, {film_pos: vertex tuple})."""
|
|
data = path.read_bytes()
|
|
assert data[:8] == b'GRLENS\x01\x00'
|
|
assert struct.unpack_from('<I', data, 8)[0] == 3
|
|
assert struct.unpack_from('<Q', data, 32)[0] == 1
|
|
nv, nt = struct.unpack_from('<QQ', data, 200)
|
|
end = 224 + nv * 108 + nt * 32
|
|
assert end + 4 == len(data)
|
|
assert zlib.crc32(data[224:end]) == struct.unpack_from('<I', data, end)[0]
|
|
vertices = [struct.unpack_from('<9dIIIQQQ', data, 224 + i * 108)
|
|
for i in range(nv)]
|
|
costs = [sum(v[j] for v in vertices) for j in (12, 13, 14)]
|
|
result = dict(vertices=nv, triangles=nt,
|
|
outcomes=dict(collections.Counter(v[10] for v in vertices)),
|
|
reasons=dict(collections.Counter(v[11] for v in vertices)),
|
|
accepted=costs[0], rejected=costs[1], rhs=costs[2])
|
|
names = ('atol_x', 'atol_Pi', 'atol_L', 'rtol', 'min_step', 'max_step',
|
|
'max_lookback_time', 'retry_lookback_increment',
|
|
'max_total_lookback_time')
|
|
result['provenance'] = dict(zip(names, struct.unpack_from('<9d', data, 100)))
|
|
result['provenance'].update(
|
|
integrator=struct.unpack_from('<I', data, 68)[0],
|
|
initial_step=struct.unpack_from('<d', data, 88)[0],
|
|
initial_max_steps=struct.unpack_from('<I', data, 96)[0],
|
|
threshold=struct.unpack_from('<d', data, 48)[0],
|
|
retry_step_increment=struct.unpack_from('<I', data, 56)[0],
|
|
max_total_steps=struct.unpack_from('<I', data, 60)[0])
|
|
for j, name in ((12, 'accepted'), (13, 'rejected'), (14, 'rhs')):
|
|
xs = sorted(v[j] for v in vertices)
|
|
result[name + '_percentiles'] = {
|
|
str(q): xs[min(len(xs) - 1, int(q * (len(xs) - 1)))]
|
|
for q in (0, .5, .9, .99, 1)}
|
|
result['sha256'] = hashlib.sha256(data).hexdigest()
|
|
result['triangle_sha256'] = hashlib.sha256(
|
|
data[224 + nv * 108:end]).hexdigest()
|
|
by_film = {(v[0], v[1]): v for v in vertices}
|
|
result['unique_film_positions'] = len(by_film)
|
|
result['duplicate_film_positions'] = len(vertices) - len(by_film)
|
|
return result, by_film
|
|
|
|
|
|
def compare(a, b):
|
|
"""Compare two {film_pos: vertex} maps by shared film-coordinate identity."""
|
|
keys = a.keys() & b.keys()
|
|
mismatches = 0
|
|
worst_sky = worst_logg = 0.
|
|
for k in keys:
|
|
x, y = a[k], b[k]
|
|
if x[9:12] != y[9:12]:
|
|
mismatches += 1
|
|
elif x[10] == 0:
|
|
n, m = x[5:8], y[5:8]
|
|
cross = (n[1] * m[2] - n[2] * m[1], n[2] * m[0] - n[0] * m[2],
|
|
n[0] * m[1] - n[1] * m[0])
|
|
angle = math.atan2(math.sqrt(sum(v * v for v in cross)),
|
|
sum(v * w for v, w in zip(n, m)))
|
|
worst_sky = max(worst_sky, angle)
|
|
worst_logg = max(worst_logg, abs(x[8] - y[8]))
|
|
return dict(shared=len(keys), only_a=len(a) - len(keys),
|
|
only_b=len(b) - len(keys), terminal_mismatches=mismatches,
|
|
max_sky_angle=worst_sky, max_logg_difference=worst_logg)
|