675 lines
35 KiB
Python
675 lines
35 KiB
Python
#!/usr/bin/env python3
|
|
"""Exercise the production DP54 default, v3 lens-map provenance/replay and the
|
|
adaptive retry/budget policy from the CLI.
|
|
|
|
The production default is now adaptive Dormand-Prince 5(4). These checks
|
|
require that a run without --integrator exports wire code 1 and is physically
|
|
and bit-for-bit equal to an explicit --integrator dp54 run, so `make test`
|
|
truly covers the production default rather than only the explicit path.
|
|
--integrator rk4 stays available for the legacy wire code and for a
|
|
fixed-step reference convergence check. Small CPU 16x8/32x16 scenes keep the
|
|
runtime short; physical comparisons use stored lens-map endpoints, not only the
|
|
rendered PNG.
|
|
"""
|
|
import math
|
|
import os
|
|
import re
|
|
import struct
|
|
import subprocess
|
|
import sys
|
|
import tempfile
|
|
import zlib
|
|
from pathlib import Path
|
|
|
|
# Keep scratch data inside the pre-approved OpenCode scratch directory instead
|
|
# of creating directories directly under the system temporary root.
|
|
TMP_ROOT = Path('/tmp/opencode')
|
|
TMP_ROOT.mkdir(parents=True, exist_ok=True)
|
|
|
|
BUILD = Path(sys.argv[1] if len(sys.argv) > 1 else 'build/Release').resolve()
|
|
TESTDIR = Path(sys.argv[2]).resolve() if len(sys.argv) > 2 else BUILD
|
|
ENV = dict(os.environ, OMP_NUM_THREADS='4')
|
|
|
|
# v3 wire offsets (see src/lens_map.c). The header is deliberately not part of
|
|
# the payload CRC.
|
|
VERSION_OFFSET = 8
|
|
FRAME_COUNT_OFFSET = 32
|
|
PROVENANCE_OFFSET = 40
|
|
PROVENANCE_V3_OFFSET = 100
|
|
VERTEX_COUNT_OFFSET = 200
|
|
TRIANGLE_COUNT_OFFSET = 208
|
|
VERTEX_START = 224
|
|
VERTEX_SIZE = 108
|
|
TRIANGLE_SIZE = 32
|
|
|
|
ATOL_FIELDS = ('atol_x', 'atol_Pi', 'atol_L', 'rtol', 'min_step', 'max_step',
|
|
'max_lookback_time', 'retry_lookback_increment',
|
|
'max_total_lookback_time')
|
|
# Machine-roundoff floor for comparing two adaptive integrations; below this a
|
|
# difference carries no convergence information.
|
|
ROUNDOFF_FLOOR = 1e-12
|
|
ENDPOINT_ASSERT = 1e-6
|
|
|
|
|
|
def run(binary, *args, ok=True, env=ENV):
|
|
result = subprocess.run([str(binary), *map(str, args)], env=env,
|
|
capture_output=True, text=True)
|
|
if (result.returncode == 0) != ok:
|
|
raise AssertionError(
|
|
f'{binary.name} {args}: rc={result.returncode}\n{result.stderr}')
|
|
return result
|
|
|
|
|
|
def image_payload(path):
|
|
data = path.read_bytes()
|
|
assert data[:8] == b'\x89PNG\r\n\x1a\n', f'not a PNG: {path}'
|
|
offset, compressed = 8, bytearray()
|
|
while offset < len(data):
|
|
count, kind = struct.unpack_from('>I4s', data, offset)
|
|
payload = data[offset + 8:offset + 8 + count]
|
|
if kind == b'IDAT':
|
|
compressed.extend(payload)
|
|
offset += count + 12
|
|
return bytes(zlib.decompress(compressed))
|
|
|
|
|
|
def map_provenance(path):
|
|
data = path.read_bytes()
|
|
assert data[:8] == b'GRLENS\x01\x00'
|
|
version = struct.unpack_from('<I', data, VERSION_OFFSET)[0]
|
|
frame_count = struct.unpack_from('<Q', data, FRAME_COUNT_OFFSET)[0]
|
|
base = PROVENANCE_OFFSET
|
|
prov = {}
|
|
prov['threshold_kind'], prov['threshold_policy_version'] = struct.unpack_from(
|
|
'<II', data, base)
|
|
prov['threshold_value'] = struct.unpack_from('<d', data, base + 8)[0]
|
|
(prov['retry_step_increment'], prov['max_total_steps'], prov['max_level'],
|
|
prov['integrator']) = struct.unpack_from('<IIII', data, base + 16)
|
|
prov['min_edge_pixels'] = struct.unpack_from('<d', data, base + 32)[0]
|
|
prov['min_area_pixels2'] = struct.unpack_from('<d', data, base + 40)[0]
|
|
prov['coordinate_time_step'] = struct.unpack_from('<d', data, base + 48)[0]
|
|
prov['initial_max_steps'] = struct.unpack_from('<I', data, base + 56)[0]
|
|
if version == 3:
|
|
off = PROVENANCE_V3_OFFSET
|
|
for name in ATOL_FIELDS:
|
|
prov[name] = struct.unpack_from('<d', data, off)[0]
|
|
off += 8
|
|
prov['max_consecutive_rejections'] = struct.unpack_from('<I', data, off)[0]
|
|
return version, frame_count, prov
|
|
|
|
|
|
def map_vertices(path):
|
|
data = path.read_bytes()
|
|
version = struct.unpack_from('<I', data, VERSION_OFFSET)[0]
|
|
assert version == 3, f'expected a v3 map, got v{version}'
|
|
vertices = struct.unpack_from('<Q', data, VERTEX_COUNT_OFFSET)[0]
|
|
triangles = struct.unpack_from('<Q', data, TRIANGLE_COUNT_OFFSET)[0]
|
|
offset = VERTEX_START
|
|
values = []
|
|
for _ in range(vertices):
|
|
values.append(struct.unpack_from('<9dIIIQQQ', data, offset))
|
|
offset += VERTEX_SIZE
|
|
return values, data[offset:offset + triangles * TRIANGLE_SIZE]
|
|
|
|
|
|
def endpoint_deviation(a, b):
|
|
"""(mismatched provenance count, max direction/log-g deviation).
|
|
|
|
Vertices whose end_id/outcome/reason differ are counted as mismatches and
|
|
excluded from the numeric deviation; the caller requires zero mismatches.
|
|
"""
|
|
mismatches = 0
|
|
worst = 0.0
|
|
assert len(a) == len(b)
|
|
for x, y in zip(a, b):
|
|
if x[9:12] != y[9:12]:
|
|
mismatches += 1
|
|
continue
|
|
for k in (3, 4, 5, 6, 7, 8):
|
|
worst = max(worst, abs(x[k] - y[k]))
|
|
return mismatches, worst
|
|
|
|
|
|
def trace_cost(stderr, label):
|
|
match = re.search(label + r' trace cost: accepted=(\d+) rejected=(\d+) '
|
|
r'rhs=(\d+)', stderr)
|
|
return None if match is None else match.groups()
|
|
|
|
|
|
def sum_rhs(vertices):
|
|
return sum(v[14] for v in vertices)
|
|
|
|
|
|
def alcubierre_budget(escape, vs, dark_threshold):
|
|
"""Reference implementation of the production Alcubierre time allowance:
|
|
B = 5*escape / max(|1-|v_s||, exp(-D)), with the DBL_MIN..DBL_MAX/4
|
|
saturation and the log-space fallback used when the ordinary division is
|
|
not finite and positive. exp(D) is never formed."""
|
|
sep = abs(1.0 - abs(vs))
|
|
floor = math.exp(-dark_threshold)
|
|
denom = max(sep, floor)
|
|
upper = sys.float_info.max / 4.0
|
|
if denom > 0.0 and math.isfinite(denom):
|
|
scaled = 5.0 * (escape / denom)
|
|
if math.isfinite(scaled) and scaled > 0.0:
|
|
return min(max(scaled, sys.float_info.min), upper)
|
|
log_sep = math.log(sep) if sep > 0.0 else -math.inf
|
|
log_budget = math.log(5.0) + math.log(escape) - max(log_sep, -dark_threshold)
|
|
if not math.isfinite(log_budget):
|
|
log_budget = math.log(upper)
|
|
log_budget = min(log_budget, math.log(upper))
|
|
log_budget = max(log_budget, math.log(sys.float_info.min))
|
|
return min(max(math.exp(log_budget), sys.float_info.min), upper)
|
|
|
|
|
|
with tempfile.TemporaryDirectory(prefix='gr-adaptive-cli-',
|
|
dir=str(TMP_ROOT)) as directory:
|
|
tmp = Path(directory)
|
|
for backend in ('minkowski', 'schwarzschild'):
|
|
binary = BUILD / f'{backend}_sky'
|
|
if not binary.exists():
|
|
print(f'{backend}: binary absent, skipping', flush=True)
|
|
continue
|
|
help_text = run(binary, '--help').stdout
|
|
for option in ('--integrator', '--ode-rtol', '--ode-atol-x',
|
|
'--ode-atol-pi', '--ode-atol-l', '--ode-initial-step',
|
|
'--ode-min-step', '--ode-max-step',
|
|
'--ode-max-rejections', '--trace-max-steps',
|
|
'--trace-lookback-time', '--retry-step-increment',
|
|
'--max-total-steps', '--retry-lookback-increment',
|
|
'--max-total-lookback-time'):
|
|
assert option in help_text, (backend, option)
|
|
ext = 'png' if '.png' in help_text else 'ppm'
|
|
hdr_available = '--hdr-output' in help_text
|
|
hdr_args = ['--hdr-output'] if hdr_available else []
|
|
common = ['--catalog', 'assets/sky_grid_5deg.csv', '--width', 32,
|
|
'--height', 16, '--fov-deg', 80, '--exposure', 1e-3,
|
|
'--coarse-cell-pixels', 8, '--refine-max-level', 0,
|
|
'--psf-relative-tail', 1e-4]
|
|
|
|
def single(name, *options, ok=True, env=ENV, use_common=common):
|
|
out = tmp / f'{backend}_{name}.{ext}'
|
|
result = run(binary, *use_common, '--output', out, *options,
|
|
ok=ok, env=env)
|
|
return out, result
|
|
|
|
# 1) The production default (no --integrator) must be DP54 and must
|
|
# match an explicit --integrator dp54 run physically and in its
|
|
# PNG/HDR output.
|
|
dflt_map = tmp / f'{backend}_dflt.grlens'
|
|
dflt_out, dflt_run = single('dflt', *hdr_args, '--verbose',
|
|
'--lens-map-output', dflt_map)
|
|
version, frame_count, dflt_prov = map_provenance(dflt_map)
|
|
assert version == 3 and frame_count == 1
|
|
assert dflt_prov['integrator'] == 1, dflt_prov
|
|
assert dflt_prov['min_step'] == 1e-12, dflt_prov
|
|
assert dflt_prov['max_step'] == {'minkowski': 16.0,
|
|
'schwarzschild': 8.0}[backend], dflt_prov
|
|
assert dflt_prov['min_step'] <= dflt_prov['coordinate_time_step'] \
|
|
<= dflt_prov['max_step']
|
|
assert dflt_prov['atol_x'] > 0 and dflt_prov['rtol'] > 0
|
|
assert dflt_prov['max_lookback_time'] > 0
|
|
assert dflt_prov['max_consecutive_rejections'] > 0
|
|
expl_map = tmp / f'{backend}_expl.grlens'
|
|
expl_out, _ = single('expl', *hdr_args, '--integrator', 'dp54',
|
|
'--lens-map-output', expl_map)
|
|
assert dflt_map.read_bytes() == expl_map.read_bytes(), \
|
|
'default map differs from explicit dp54'
|
|
assert image_payload(dflt_out) == image_payload(expl_out)
|
|
if hdr_available:
|
|
dflt_hdr = dflt_out.with_name(dflt_out.stem + '_HDR.fits')
|
|
expl_hdr = expl_out.with_name(expl_out.stem + '_HDR.fits')
|
|
assert dflt_hdr.read_bytes() == expl_hdr.read_bytes()
|
|
dflt_vertices, _ = map_vertices(dflt_map)
|
|
|
|
# 2) The legacy RK4 wire code must be 0 and its cost counters must be
|
|
# real (nonzero RHS evaluations), not legacy zeros.
|
|
rk4_map = tmp / f'{backend}_rk4.grlens'
|
|
single('rk4', '--integrator', 'rk4', '--lens-map-output', rk4_map)
|
|
_, _, rk4_prov = map_provenance(rk4_map)
|
|
assert rk4_prov['integrator'] == 0, rk4_prov
|
|
rk4_vertices, _ = map_vertices(rk4_map)
|
|
assert sum_rhs(rk4_vertices) > 0, 'RK4 RHS cost counters are not real'
|
|
|
|
# 3) Same-camera tolerance convergence with three levels. Outcome,
|
|
# reason and end must not change between levels, and the deviation
|
|
# from the tightest reference must shrink as the tolerance tightens.
|
|
# These are local ODE tolerances, not a global sky-error bound.
|
|
tol_maps = {}
|
|
for tol in ('1e-7', '1e-9', '1e-12'):
|
|
m = tmp / f'{backend}_tol_{tol}.grlens'
|
|
single(f'tol_{tol}', '--integrator', 'dp54', '--ode-rtol', tol,
|
|
'--ode-atol-x', tol, '--ode-atol-pi', tol,
|
|
'--ode-atol-l', tol, '--lens-map-output', m)
|
|
tol_maps[tol] = map_vertices(m)[0]
|
|
mism7, err7 = endpoint_deviation(tol_maps['1e-7'], tol_maps['1e-12'])
|
|
mism9, err9 = endpoint_deviation(tol_maps['1e-9'], tol_maps['1e-12'])
|
|
assert mism7 == 0 and mism9 == 0, \
|
|
f'{backend}: tolerance levels disagree on outcome/end'
|
|
assert err9 < ENDPOINT_ASSERT, (backend, 'default vs tight', err9)
|
|
assert err7 + ROUNDOFF_FLOOR >= err9, \
|
|
f'{backend}: tightening tolerance did not reduce error ' \
|
|
f'({err7} -> {err9})'
|
|
print(f'{backend}: default==dp54, tol errors 1e-7={err7:.3g} '
|
|
f'1e-9={err9:.3g}', flush=True)
|
|
|
|
# 4) Render-only replay of the default map must be bit-identical and
|
|
# its stored statistics must equal the live trace cost.
|
|
replay_out = tmp / f'{backend}_replay.{ext}'
|
|
replay_run = run(binary, *common, *hdr_args, '--lens-map-input',
|
|
dflt_map, '--verbose', '--output', replay_out)
|
|
assert image_payload(replay_out) == image_payload(dflt_out)
|
|
if hdr_available:
|
|
replay_hdr = replay_out.with_name(replay_out.stem + '_HDR.fits')
|
|
assert dflt_out.with_name(dflt_out.stem + '_HDR.fits').read_bytes() \
|
|
== replay_hdr.read_bytes()
|
|
live_cost = trace_cost(dflt_run.stderr, 'Frame 0')
|
|
replay_cost = trace_cost(replay_run.stderr, 'Imported map')
|
|
assert live_cost is not None and replay_cost is not None
|
|
assert live_cost == replay_cost, (live_cost, replay_cost)
|
|
|
|
# 5) Explicit zeros in the retry policy must survive, not be filled in
|
|
# by the derived defaults.
|
|
zero_step_map = tmp / f'{backend}_zero_step.grlens'
|
|
single('zero_step', '--integrator', 'dp54', '--retry-step-increment',
|
|
'0', '--lens-map-output', zero_step_map)
|
|
_, _, zero_step = map_provenance(zero_step_map)
|
|
assert zero_step['retry_step_increment'] == 0, zero_step
|
|
assert zero_step['max_total_steps'] > 0, zero_step
|
|
zero_time_map = tmp / f'{backend}_zero_time.grlens'
|
|
single('zero_time', '--integrator', 'dp54',
|
|
'--retry-lookback-increment', '0', '--lens-map-output',
|
|
zero_time_map)
|
|
_, _, zero_time = map_provenance(zero_time_map)
|
|
assert zero_time['retry_lookback_increment'] == 0, zero_time
|
|
assert zero_time['max_total_lookback_time'] >= \
|
|
zero_time['max_lookback_time'], zero_time
|
|
|
|
# 6) RK4 rejects every DP-only option rather than silently ignoring it.
|
|
rk4_errors = [
|
|
(['--integrator', 'rk4', '--ode-rtol', 1e-9], 'applies only'),
|
|
(['--integrator', 'rk4', '--ode-min-step', 1e-9], 'applies only'),
|
|
(['--integrator', 'rk4', '--ode-max-step', 1e-3], 'applies only'),
|
|
(['--integrator', 'rk4', '--ode-max-rejections', 4], 'applies only'),
|
|
(['--integrator', 'rk4', '--trace-lookback-time', 1], 'applies only'),
|
|
(['--integrator', 'rk4', '--retry-lookback-increment', 1], 'applies only'),
|
|
(['--integrator', 'rk4', '--max-total-lookback-time', 2], 'applies only'),
|
|
]
|
|
# DP cross-field validation and malformed values fail immediately.
|
|
invalid = [
|
|
(['--integrator', 'bogus'], None),
|
|
(['--integrator'], None),
|
|
(['--integrator', 'dp54', '--ode-rtol', 0], None),
|
|
(['--integrator', 'dp54', '--ode-rtol', -1], None),
|
|
(['--integrator', 'dp54', '--ode-atol-x', 'nan'], None),
|
|
(['--integrator', 'dp54', '--ode-max-rejections', 0], None),
|
|
(['--integrator', 'dp54', '--trace-max-steps', 0], None),
|
|
(['--integrator', 'dp54', '--ode-initial-step', 5,
|
|
'--ode-max-step', 1], 'min <= initial <= max'),
|
|
(['--integrator', 'dp54', '--ode-min-step', 4,
|
|
'--ode-initial-step', 2], 'min <= initial <= max'),
|
|
(['--integrator', 'dp54', '--trace-lookback-time', 0], None),
|
|
(['--integrator', 'dp54', '--retry-step-increment', -1], None),
|
|
]
|
|
for options, message in rk4_errors + invalid:
|
|
missing = tmp / 'adaptive_should_not_exist.csv'
|
|
result = run(binary, '--catalog', missing, *options, ok=False)
|
|
if message:
|
|
assert message in result.stderr, (options, result.stderr)
|
|
assert not missing.exists(), result.stderr
|
|
|
|
# 7) Replay must consume the stored policy: explicit tracing options
|
|
# cannot be layered on top of --lens-map-input.
|
|
conflict = run(binary, *common, '--lens-map-input', dflt_map,
|
|
'--integrator', 'rk4', '--output',
|
|
tmp / 'conflict.png', ok=False)
|
|
assert 'cannot be combined with --lens-map-input' in conflict.stderr, \
|
|
conflict.stderr
|
|
conflict2 = run(binary, *common, '--lens-map-input', dflt_map,
|
|
'--trace-max-steps', 100, '--output',
|
|
tmp / 'conflict2.png', ok=False)
|
|
assert 'cannot be combined with --lens-map-input' in conflict2.stderr, \
|
|
conflict2.stderr
|
|
conflict3 = run(binary, *common, '--lens-map-input', dflt_map,
|
|
'--integrator', 'dp54', '--ode-rtol', 1e-9,
|
|
'--output', tmp / 'conflict3.png', ok=False)
|
|
assert 'cannot be combined with --lens-map-input' in conflict3.stderr, \
|
|
conflict3.stderr
|
|
|
|
# 8) Single vs movie: the same physical observer event must agree, and
|
|
# threads/slabs must not change the stored adaptive endpoints.
|
|
track = tmp / f'{backend}.csv'
|
|
observer_test = TESTDIR / f'test_observer_{backend}'
|
|
if observer_test.exists():
|
|
run(observer_test, track)
|
|
moving_single = tmp / f'{backend}_moving.grlens'
|
|
single('moving', '--observer-position', 3, -4, 5,
|
|
'--observer-velocity', 0.2, -0.1, 0.3,
|
|
'--look-ra-deg', 37, '--look-dec-deg', -23,
|
|
'--camera-roll-deg', 19, '--lens-map-output', moving_single)
|
|
movie_map = tmp / f'{backend}_movie.grlens'
|
|
run(binary, *common, '--observer-track', track, '--frames-dir', tmp,
|
|
'--frames-prefix', f'{backend}_movie', '--duration', 0,
|
|
'--fps', 1, '--lens-map-output', movie_map)
|
|
single_v, _ = map_vertices(moving_single)
|
|
movie_v, _ = map_vertices(movie_map)
|
|
mism, dev = endpoint_deviation(single_v, movie_v)
|
|
assert mism == 0 and dev < ENDPOINT_ASSERT, (mism, dev)
|
|
|
|
thread_maps = []
|
|
for threads in (1, 2, 4):
|
|
m = tmp / f'{backend}_movie_{threads}thr.grlens'
|
|
run(binary, *common, '--observer-track', track, '--frames-dir',
|
|
tmp, '--frames-prefix', f'{backend}_m{threads}',
|
|
'--duration', 0, '--fps', 1, '--lens-map-output', m,
|
|
env=dict(ENV, OMP_NUM_THREADS=str(threads)))
|
|
thread_maps.append(m)
|
|
reference = thread_maps[0].read_bytes()
|
|
for m in thread_maps[1:]:
|
|
assert m.read_bytes() == reference, \
|
|
f'{backend}: DP movie map changed across threads'
|
|
|
|
slab_maps = []
|
|
for slab in (2, 8, 64):
|
|
m = tmp / f'{backend}_slab_{slab}.grlens'
|
|
run(binary, *common, '--observer-track', track, '--frames-dir',
|
|
tmp, '--frames-prefix', f'{backend}_s{slab}',
|
|
'--slab-duration', slab, '--duration', 0, '--fps', 1,
|
|
'--lens-map-output', m)
|
|
slab_maps.append(m)
|
|
base_v, _ = map_vertices(slab_maps[0])
|
|
for m in slab_maps[1:]:
|
|
other_v, _ = map_vertices(m)
|
|
mism, dev = endpoint_deviation(base_v, other_v)
|
|
assert mism == 0 and dev < ENDPOINT_ASSERT, (mism, dev)
|
|
print(f'{backend}: DP movie threads/slabs agree', flush=True)
|
|
|
|
# 9) Budget: a tiny initial coordinate-time budget leaves UNRESOLVED
|
|
# rays; with no room to grow the publication gate refuses the frame,
|
|
# while retry increments that can grow resolve it.
|
|
if backend == 'schwarzschild':
|
|
refused = tmp / f'{backend}_refused.{ext}'
|
|
small = ['--integrator', 'dp54', '--trace-lookback-time', 1e-6,
|
|
'--retry-lookback-increment', 0,
|
|
'--max-total-lookback-time', 1e-6]
|
|
result = run(binary, *common, '--output', refused, *small, ok=False)
|
|
assert 'Incomplete render refused' in result.stderr, result.stderr
|
|
assert not refused.exists()
|
|
allow = tmp / f'{backend}_allow.{ext}'
|
|
run(binary, *common, '--output', allow, '--allow-incomplete', *small)
|
|
assert allow.exists()
|
|
|
|
direct = tmp / f'{backend}_direct.{ext}'
|
|
direct_result = subprocess.run(
|
|
[str(binary), *map(str, common), '--output', str(direct),
|
|
'--integrator', 'dp54', '--trace-lookback-time', '4000'],
|
|
env=ENV, capture_output=True, text=True)
|
|
if direct_result.returncode == 0:
|
|
retried = tmp / f'{backend}_retried.{ext}'
|
|
run(binary, *common, '--output', retried, '--integrator', 'dp54',
|
|
'--trace-lookback-time', 1e-6,
|
|
'--retry-lookback-increment', 25,
|
|
'--max-total-lookback-time', 4000)
|
|
assert image_payload(retried)
|
|
print('schwarzschild: retry budget resolves previously '
|
|
'unresolved rays', flush=True)
|
|
else:
|
|
print('schwarzschild: direct budget scene still unresolved; '
|
|
'retry-resolution case skipped', flush=True)
|
|
|
|
# 10) Fixed-step reference convergence on a small scene. The
|
|
# accepted budget is explicit and large so the finer step does
|
|
# not silently shorten the traced history. Both fixed-step
|
|
# levels and the DP54 default must agree below ENDPOINT_ASSERT.
|
|
# If the reference itself does not converge this must FAIL and
|
|
# be reported, never loosened into a false zero.
|
|
ref_common = ['--catalog', 'assets/sky_grid_5deg.csv', '--width',
|
|
16, '--height', 8, '--fov-deg', 80, '--exposure',
|
|
1e-3, '--coarse-cell-pixels', 8, '--refine-max-level',
|
|
0, '--psf-relative-tail', 1e-4]
|
|
|
|
def ref_map(tag, *options):
|
|
m = tmp / f'{backend}_ref_{tag}.grlens'
|
|
run(binary, *ref_common, '--output',
|
|
tmp / f'{backend}_ref_{tag}.{ext}', '--lens-map-output', m,
|
|
*options)
|
|
return map_vertices(m)[0]
|
|
|
|
rk4_04 = ref_map('rk4_04', '--integrator', 'rk4',
|
|
'--ode-initial-step', 0.04, '--trace-max-steps',
|
|
'262144')
|
|
rk4_02 = ref_map('rk4_02', '--integrator', 'rk4',
|
|
'--ode-initial-step', 0.02, '--trace-max-steps',
|
|
'262144')
|
|
dp_tight = ref_map('dp_tight', '--ode-rtol', '1e-12',
|
|
'--ode-atol-x', '1e-12', '--ode-atol-pi',
|
|
'1e-12', '--ode-atol-l', '1e-12')
|
|
assert sum(1 for v in rk4_02 if v[10] == 0) > 0, \
|
|
'reference scene has no escaped ray'
|
|
mism, ref_dev = endpoint_deviation(rk4_04, rk4_02)
|
|
assert mism == 0, 'RK4 reference levels disagree on outcome/end'
|
|
assert ref_dev < ENDPOINT_ASSERT, \
|
|
f'RK4 reference not converged: {ref_dev}; report to parent'
|
|
mism, dp_dev = endpoint_deviation(dp_tight, rk4_02)
|
|
assert mism == 0, 'DP54 default disagrees with RK4 reference'
|
|
assert dp_dev < ENDPOINT_ASSERT, (dp_dev,)
|
|
print(f'schwarzschild: RK4 ref convergence {ref_dev:.3g}, '
|
|
f'DP default vs fine RK4 {dp_dev:.3g}', flush=True)
|
|
|
|
print(f'{backend}: adaptive CLI checks passed', flush=True)
|
|
|
|
# Alcubierre production policy: sub- and super-luminal velocities share one
|
|
# finite resource allowance B = 5*escape / max(|1-|v_s||, exp(-D)); the DP54
|
|
# default must validate, actually integrate the warp feature, and reach both
|
|
# the shared DARK terminal and escapes. All images are tiny 8x4/16x8.
|
|
alc = BUILD / 'alcubierre_sky'
|
|
if not alc.exists():
|
|
print('alcubierre: binary absent, skipping', flush=True)
|
|
else:
|
|
alc_help = run(alc, '--help').stdout
|
|
ext = 'png' if '.png' in alc_help else 'ppm'
|
|
assert '|v_s| < 1' not in alc_help, 'help still claims |v_s| < 1'
|
|
escape1 = 1.0 + 20.0 / 1.0 # R = sigma = 1
|
|
|
|
def alc_camera(vs, ra_deg, dec_deg=0.0):
|
|
return ['--observer-position', '0', '0', '0',
|
|
'--observer-velocity', repr(vs), '0', '0',
|
|
'--look-ra-deg', repr(ra_deg),
|
|
'--look-dec-deg', repr(dec_deg)]
|
|
|
|
def alc_map_prov(tag, vs, D=8.0, radius=1.0, extra=(), allow=True,
|
|
width=8, height=4, cell=4, camera=None,
|
|
catalog=True):
|
|
m = tmp / f'alc_{tag}.grlens'
|
|
args = ['--alcubierre-vs', repr(vs), '--alcubierre-radius',
|
|
repr(radius), '--alcubierre-sigma', '1',
|
|
'--dark-threshold', repr(D), '--width', str(width),
|
|
'--height', str(height), '--fov-deg', 80, '--exposure',
|
|
'1e-3', '--coarse-cell-pixels', str(cell),
|
|
'--refine-max-level', 0, '--psf-relative-tail', 1e-4]
|
|
if catalog:
|
|
args += ['--catalog', 'assets/sky_grid_5deg.csv']
|
|
if camera is not None:
|
|
args += camera
|
|
if allow:
|
|
args.append('--allow-incomplete')
|
|
args += ['--output', str(tmp / f'alc_{tag}.{ext}'),
|
|
'--lens-map-output', str(m), *extra]
|
|
run(alc, *args)
|
|
return map_provenance(m)[2], m
|
|
|
|
# 1) Ordinary sub-luminal separation: B is exactly 5*escape/sep and is
|
|
# independent of the step count and of the initial step.
|
|
p_sub, _ = alc_map_prov('sub', 0.3, camera=alc_camera(0.3, 0.0))
|
|
assert p_sub['integrator'] == 1, p_sub
|
|
assert math.isclose(p_sub['max_lookback_time'],
|
|
alcubierre_budget(escape1, 0.3, 8.0), rel_tol=1e-12)
|
|
assert p_sub['min_step'] <= p_sub['coordinate_time_step'] \
|
|
<= p_sub['max_step']
|
|
p_steps, _ = alc_map_prov('sub_steps', 0.3,
|
|
extra=('--trace-max-steps', '8'),
|
|
camera=alc_camera(0.3, 0.0))
|
|
assert p_steps['initial_max_steps'] == 8
|
|
assert math.isclose(p_steps['max_lookback_time'],
|
|
alcubierre_budget(escape1, 0.3, 8.0), rel_tol=1e-12)
|
|
p_step, _ = alc_map_prov('sub_step', 0.3,
|
|
extra=('--ode-initial-step', '0.02'),
|
|
camera=alc_camera(0.3, 0.0))
|
|
assert abs(p_step['coordinate_time_step'] - 0.02) < 1e-15
|
|
assert math.isclose(p_step['max_lookback_time'],
|
|
alcubierre_budget(escape1, 0.3, 8.0), rel_tol=1e-12)
|
|
|
|
# 2) Near-luminal (vs = 1, separation 0): the exp(-D) floor makes the
|
|
# allowance finite, and it grows with the dark threshold D.
|
|
p8, _ = alc_map_prov('near8', 1.0, D=8.0, camera=alc_camera(1.0, 0.0))
|
|
p12, _ = alc_map_prov('near12', 1.0, D=12.0, camera=alc_camera(1.0, 0.0))
|
|
assert math.isclose(p8['max_lookback_time'],
|
|
alcubierre_budget(escape1, 1.0, 8.0), rel_tol=1e-12)
|
|
assert math.isclose(p12['max_lookback_time'],
|
|
alcubierre_budget(escape1, 1.0, 12.0), rel_tol=1e-12)
|
|
assert p12['max_lookback_time'] > p8['max_lookback_time']
|
|
|
|
# 3) Both sides of the threshold-derived vcut and exactly luminal
|
|
# values use either separation or the finite floor; a tiny film traces
|
|
# and resolves without INCOMPLETE outcomes.
|
|
vcut = 1.0 - math.exp(-8.0)
|
|
for vs in (0.999, math.nextafter(vcut, 0.0),
|
|
math.nextafter(vcut, 1.0), 0.9999, math.nextafter(1.0, 0.0),
|
|
math.nextafter(1.0, 2.0)):
|
|
pv, mv = alc_map_prov(f'vcut_{vs!r}', vs, camera=alc_camera(vs, 0.0))
|
|
assert math.isclose(pv['max_lookback_time'],
|
|
alcubierre_budget(escape1, vs, 8.0),
|
|
rel_tol=1e-12), (vs, pv['max_lookback_time'])
|
|
vv, _ = map_vertices(mv)
|
|
assert 3 not in {v[10] for v in vv}, (vs, 'INCOMPLETE ray')
|
|
assert sum_rhs(vv) > 0, vs
|
|
|
|
# 4) The default camera (0,0,15 for R=5) must work with a superluminal
|
|
# bubble: the generic radius-15 camera lies inside escape radius 25.
|
|
pd, md = alc_map_prov('defaultcam', 2.0, radius=5.0)
|
|
dv, _ = map_vertices(md)
|
|
assert 3 not in {v[10] for v in dv}, 'default camera left INCOMPLETE rays'
|
|
assert sum_rhs(dv) > 0
|
|
|
|
# 5) Real small images at the critical velocities: the direction along
|
|
# the bubble motion is the DARK direction (Pi_x = sign(vs)); the
|
|
# central vertex must be DARK and some edge ray must escape. The
|
|
# 16x8 film with a 4-pixel coarse cell places a vertex exactly at the
|
|
# image center.
|
|
for vs in (1.0, -1.0, 2.0, -2.0):
|
|
ra = 0.0 if vs < 0.0 else 180.0
|
|
pimg, mimg = alc_map_prov(f'img_{vs!r}', vs, width=16, height=8,
|
|
cell=4, camera=alc_camera(vs, ra))
|
|
verts, _ = map_vertices(mimg)
|
|
outcomes = {v[10] for v in verts}
|
|
assert 3 not in outcomes, (vs, 'INCOMPLETE outcome', outcomes)
|
|
assert 0 in outcomes and 1 in outcomes, (vs, outcomes)
|
|
assert sum_rhs(verts) > 0, vs
|
|
central = min(verts, key=lambda v: (v[0] - 8.0) ** 2
|
|
+ (v[1] - 4.0) ** 2)
|
|
assert central[10] == 1, \
|
|
(vs, 'central vertex is not DARK', central[0], central[1],
|
|
central[10])
|
|
|
|
# 6) Extreme tiny-budget quota path: an explicit 4-step, 1-time
|
|
# allowance from an inside camera must stop as UNRESOLVED, never as
|
|
# a fabricated escape from a trace that took no steps.
|
|
extreme = (alc_camera(0.99999999, 0.0)
|
|
+ ['--alcubierre-vs', '0.99999999', '--alcubierre-radius',
|
|
'1', '--alcubierre-sigma', '1', '--catalog',
|
|
'assets/sky_grid_5deg.csv', '--width', '8', '--height',
|
|
'4', '--fov-deg', '80', '--exposure', '1e-3',
|
|
'--coarse-cell-pixels', '4', '--refine-max-level', 0,
|
|
'--psf-relative-tail', '1e-4'])
|
|
ext_map = tmp / 'alc_extreme.grlens'
|
|
run(alc, *extreme, '--integrator', 'dp54',
|
|
'--trace-lookback-time', '1', '--trace-max-steps', '4',
|
|
'--max-total-steps', '4', '--retry-step-increment', '0',
|
|
'--max-total-lookback-time', '1', '--retry-lookback-increment', '0',
|
|
'--allow-incomplete', '--output', tmp / f'alc_extreme.{ext}',
|
|
'--lens-map-output', ext_map)
|
|
_, _, ext_prov = map_provenance(ext_map)
|
|
assert ext_prov['integrator'] == 1
|
|
assert abs(ext_prov['max_lookback_time'] - 1.0) < 1e-15, ext_prov
|
|
assert ext_prov['initial_max_steps'] == 4, ext_prov
|
|
ext_vertices, _ = map_vertices(ext_map)
|
|
assert all(v[10] == 2 for v in ext_vertices), \
|
|
('a 1-time/4-step trace fabricated a non-UNRESOLVED outcome',
|
|
[v[10] for v in ext_vertices])
|
|
|
|
# 7) Large dark threshold (D=1000) with vs=1: exp(-1000) underflows to
|
|
# 0 and the separation is exactly 0, so the log fallback saturates
|
|
# the default allowance to DBL_MAX/4. A 4-step/4-total cap with
|
|
# retry disabled keeps the trace short while the provenance records
|
|
# the non-masked saturated default. An explicit --trace-lookback
|
|
# overrides it without masking the independent step allowance.
|
|
psat, _ = alc_map_prov('sat1000', 1.0, D=1000.0,
|
|
extra=('--trace-max-steps', '4',
|
|
'--max-total-steps', '4',
|
|
'--retry-step-increment', '0'),
|
|
camera=alc_camera(1.0, 0.0))
|
|
assert math.isclose(psat['max_lookback_time'],
|
|
sys.float_info.max / 4.0, rel_tol=1e-12), \
|
|
psat['max_lookback_time']
|
|
assert psat['initial_max_steps'] == 4
|
|
psat_ov, _ = alc_map_prov('sat1000_ov', 1.0, D=1000.0,
|
|
extra=('--trace-lookback-time', '1',
|
|
'--trace-max-steps', '4',
|
|
'--max-total-steps', '4',
|
|
'--retry-step-increment', '0'),
|
|
camera=alc_camera(1.0, 0.0))
|
|
assert psat_ov['max_lookback_time'] == 1.0, psat_ov
|
|
assert psat_ov['initial_max_steps'] == 4, psat_ov
|
|
# Explicit time override must not mask the (default) step allowance.
|
|
pind, _ = alc_map_prov('indep', 1.0, D=8.0,
|
|
extra=('--trace-lookback-time', '1'),
|
|
camera=alc_camera(1.0, 0.0))
|
|
assert pind['max_lookback_time'] == 1.0, pind
|
|
assert pind['initial_max_steps'] > 1, pind
|
|
|
|
# 8) RK4 guard: a super-luminal default (vs=2, small estimate) is no
|
|
# longer rejected; the guard triggers only when the *derived default*
|
|
# estimate exceeds the cap (D=12, vs=1), and an explicit
|
|
# --trace-max-steps is a user allowance that bypasses it.
|
|
run(alc, '--integrator', 'rk4', '--alcubierre-vs', '2',
|
|
'--alcubierre-radius', '1', '--alcubierre-sigma', '1',
|
|
'--catalog', 'assets/sky_grid_5deg.csv', '--width', '8',
|
|
'--height', '4', '--fov-deg', '80', '--exposure', '1e-3',
|
|
'--coarse-cell-pixels', '4', '--refine-max-level', 0,
|
|
'--psf-relative-tail', '1e-4', '--allow-incomplete',
|
|
'--output', tmp / f'alc_rk4_v2.{ext}')
|
|
rk4_guard = run(
|
|
alc, '--integrator', 'rk4', '--alcubierre-vs', '1',
|
|
'--alcubierre-radius', '1', '--alcubierre-sigma', '1',
|
|
'--dark-threshold', '12', '--catalog', 'assets/sky_grid_5deg.csv',
|
|
'--width', '8', '--height', '4', '--fov-deg', '80',
|
|
'--exposure', '1e-3', '--coarse-cell-pixels', '4',
|
|
'--refine-max-level', 0, '--psf-relative-tail', '1e-4',
|
|
'--output', tmp / f'alc_rk4_guard.{ext}', ok=False)
|
|
assert rk4_guard.returncode == 2, rk4_guard.stderr
|
|
assert 'cap' in rk4_guard.stderr, rk4_guard.stderr
|
|
run(alc, '--integrator', 'rk4', '--alcubierre-vs', '1',
|
|
'--alcubierre-radius', '1', '--alcubierre-sigma', '1',
|
|
'--dark-threshold', '12', '--trace-max-steps', '4',
|
|
'--catalog', 'assets/sky_grid_5deg.csv', '--width', '8',
|
|
'--height', '4', '--fov-deg', '80', '--exposure', '1e-3',
|
|
'--coarse-cell-pixels', '4', '--refine-max-level', 0,
|
|
'--psf-relative-tail', '1e-4', '--allow-incomplete',
|
|
'--output', tmp / f'alc_rk4_expl.{ext}')
|
|
|
|
# 9) Impossible parameters stay clearly rejected.
|
|
for bad in (['--alcubierre-vs', 'nan'], ['--alcubierre-vs', 'inf'],
|
|
['--alcubierre-radius', '0'], ['--alcubierre-sigma', '0']):
|
|
bad_out = tmp / f'alc_bad.{ext}'
|
|
rejected = run(alc, *bad, '--catalog', 'assets/sky_grid_5deg.csv',
|
|
'--width', '8', '--height', '4', '--fov-deg', '80',
|
|
'--exposure', '1e-3', '--coarse-cell-pixels', '4',
|
|
'--refine-max-level', 0, '--psf-relative-tail', '1e-4',
|
|
'--output', bad_out, ok=False)
|
|
assert rejected.returncode != 0, rejected.stderr
|
|
assert not bad_out.exists()
|
|
|
|
print('alcubierre: budget/provenance, vcut scans, vcut images with '
|
|
'central DARK, default camera, large-D saturation, explicit '
|
|
'overrides and RK4 guard all passed', flush=True)
|