From f7380cbf75b76ad2f2b46a80f7af8d8fe02a1c57 Mon Sep 17 00:00:00 2001 From: Yingjie Wang Date: Mon, 5 Oct 2026 06:22:47 -0400 Subject: [PATCH] Feat: Rework ray termination into escaped/dark/unresolved/incomplete Replace the position capture cutoff with a camera-relative dark threshold shared by every backend, and carry explicit outcome/reason provenance through the ray, RayPool, adaptive mesh, lens-map and replay paths. - eval/eval_slab return SpacetimePointStatus; remove SPACETIME_RAY_CAPTURED and the Schwarzschild capture radius; decouple observer construction from ray position. - RayEndpoint stores RayOutcome/RayReason plus the last trusted state; budget exhaustion is retryable UNRESOLVED, data/integration failures are INCOMPLETE. - Normal dark terminal is L - L0 >= --dark-threshold (default 8), with L0 taken at the camera event and kept distinct from the worldtube entry energy; photon energy and frequency ratio are never reset. - Implement E/D/U triangle decisions with merged budget retries, persistent probe witnesses promoted in place by vertex identity, conformity settling, and approximate-black boundary provenance with achieved-scale statistics. - Add RayPool continuation state and per-ray step budgets. - Bump lens-map to v2 with explicit end/outcome/reason, approx_black, threshold/retry/geometry provenance and per-frame retry counts; reject v1. - Gate production output on incomplete/error results, overridable with --allow-incomplete. - Update AGENTS.md, the design document and usage docs; add the termination oracle and regression coverage. make -B -j4 BUILD_TYPE=Debug test passes with bit-identical reference HDRs. --- AGENTS.md | 13 +- Makefile | 9 +- README.md | 9 +- nr_spacetime_movie_renderer_design.md | 132 ++++-- src/asymptotic.c | 20 +- src/asymptotic.h | 4 + src/frame.c | 609 +++++++++++++++++++++++--- src/frame.h | 75 +++- src/geodesic.c | 303 ++++++++++--- src/geodesic.h | 82 +++- src/lens_map.c | 96 +++- src/lens_map.h | 26 +- src/main.c | 253 +++++++++-- src/ray.c | 93 +++- src/ray.h | 15 +- src/spacetime.h | 36 +- src/spacetime_alcubierre.c | 9 +- src/spacetime_common.c | 31 +- src/spacetime_minkowski.c | 7 +- src/spacetime_schwarzschild.c | 40 +- tests/capture_psf.c | 2 +- tests/test_alcubierre.c | 10 +- tests/test_asymptotic.c | 25 +- tests/test_asymptotic_schwarzschild.c | 66 ++- tests/test_camera_cli.py | 49 ++- tests/test_frame.c | 297 ++++++++++++- tests/test_geodesic.c | 4 +- tests/test_observer.c | 8 +- tests/test_schwarzschild.c | 94 +++- tests/test_termination_oracle.c | 293 +++++++++++++ usage.md | 41 +- 31 files changed, 2349 insertions(+), 402 deletions(-) create mode 100644 tests/test_termination_oracle.c diff --git a/AGENTS.md b/AGENTS.md index 3c0e72e..136d20c 100644 --- a/AGENTS.md +++ b/AGENTS.md @@ -59,9 +59,18 @@ catalog 内部数据保留 `(direction, temperature, amplitude)`，而非 RGB。 优先从 3+1 identities、已知 gauge RHS 或 temporal interpolant 的解析导数获得时间导数；不要为已有插值量另行做低阶 finite difference，也不要在 geodesic RHS 中计算不会使用的量。 -## 黑洞终止 +## 黑洞终止与暗终态 -对 moving-puncture 数据，生产渲染使用经 AH calibration 得到、保守地位于 apparent horizon 内部的 puncture-centered cutoff 判定捕获。不要假定每次生产演化都会运行昂贵的 AH finder。可在未来加入 common-horizon 终止优化，但不得改变物理分类。 +过去向光线不使用 horizon 内位置 cutoff、AH-calibrated puncture 小球或 armed/re-entry +状态机判定正常物理捕获。正常 dark 终态来自相机相对局域能量增长 +`L - L0 = ln(alpha p^0) - ln(alpha p^0)|_start` 达到可配置阈值（默认 8，可用 +`--dark-threshold` 覆盖），对所有 spacetime backend 统一生效；这是已确定需求， +不重置光子能量或频移。不同 dark reason 不制造 mesh seam。无法可靠推进的 +积分必须报告具体数值失败，不得改写成 capture。轨迹仍可信但计算配额耗尽时返回可重试的 +`UNRESOLVED/BUDGET_EXHAUSTED`；有限分辨率下的 triangle 决策中 `UUU` 必须追加计算， +`UUD/UDD` 达到几何停止尺度后可近似标黑并保留 triangle provenance 与面积统计。 +不假定每次生产演化都会运行昂贵的 AH finder，也不依赖 capture sidecar。跨 chart、 +跨 region 或穿越视界本身不是暗终态。 ## 开发与验证顺序 diff --git a/Makefile b/Makefile index 7fb16f3..480992f 100644 --- a/Makefile +++ b/Makefile @@ -108,6 +108,7 @@ TEST_OUT_DIR := $(OBJECT_DIR)/$(HDR_BUILD_TAG) TEST_TARGET := $(TEST_OUT_DIR)/test_geodesic ASYMPTOTIC_TEST_TARGET := $(TEST_OUT_DIR)/test_asymptotic ASYMPTOTIC_SCHWARZSCHILD_TEST_TARGET := $(TEST_OUT_DIR)/test_asymptotic_schwarzschild +TERMINATION_ORACLE_TEST_TARGET := $(TEST_OUT_DIR)/test_termination_oracle FRAME_TEST_TARGET := $(TEST_OUT_DIR)/test_frame SCHWARZSCHILD_TEST_TARGET := $(TEST_OUT_DIR)/test_schwarzschild ALCUBIERRE_TEST_TARGET := $(TEST_OUT_DIR)/test_alcubierre @@ -218,6 +219,11 @@ $(FRAME_TEST_TARGET): tests/test_frame.c $(CORE_MINKOWSKI_SOURCES) $(CPU_FFTW_SO $(SCHWARZSCHILD_TEST_TARGET): tests/test_schwarzschild.c $(COMMON_SOURCES) src/spacetime_schwarzschild.c $(CPU_FFTW_SOURCES) | $(TEST_OUT_DIR) $(CC) $(CPPFLAGS) $(BUILD_CPPFLAGS) $(CFLAGS) $(BUILD_CFLAGS) $(OPENMP_FLAGS) -Isrc $^ $(LDLIBS) -o $@ +# Independent physics oracle for the termination policy (plan P0); links the +# analytic Schwarzschild backend and its exterior module. +$(TERMINATION_ORACLE_TEST_TARGET): tests/test_termination_oracle.c $(COMMON_SOURCES) src/spacetime_schwarzschild.c $(CPU_FFTW_SOURCES) | $(TEST_OUT_DIR) + $(CC) $(CPPFLAGS) $(BUILD_CPPFLAGS) $(CFLAGS) $(BUILD_CFLAGS) $(OPENMP_FLAGS) -DSPACETIME_SCHWARZSCHILD -Isrc $^ $(LDLIBS) -o $@ + $(ALCUBIERRE_TEST_TARGET): tests/test_alcubierre.c $(COMMON_SOURCES) src/spacetime_alcubierre.c $(CPU_FFTW_SOURCES) | $(TEST_OUT_DIR) $(CC) $(CPPFLAGS) $(BUILD_CPPFLAGS) $(CFLAGS) $(BUILD_CFLAGS) $(OPENMP_FLAGS) -Isrc $^ $(LDLIBS) -o $@ @@ -267,12 +273,13 @@ FAST_PSF_FFTW_TEST_DEP := FAST_PSF_FFTW_TEST_RUN := endif -test: $(CAMERA_TEST_TARGETS) $(TEST_TARGET) $(ASYMPTOTIC_TEST_TARGET) $(ASYMPTOTIC_SCHWARZSCHILD_TEST_TARGET) $(FRAME_TEST_TARGET) $(SCHWARZSCHILD_TEST_TARGET) $(ALCUBIERRE_TEST_TARGET) $(OBSERVER_TRACK_TEST_TARGET) $(CATALOG_PREFETCH_TEST_TARGET) $(FAST_PSF_FFTW_TEST_DEP) $(TONE_MAP_TEST_TARGET) $(MOVIE_OUTPUT_TEST_TARGET) $(SENSOR_BLOOM_TEST_TARGET) +test: $(CAMERA_TEST_TARGETS) $(TEST_TARGET) $(ASYMPTOTIC_TEST_TARGET) $(ASYMPTOTIC_SCHWARZSCHILD_TEST_TARGET) $(TERMINATION_ORACLE_TEST_TARGET) $(FRAME_TEST_TARGET) $(SCHWARZSCHILD_TEST_TARGET) $(ALCUBIERRE_TEST_TARGET) $(OBSERVER_TRACK_TEST_TARGET) $(CATALOG_PREFETCH_TEST_TARGET) $(FAST_PSF_FFTW_TEST_DEP) $(TONE_MAP_TEST_TARGET) $(MOVIE_OUTPUT_TEST_TARGET) $(SENSOR_BLOOM_TEST_TARGET) $(TEST_OUT_DIR)/test_observer_minkowski $(TEST_OUT_DIR)/test_observer_schwarzschild $(TEST_TARGET) $(ASYMPTOTIC_TEST_TARGET) $(ASYMPTOTIC_SCHWARZSCHILD_TEST_TARGET) + $(TERMINATION_ORACLE_TEST_TARGET) $(FRAME_TEST_TARGET) $(SCHWARZSCHILD_TEST_TARGET) $(ALCUBIERRE_TEST_TARGET) diff --git a/README.md b/README.md index 263eb9d..ba75669 100644 --- a/README.md +++ b/README.md @@ -39,7 +39,7 @@ HIP retains parallel CPU catalog mapping and uses bounded, completion-protected event uploads. See [HIP configuration and bounded performance checks](build.md#optional-hip-psf-acceleration). The [Nmesh](https://github.com/nmeshsource/nmesh) numerical-spacetime backend and BBH rendering are still planned. -The current scope is black-hole capture and distant stellar backgrounds; +The current scope is black-hole shadows and distant stellar backgrounds; local matter emission, accretion disks, and plasma are outside this stage. See the [design document](nr_spacetime_movie_renderer_design.md) for the architecture and development roadmap. @@ -254,8 +254,11 @@ derivatives. Defaults: `--rtol 1e-10 --atol 1e-12 --stop-radius 0.001`. Integration crosses the horizon and stops at this numerical guard before the singularity, reporting its proper time and retaining only regular cadence samples. The guard is not the exact singularity; reduce it and tolerances to check convergence. -The renderer's independent ray capture cutoff remains `r=1.5M`: rows inside it are -valid trajectory data but the current renderer captures those rays immediately. +The renderer no longer uses a position capture cutoff: cameras at and inside the +old `r=1.5M` guard are valid targets, and a normal dark pixel comes from the +redshift-threshold truncation `log(alpha p^0) >= 8`. Budget-exhausted and +data/integration failures are separate unresolved/incomplete outcomes and are not +silently rendered as dark. The script reports maximum tetrad drift and rejects errors above `1e-6` rather than silently repairing the transported frame. Run the orbit, transport and CSV render regressions after building with `python3 tests/test_schwarzschild_camera_track.py`. diff --git a/nr_spacetime_movie_renderer_design.md b/nr_spacetime_movie_renderer_design.md index 2e5efa9..abcfe3f 100644 --- a/nr_spacetime_movie_renderer_design.md +++ b/nr_spacetime_movie_renderer_design.md @@ -14,10 +14,9 @@ - 让平直时空、解析时空、数值时空在同一渲染框架中作为可替换 backend； - 最终能够“看到”每次 NR 代码实际跑出来的时空，而不是只看 waveform 或标量诊断。 -第一阶段不考虑物质辐射、吸积盘、流体、等离子体等局域发射源。每条 ray 的终点暂时只有两类： - -1. 被黑洞捕获； -2. 到达无穷远天球。 +第一阶段不考虑物质辐射、吸积盘、流体、等离子体等局域发射源。每条 ray 的正常终点 +为：逃逸到某个无穷远天球，或达到红移暗阈值。预算耗尽与数据/积分失败是单独的 +`UNRESOLVED` / `INCOMPLETE` 类别，不与物理暗终态混同（见 §18）。 --- @@ -79,7 +78,7 @@ C++ 并不是当前项目的必要条件。需要的抽象主要可以通过： 4. 将这些 rays 组成一个全局 `RayPool`； 5. 从视频结束时刻向过去，按 time slab 顺序加载数值时空； 6. 在每个 slab 内，把所有 active rays 一起推进到 slab 左边界； -7. ray 若到达无穷远或进入黑洞，则立即终止； +7. ray 若在某个渐近端逃逸、达到红移暗阈值、耗尽预算或遇到数据/积分失败，则按类别终止； 8. 一轮 ray tracing 完成后，把 endpoint 数据回填到各帧 image mesh； 9. 根据局部 lens mapping 误差判断哪些 image-plane triangles 需要进一步细分； 10. 生成下一批新增 rays； @@ -116,7 +115,7 @@ Spacetime time-slab stream │ ▼ ray endpoint: -n∞, frequency shift, captured/escaped +n∞, frequency shift, end/outcome/reason │ ├──────────► next refinement pass │ @@ -222,8 +221,8 @@ F^{-1}:\ \hat n_\infty \to (x,y)_\text{image} 每个 image-plane triangle 的三个顶点都保存： - image-plane 坐标 `(x,y)`； -- ray 是否 escaped/captured； -- 若 escaped：无穷远方向 `n_inf`； +- ray 的终态类别 `ESCAPED`/`DARK`/`UNRESOLVED`/`INCOMPLETE` 及其 reason； +- 若 escaped：无穷远方向 `n_inf` 与所属 end； - frequency shift / redshift accumulator。 示意： @@ -235,7 +234,9 @@ typedef struct { double n_inf[3]; double log_g; - uint8_t ray_status; + uint8_t outcome; /* ESCAPED / DARK / UNRESOLVED / INCOMPLETE */ + uint8_t reason; /* REDSHIFT_LIMIT / BUDGET_EXHAUSTED / ... */ + uint32_t end_id; } LensVertex; typedef struct { @@ -417,8 +418,9 @@ $e'_3=-\sin\rho\,e_2+\cos\rho\,e_3$ 定义。 observer 构造只接收当地 metric 和已补全的参数，不加载 slab、不分类 ray。 调用方在昂贵的 catalog/PSF 初始化前验证相机及 backend 数据域。 -现有 Schwarzschild cutoff 为 $r=1.5M$；相机必须在 cutoff 外，但允许在视界内。 -此功能不改变捕获 cutoff 或向过去追踪的高红移终止条件。 +相机合法性只由 metric 可用性、四速度 timelike、时间定向和 tetrad 正交归一决定； +视界内、旧 cutoff 内的相机都是正常渲染目标，位置本身不决定 ray 终态。 +向过去追踪的高红移阈值截断仍正常生效（见 §18）。 单张相机参数与轨迹输入、lens-map 导入互斥；导入仍跳过 metric 与 observer 初始化。 验证包括 tetrad 正交归一和 null 初始化、平直时空平移不变性与解析光行差/多普勒、 @@ -708,27 +710,72 @@ slab 边界需要少量 overlapping temporal ghost slices。 --- -# 18. 黑洞捕获判据 +# 18. 黑洞终止与暗终态 -目标使用 moving-puncture BBH，而不是 excision。 +目标使用 moving-puncture BBH，而不是 excision。正常终态**不使用** horizon 内位置 +cutoff、AH-calibrated puncture 小球或 armed/re-entry 状态机判定物理捕获。过去向光线 +围绕渐近端逃逸、能量阈值截断和经可靠识别的渐近轨道组织；达到红移阈值后停止、渲染 +为黑，是已确定需求。 -production renderer 不希望依赖每次 NR run 都开启昂贵的 AH finder。 - -计划： - -1. 用低分辨率 single-BH / BBH calibration run 开 AH finder； -2. 测量 horizon 相对于 puncture 的最小 coordinate radius； -3. 选择明显保守、始终位于 AH 内部的 puncture-centered cutoff； -4. 正式 renderer 只根据 puncture trajectory 做判断。 - -形式： +形式（相机相对局域能量增长，对全部 backend 统一；具体阈值通过小型 +oracle/convergence test 标定，不宣称由论文给定）： \[ -|\mathbf x-\mathbf x_p(t)|= L_dark`（默认 `L_dark=8`，可用 + `--dark-threshold` 覆盖，对全部 spacetime backend 统一生效）。计算配额耗尽 + 返回可重试的 `UNRESOLVED/BUDGET_EXHAUSTED`， + 保留最后可信连续状态；数据/积分/历史耗尽返回具体 `INCOMPLETE` reason。 +- `eval`/`eval_slab` 返回 `SpacetimePointStatus`（时间不足、域外、invalid + metric、内部错误）；observer 合法性只由 metric 可用性、timelike 四速度、 + 时间定向和 tetrad 正交归一决定，不再调用位置分类。 +- 三角形决策实现 $E/D/U$ 账目：`UUU` 与含逃逸顶点的未决组合强制追加预算重试； + `UUD/UDD` 在几何停止尺度处近似标黑并记录 image-plane 面积与 triangle + provenance；含 `INCOMPLETE` 的三角形计为错误，不参与近似标黑。总资源上限 + 耗尽而未决时报告 incomplete，诊断可用 `--allow-incomplete` 覆盖。 +- lens-map 文件格式升级为 v2：显式存储 `end_id/outcome/reason`、triangle + `approx_black` 以及 policy/retry/几何阈值 provenance；v1 文件被明确拒绝。 ## 18A.11 $M>0$ 解析外区：闭式约化与验证 @@ -1707,7 +1769,9 @@ void rays_trace_generation( - metric interpolation； - spatial derivatives； - 必要的 temporal derivatives； -- capture/infinity classification。 +- metric/data 状态（成功、时间不足、空间域不足、invalid metric、内部错误）； +- 渐近端声明与 escape worldtube 能力。物理暗终态由能量阈值 policy 判定，不由 + 位置分类决定。 --- @@ -1811,7 +1875,7 @@ map 共用同一入口。每个 RGB 通道独立、各向同性地把超过有 验证： -- capture； +- 红移暗阈值截断与 shadow； - Einstein ring； - multiple images； - adaptive refinement； @@ -1839,7 +1903,7 @@ map 共用同一入口。每个 RGB 通道独立、各向同性地把超过有 \quad g, \quad -\text{captured/escaped classification} +\text{end/outcome classification} \] --- @@ -1868,7 +1932,7 @@ renderer 顶层架构原则上不应为 BBH 重新设计。 - time slab 最佳内存大小； - BBH production node 数； - 4D 输出总数据量； -- AH calibration 后 puncture cutoff； +- 相机相对能量阈值 `L-L0` 的默认值标定（当前 `L_dark=8`，可 CLI 覆盖）； - adaptive triangle refinement criterion； - critical curve 附近最大 refinement level； - Gaia 与 2MASS 的最终组合； @@ -1879,7 +1943,9 @@ renderer 顶层架构原则上不应为 BBH 重新设计。 Reinhard 作为兼容模式保留；最终 production color management、传感器模型、 曝光标定和 HDR 视频编码规则仍未决定； - 是否需要 diffuse Milky Way background； -- 是否将 ray redshift 变量定义为 `log(alpha p^0)` 或其他更方便的量。 +- 阈值监测量当前统一采用相机相对增长 `L-L0`（对全部 backend 生效），不再使用 + 绝对 `L`、`ln(p^0)` 或 Killing 相对量；这些量均不得与真正的 infinity + `g=E_camera/E_source` 混同。 --- @@ -1917,8 +1983,8 @@ $e_{(0)}=u$ 同时满足自由落体方程；四加速度为零时费米–沃 DOP853 默认 rtol=$10^{-10}$、atol=$10^{-12}$，可配置并通过解析径向自由落体、 圆轨道和圆轨道平行输运的收敛回归验证。视界不终止相机；默认 $r=10^{-3}M$ -只是可配置的奇点数值保护边界，不等于精确撞击奇点。它独立于光线的 $1.5M$ -捕获 cutoff；该 cutoff 内的轨迹可输出，但当前 renderer 的光线会立即被捕获。 +只是可配置的奇点数值保护边界，不等于精确撞击奇点。相机轨迹与光线终态解耦： +相机合法性与位置 cutoff 无关，光线正常暗终态由 §18 的能量阈值 policy 决定。 CSV 保持 21 列不变，记录 $\tau=k/\mathrm{fps}$ 与积分得到的真实坐标时间。 只保留不超过请求持续本征时或提前终止时刻的规则采样，包含 $\tau=0$。 diff --git a/src/asymptotic.c b/src/asymptotic.c index c7e372c..8bc9b6a 100644 --- a/src/asymptotic.c +++ b/src/asymptotic.c @@ -539,6 +539,7 @@ AsymptoticStatus asymptotic_route_camera(const SpacetimeSource *source, route->Pi[i] = state.Pi[i]; } route->log_alpha_p0 = state.log_alpha_p0; + route->log_alpha_p0_camera = state.log_alpha_p0; return ASYMPTOTIC_OK; } /* A backend that declares ends must describe them consistently and use a @@ -587,6 +588,7 @@ AsymptoticStatus asymptotic_route_camera(const SpacetimeSource *source, route->Pi[k] = state.Pi[k]; } route->log_alpha_p0 = state.log_alpha_p0; + route->log_alpha_p0_camera = state.log_alpha_p0; return ASYMPTOTIC_OK; } } @@ -603,8 +605,13 @@ AsymptoticStatus asymptotic_route_camera(const SpacetimeSource *source, return ASYMPTOTIC_INVALID; if (first_end == SPACETIME_END_NONE) first_end = end.end_id; - if (end.exterior_kind == ASYMPTOTIC_EXTERIOR_SCHWARZSCHILD_MONOPOLE) - return schwarzschild_route(source, &end, &metric, &state, route); + if (end.exterior_kind == ASYMPTOTIC_EXTERIOR_SCHWARZSCHILD_MONOPOLE) { + const AsymptoticStatus status = + schwarzschild_route(source, &end, &metric, &state, route); + if (status == ASYMPTOTIC_OK) + route->log_alpha_p0_camera = state.log_alpha_p0; + return status; + } if (end.exterior_kind != ASYMPTOTIC_EXTERIOR_MINKOWSKI) return ASYMPTOTIC_UNSUPPORTED; if (asymptotic_canonical_from_backend(source, end.end_id, &metric, @@ -639,6 +646,7 @@ AsymptoticStatus asymptotic_route_camera(const SpacetimeSource *source, if (have_entry) { *route = best; route->log_alpha_p0 = state.log_alpha_p0; + route->log_alpha_p0_camera = state.log_alpha_p0; return ASYMPTOTIC_OK; } if (!have_miss) { @@ -678,8 +686,10 @@ AsymptoticStatus asymptotic_finish_escape(const SpacetimeSource *source, endpoint->n_infinity[i] = n_inf[i]; endpoint->frequency_ratio = frequency; endpoint->end_id = end_id; - endpoint->status = RAY_ENDPOINT_ESCAPED; + endpoint->outcome = RAY_OUTCOME_ESCAPED; + endpoint->reason = RAY_REASON_NONE; endpoint->magnification = 1.0; + endpoint->threshold_value = NAN; return ASYMPTOTIC_OK; } if (end.exterior_kind != ASYMPTOTIC_EXTERIOR_MINKOWSKI) @@ -704,7 +714,9 @@ AsymptoticStatus asymptotic_finish_escape(const SpacetimeSource *source, endpoint->n_infinity[i] = n[i]; endpoint->frequency_ratio = 1.0 / energy; endpoint->end_id = end_id; - endpoint->status = RAY_ENDPOINT_ESCAPED; + endpoint->outcome = RAY_OUTCOME_ESCAPED; + endpoint->reason = RAY_REASON_NONE; endpoint->magnification = 1.0; + endpoint->threshold_value = NAN; return ASYMPTOTIC_OK; } diff --git a/src/asymptotic.h b/src/asymptotic.h index c63ca39..8ab6e4f 100644 --- a/src/asymptotic.h +++ b/src/asymptotic.h @@ -40,7 +40,11 @@ typedef struct { double activate_t; double x[3]; double Pi[3]; + /* Current L at the activation event (camera when inside, entry event when + * externing). `log_alpha_p0_camera` is the reference L at the camera event + * used by the camera-relative dark threshold, and must be kept separate. */ double log_alpha_p0; + double log_alpha_p0_camera; /* Terminal infinity endpoint for ESCAPED. */ double n_infinity[3]; double frequency_ratio; diff --git a/src/frame.c b/src/frame.c index d1706a8..46a3304 100644 --- a/src/frame.c +++ b/src/frame.c @@ -89,9 +89,9 @@ int frame_lens_mesh_build_coarse(FrameLensMesh *mesh, int width, int height, const size_t bottom_left = vertex_index(column, row + 1, columns); const size_t bottom_right = vertex_index(column + 1, row + 1, columns); triangles[next_triangle++] = - (LensTriangle){{top_left, bottom_left, bottom_right}, 0, 0}; + (LensTriangle){{top_left, bottom_left, bottom_right}, 0, 0, 0}; triangles[next_triangle++] = - (LensTriangle){{top_left, bottom_right, top_right}, 0, 0}; + (LensTriangle){{top_left, bottom_right, top_right}, 0, 0, 0}; } *mesh = (FrameLensMesh){.vertices = vertices, .triangles = triangles, @@ -102,6 +102,38 @@ int frame_lens_mesh_build_coarse(FrameLensMesh *mesh, int width, int height, return 0; } +/* Store an endpoint into a lens vertex. For an UNRESOLVED result the last + * accepted continuous state is kept so the ray can be resumed; `granted_limit` + * is the total accepted-step budget that produced this result (0 when it is + * the base trace config). */ +static void store_endpoint(LensVertex *vertex, const RayEndpoint *endpoint, + unsigned int granted_limit) { + vertex->outcome = endpoint->outcome; + vertex->reason = endpoint->reason; + vertex->end_id = endpoint->end_id; + vertex->traced = 1; + if (endpoint->outcome == RAY_OUTCOME_ESCAPED) { + for (int axis = 0; axis < 3; ++axis) + vertex->n_infinity[axis] = endpoint->n_infinity[axis]; + vertex->log_frequency_ratio = log(endpoint->frequency_ratio); + return; + } + if (endpoint->outcome == RAY_OUTCOME_UNRESOLVED) { + vertex->continuation_t = endpoint->stop_coordinate_time; + for (int axis = 0; axis < 3; ++axis) { + vertex->continuation_x[axis] = endpoint->final_x[axis]; + vertex->continuation_Pi[axis] = endpoint->final_Pi[axis]; + } + vertex->continuation_log_alpha_p0 = endpoint->final_log_alpha_p0; + vertex->continuation_log_alpha_p0_0 = endpoint->final_log_alpha_p0_0; + vertex->continuation_steps = endpoint->accepted_steps; + const unsigned int used = + granted_limit != 0 ? granted_limit : endpoint->accepted_steps; + if (used > vertex->continuation_limit) + vertex->continuation_limit = used; + } +} + int frame_lens_mesh_trace(FrameLensMesh *mesh, const SpacetimeSource *spacetime, const ObserverState *observer, const GeodesicTraceConfig *trace) { @@ -115,14 +147,7 @@ int frame_lens_mesh_trace(FrameLensMesh *mesh, const SpacetimeSource *spacetime, LensVertex *vertex = &mesh->vertices[i]; RayEndpoint endpoint = geodesic_trace_past(spacetime, observer, vertex->camera_direction, trace); - vertex->status = endpoint.status; - vertex->end_id = endpoint.end_id; - vertex->traced = 1; - if (endpoint.status == RAY_ENDPOINT_ESCAPED) { - for (int axis = 0; axis < 3; ++axis) - vertex->n_infinity[axis] = endpoint.n_infinity[axis]; - vertex->log_frequency_ratio = log(endpoint.frequency_ratio); - } + store_endpoint(vertex, &endpoint, 0); } return 0; } @@ -436,24 +461,54 @@ static int all_vertices_traced(const FrameLensMesh *mesh) { return 1; } +typedef struct { + int e, d, u, bad; +} VertexMix; + +static VertexMix triangle_mix(const FrameLensMesh *mesh, + const LensTriangle *triangle) { + VertexMix mix = {0, 0, 0, 0}; + for (int i = 0; i < 3; ++i) { + switch (mesh->vertices[triangle->vertex[i]].outcome) { + case RAY_OUTCOME_ESCAPED: ++mix.e; break; + case RAY_OUTCOME_DARK: ++mix.d; break; + case RAY_OUTCOME_UNRESOLVED: ++mix.u; break; + default: ++mix.bad; break; + } + } + return mix; +} + +/* A discontinuous boundary that must be red-refined. Escape/dark and + * unresolved/dark (with no escape vertex) are boundaries; errors are not. */ static int terminal_mismatch(const LensVertex *a, const LensVertex *b, const LensVertex *c) { - int escaped = 0, captured = 0, have_end = 0; + int escaped = 0, dark = 0, unresolved = 0, bad = 0, have_end = 0; SpacetimeEndId end = SPACETIME_END_NONE; const LensVertex *vertices[] = {a, b, c}; for (size_t i = 0; i < 3; ++i) { - escaped |= vertices[i]->status == RAY_ENDPOINT_ESCAPED; - captured |= vertices[i]->status == RAY_ENDPOINT_CAPTURED; - if (vertices[i]->status == RAY_ENDPOINT_ESCAPED) { + switch (vertices[i]->outcome) { + case RAY_OUTCOME_ESCAPED: + ++escaped; if (!have_end) { end = vertices[i]->end_id; have_end = 1; } else if (vertices[i]->end_id != end) { return 1; /* two different infinity ends must not be interpolated */ } + break; + case RAY_OUTCOME_DARK: ++dark; break; + case RAY_OUTCOME_UNRESOLVED: ++unresolved; break; + default: ++bad; break; } } - return escaped && captured; + if (bad > 0) + return 0; + if (escaped > 0 && dark > 0) + return 1; + if (unresolved > 0 && dark > 0 && escaped == 0) + return 1; + return 0; } static int add_sample(FrameLensMesh *mesh, const FrameSample *sample) { @@ -512,6 +567,94 @@ static void index_probe(FrameLensMesh *mesh, size_t sample_id) { mesh->probe_slots[slot] = sample_id + 1; } +/* Rebuild the persistent witness index from the mesh vertices: a witness is a + * vertex flagged diagnostic_probe that no triangle references. Consumed + * witnesses (now formal midpoints) are un-flagged here. */ +static int witness_rebuild(FrameLensMesh *mesh) { + unsigned char *used = calloc(mesh->vertex_count ? mesh->vertex_count : 1, 1); + if (used == NULL) + return -1; + for (size_t t = 0; t < mesh->triangle_count; ++t) + for (int j = 0; j < 3; ++j) + if (mesh->triangles[t].vertex[j] < mesh->vertex_count) + used[mesh->triangles[t].vertex[j]] = 1; + size_t count = 0; + for (size_t v = 0; v < mesh->vertex_count; ++v) { + if (!mesh->vertices[v].diagnostic_probe) + continue; + if (used[v]) { + mesh->vertices[v].diagnostic_probe = 0; /* promoted to a midpoint */ + continue; + } + if (count == mesh->witness_capacity) { + size_t cap = mesh->witness_capacity ? mesh->witness_capacity * 2 : 8; + size_t *list = realloc(mesh->witness_vertices, cap * sizeof *list); + if (list == NULL) { + free(used); + return -1; + } + mesh->witness_vertices = list; + mesh->witness_capacity = cap; + } + mesh->witness_vertices[count++] = v; + } + mesh->witness_count = count; + mesh->diagnostic_probe_count = count; + size_t cap = 16; + while (cap < count * 2) + cap *= 2; + if (cap > mesh->witness_slot_capacity) { + size_t *slots = realloc(mesh->witness_slots, cap * sizeof *slots); + if (slots == NULL) { + free(used); + return -1; + } + mesh->witness_slots = slots; + mesh->witness_slot_capacity = cap; + } + if (mesh->witness_slot_capacity != 0) + memset(mesh->witness_slots, 0, + mesh->witness_slot_capacity * sizeof *mesh->witness_slots); + for (size_t i = 0; i < count; ++i) { + const size_t id = mesh->witness_vertices[i]; + const size_t a = mesh->vertices[id].probe_edge[0]; + const size_t b = mesh->vertices[id].probe_edge[1]; + size_t slot = probe_hash(a, b) & (mesh->witness_slot_capacity - 1); + while (mesh->witness_slots[slot] != 0) + slot = (slot + 1) & (mesh->witness_slot_capacity - 1); + mesh->witness_slots[slot] = i + 1; + } + free(used); + return 0; +} + +static size_t witness_find(const FrameLensMesh *mesh, size_t a, size_t b) { + if (a > b) { const size_t swap = a; a = b; b = swap; } + if (mesh->witness_slot_capacity == 0) + return SIZE_MAX; + size_t slot = probe_hash(a, b) & (mesh->witness_slot_capacity - 1); + while (mesh->witness_slots[slot] != 0) { + const size_t id = mesh->witness_vertices[mesh->witness_slots[slot] - 1]; + if (mesh->vertices[id].probe_edge[0] == a && + mesh->vertices[id].probe_edge[1] == b) + return id; + slot = (slot + 1) & (mesh->witness_slot_capacity - 1); + } + return SIZE_MAX; +} + +static unsigned int retry_limit(unsigned int current, + const RefinementConfig *config) { + const unsigned int remaining = config->max_total_steps - current; + return current + (config->retry_step_increment < remaining + ? config->retry_step_increment + : remaining); +} + +static int triangle_allows_children(const FrameLensMesh *mesh, + const LensTriangle *triangle, + const RefinementConfig *config); + int frame_lens_mesh_prepare_generation(FrameLensMesh *mesh, const RefinementConfig *config) { if (mesh == NULL || config == NULL || mesh->sample_count != 0) @@ -528,16 +671,96 @@ int frame_lens_mesh_prepare_generation(FrameLensMesh *mesh, .vertex = mesh->vertices[i]})) return -1; } + int emitted_probes = 0; + /* Retry pass: merge one request per physical sample id. A triangle vertex + * that is UNRESOLVED with an escape side (or an all-U triangle) is retried, + * as is an off-mesh unresolved witness. A witness retry carries its edge so + * the refinement decision sees the updated state; both passes share + * retry_seen so a promoted witness is never requested twice. */ + if (config->retry_step_increment > 0 && mesh->vertex_count > 0) { + unsigned char *retry_seen = calloc(mesh->vertex_count, 1); + if (retry_seen == NULL) + return -1; + for (size_t i = 0; i < mesh->triangle_count; ++i) { + const LensTriangle *triangle = &mesh->triangles[i]; + const VertexMix mix = triangle_mix(mesh, triangle); + if (mix.bad > 0) + continue; + if (!((mix.u > 0 && mix.e > 0) || mix.u == 3)) + continue; + for (int corner = 0; corner < 3; ++corner) { + const size_t v = triangle->vertex[corner]; + LensVertex *vertex = &mesh->vertices[v]; + if (vertex->outcome != RAY_OUTCOME_UNRESOLVED || retry_seen[v]) + continue; + retry_seen[v] = 1; + if (vertex->continuation_limit >= config->max_total_steps) + continue; /* capped: reported as budget-incomplete, not retried */ + const unsigned int limit = retry_limit(vertex->continuation_limit, config); + if (limit <= vertex->continuation_limit) + continue; + FrameSample retry = {.kind = FRAME_SAMPLE_RETRY, + .vertex_id = v, + .step_limit = limit, + .vertex = *vertex}; + if (add_sample(mesh, &retry)) { + free(retry_seen); + return -1; + } + ++mesh->retry_requests; + } + } + for (size_t i = 0; i < mesh->witness_count; ++i) { + const size_t v = mesh->witness_vertices[i]; + LensVertex *vertex = &mesh->vertices[v]; + if (retry_seen[v] || vertex->outcome != RAY_OUTCOME_UNRESOLVED) + continue; + retry_seen[v] = 1; + if (vertex->continuation_limit >= config->max_total_steps) + continue; + const unsigned int limit = retry_limit(vertex->continuation_limit, config); + if (limit <= vertex->continuation_limit) + continue; + FrameSample retry = {.kind = FRAME_SAMPLE_RETRY, .vertex_id = v, + .edge_vertex = {vertex->probe_edge[0], + vertex->probe_edge[1]}, + .step_limit = limit, .vertex = *vertex}; + if (add_sample(mesh, &retry)) { + free(retry_seen); + return -1; + } + if (mesh->probe_slot_capacity != 0) + index_probe(mesh, mesh->sample_count - 1); + emitted_probes = 1; + ++mesh->retry_requests; + } + free(retry_seen); + } if (config->max_level == 0) /* Coarse vertices still need tracing when refinement is disabled. */ return (int)mesh->sample_count; /* Every generation may batch newly inserted vertices with probes for its * new leaves: probe positions depend only on image-plane geometry. Their - * endpoints are considered only after this complete generation finishes. */ + * endpoints are considered only after this complete generation finishes. + * Unresolved/error triangles are handled by the retry pass or the boundary + * accounting, so they request no probes here. */ for (size_t i = 0; i < mesh->triangle_count; ++i) { const LensTriangle *triangle = &mesh->triangles[i]; if (triangle->level >= config->max_level || triangle->evaluated) continue; + const VertexMix mix = triangle_mix(mesh, triangle); + if (mix.bad > 0) + continue; + if (mix.u > 0) { + /* UUD/UDD may still be red-refined to locate the boundary, but only + * while the geometry can support children; otherwise it is an + * approximate-black boundary decision. U+escape and UUU are handled by + * the retry pass instead. */ + const int dark_side_only = mix.e == 0 && mix.d > 0; + if (!dark_side_only || + !triangle_allows_children(mesh, triangle, config)) + continue; + } const unsigned int first_side = longest_side(mesh, triangle); const unsigned int side_count = terminal_mismatch( &mesh->vertices[triangle->vertex[0]], @@ -551,6 +774,21 @@ int frame_lens_mesh_prepare_generation(FrameLensMesh *mesh, const size_t b = triangle->vertex[(side + 1) % 3]; if (find_probe(mesh, a, b) != SIZE_MAX) continue; + /* A persistent witness (terminal or capped unresolved) is reused in + * place: no retrace and, if its edge is later split, no duplicate. */ + const size_t wid = witness_find(mesh, a, b); + if (wid != SIZE_MAX) { + FrameSample cached = {.kind = FRAME_SAMPLE_PROBE, + .vertex_id = wid, + .edge_vertex = {a, b}, + .cached = 1, + .vertex = mesh->vertices[wid]}; + if (add_sample(mesh, &cached)) + return -1; + index_probe(mesh, mesh->sample_count - 1); + emitted_probes = 1; + continue; + } FrameSample probe = {.kind = FRAME_SAMPLE_PROBE, .edge_vertex = {a, b}}; const LensVertex *left = &mesh->vertices[a]; const LensVertex *right = &mesh->vertices[b]; @@ -564,9 +802,10 @@ int frame_lens_mesh_prepare_generation(FrameLensMesh *mesh, if (add_sample(mesh, &probe)) return -1; index_probe(mesh, mesh->sample_count - 1); + emitted_probes = 1; } } - mesh->samples_include_probes = mesh->sample_count != 0; + mesh->samples_include_probes = emitted_probes; return (int)mesh->sample_count; } @@ -582,17 +821,12 @@ int frame_lens_mesh_install_sample(FrameLensMesh *mesh, size_t sample_id, if (mesh == NULL || endpoint == NULL || sample_id >= mesh->sample_count) return -1; FrameSample *sample = &mesh->samples[sample_id]; - LensVertex *vertex = sample->kind == FRAME_SAMPLE_VERTEX - ? &mesh->vertices[sample->vertex_id] - : &sample->vertex; - vertex->status = endpoint->status; - vertex->end_id = endpoint->end_id; - vertex->traced = 1; - if (endpoint->status == RAY_ENDPOINT_ESCAPED) { - for (int axis = 0; axis < 3; ++axis) - vertex->n_infinity[axis] = endpoint->n_infinity[axis]; - vertex->log_frequency_ratio = log(endpoint->frequency_ratio); - } + LensVertex *vertex = sample->kind == FRAME_SAMPLE_PROBE + ? &sample->vertex + : &mesh->vertices[sample->vertex_id]; + store_endpoint(vertex, endpoint, sample->step_limit); + if (sample->kind == FRAME_SAMPLE_RETRY) + sample->vertex = *vertex; /* witness retries are read back by find_probe */ return 0; } @@ -618,8 +852,8 @@ static int discrete_jacobian(const FrameLensMesh *mesh, const LensVertex *a = &mesh->vertices[triangle->vertex[0]]; const LensVertex *b = &mesh->vertices[triangle->vertex[1]]; const LensVertex *c = &mesh->vertices[triangle->vertex[2]]; - if (a->status != RAY_ENDPOINT_ESCAPED || b->status != RAY_ENDPOINT_ESCAPED || - c->status != RAY_ENDPOINT_ESCAPED) + if (a->outcome != RAY_OUTCOME_ESCAPED || b->outcome != RAY_OUTCOME_ESCAPED || + c->outcome != RAY_OUTCOME_ESCAPED) return 0; if (a->end_id != b->end_id || a->end_id != c->end_id) return 0; @@ -651,16 +885,22 @@ static int probe_requires_split(const FrameLensMesh *mesh, if (probe_id == SIZE_MAX) return 0; const LensVertex *probe = &mesh->samples[probe_id].vertex; - if ((probe->status == RAY_ENDPOINT_ESCAPED) != - (a->status == RAY_ENDPOINT_ESCAPED) || - (probe->status == RAY_ENDPOINT_ESCAPED) != - (b->status == RAY_ENDPOINT_ESCAPED)) - return (probe->status == RAY_ENDPOINT_ESCAPED || - a->status == RAY_ENDPOINT_ESCAPED || b->status == RAY_ENDPOINT_ESCAPED) && - (probe->status == RAY_ENDPOINT_CAPTURED || - a->status == RAY_ENDPOINT_CAPTURED || b->status == RAY_ENDPOINT_CAPTURED); - if (a->status != RAY_ENDPOINT_ESCAPED || b->status != RAY_ENDPOINT_ESCAPED || - probe->status != RAY_ENDPOINT_ESCAPED) + if (probe->outcome == RAY_OUTCOME_INCOMPLETE || + probe->outcome == RAY_OUTCOME_UNRESOLVED) + return 0; /* retained as a witness and retried, not a mapping estimate */ + if ((probe->outcome == RAY_OUTCOME_ESCAPED) != + (a->outcome == RAY_OUTCOME_ESCAPED) || + (probe->outcome == RAY_OUTCOME_ESCAPED) != + (b->outcome == RAY_OUTCOME_ESCAPED)) + return (probe->outcome == RAY_OUTCOME_ESCAPED || + a->outcome == RAY_OUTCOME_ESCAPED || + b->outcome == RAY_OUTCOME_ESCAPED) && + (probe->outcome == RAY_OUTCOME_DARK || + a->outcome == RAY_OUTCOME_DARK || + b->outcome == RAY_OUTCOME_DARK); + if (a->outcome != RAY_OUTCOME_ESCAPED || + b->outcome != RAY_OUTCOME_ESCAPED || + probe->outcome != RAY_OUTCOME_ESCAPED) return 0; double predicted[3] = {a->n_infinity[0] + b->n_infinity[0], a->n_infinity[1] + b->n_infinity[1], @@ -686,7 +926,72 @@ static int triangle_allows_children(const FrameLensMesh *mesh, const double edge = fmax(image_edge_length(a, b), fmax(image_edge_length(b, c), image_edge_length(c, a))); const double area = image_triangle_area(a, b, c); - return edge > config->min_edge_pixels && area > config->min_area_pixels2; + return triangle->level < config->max_level && + edge > config->min_edge_pixels && area > config->min_area_pixels2; +} + +void frame_lens_mesh_boundary_stats(FrameLensMesh *mesh, + const RefinementConfig *config, + FrameBoundaryStats *stats) { + if (stats == NULL) + return; + *stats = (FrameBoundaryStats){0}; + if (mesh == NULL || config == NULL) + return; + /* Off-mesh probes are samples too. Their failures cannot disappear just + * because no inverse patch uses them. Infer orphanhood for replay as well. */ + unsigned char *used = calloc(mesh->vertex_count, 1); + if (used == NULL) { ++stats->error; return; } + for (size_t i = 0; i < mesh->triangle_count; ++i) + for (int j = 0; j < 3; ++j) used[mesh->triangles[i].vertex[j]] = 1; + for (size_t i = 0; i < mesh->vertex_count; ++i) { + if (!used[i] && mesh->vertices[i].outcome == RAY_OUTCOME_INCOMPLETE) + ++stats->error; + if (!used[i] && mesh->vertices[i].outcome == RAY_OUTCOME_UNRESOLVED) + ++stats->budget_incomplete_triangles; + } + free(used); + for (size_t i = 0; i < mesh->triangle_count; ++i) { + LensTriangle *triangle = &mesh->triangles[i]; + triangle->approx_black = 0; + const VertexMix mix = triangle_mix(mesh, triangle); + if (mix.bad > 0) { + ++stats->error; + continue; + } + const int unresolved_with_escape = mix.u > 0 && mix.e > 0; + const int all_unresolved = mix.u == 3; + if (mix.u == 0) { + if (mix.d == 0) + ++stats->escaped_only; + else if (mix.e == 0) + ++stats->dark_only; + else + ++stats->eed_edd; + } else if (unresolved_with_escape) { + ++stats->u_with_escape; + } else if (all_unresolved) { + ++stats->uuu; + } else { + ++stats->uud_udd; /* UUD / UDD */ + if (!triangle_allows_children(mesh, triangle, config)) { + triangle->approx_black = 1; + ++stats->approx_black_triangles; + const LensVertex *a = &mesh->vertices[triangle->vertex[0]]; + const LensVertex *b = &mesh->vertices[triangle->vertex[1]]; + const LensVertex *c = &mesh->vertices[triangle->vertex[2]]; + stats->approx_black_area_pixels2 += image_triangle_area(a, b, c); + stats->approx_black_max_edge_pixels = fmax(stats->approx_black_max_edge_pixels, + fmax(image_edge_length(a,b), fmax(image_edge_length(b,c), image_edge_length(c,a)))); + stats->approx_black_max_area_pixels2 = fmax(stats->approx_black_max_area_pixels2, + image_triangle_area(a,b,c)); + if (triangle->level >= config->max_level) ++stats->approx_black_level_stops; + } else ++stats->budget_incomplete_triangles; + } + if (unresolved_with_escape || all_unresolved) + ++stats->budget_incomplete_triangles; + } + stats->retry_requests = mesh->retry_requests; } static LensVertex midpoint_vertex(const LensVertex *a, const LensVertex *b) { @@ -727,7 +1032,8 @@ static int append_triangle(LensTriangle *triangles, size_t *count, unsigned int level, int evaluated) { if (*count >= capacity) return -1; - triangles[(*count)++] = (LensTriangle){{a, b, c}, level, evaluated}; + triangles[(*count)++] = + (LensTriangle){{a, b, c}, level, evaluated, 0}; return 0; } @@ -756,10 +1062,75 @@ static int append_triangle_with_parent_winding( evaluated); } +/* Append off-mesh failed/unresolved probes as persistent witnesses and + * rebuild the witness index. `consumed[i]` marks a probe sample that became + * a formal midpoint this generation. */ +static int promote_witnesses(FrameLensMesh *mesh, + const unsigned char *consumed) { + int added = 0; + for (size_t i = 0; i < mesh->sample_count; ++i) { + const FrameSample *s = &mesh->samples[i]; + if (s->cached || s->kind != FRAME_SAMPLE_PROBE) + continue; + if (consumed != NULL && consumed[i]) + continue; + if (s->vertex.outcome != RAY_OUTCOME_INCOMPLETE && + s->vertex.outcome != RAY_OUTCOME_UNRESOLVED) + continue; + if (ensure_vertices(mesh, mesh->vertex_count + 1)) + return -1; + LensVertex witness = s->vertex; + witness.diagnostic_probe = 1; + size_t a = s->edge_vertex[0], b = s->edge_vertex[1]; + if (a > b) { const size_t swap = a; a = b; b = swap; } + witness.probe_edge[0] = a; + witness.probe_edge[1] = b; + mesh->vertices[mesh->vertex_count++] = witness; + ++added; + } + /* Rebuild only when the witness set can have changed; the common no-witness + * refinement path stays O(T) without an O(V) scan. */ + if (added == 0 && mesh->witness_count == 0) + return 0; + return witness_rebuild(mesh); +} + int frame_lens_mesh_finish_generation(FrameLensMesh *mesh, const RefinementConfig *config) { if (mesh == NULL || config == NULL || mesh->sample_count == 0) return -1; + /* Serial decision invalidation: installation is an OpenMP bulk loop and + * must not write shared triangle flags from individual ray workers. Mark + * retried vertex ids once, then make a single pass over the mesh; witness + * edges are resolved through the O(1) witness index rather than scanning all + * retries for every triangle. */ + unsigned char *vertex_retried = + calloc(mesh->vertex_count ? mesh->vertex_count : 1, 1); + if (vertex_retried == NULL) + return -1; + for (size_t j = 0; j < mesh->sample_count; ++j) { + const FrameSample *s = &mesh->samples[j]; + if (s->kind == FRAME_SAMPLE_RETRY && s->vertex_id < mesh->vertex_count) + vertex_retried[s->vertex_id] = 1; + } + for (size_t i = 0; i < mesh->triangle_count; ++i) { + LensTriangle *t = &mesh->triangles[i]; + int touches = vertex_retried[t->vertex[0]] || + vertex_retried[t->vertex[1]] || + vertex_retried[t->vertex[2]]; + for (unsigned int side = 0; side < 3 && !touches; ++side) { + const size_t wid = witness_find(mesh, t->vertex[side], + t->vertex[(side + 1) % 3]); + if (wid != SIZE_MAX && wid < mesh->vertex_count && + vertex_retried[wid]) + touches = 1; + } + if (touches) { + t->evaluated = 0; + t->approx_black = 0; + } + } + free(vertex_retried); for (size_t i = 0; i < mesh->sample_count; ++i) if (!mesh->samples[i].vertex.traced && mesh->samples[i].kind == FRAME_SAMPLE_PROBE) @@ -769,12 +1140,31 @@ int frame_lens_mesh_finish_generation(FrameLensMesh *mesh, mesh->samples_include_probes = 0; return 0; } - for (size_t i = 0; i < mesh->triangle_count; ++i) - if (mesh->triangles[i].level < config->max_level) - mesh->triangles[i].evaluated = 1; + /* A triangle with a completed (or capped) probe can settle; one whose probe + * is still being retried must stay pending so the next generation sees the + * updated state. */ + for (size_t i = 0; i < mesh->triangle_count; ++i) { + LensTriangle *t = &mesh->triangles[i]; + const VertexMix mix = triangle_mix(mesh, t); + if (mix.u || mix.bad) + continue; + const unsigned side = longest_side(mesh, t); + const size_t probe = + find_probe(mesh, t->vertex[side], t->vertex[(side + 1) % 3]); + if (probe == SIZE_MAX) + continue; + const FrameSample *ps = &mesh->samples[probe]; + if (ps->vertex.outcome == RAY_OUTCOME_UNRESOLVED && !ps->cached) + continue; /* retry in flight */ + t->evaluated = 1; + } const size_t edge_count = mesh->triangle_count * 3; MeshEdge *edges = calloc(edge_count, sizeof *edges); unsigned char *requested = calloc(edge_count, sizeof *requested); + /* `wanted` records that a triangle asked to split before conformity may + * cancel its edges, so a fully blocked triangle can settle instead of + * re-requesting the same probes forever. */ + unsigned char *wanted = calloc(mesh->triangle_count, sizeof *wanted); unsigned char *allowed = calloc(mesh->triangle_count, sizeof *allowed); signed char *parity = calloc(mesh->triangle_count, sizeof *parity); double *jacobians = calloc(mesh->triangle_count, sizeof *jacobians); @@ -782,10 +1172,12 @@ int frame_lens_mesh_finish_generation(FrameLensMesh *mesh, * Keep that relation in the original triangle-side order so child emission * stays O(T), rather than scanning every sorted edge for every child side. */ size_t *side_midpoints = malloc(edge_count * sizeof *side_midpoints); - if (edges == NULL || requested == NULL || allowed == NULL || parity == NULL || - jacobians == NULL || side_midpoints == NULL) { - free(edges); free(requested); free(allowed); free(parity); free(jacobians); - free(side_midpoints); + unsigned char *consumed = calloc(mesh->sample_count ? mesh->sample_count : 1, 1); + if (edges == NULL || requested == NULL || wanted == NULL || allowed == NULL || + parity == NULL || jacobians == NULL || side_midpoints == NULL || + consumed == NULL) { + free(edges); free(requested); free(wanted); free(allowed); free(parity); + free(jacobians); free(side_midpoints); free(consumed); return -1; } for (size_t i = 0; i < edge_count; ++i) @@ -814,6 +1206,8 @@ int frame_lens_mesh_finish_generation(FrameLensMesh *mesh, probe_requires_split(mesh, triangle, side, config); } (void)discrete_jacobian(mesh, triangle, &jacobians[i], &parity[i]); + for (unsigned int side = 0; side < 3; ++side) + wanted[i] |= requested[3 * i + side] != 0; } } qsort(edges, edge_count, sizeof *edges, compare_mesh_edge); @@ -834,6 +1228,7 @@ int frame_lens_mesh_finish_generation(FrameLensMesh *mesh, if (allowed[left] && allowed[right]) { requested[3 * left + edges[first].side] = 1; requested[3 * right + edges[first + 1].side] = 1; + wanted[left] = wanted[right] = 1; } } } @@ -864,11 +1259,32 @@ int frame_lens_mesh_finish_generation(FrameLensMesh *mesh, size_t split_edges = 0; for (size_t i = 0; i < edge_count; ++i) split_edges += requested[3 * edges[i].triangle + edges[i].side] != 0; + /* A triangle whose requested split was fully cancelled by conformity or a + * geometric limit must settle here, or it re-requests the same probes every + * generation forever. A UUD/UDD blocked while geometry still allows is + * counted budget-incomplete by boundary_stats; at the stop scale it is an + * approximate-black boundary. */ + for (size_t i = 0; i < mesh->triangle_count; ++i) { + if (!wanted[i] || !allowed[i]) + continue; + int still = 0; + for (unsigned int side = 0; side < 3; ++side) + still |= requested[3 * i + side] != 0; + if (!still) { + mesh->triangles[i].evaluated = 1; + mesh->triangles[i].approx_black = 0; + } + } if (split_edges == 0) { + if (promote_witnesses(mesh, NULL)) { + free(edges); free(requested); free(wanted); free(allowed); free(parity); + free(jacobians); free(side_midpoints); free(consumed); + return -1; + } mesh->sample_count = 0; mesh->samples_include_probes = 0; - free(edges); free(requested); free(allowed); free(parity); free(jacobians); - free(side_midpoints); + free(edges); free(requested); free(wanted); free(allowed); free(parity); + free(jacobians); free(side_midpoints); free(consumed); return 0; } /* Allocate a single stable midpoint vertex for each requested edge group. */ @@ -885,8 +1301,8 @@ int frame_lens_mesh_finish_generation(FrameLensMesh *mesh, first = last; } if (ensure_vertices(mesh, mesh->vertex_count + midpoint_count)) { - free(edges); free(requested); free(allowed); free(parity); free(jacobians); - free(side_midpoints); + free(edges); free(requested); free(wanted); free(allowed); free(parity); + free(jacobians); free(side_midpoints); free(consumed); return -1; } size_t next_vertex = mesh->vertex_count; @@ -900,21 +1316,40 @@ int frame_lens_mesh_finish_generation(FrameLensMesh *mesh, any |= requested[3 * edges[i].triangle + edges[i].side] != 0; if (any) { const size_t probe_id = find_probe(mesh, edges[first].a, edges[first].b); - LensVertex midpoint = midpoint_vertex(&mesh->vertices[edges[first].a], - &mesh->vertices[edges[first].b]); - if (probe_id != SIZE_MAX) - midpoint = mesh->samples[probe_id].vertex; - mesh->vertices[next_vertex++] = midpoint; + size_t midpoint_id; + const size_t wid = witness_find(mesh, edges[first].a, edges[first].b); + if (wid != SIZE_MAX && wid < mesh->vertex_count && + mesh->vertices[wid].diagnostic_probe) { + /* Promote the existing off-mesh witness in place: one physical sample + * keeps a single stable vertex id. This must key on the persistent + * vertex identity, not on the generating sample being a cached PROBE: + * a witness retry is a FRAME_SAMPLE_RETRY and can trigger the split in + * the same generation. */ + midpoint_id = wid; + mesh->vertices[wid].diagnostic_probe = 0; + if (probe_id != SIZE_MAX) + consumed[probe_id] = 1; + } else { + LensVertex midpoint = midpoint_vertex(&mesh->vertices[edges[first].a], + &mesh->vertices[edges[first].b]); + if (probe_id != SIZE_MAX) { + midpoint = mesh->samples[probe_id].vertex; + consumed[probe_id] = 1; + } + midpoint.diagnostic_probe = 0; + midpoint_id = next_vertex; + mesh->vertices[next_vertex++] = midpoint; + } for (size_t i = first; i < last; ++i) - side_midpoints[3 * edges[i].triangle + edges[i].side] = next_vertex - 1; + side_midpoints[3 * edges[i].triangle + edges[i].side] = midpoint_id; } first = last; } const size_t old_count = mesh->triangle_count; LensTriangle *children = calloc(old_count * 4, sizeof *children); if (children == NULL) { - free(edges); free(requested); free(allowed); free(parity); free(jacobians); - free(side_midpoints); + free(edges); free(requested); free(wanted); free(allowed); free(parity); + free(jacobians); free(side_midpoints); free(consumed); return -1; } size_t child_count = 0; @@ -985,11 +1420,24 @@ int frame_lens_mesh_finish_generation(FrameLensMesh *mesh, mesh->triangle_count = child_count; mesh->triangle_capacity = old_count * 4; mesh->vertex_count = next_vertex; + const int promote_failed = promote_witnesses(mesh, consumed); mesh->sample_count = 0; mesh->samples_include_probes = 0; - free(edges); free(requested); free(allowed); free(parity); free(jacobians); - free(side_midpoints); - return (int)midpoint_count; + free(edges); free(requested); free(wanted); free(allowed); free(parity); + free(jacobians); free(side_midpoints); free(consumed); + return promote_failed ? -1 : (int)midpoint_count; +} + +void frame_retry_config_defaults(RefinementConfig *config, + const GeodesicTraceConfig *trace) { + if (config == NULL || trace == NULL) + return; + if (config->retry_step_increment == 0 && config->max_total_steps == 0) { + config->retry_step_increment = trace->max_steps; + config->max_total_steps = + trace->max_steps > (UINT_MAX / 4u) ? trace->max_steps + : trace->max_steps * 4u; + } } int frame_lens_mesh_refine(FrameLensMesh *mesh, @@ -1009,8 +1457,10 @@ int frame_lens_mesh_refine_with_progress( if (mesh == NULL || spacetime == NULL || observer == NULL || trace == NULL || config == NULL) return -1; + RefinementConfig effective = *config; + frame_retry_config_defaults(&effective, trace); for (size_t generation = 0;; ++generation) { - const int requested = frame_lens_mesh_prepare_generation(mesh, config); + const int requested = frame_lens_mesh_prepare_generation(mesh, &effective); if (requested < 0) return -1; if (requested == 0) return 0; if (callback != NULL) @@ -1019,12 +1469,31 @@ int frame_lens_mesh_refine_with_progress( #pragma omp parallel for schedule(static) for (size_t i = 0; i < mesh->sample_count; ++i) { const FrameSample *sample = &mesh->samples[i]; - const RayEndpoint endpoint = geodesic_trace_past( - spacetime, observer, sample->vertex.camera_direction, trace); + if (sample->cached) continue; + RayEndpoint endpoint; + if (sample->kind == FRAME_SAMPLE_RETRY) { + const LensVertex *v = &sample->vertex; + GeodesicRayState state = { + .coordinate_time = v->continuation_t, + .x = {v->continuation_x[0], v->continuation_x[1], + v->continuation_x[2]}, + .Pi = {v->continuation_Pi[0], v->continuation_Pi[1], + v->continuation_Pi[2]}, + .log_alpha_p0 = v->continuation_log_alpha_p0, + .log_alpha_p0_0 = v->continuation_log_alpha_p0_0, + .steps = v->continuation_steps}; + GeodesicTraceConfig retry_config = *trace; + retry_config.max_steps = sample->step_limit; + endpoint = geodesic_trace_past_from_state(spacetime, &state, + &retry_config); + } else { + endpoint = geodesic_trace_past(spacetime, observer, + sample->vertex.camera_direction, trace); + } /* Each request has a distinct destination vertex or probe slot. */ (void)frame_lens_mesh_install_sample(mesh, i, &endpoint); } - const int added = frame_lens_mesh_finish_generation(mesh, config); + const int added = frame_lens_mesh_finish_generation(mesh, &effective); if (added < 0) return -1; if (callback != NULL) callback(context, generation, 0, mesh->vertex_count, mesh->triangle_count, @@ -1088,7 +1557,7 @@ static int usable_triangle(const FrameLensMesh *mesh, const LensVertex *vertices[3]) { for (int i = 0; i < 3; ++i) { vertices[i] = &mesh->vertices[triangle->vertex[i]]; - if (vertices[i]->status != RAY_ENDPOINT_ESCAPED) + if (vertices[i]->outcome != RAY_OUTCOME_ESCAPED) return 0; } return spherical_area(vertices[0]->n_infinity, vertices[1]->n_infinity, @@ -1989,5 +2458,7 @@ void frame_lens_mesh_destroy(FrameLensMesh *mesh) { free(mesh->triangles); free(mesh->samples); free(mesh->probe_slots); + free(mesh->witness_slots); + free(mesh->witness_vertices); *mesh = (FrameLensMesh){0}; } diff --git a/src/frame.h b/src/frame.h index 2d0b21c..e883329 100644 --- a/src/frame.h +++ b/src/frame.h @@ -14,17 +14,34 @@ typedef struct { double camera_direction[3]; double n_infinity[3]; double log_frequency_ratio; - RayEndpointStatus status; + RayOutcome outcome; + RayReason reason; /* Asymptotic end this escaped vertex belongs to; a triangle must not * interpolate across two different ends. */ SpacetimeEndId end_id; int traced; + /* Retry continuation, valid when outcome == RAY_OUTCOME_UNRESOLVED: resume + * from this last accepted state instead of replaying the ray. */ + double continuation_t; + double continuation_x[3]; + double continuation_Pi[3]; + double continuation_log_alpha_p0; + double continuation_log_alpha_p0_0; + unsigned int continuation_steps; + unsigned int continuation_limit; + /* Persistent off-mesh probe witness; also participates in completion checks. + * When set, probe_edge holds the sorted edge (a,b) this witness samples. */ + int diagnostic_probe; + size_t probe_edge[2]; } LensVertex; typedef struct { size_t vertex[3]; unsigned int level; int evaluated; + /* Set when a boundary triangle containing UNRESOLVED vertices was blackened + * as a finite-resolution approximation rather than resolved. */ + int approx_black; } LensTriangle; /* Per-frame staged wall-clock breakdown for one movie frame. All fields are @@ -54,20 +71,51 @@ typedef struct { double jacobian_minimum; double min_edge_pixels; double min_area_pixels2; + /* Retry budget for UNRESOLVED vertices. retry_step_increment == 0 disables + * retry. max_total_steps is the per-ray hard cap on accepted steps; when a + * UUU / escape-containing triangle reaches it, the frame is reported as + * budget-incomplete instead of silently blackened. */ + unsigned int retry_step_increment; + unsigned int max_total_steps; } RefinementConfig; typedef enum { FRAME_SAMPLE_VERTEX, - FRAME_SAMPLE_PROBE + FRAME_SAMPLE_PROBE, + FRAME_SAMPLE_RETRY } FrameSampleKind; typedef struct { FrameSampleKind kind; size_t vertex_id; size_t edge_vertex[2]; + /* For FRAME_SAMPLE_RETRY: the new total accepted-step budget. */ + unsigned int step_limit; + int cached; LensVertex vertex; } FrameSample; +/* E/D/U triangle accounting for one finalized mesh. Counts use the + * rendering categories: E=ESCAPED, D=DARK, U=UNRESOLVED; triangles + * containing an INCOMPLETE vertex are counted as errors and are never + * blackened. */ +typedef struct { + size_t escaped_only; /* EEE */ + size_t dark_only; /* DDD */ + size_t eed_edd; /* EED / EDD boundary */ + size_t uud_udd; /* UUD / UDD */ + size_t u_with_escape; /* UEE / UED / UUE */ + size_t uuu; /* UUU */ + size_t error; /* any INCOMPLETE vertex */ + size_t approx_black_triangles; + double approx_black_area_pixels2; + double approx_black_max_edge_pixels; + double approx_black_max_area_pixels2; + size_t approx_black_level_stops; + size_t retry_requests; + size_t budget_incomplete_triangles; +} FrameBoundaryStats; + typedef struct { LensVertex *vertices; LensTriangle *triangles; @@ -80,6 +128,18 @@ typedef struct { int samples_include_probes; size_t *probe_slots; size_t probe_slot_capacity; + /* Cumulative count of retry rays requested across all generations. */ + size_t retry_requests; + /* Persistent off-mesh probe witnesses, keyed by their edge (sorted vertex + * ids). Slot value is witness_vertex_id + 1; 0 is empty. Witnesses are + * promoted in place to a formal midpoint when their edge is later split, so + * one physical sample always has one stable vertex id. */ + size_t *witness_slots; + size_t witness_slot_capacity; + /* Compact list of live witness vertex ids for rehashing and accounting. */ + size_t *witness_vertices; + size_t witness_count, witness_capacity; + size_t diagnostic_probe_count; } FrameLensMesh; typedef enum { @@ -143,6 +203,17 @@ int frame_lens_mesh_refine_with_progress( const ObserverState *observer, const GeodesicTraceConfig *trace, const RefinementConfig *config, FrameRefinementProgressCallback callback, void *context); +/* Recompute per-triangle approximate-black provenance and E/D/U accounting + * from the finalized mesh. Call after refinement has converged. Leaves the + * shared UNRESOLVED vertices untouched. */ +void frame_lens_mesh_boundary_stats(FrameLensMesh *mesh, + const RefinementConfig *config, + FrameBoundaryStats *stats); +/* Fill retry defaults (derived from the trace step budget) when the caller + * did not configure them explicitly. A zero max_total_steps disables retry + * only if retry_step_increment is also zero. */ +void frame_retry_config_defaults(RefinementConfig *config, + const GeodesicTraceConfig *trace); /* Whether frame_splat_catalog() should run the per-frame catalog prefetch or * rely on a movie-level union prefetch that already completed. */ diff --git a/src/geodesic.c b/src/geodesic.c index 4262693..0ec9327 100644 --- a/src/geodesic.c +++ b/src/geodesic.c @@ -30,14 +30,19 @@ static int invert(double a[3][3], double b[3][3]) { return 0; } -/* Equation (4) and (5) of Bohn et al., arXiv:1410.7775. */ -static int rhs(const MetricSlab *slab, double t, const State *s, - Derivative *out) { +/* Equation (4) and (5) of Bohn et al., arXiv:1410.7775. Returns a metric/data + * status so the integrator can report why a step failed instead of collapsing + * every failure into one generic error. */ +static SpacetimePointStatus rhs(const MetricSlab *slab, double t, + const State *s, Derivative *out) { MetricData m; double inv[3][3], up[3] = {0}, da_pi = 0, k_pi_pi = 0; - if (spacetime_slab_eval(slab, t, s->x, &m) || m.alpha <= 0 || - invert(m.gamma, inv)) - return -1; + const SpacetimePointStatus metric_status = + spacetime_slab_eval(slab, t, s->x, &m); + if (metric_status != SPACETIME_POINT_OK) + return metric_status; + if (m.alpha <= 0 || invert(m.gamma, inv)) + return SPACETIME_POINT_INVALID_METRIC; for (int i = 0; i < 3; i++) for (int j = 0; j < 3; j++) up[i] += inv[i][j] * s->Pi[j]; @@ -68,7 +73,7 @@ static int rhs(const MetricSlab *slab, double t, const State *s, db_pi - 0.5 * m.alpha * dg_pi_pi; } out->log_alpha_p0 = -da_pi + m.alpha * k_pi_pi; - return 0; + return SPACETIME_POINT_OK; } static State add(const State *s, const Derivative *d, double h) { @@ -81,20 +86,25 @@ static State add(const State *s, const Derivative *d, double h) { return r; } -static int rk4(const MetricSlab *slab, double t, double h, State *s) { +static SpacetimePointStatus rk4(const MetricSlab *slab, double t, double h, + State *s) { Derivative a, b, c, d; State q; - if (rhs(slab, t, s, &a)) - return -1; + SpacetimePointStatus status = rhs(slab, t, s, &a); + if (status != SPACETIME_POINT_OK) + return status; q = add(s, &a, h / 2); - if (rhs(slab, t + h / 2, &q, &b)) - return -1; + status = rhs(slab, t + h / 2, &q, &b); + if (status != SPACETIME_POINT_OK) + return status; q = add(s, &b, h / 2); - if (rhs(slab, t + h / 2, &q, &c)) - return -1; + status = rhs(slab, t + h / 2, &q, &c); + if (status != SPACETIME_POINT_OK) + return status; q = add(s, &c, h); - if (rhs(slab, t + h, &q, &d)) - return -1; + status = rhs(slab, t + h, &q, &d); + if (status != SPACETIME_POINT_OK) + return status; for (int i = 0; i < 3; i++) { s->x[i] += h * (a.x[i] + 2 * b.x[i] + 2 * c.x[i] + d.x[i]) / 6; s->Pi[i] += h * (a.Pi[i] + 2 * b.Pi[i] + 2 * c.Pi[i] + d.Pi[i]) / 6; @@ -103,7 +113,7 @@ static int rk4(const MetricSlab *slab, double t, double h, State *s) { (a.log_alpha_p0 + 2 * b.log_alpha_p0 + 2 * c.log_alpha_p0 + d.log_alpha_p0) / 6; - return 0; + return SPACETIME_POINT_OK; } int geodesic_initialize_past_ray_metric(const MetricData *metric, @@ -127,6 +137,7 @@ int geodesic_initialize_past_ray_metric(const MetricData *metric, s->Pi[i] /= m->alpha * k[0]; } s->log_alpha_p0 = log(m->alpha * k[0]); + s->log_alpha_p0_0 = s->log_alpha_p0; s->coordinate_time = o->coordinate_time; s->steps = 0; return isfinite(s->log_alpha_p0) ? 0 : -1; @@ -218,22 +229,104 @@ static AsymptoticStatus localize_worldtube_crossing(const MetricSlab *slab, return ASYMPTOTIC_OK; } -static GeodesicAdvanceResult legacy_escape_or_capture(const MetricSlab *slab, - const State *s, - SpacetimeRayStatus status, - RayEndpoint *out) { - out->status = - status == SPACETIME_RAY_ESCAPED ? RAY_ENDPOINT_ESCAPED - : RAY_ENDPOINT_CAPTURED; - if (out->status == RAY_ENDPOINT_ESCAPED) { - if (escaped_direction(slab, s->coordinate_time, s, out->n_infinity) == 0) - out->frequency_ratio = exp(-s->log_alpha_p0); - else - out->status = RAY_ENDPOINT_INTEGRATION_FAILURE; +static RayReason reason_from_point_status(SpacetimePointStatus status) { + switch (status) { + case SPACETIME_POINT_TIME_UNAVAILABLE: + return RAY_REASON_TIME_RANGE_EXHAUSTED; + case SPACETIME_POINT_OUT_OF_DOMAIN: + return RAY_REASON_OUT_OF_DOMAIN; + case SPACETIME_POINT_INVALID_METRIC: + return RAY_REASON_INVALID_METRIC; + case SPACETIME_POINT_INTERNAL_ERROR: + return RAY_REASON_PROTOCOL_ERROR; + default: + return RAY_REASON_INTEGRATION_ERROR; } - return out->status == RAY_ENDPOINT_INTEGRATION_FAILURE - ? GEODESIC_ADVANCE_FAILED - : GEODESIC_ADVANCE_TERMINATED; +} + +static void record_final_state(RayEndpoint *out, const State *s) { + if (s == NULL) { + out->stop_coordinate_time = NAN; + out->accepted_steps = 0; + for (int i = 0; i < 3; ++i) { + out->final_x[i] = NAN; + out->final_Pi[i] = NAN; + } + out->final_log_alpha_p0 = NAN; + out->final_log_alpha_p0_0 = NAN; + return; + } + out->stop_coordinate_time = s->coordinate_time; + out->accepted_steps = s->steps; + for (int i = 0; i < 3; ++i) { + out->final_x[i] = s->x[i]; + out->final_Pi[i] = s->Pi[i]; + } + out->final_log_alpha_p0 = s->log_alpha_p0; + out->final_log_alpha_p0_0 = s->log_alpha_p0_0; +} + +static void set_incomplete(RayEndpoint *out, RayReason reason, + const State *last) { + out->outcome = RAY_OUTCOME_INCOMPLETE; + out->reason = reason; + out->end_id = SPACETIME_END_NONE; + record_final_state(out, last); + out->threshold_value = NAN; +} + +/* Value of the monitored dark-threshold quantity at a trusted state. */ +static double monitored_threshold_value(const MetricSlab *slab, + ThresholdKind kind, const State *s) { + if (kind == THRESHOLD_LOG_ALPHA_P0) + return s->log_alpha_p0; + if (kind == THRESHOLD_LOG_ENERGY_GROWTH) + return s->log_alpha_p0 - s->log_alpha_p0_0; + if (kind == THRESHOLD_LOG_P0) { + MetricData m; + if (spacetime_slab_eval(slab, s->coordinate_time, s->x, &m) != + SPACETIME_POINT_OK || + m.alpha <= 0.0) + return NAN; + return s->log_alpha_p0 - log(m.alpha); + } + return NAN; +} + +/* Check the dark threshold on one trusted state. Returns nonzero and fills a + * DARK endpoint when the monitored quantity has reached the threshold. */ +static int threshold_reached(const MetricSlab *slab, + const GeodesicTraceConfig *config, const State *s, + RayEndpoint *out) { + if (config->threshold.kind == THRESHOLD_DISABLED) + return 0; + const double value = + monitored_threshold_value(slab, config->threshold.kind, s); + if (!isfinite(value) || value < config->threshold.value) + return 0; + out->outcome = RAY_OUTCOME_DARK; + out->reason = RAY_REASON_REDSHIFT_LIMIT; + out->end_id = SPACETIME_END_NONE; + record_final_state(out, s); + out->threshold_value = value; + return 1; +} + +/* Legacy no-declared-ends backends may still report a region escape. There + * is no position-based physical capture; a failed sky direction is reported + * as an integration failure, never as capture. */ +static GeodesicAdvanceResult legacy_escape(const MetricSlab *slab, + const State *s, RayEndpoint *out) { + if (escaped_direction(slab, s->coordinate_time, s, out->n_infinity) != 0) { + set_incomplete(out, RAY_REASON_INTEGRATION_ERROR, s); + return GEODESIC_ADVANCE_FAILED; + } + out->outcome = RAY_OUTCOME_ESCAPED; + out->reason = RAY_REASON_NONE; + out->frequency_ratio = exp(-s->log_alpha_p0); + record_final_state(out, s); + out->threshold_value = NAN; + return GEODESIC_ADVANCE_TERMINATED; } GeodesicAdvanceResult geodesic_advance_past_ray( @@ -241,53 +334,58 @@ GeodesicAdvanceResult geodesic_advance_past_ray( const GeodesicTraceConfig *config, RayEndpoint *out) { if (!slab || !s || !config || !out || config->coordinate_time_step <= 0 || !config->max_steps || !isfinite(slab_left_time) || - slab_left_time > s->coordinate_time) + slab_left_time > s->coordinate_time) { + if (out != NULL) { + out->outcome = RAY_OUTCOME_INCOMPLETE; + out->reason = RAY_REASON_PROTOCOL_ERROR; + record_final_state(out, s); + out->threshold_value = NAN; + } return GEODESIC_ADVANCE_FAILED; + } const AsymLifecycleMode mode = asym_lifecycle_mode(slab->source); if (mode == ASYM_LIFECYCLE_PROTOCOL_ERROR) { - out->status = RAY_ENDPOINT_INVALID; + set_incomplete(out, RAY_REASON_PROTOCOL_ERROR, s); return GEODESIC_ADVANCE_FAILED; } const int directed = mode == ASYM_LIFECYCLE_READY; const size_t end_count = directed ? spacetime_asymptotic_end_count(slab->source) : 0; while (s->coordinate_time > slab_left_time) { - if (config->capture_log_alpha_p0 > 0.0 && - s->log_alpha_p0 >= config->capture_log_alpha_p0) { - out->status = RAY_ENDPOINT_CAPTURED; + /* The threshold is checked only on trusted initial/accepted states. A + * trial stage that crosses it does not by itself produce DARK. */ + if (threshold_reached(slab, config, s, out)) return GEODESIC_ADVANCE_TERMINATED; - } - const SpacetimeRayStatus status = - spacetime_slab_classify(slab, s->coordinate_time, s->x); - if (status == SPACETIME_RAY_CAPTURED) { - out->status = RAY_ENDPOINT_CAPTURED; - return GEODESIC_ADVANCE_TERMINATED; - } - if (!directed && status != SPACETIME_RAY_ACTIVE) - return legacy_escape_or_capture(slab, s, status, out); + if (!directed && + spacetime_slab_classify(slab, s->coordinate_time, s->x) == + SPACETIME_RAY_ESCAPED) + return legacy_escape(slab, s, out); if (s->steps >= config->max_steps) { - out->status = RAY_ENDPOINT_MAX_STEPS; + /* Trustworthy trajectory, compute budget exhausted: retryable. */ + out->outcome = RAY_OUTCOME_UNRESOLVED; + out->reason = RAY_REASON_BUDGET_EXHAUSTED; + out->end_id = SPACETIME_END_NONE; + record_final_state(out, s); + out->threshold_value = NAN; return GEODESIC_ADVANCE_TERMINATED; } const double h = -fmin(config->coordinate_time_step, s->coordinate_time - slab_left_time); const State before = *s; - if (rk4(slab, s->coordinate_time, h, s)) + const SpacetimePointStatus step_status = rk4(slab, s->coordinate_time, h, s); + if (step_status != SPACETIME_POINT_OK) { + set_incomplete(out, reason_from_point_status(step_status), &before); return GEODESIC_ADVANCE_FAILED; + } s->coordinate_time += h; ++s->steps; if (!directed) continue; - if (spacetime_slab_classify(slab, s->coordinate_time, s->x) == - SPACETIME_RAY_CAPTURED) { - out->status = RAY_ENDPOINT_CAPTURED; - return GEODESIC_ADVANCE_TERMINATED; - } for (size_t i = 0; i < end_count; ++i) { SpacetimeAsymptoticEnd end; if (spacetime_asymptotic_end(slab->source, i, &end)) { - out->status = RAY_ENDPOINT_INVALID; + set_incomplete(out, RAY_REASON_PROTOCOL_ERROR, &before); return GEODESIC_ADVANCE_FAILED; } double f_before, f_after; @@ -300,14 +398,14 @@ GeodesicAdvanceResult geodesic_advance_past_ray( after_status == ASYMPTOTIC_TIME_RANGE_EXHAUSTED) { /* A first-class terminal reason, matching pre-route exhaustion: * preserve the end id and install it as terminated provenance. */ - out->status = RAY_ENDPOINT_TIME_RANGE_EXHAUSTED; + set_incomplete(out, RAY_REASON_TIME_RANGE_EXHAUSTED, &before); out->end_id = end.end_id; return GEODESIC_ADVANCE_TERMINATED; } if (before_status != ASYMPTOTIC_OK || after_status != ASYMPTOTIC_OK) { /* The backend cannot describe its own worldtube; this is an explicit * failure, not a physical escape. */ - out->status = RAY_ENDPOINT_INVALID; + set_incomplete(out, RAY_REASON_PROTOCOL_ERROR, &before); return GEODESIC_ADVANCE_FAILED; } /* Strict inside->outside: the step must end strictly outside, so a @@ -320,29 +418,37 @@ GeodesicAdvanceResult geodesic_advance_past_ray( const AsymptoticStatus localized = localize_worldtube_crossing( slab, end.end_id, &before, h, &crossing); if (localized == ASYMPTOTIC_TIME_RANGE_EXHAUSTED) { - out->status = RAY_ENDPOINT_TIME_RANGE_EXHAUSTED; + set_incomplete(out, RAY_REASON_TIME_RANGE_EXHAUSTED, &before); out->end_id = end.end_id; return GEODESIC_ADVANCE_TERMINATED; } if (localized != ASYMPTOTIC_OK) { - out->status = RAY_ENDPOINT_INVALID; + set_incomplete(out, RAY_REASON_PROTOCOL_ERROR, &before); return GEODESIC_ADVANCE_FAILED; } const AsymptoticStatus transfer = asymptotic_finish_escape( slab->source, end.end_id, crossing.coordinate_time, crossing.x, crossing.Pi, crossing.log_alpha_p0, out); - if (transfer == ASYMPTOTIC_OK) + if (transfer == ASYMPTOTIC_OK) { + record_final_state(out, &crossing); + out->threshold_value = NAN; return GEODESIC_ADVANCE_TERMINATED; - out->status = transfer == ASYMPTOTIC_TIME_RANGE_EXHAUSTED - ? RAY_ENDPOINT_TIME_RANGE_EXHAUSTED - : RAY_ENDPOINT_INVALID; + } if (transfer == ASYMPTOTIC_TIME_RANGE_EXHAUSTED) { + set_incomplete(out, RAY_REASON_TIME_RANGE_EXHAUSTED, &before); out->end_id = end.end_id; return GEODESIC_ADVANCE_TERMINATED; } + set_incomplete(out, RAY_REASON_PROTOCOL_ERROR, &before); return GEODESIC_ADVANCE_FAILED; } } + /* The loop can stop exactly at the slab's left boundary, so the final + * accepted step must also be checked before reporting ACTIVE; otherwise a + * ray that crossed the dark threshold on its last step would be settled as + * budget-unresolved and require another slab/retry. */ + if (threshold_reached(slab, config, s, out)) + return GEODESIC_ADVANCE_TERMINATED; return GEODESIC_ADVANCE_ACTIVE; } @@ -353,7 +459,12 @@ RayEndpoint geodesic_trace_past(const SpacetimeSource *source, RayEndpoint out = {.frequency_ratio = 0, .magnification = 1, .end_id = SPACETIME_END_NONE, - .status = RAY_ENDPOINT_INTEGRATION_FAILURE}; + .outcome = RAY_OUTCOME_INCOMPLETE, + .reason = RAY_REASON_INTEGRATION_ERROR, + .stop_coordinate_time = NAN, + .accepted_steps = 0, + .threshold_value = NAN}; + record_final_state(&out, NULL); if (!source || !observer || !config || config->coordinate_time_step <= 0 || !config->max_steps || fabs(dot(n, n) - 1) > 1e-10) return out; @@ -361,16 +472,19 @@ RayEndpoint geodesic_trace_past(const SpacetimeSource *source, const AsymptoticStatus route_status = asymptotic_route_camera(source, observer, n, &route); if (route_status == ASYMPTOTIC_UNSUPPORTED) { - out.status = RAY_ENDPOINT_INVALID; + out.outcome = RAY_OUTCOME_INCOMPLETE; + out.reason = RAY_REASON_UNSUPPORTED; return out; } if (route_status == ASYMPTOTIC_TIME_RANGE_EXHAUSTED) { - out.status = RAY_ENDPOINT_TIME_RANGE_EXHAUSTED; + out.outcome = RAY_OUTCOME_INCOMPLETE; + out.reason = RAY_REASON_TIME_RANGE_EXHAUSTED; out.end_id = route.end_id; return out; } if (route_status != ASYMPTOTIC_OK) { - out.status = RAY_ENDPOINT_INVALID; + out.outcome = RAY_OUTCOME_INCOMPLETE; + out.reason = RAY_REASON_PROTOCOL_ERROR; return out; } if (route.kind == ASYMPTOTIC_ROUTE_ESCAPED) { @@ -378,34 +492,79 @@ RayEndpoint geodesic_trace_past(const SpacetimeSource *source, out.n_infinity[i] = route.n_infinity[i]; out.frequency_ratio = route.frequency_ratio; out.end_id = route.end_id; - out.status = RAY_ENDPOINT_ESCAPED; + out.outcome = RAY_OUTCOME_ESCAPED; + out.reason = RAY_REASON_NONE; return out; } if (route.kind == ASYMPTOTIC_ROUTE_TIME_RANGE_EXHAUSTED) { - out.status = RAY_ENDPOINT_TIME_RANGE_EXHAUSTED; + out.outcome = RAY_OUTCOME_INCOMPLETE; + out.reason = RAY_REASON_TIME_RANGE_EXHAUSTED; out.end_id = route.end_id; return out; } if (route.kind != ASYMPTOTIC_ROUTE_INSIDE && route.kind != ASYMPTOTIC_ROUTE_ENTRY) { - out.status = RAY_ENDPOINT_INVALID; + out.outcome = RAY_OUTCOME_INCOMPLETE; + out.reason = RAY_REASON_PROTOCOL_ERROR; return out; } State state = {.coordinate_time = route.activate_t, .x = {route.x[0], route.x[1], route.x[2]}, .Pi = {route.Pi[0], route.Pi[1], route.Pi[2]}, .log_alpha_p0 = route.log_alpha_p0, + /* Camera-event reference, distinct from the entry-state L for + * an external camera. */ + .log_alpha_p0_0 = route.log_alpha_p0_camera, .steps = 0}; const double last_time = route.activate_t - config->coordinate_time_step * config->max_steps; MetricSlab *slab = NULL; if (spacetime_load_slab(source, route.activate_t, last_time - 1.0, &slab)) { - out.status = RAY_ENDPOINT_INVALID; + out.outcome = RAY_OUTCOME_INCOMPLETE; + out.reason = RAY_REASON_IO_ERROR; return out; } if (geodesic_advance_past_ray(slab, &state, last_time, config, &out) == - GEODESIC_ADVANCE_ACTIVE) - out.status = RAY_ENDPOINT_MAX_STEPS; + GEODESIC_ADVANCE_ACTIVE) { + out.outcome = RAY_OUTCOME_UNRESOLVED; + out.reason = RAY_REASON_BUDGET_EXHAUSTED; + record_final_state(&out, &state); + } + spacetime_free_slab(slab); + return out; +} + +RayEndpoint geodesic_trace_past_from_state(const SpacetimeSource *source, + const GeodesicRayState *state_in, + const GeodesicTraceConfig *config) { + RayEndpoint out = {.frequency_ratio = 0, + .magnification = 1, + .end_id = SPACETIME_END_NONE, + .outcome = RAY_OUTCOME_INCOMPLETE, + .reason = RAY_REASON_INTEGRATION_ERROR, + .stop_coordinate_time = NAN, + .accepted_steps = 0, + .threshold_value = NAN}; + record_final_state(&out, NULL); + if (!source || !state_in || !config || config->coordinate_time_step <= 0 || + !config->max_steps || state_in->steps >= config->max_steps) + return out; + State state = *state_in; + const double remaining = + config->coordinate_time_step * (double)(config->max_steps - state.steps); + const double last_time = state.coordinate_time - remaining - 1.0; + MetricSlab *slab = NULL; + if (spacetime_load_slab(source, state.coordinate_time, last_time, &slab)) { + out.outcome = RAY_OUTCOME_INCOMPLETE; + out.reason = RAY_REASON_IO_ERROR; + return out; + } + if (geodesic_advance_past_ray(slab, &state, last_time, config, &out) == + GEODESIC_ADVANCE_ACTIVE) { + out.outcome = RAY_OUTCOME_UNRESOLVED; + out.reason = RAY_REASON_BUDGET_EXHAUSTED; + record_final_state(&out, &state); + } spacetime_free_slab(slab); return out; } diff --git a/src/geodesic.h b/src/geodesic.h index 7d8bfb9..20b97fb 100644 --- a/src/geodesic.h +++ b/src/geodesic.h @@ -4,33 +4,75 @@ #include "observer.h" #include "spacetime.h" +/* Rendering/completion category. This is deliberately separate from the + * diagnostic reason below, and from the ray-pool lifecycle. */ typedef enum { - RAY_ENDPOINT_ESCAPED, - RAY_ENDPOINT_CAPTURED, - RAY_ENDPOINT_MAX_STEPS, - RAY_ENDPOINT_INTEGRATION_FAILURE, - RAY_ENDPOINT_TIME_RANGE_EXHAUSTED, - RAY_ENDPOINT_INVALID -} RayEndpointStatus; + RAY_OUTCOME_ESCAPED = 0, /* reached an infinity end; carries a payload */ + RAY_OUTCOME_DARK, /* normal dark terminal (currently redshift limit) */ + RAY_OUTCOME_UNRESOLVED, /* trustworthy trajectory, compute budget exhausted */ + RAY_OUTCOME_INCOMPLETE /* history/domain/metric/integration/protocol error */ +} RayOutcome; + +/* Diagnostic reason. Different DARK reasons must not create a mesh seam; the + * reason is for accounting and provenance only. */ +typedef enum { + RAY_REASON_NONE = 0, + RAY_REASON_REDSHIFT_LIMIT, + RAY_REASON_BUDGET_EXHAUSTED, + RAY_REASON_TIME_RANGE_EXHAUSTED, + RAY_REASON_OUT_OF_DOMAIN, + RAY_REASON_INVALID_METRIC, + RAY_REASON_INTEGRATION_ERROR, + RAY_REASON_UNSUPPORTED, + RAY_REASON_PROTOCOL_ERROR, + RAY_REASON_IO_ERROR +} RayReason; + +/* Monitored quantity used by the dark-redshift termination policy. `LOG_P0` + * is ln(p^0) = L - ln(alpha); it differs from `LOG_ALPHA_P0` by a local + * function of position and must not be confused with the true infinity + * frequency ratio g. */ +typedef enum { + THRESHOLD_DISABLED = 0, + THRESHOLD_LOG_ALPHA_P0, /* absolute L = ln(alpha p^0) */ + THRESHOLD_LOG_P0, /* ln(p^0) = L - ln(alpha) */ + THRESHOLD_LOG_ENERGY_GROWTH /* L - L0, local energy growth since the start */ +} ThresholdKind; + +typedef struct { + ThresholdKind kind; + double value; /* terminate when the monitored quantity reaches this */ + unsigned int policy_version; +} ThresholdPolicy; typedef struct { double n_infinity[3]; double frequency_ratio; /* E_camera / E_infinity */ double magnification; /* Filled by the future local inverse lens map. */ - /* Meaningful for RAY_ENDPOINT_ESCAPED and for - * RAY_ENDPOINT_TIME_RANGE_EXHAUSTED; SPACETIME_END_NONE otherwise. */ + /* End this escape belongs to; SPACETIME_END_NONE when not applicable. */ SpacetimeEndId end_id; - RayEndpointStatus status; + RayOutcome outcome; + RayReason reason; + /* Last trusted state at termination. For an ESCAPED endpoint this is the + * (finite) numerical truncation position, not the true parameter end. */ + double stop_coordinate_time; + unsigned int accepted_steps; + /* Last trusted continuous state, used to resume an UNRESOLVED ray from its + * last accepted step instead of replaying it from the camera. */ + double final_x[3]; + double final_Pi[3]; + double final_log_alpha_p0; + double final_log_alpha_p0_0; /* original reference L0 for retries */ + /* Monitored threshold value at termination, NAN when not applicable. */ + double threshold_value; } RayEndpoint; typedef struct { double coordinate_time_step; unsigned int max_steps; - /* A positive value terminates a backwards ray whose horizon redshift has - * made log(alpha p^0) reach this value. Zero disables this analytic/demo - * criterion; numerical moving-puncture backends use their AH-calibrated - * spatial cutoff instead. */ - double capture_log_alpha_p0; + /* Normal dark terminal for every backend. No backend may substitute a + * position/horizon cutoff for physical capture. */ + ThresholdPolicy threshold; } GeodesicTraceConfig; typedef struct { @@ -38,6 +80,9 @@ typedef struct { double x[3]; double Pi[3]; double log_alpha_p0; + /* Reference L at the start of this ray's integration, carried unchanged + * through retries so THRESHOLD_LOG_ENERGY_GROWTH stays camera-relative. */ + double log_alpha_p0_0; unsigned int steps; } GeodesicRayState; @@ -67,4 +112,11 @@ GeodesicAdvanceResult geodesic_advance_past_ray( const MetricSlab *slab, GeodesicRayState *state, double slab_left_time, const GeodesicTraceConfig *config, RayEndpoint *endpoint); +/* Resume a past ray from its last trusted state and integrate to the total + * step budget in `config->max_steps` (state->steps counts steps already + * consumed). Used to retry UNRESOLVED rays without replaying them from the + * camera. */ +RayEndpoint geodesic_trace_past_from_state(const SpacetimeSource *source, + const GeodesicRayState *state, + const GeodesicTraceConfig *config); #endif diff --git a/src/lens_map.c b/src/lens_map.c index fea9b73..f407c50 100644 --- a/src/lens_map.c +++ b/src/lens_map.c @@ -10,7 +10,10 @@ /* All scalar fields are explicitly little-endian; never serialize C structs * because their padding and size_t width are ABI-dependent. */ static const unsigned char lens_map_magic[8] = {'G', 'R', 'L', 'E', 'N', 'S', 1, 0}; -enum { LENS_MAP_VERSION = 1, LENS_MAP_ENDIAN = 0x01020304u }; +/* Version 2 stores the two-level RayOutcome instead of the removed + * RayEndpointStatus. Version 1 files are rejected: their old captured bit + * cannot be upgraded into the new dark/unresolved/error provenance. */ +enum { LENS_MAP_VERSION = 2, LENS_MAP_ENDIAN = 0x01020304u }; static uint32_t crc32_update(uint32_t crc, const void *data, size_t size) { const unsigned char *bytes = data; @@ -64,11 +67,14 @@ static int valid_mesh(const FrameLensMesh *m) { if (m == NULL || m->vertex_count == 0 || m->triangle_count == 0) return 0; for (size_t i = 0; i < m->vertex_count; ++i) { const LensVertex *v = &m->vertices[i]; - if (!v->traced || v->status < RAY_ENDPOINT_ESCAPED || - v->status > RAY_ENDPOINT_INTEGRATION_FAILURE || !isfinite(v->image_x) || + if (!v->traced || v->outcome > RAY_OUTCOME_INCOMPLETE || + v->reason > RAY_REASON_IO_ERROR || !isfinite(v->image_x) || !isfinite(v->image_y) || !isfinite(v->log_frequency_ratio) || - !unit_vector(v->camera_direction) || - (v->status == RAY_ENDPOINT_ESCAPED && !unit_vector(v->n_infinity))) return 0; + !unit_vector(v->camera_direction)) + return 0; + if (v->outcome == RAY_OUTCOME_ESCAPED && + (!unit_vector(v->n_infinity) || v->end_id == SPACETIME_END_NONE)) + return 0; } for (size_t i = 0; i < m->triangle_count; ++i) for (int j = 0; j < 3; ++j) @@ -77,34 +83,51 @@ static int valid_mesh(const FrameLensMesh *m) { } int lens_map_write(const char *path, int width, int height, double fov, + const LensMapProvenance *provenance, const LensMapFrame *frames, size_t frame_count) { - if (path == NULL || frames == NULL || width <= 0 || height <= 0 || - !isfinite(fov) || fov <= 0.0 || fov >= 179.0 || frame_count == 0 || - frame_count > UINT64_MAX) return -1; + if (path == NULL || provenance == NULL || frames == NULL || width <= 0 || + height <= 0 || !isfinite(fov) || fov <= 0.0 || fov >= 179.0 || + frame_count == 0 || frame_count > UINT64_MAX) + return -1; for (size_t i = 0; i < frame_count; ++i) if (!valid_mesh(&frames[i].mesh)) return -1; FILE *file = fopen(path, "wb"); if (file == NULL) return -1; int failed = write_bytes(file, lens_map_magic, sizeof lens_map_magic, NULL) || write_u32(file, LENS_MAP_VERSION, NULL) || write_u32(file, LENS_MAP_ENDIAN, NULL) || write_u32(file, (uint32_t)width, NULL) || write_u32(file, (uint32_t)height, NULL) || - write_double(file, fov, NULL) || write_u64(file, (uint64_t)frame_count, NULL); + write_double(file, fov, NULL) || write_u64(file, (uint64_t)frame_count, NULL) || + write_u32(file, provenance->threshold_kind, NULL) || + write_u32(file, provenance->threshold_policy_version, NULL) || + write_double(file, provenance->threshold_value, NULL) || + write_u32(file, provenance->retry_step_increment, NULL) || + write_u32(file, provenance->max_total_steps, NULL) || + write_u32(file, provenance->max_level, NULL) || + write_u32(file, provenance->integrator, NULL) || + write_double(file, provenance->min_edge_pixels, NULL) || + write_double(file, provenance->min_area_pixels2, NULL) || + write_double(file, provenance->coordinate_time_step, NULL) || + write_u32(file, provenance->initial_max_steps, NULL); for (size_t f = 0; !failed && f < frame_count; ++f) { const FrameLensMesh *m = &frames[f].mesh; uint32_t crc = UINT32_MAX; failed = write_u64(file, frames[f].frame_id, NULL) || write_double(file, frames[f].coordinate_time, NULL) || write_double(file, frames[f].proper_time, NULL) || write_u64(file, (uint64_t)m->vertex_count, NULL) || - write_u64(file, (uint64_t)m->triangle_count, NULL); + write_u64(file, (uint64_t)m->triangle_count, NULL) || + write_u64(file, (uint64_t)m->retry_requests, NULL); for (size_t i = 0; !failed && i < m->vertex_count; ++i) { const LensVertex *v = &m->vertices[i]; failed = write_double(file, v->image_x, &crc) || write_double(file, v->image_y, &crc); for (int j = 0; !failed && j < 3; ++j) failed = write_double(file, v->camera_direction[j], &crc); for (int j = 0; !failed && j < 3; ++j) failed = write_double(file, v->n_infinity[j], &crc); failed = failed || write_double(file, v->log_frequency_ratio, &crc) || - write_u32(file, (uint32_t)v->status, &crc); + write_u32(file, (uint32_t)v->end_id, &crc) || + write_u32(file, (uint32_t)v->outcome, &crc) || + write_u32(file, (uint32_t)v->reason, &crc); } for (size_t i = 0; !failed && i < m->triangle_count; ++i) { for (int j = 0; j < 3; ++j) failed = failed || write_u64(file, m->triangles[i].vertex[j], &crc); - failed = failed || write_u32(file, m->triangles[i].level, &crc); + failed = failed || write_u32(file, m->triangles[i].level, &crc) || + write_u32(file, (uint32_t)m->triangles[i].approx_black, &crc); } failed = failed || write_u32(file, crc ^ UINT32_MAX, NULL); } @@ -118,41 +141,69 @@ void lens_map_destroy(LensMap *map) { free(map->frames); *map = (LensMap){0}; } -int lens_map_read(const char *path, LensMap *map) { +int lens_map_read(const char *path, LensMapProvenance *provenance, + LensMap *map) { if (path == NULL || map == NULL) return -1; *map = (LensMap){0}; FILE *file = fopen(path, "rb"); if (file == NULL) return -1; unsigned char magic[8]; uint32_t version, endian, width, height; uint64_t count; + LensMapProvenance prov = {0}; int failed = read_bytes(file, magic, sizeof magic, NULL) || memcmp(magic, lens_map_magic, sizeof magic) || read_u32(file, &version, NULL) || read_u32(file, &endian, NULL) || read_u32(file, &width, NULL) || read_u32(file, &height, NULL) || read_double(file, &map->horizontal_fov_deg, NULL) || read_u64(file, &count, NULL) || + read_u32(file, &prov.threshold_kind, NULL) || + read_u32(file, &prov.threshold_policy_version, NULL) || + read_double(file, &prov.threshold_value, NULL) || + read_u32(file, &prov.retry_step_increment, NULL) || + read_u32(file, &prov.max_total_steps, NULL) || + read_u32(file, &prov.max_level, NULL) || + read_u32(file, &prov.integrator, NULL) || + read_double(file, &prov.min_edge_pixels, NULL) || + read_double(file, &prov.min_area_pixels2, NULL) || + read_double(file, &prov.coordinate_time_step, NULL) || + read_u32(file, &prov.initial_max_steps, NULL) || version != LENS_MAP_VERSION || endian != LENS_MAP_ENDIAN || width == 0 || height == 0 || width > INT32_MAX || height > INT32_MAX || !isfinite(map->horizontal_fov_deg) || map->horizontal_fov_deg <= 0.0 || map->horizontal_fov_deg >= 179.0 || count == 0 || + prov.threshold_kind > THRESHOLD_LOG_ENERGY_GROWTH || + !isfinite(prov.threshold_value) || + prov.coordinate_time_step < 0.0 || !isfinite(prov.coordinate_time_step) || count > SIZE_MAX / sizeof *map->frames; if (failed) goto done; + map->provenance = prov; + if (provenance != NULL) + *provenance = prov; map->width = (int)width; map->height = (int)height; map->frame_count = (size_t)count; map->frames = calloc(map->frame_count, sizeof *map->frames); if (map->frames == NULL) { failed = 1; goto done; } for (size_t f = 0; !failed && f < map->frame_count; ++f) { - LensMapFrame *frame = &map->frames[f]; uint64_t vertices, triangles; uint32_t stored_crc, crc = UINT32_MAX; + LensMapFrame *frame = &map->frames[f]; uint64_t vertices, triangles, retry_requests; + uint32_t stored_crc, crc = UINT32_MAX; failed = read_u64(file, &frame->frame_id, NULL) || read_double(file, &frame->coordinate_time, NULL) || read_double(file, &frame->proper_time, NULL) || read_u64(file, &vertices, NULL) || read_u64(file, &triangles, NULL) || + read_u64(file, &retry_requests, NULL) || !isfinite(frame->coordinate_time) || !isfinite(frame->proper_time) || vertices == 0 || triangles == 0 || - vertices > SIZE_MAX / sizeof *frame->mesh.vertices || triangles > SIZE_MAX / sizeof *frame->mesh.triangles; + vertices > SIZE_MAX / sizeof *frame->mesh.vertices || triangles > SIZE_MAX / sizeof *frame->mesh.triangles || + retry_requests > SIZE_MAX; if (failed) break; frame->mesh.vertices = calloc((size_t)vertices, sizeof *frame->mesh.vertices); frame->mesh.triangles = calloc((size_t)triangles, sizeof *frame->mesh.triangles); if (frame->mesh.vertices == NULL || frame->mesh.triangles == NULL) { failed = 1; break; } frame->mesh.vertex_count = frame->mesh.vertex_capacity = (size_t)vertices; frame->mesh.triangle_count = frame->mesh.triangle_capacity = (size_t)triangles; + frame->mesh.retry_requests = (size_t)retry_requests; for (size_t i = 0; !failed && i < frame->mesh.vertex_count; ++i) { - LensVertex *v = &frame->mesh.vertices[i]; uint32_t status; + LensVertex *v = &frame->mesh.vertices[i]; uint32_t end_id, outcome, reason; failed = read_double(file, &v->image_x, &crc) || read_double(file, &v->image_y, &crc); for (int j = 0; !failed && j < 3; ++j) failed = read_double(file, &v->camera_direction[j], &crc); for (int j = 0; !failed && j < 3; ++j) failed = read_double(file, &v->n_infinity[j], &crc); - failed = failed || read_double(file, &v->log_frequency_ratio, &crc) || read_u32(file, &status, &crc) || - status > RAY_ENDPOINT_INTEGRATION_FAILURE; - v->status = (RayEndpointStatus)status; v->traced = 1; + failed = failed || read_double(file, &v->log_frequency_ratio, &crc) || + read_u32(file, &end_id, &crc) || read_u32(file, &outcome, &crc) || + read_u32(file, &reason, &crc) || outcome > RAY_OUTCOME_INCOMPLETE || + reason > RAY_REASON_IO_ERROR; + v->end_id = (SpacetimeEndId)end_id; + v->outcome = (RayOutcome)outcome; + v->reason = (RayReason)reason; + v->traced = 1; } for (size_t i = 0; !failed && i < frame->mesh.triangle_count; ++i) { for (int j = 0; j < 3; ++j) { @@ -163,7 +214,12 @@ int lens_map_read(const char *path, LensMap *map) { } frame->mesh.triangles[i].vertex[j] = (size_t)index; } - failed = failed || read_u32(file, &frame->mesh.triangles[i].level, &crc); frame->mesh.triangles[i].evaluated = 1; + uint32_t approx_black = 0; + failed = failed || + read_u32(file, &frame->mesh.triangles[i].level, &crc) || + read_u32(file, &approx_black, &crc); + frame->mesh.triangles[i].approx_black = approx_black != 0; + frame->mesh.triangles[i].evaluated = 1; } failed = failed || read_u32(file, &stored_crc, NULL) || stored_crc != (crc ^ UINT32_MAX) || !valid_mesh(&frame->mesh); } diff --git a/src/lens_map.h b/src/lens_map.h index e269922..5f5dbd6 100644 --- a/src/lens_map.h +++ b/src/lens_map.h @@ -16,17 +16,37 @@ typedef struct { FrameLensMesh mesh; } LensMapFrame; +/* File-level provenance. Stored explicitly so a replay can be attributed to + * the terminal policy and integration settings that produced it. */ +typedef struct { + uint32_t threshold_kind; /* ThresholdKind */ + uint32_t threshold_policy_version; + double threshold_value; + uint32_t retry_step_increment; + uint32_t max_total_steps; + uint32_t max_level; + uint32_t integrator; /* 0 = fixed-step RK4 (transitional) */ + double min_edge_pixels; + double min_area_pixels2; + /* Integration source: the coordinate-time step and the initial per-ray + * accepted-step budget used for the first trace. */ + double coordinate_time_step; + uint32_t initial_max_steps; +} LensMapProvenance; + typedef struct { int width, height; double horizontal_fov_deg; + LensMapProvenance provenance; LensMapFrame *frames; size_t frame_count; } LensMap; int lens_map_write(const char *path, int width, int height, - double horizontal_fov_deg, const LensMapFrame *frames, - size_t frame_count); -int lens_map_read(const char *path, LensMap *map); + double horizontal_fov_deg, + const LensMapProvenance *provenance, + const LensMapFrame *frames, size_t frame_count); +int lens_map_read(const char *path, LensMapProvenance *provenance, LensMap *map); void lens_map_destroy(LensMap *map); #endif diff --git a/src/main.c b/src/main.c index 4cdae04..956fc88 100644 --- a/src/main.c +++ b/src/main.c @@ -23,9 +23,11 @@ typedef struct { int width, height; int coarse_cell_pixels; int draw_mesh; + int allow_incomplete; int psf_direct; int verbose; double horizontal_fov_deg, look_ra_deg, look_dec_deg, exposure; + double dark_threshold; double max_magnification; double max_cache_psf_flux; double psf_relative_tail; @@ -98,7 +100,10 @@ static int parse_double(const char *text, double *value) { char *end; errno = 0; *value = strtod(text, &end); - return errno || *end || *value <= 0.0 || *value >= 179.0 ? -1 : 0; + return errno || *end || !isfinite(*value) || *value <= 0.0 || + *value >= 179.0 + ? -1 + : 0; } static int parse_ra_deg(const char *text, double *value) { @@ -119,7 +124,7 @@ static int parse_positive(const char *text, double *value) { char *end; errno = 0; *value = strtod(text, &end); - return errno || *end || *value <= 0.0 ? -1 : 0; + return errno || *end || !isfinite(*value) || *value <= 0.0 ? -1 : 0; } static int parse_nonnegative(const char *text, double *value) { @@ -379,6 +384,7 @@ static int parse_args(int argc, char **argv, Settings *s, .look_ra_deg = 90.0, .look_dec_deg = -90.0, .exposure = 1e-3, + .dark_threshold = 8.0, .max_magnification = INFINITY, .max_cache_psf_flux = 1.0, .psf_relative_tail = 1e-8, @@ -453,6 +459,8 @@ static int parse_args(int argc, char **argv, Settings *s, !parse_positive(argv[++i], &s->refinement.min_area_pixels2)) { } else if (!strcmp(argv[i], "--draw-mesh")) { s->draw_mesh = 1; + } else if (!strcmp(argv[i], "--allow-incomplete")) { + s->allow_incomplete = 1; } else if (!strcmp(argv[i], "--fov-deg") && i + 1 < argc && !parse_double(argv[++i], &s->horizontal_fov_deg)) { s->fov_specified = 1; @@ -464,6 +472,8 @@ static int parse_args(int argc, char **argv, Settings *s, s->look_specified = 1; } else if (!strcmp(argv[i], "--exposure") && i + 1 < argc && !parse_positive(argv[++i], &s->exposure)) { + } else if (!strcmp(argv[i], "--dark-threshold") && i + 1 < argc && + !parse_positive(argv[++i], &s->dark_threshold)) { } else if (!strcmp(argv[i], "--tone-map") && i + 1 < argc) { const char *mode = argv[++i]; if (!strcmp(mode, "softclip")) @@ -636,6 +646,7 @@ static void print_help(const char *program) { " --look-ra-deg D ICRS look direction right ascension in degrees (default: 90)\n" " --look-dec-deg D ICRS look direction declination in degrees (default: -90)\n" " --exposure E Linear exposure multiplier (default: 1e-3)\n" + " --dark-threshold T Camera-relative dark cutoff L-L0 (default: 8)\n" " --tone-map MODE Display transform: softclip or reinhard\n" " (default: softclip)\n" " --tone-map-p P Softclip hardness P >= 1 (default: 2)\n" @@ -706,6 +717,8 @@ static void print_help(const char *program) { fputs(" --draw-mesh Also write the final lens-mesh overlay as _mesh.ppm\n", stdout); #endif + fputs(" --allow-incomplete Publish even when unresolved/error rays remain (diagnostic; output is marked incomplete)\n", + stdout); #ifdef SPACETIME_ALCUBIERRE fputs( "\nAlcubierre warp bubble (moving x_s(t)=v_s*t; no capture):\n" @@ -806,13 +819,14 @@ static void report_splat_worker_progress(void *context, size_t worker_id, static void ray_pool_status_counts(const RayPool *rays, size_t *pending, size_t *active, size_t *terminated, - size_t *failed) { - *pending = *active = *terminated = *failed = 0; + size_t *unresolved, size_t *failed) { + *pending = *active = *terminated = *unresolved = *failed = 0; for (size_t i = 0; i < rays->count; ++i) switch (rays->status[i]) { case RAY_POOL_PENDING: ++*pending; break; case RAY_POOL_ACTIVE: ++*active; break; case RAY_POOL_TERMINATED: ++*terminated; break; + case RAY_POOL_UNRESOLVED: ++*unresolved; break; case RAY_POOL_FAILED: ++*failed; break; } } @@ -839,7 +853,7 @@ static void report_frame_refinement(void *context, size_t generation, #ifdef SPACETIME_ALCUBIERRE /* Upper bound on the per-ray step budget. Legal parameters whose worst-case * near-comoving ray could need more than this are rejected at startup rather - * than silently terminating as RAY_ENDPOINT_MAX_STEPS. */ + * than silently terminating as UNRESOLVED/BUDGET_EXHAUSTED. */ #define ALCUBIERRE_MAX_TRACE_STEPS (1u << 24) /* Safety margin over the straight-line worst case: wall-region deflection can @@ -867,16 +881,109 @@ static double alcubierre_step_budget(const Settings *s) { } #endif +/* Recompute approximate-black provenance and report E/D/U accounting. The + * result is diagnostic for now; production failure gating is layered on top + * of the same counters. */ +static void report_boundary_stats(const Settings *s, + FrameLensMesh *const *meshes, + size_t frame_count, + const RefinementConfig *config, + FrameBoundaryStats *out_total) { + FrameBoundaryStats total = {0}; + for (size_t i = 0; i < frame_count; ++i) { + FrameBoundaryStats frame_stats; + frame_lens_mesh_boundary_stats(meshes[i], config, &frame_stats); + total.escaped_only += frame_stats.escaped_only; + total.dark_only += frame_stats.dark_only; + total.eed_edd += frame_stats.eed_edd; + total.uud_udd += frame_stats.uud_udd; + total.u_with_escape += frame_stats.u_with_escape; + total.uuu += frame_stats.uuu; + total.error += frame_stats.error; + total.approx_black_triangles += frame_stats.approx_black_triangles; + total.approx_black_area_pixels2 += frame_stats.approx_black_area_pixels2; + total.approx_black_max_edge_pixels = fmax(total.approx_black_max_edge_pixels, + frame_stats.approx_black_max_edge_pixels); + total.approx_black_max_area_pixels2 = fmax(total.approx_black_max_area_pixels2, + frame_stats.approx_black_max_area_pixels2); + total.approx_black_level_stops += frame_stats.approx_black_level_stops; + total.retry_requests += frame_stats.retry_requests; + total.budget_incomplete_triangles += frame_stats.budget_incomplete_triangles; + } + if (out_total != NULL) + *out_total = total; + if (s->verbose && frame_count > 0) + fprintf(stderr, + "Boundary accounting: EEE=%zu DDD=%zu EED/EDD=%zu UUD/UDD=%zu " + "U+E=%zu UUU=%zu error=%zu; approx-black=%zu (%.6g px^2), " + "retries=%zu, budget-incomplete=%zu.\n", + total.escaped_only, total.dark_only, total.eed_edd, total.uud_udd, + total.u_with_escape, total.uuu, total.error, + total.approx_black_triangles, total.approx_black_area_pixels2, + total.retry_requests, total.budget_incomplete_triangles); + if (s->verbose && total.approx_black_triangles) + fprintf(stderr, "Approx-black achieved scale: max edge=%.6g px, max area=%.6g px^2, level stops=%zu.\n", + total.approx_black_max_edge_pixels, total.approx_black_max_area_pixels2, + total.approx_black_level_stops); +} + +static LensMapProvenance lens_map_provenance(const Settings *s, + const GeodesicTraceConfig *trace) { + RefinementConfig rc = s->refinement; + frame_retry_config_defaults(&rc, trace); + return (LensMapProvenance){.threshold_kind = (uint32_t)trace->threshold.kind, + .threshold_policy_version = + trace->threshold.policy_version, + .threshold_value = trace->threshold.value, + .retry_step_increment = rc.retry_step_increment, + .max_total_steps = rc.max_total_steps, + .max_level = rc.max_level, + .integrator = 0, + .min_edge_pixels = rc.min_edge_pixels, + .min_area_pixels2 = rc.min_area_pixels2, + .coordinate_time_step = trace->coordinate_time_step, + .initial_max_steps = trace->max_steps}; +} + +/* Refuse to publish silently on true errors or on unresolved triangles that + * exhausted the configured total budget. Approximate-black UUD/UDD boundary + * triangles are an accepted finite-resolution error and do not block output. + * A diagnostic run may override this, but the incompleteness is reported. */ +static int boundary_allows_publish(const Settings *s, + const FrameBoundaryStats *total) { + const size_t blocking = total->error + total->budget_incomplete_triangles; + if (blocking == 0) + return 1; + if (s->allow_incomplete) { + fprintf(stderr, + "WARNING: publishing incomplete render (error triangles=%zu, " + "budget-incomplete unresolved triangles=%zu).\n", + total->error, total->budget_incomplete_triangles); + return 1; + } + fprintf(stderr, + "Incomplete render refused: %zu error triangle(s), %zu " + "budget-incomplete unresolved triangle(s). Raise the retry budget or " + "pass --allow-incomplete for a diagnostic output.\n", + total->error, total->budget_incomplete_triangles); + return 0; +} + static GeodesicTraceConfig trace_config(const Settings *s) { + /* The normal dark terminal is the camera-relative local energy growth + * L - L0, independent of the spacetime backend. A constant camera boost + * cancels; the photon energy and frequency ratio are never reset. */ + const ThresholdPolicy threshold = {.kind = THRESHOLD_LOG_ENERGY_GROWTH, + .value = s->dark_threshold, + .policy_version = 3}; #ifdef SPACETIME_SCHWARZSCHILD - (void)s; /* The directed worldtube crossing makes an escaping ray traverse the * interior as a round trip from the entry sphere (in, turn, back out), * rather than the old one-way stop at the first outside sample. The step * budget must cover roughly twice the escape sphere plus margin. */ return (GeodesicTraceConfig){.coordinate_time_step = 0.1, .max_steps = 65536, - .capture_log_alpha_p0 = 8.0}; + .threshold = threshold}; #elif defined(SPACETIME_ALCUBIERRE) const double step = alcubierre_time_step(s); const double budget = alcubierre_step_budget(s); @@ -886,11 +993,12 @@ static GeodesicTraceConfig trace_config(const Settings *s) { if (max_steps < 1024u) max_steps = 1024u; return (GeodesicTraceConfig){.coordinate_time_step = step, - .max_steps = max_steps}; + .max_steps = max_steps, + .threshold = threshold}; #else - (void)s; return (GeodesicTraceConfig){.coordinate_time_step = 1.0, - .max_steps = 2048}; + .max_steps = 2048, + .threshold = threshold}; #endif } @@ -947,17 +1055,13 @@ static int build_observer(const Settings *s, const SpacetimeSource *spacetime, camera.position[i] = s->observer_position[i]; camera.velocity[i] = s->observer_velocity[i]; } - const SpacetimeRayStatus camera_status = - spacetime_classify(spacetime, camera.coordinate_time, camera.position); - if (camera_status == SPACETIME_RAY_CAPTURED) { - fputs("Camera position is inside the backend capture cutoff or invalid.\n", stderr); - return -1; - } - /* A camera outside the escape sphere is supported by the asymptotic - * exterior module for every declared end kind; unsupported exteriors are - * reported through the ray endpoints instead. */ + /* Camera legality depends only on metric availability, a timelike + * four-velocity, time orientation, and an orthonormal tetrad. A camera + * inside a horizon or an old spatial cutoff is a normal rendering target; + * its position never decides a ray's terminal category. */ MetricData metric; - if (spacetime_eval(spacetime, camera.coordinate_time, camera.position, &metric)) { + if (spacetime_eval(spacetime, camera.coordinate_time, camera.position, + &metric) != SPACETIME_POINT_OK) { fputs("Could not evaluate metric at the camera event.\n", stderr); return -1; } @@ -1032,13 +1136,26 @@ static int render_observer_frame(const Settings *s, StarCatalog *catalog, "%zu vertices, %zu triangles.\n", omp_get_wtime() - refinement_start, mesh.vertex_count, mesh.triangle_count); + FrameBoundaryStats boundary_totals; + { + RefinementConfig boundary_config = s->refinement; + frame_retry_config_defaults(&boundary_config, &trace); + FrameLensMesh *meshes[1] = {&mesh}; + report_boundary_stats(s, meshes, 1, &boundary_config, &boundary_totals); + } + if (!boundary_allows_publish(s, &boundary_totals)) { + frame_lens_mesh_destroy(&mesh); + free(hdr); + return -1; + } if (s->lens_map_output_path != NULL) { const LensMapFrame map_frame = {.frame_id = 0, .coordinate_time = 0.0, .proper_time = 0.0, .mesh = mesh}; + const LensMapProvenance provenance = lens_map_provenance(s, &trace); if (lens_map_write(s->lens_map_output_path, s->width, s->height, - s->horizontal_fov_deg, &map_frame, 1)) { + s->horizontal_fov_deg, &provenance, &map_frame, 1)) { fprintf(stderr, "Failed to write lens map: %s\n", s->lens_map_output_path); frame_lens_mesh_destroy(&mesh); free(hdr); @@ -1113,9 +1230,11 @@ static int trace_movie_generation(Movie *movie, const Settings *s, RayPool rays = {0}; size_t ray_count = 0; size_t total_added = 0; + RefinementConfig effective = s->refinement; + frame_retry_config_defaults(&effective, trace); for (size_t f = 0; f < movie->frame_count; ++f) { const int prepared = frame_lens_mesh_prepare_generation( - &movie->frames[f].mesh, &s->refinement); + &movie->frames[f].mesh, &effective); if (prepared < 0) return -1; ray_count += (size_t)prepared; } @@ -1128,12 +1247,26 @@ static int trace_movie_generation(Movie *movie, const Settings *s, for (size_t f = 0; f < movie->frame_count; ++f) { size_t count = 0; const FrameSample *samples = frame_lens_mesh_samples(&movie->frames[f].mesh, &count); - for (size_t sample = 0; sample < count; ++sample) - if (ray_pool_append(&rays, &movie->frames[f].observer, - samples[sample].vertex.camera_direction, f, sample)) { + for (size_t sample = 0; sample < count; ++sample) { + const FrameSample *fs = &samples[sample]; + if (fs->cached) continue; + int rc; + if (fs->kind == FRAME_SAMPLE_RETRY) { + const LensVertex *v = &fs->vertex; + rc = ray_pool_append_continuation( + &rays, f, sample, v->continuation_t, v->continuation_x, + v->continuation_Pi, v->continuation_log_alpha_p0, + v->continuation_log_alpha_p0_0, v->continuation_steps, + fs->step_limit); + } else { + rc = ray_pool_append(&rays, &movie->frames[f].observer, + fs->vertex.camera_direction, f, sample); + } + if (rc) { ray_pool_destroy(&rays); return -1; } + } } ray_pool_preroute(&rays, spacetime); double slab_hi = movie->frames[movie->frame_count - 1].coordinate_time; @@ -1141,29 +1274,35 @@ static int trace_movie_generation(Movie *movie, const Settings *s, while (ray_pool_has_live(&rays)) { const double slab_lo = slab_hi - s->slab_duration; MetricSlab *slab = NULL; - size_t pending_before, active_before, terminated_before, failed_before; + size_t pending_before, active_before, terminated_before, unresolved_before, + failed_before; ray_pool_status_counts(&rays, &pending_before, &active_before, - &terminated_before, &failed_before); + &terminated_before, &unresolved_before, + &failed_before); if (spacetime_load_slab(spacetime, slab_hi, slab_lo, &slab)) { ray_pool_destroy(&rays); return -1; } ray_pool_activate_in_time_range(&rays, slab); - size_t pending_active, active_active, terminated_active, failed_active; + size_t pending_active, active_active, terminated_active, unresolved_active, + failed_active; ray_pool_status_counts(&rays, &pending_active, &active_active, - &terminated_active, &failed_active); + &terminated_active, &unresolved_active, + &failed_active); ray_pool_advance_active(&rays, slab, trace); spacetime_free_slab(slab); - size_t pending_after, active_after, terminated_after, failed_after; + size_t pending_after, active_after, terminated_after, unresolved_after, + failed_after; ray_pool_status_counts(&rays, &pending_after, &active_after, - &terminated_after, &failed_after); + &terminated_after, &unresolved_after, &failed_after); if (s->verbose) fprintf(stderr, "Ray trace generation %zu, slab %zu [%.6g, %.6g]: activated %zu; " - "live %zu -> %zu, terminated %zu, failed %zu.\n", + "live %zu -> %zu, terminated %zu, unresolved %zu, failed %zu.\n", generation, ++slab_id, slab_hi, slab_lo, active_active - active_before, pending_before + active_before, - pending_after + active_after, terminated_after, failed_after); + pending_after + active_after, terminated_after, unresolved_after, + failed_after); slab_hi = slab_lo; } for (size_t i = 0; i < rays.count; ++i) @@ -1181,7 +1320,7 @@ static int trace_movie_generation(Movie *movie, const Settings *s, if (movie->frames[f].mesh.sample_count == 0) continue; const int added = frame_lens_mesh_finish_generation(&movie->frames[f].mesh, - &s->refinement); + &effective); if (added < 0) { fprintf(stderr, "Ray trace generation %zu: frame %zu refinement failed.\n", generation, movie->frames[f].frame_id); @@ -1338,6 +1477,21 @@ static int render_movie(const Settings *s, StarCatalog *catalog, if (traced == 0) break; } + { + RefinementConfig boundary_config = s->refinement; + frame_retry_config_defaults(&boundary_config, &trace); + FrameLensMesh **meshes = malloc(movie.frame_count * sizeof *meshes); + if (meshes == NULL) + goto done; + for (size_t i = 0; i < movie.frame_count; ++i) + meshes[i] = &movie.frames[i].mesh; + FrameBoundaryStats boundary_totals; + report_boundary_stats(s, meshes, movie.frame_count, &boundary_config, + &boundary_totals); + free(meshes); + if (!boundary_allows_publish(s, &boundary_totals)) + goto done; + } if (s->lens_map_output_path != NULL) { LensMapFrame *map_frames = calloc(movie.frame_count, sizeof *map_frames); if (map_frames == NULL) goto done; @@ -1346,9 +1500,11 @@ static int render_movie(const Settings *s, StarCatalog *catalog, .coordinate_time = movie.frames[i].coordinate_time, .proper_time = movie.frames[i].proper_time, .mesh = movie.frames[i].mesh}; + const LensMapProvenance provenance = lens_map_provenance(s, &trace); const int write_failed = lens_map_write(s->lens_map_output_path, s->width, s->height, s->horizontal_fov_deg, - map_frames, movie.frame_count); + &provenance, map_frames, + movie.frame_count); free(map_frames); if (write_failed) { fprintf(stderr, "Failed to write lens map: %s\n", s->lens_map_output_path); @@ -1468,10 +1624,26 @@ done: static int render_lens_map(const Settings *s, StarCatalog *catalog) { LensMap map = {0}; - if (lens_map_read(s->lens_map_input_path, &map)) { + if (lens_map_read(s->lens_map_input_path, NULL, &map)) { fprintf(stderr, "Failed to read or validate lens map: %s\n", s->lens_map_input_path); return -1; } + RefinementConfig replay_config = { + .max_level = map.provenance.max_level, + .min_edge_pixels = map.provenance.min_edge_pixels, + .min_area_pixels2 = map.provenance.min_area_pixels2, + .retry_step_increment = map.provenance.retry_step_increment, + .max_total_steps = map.provenance.max_total_steps}; + int replay_incomplete = 0; + /* Replay consumes stored decisions, not current CLI refinement defaults. */ + for (size_t f = 0; f < map.frame_count; ++f) { + FrameBoundaryStats completion; + frame_lens_mesh_boundary_stats(&map.frames[f].mesh, &replay_config, &completion); + replay_incomplete |= completion.error != 0 || completion.budget_incomplete_triangles != 0; + if (!boundary_allows_publish(s, &completion)) { + lens_map_destroy(&map); return -1; + } + } if ((s->width_specified && s->width != map.width) || (s->height_specified && s->height != map.height) || (s->fov_specified && fabs(s->horizontal_fov_deg - map.horizontal_fov_deg) > 1e-12)) { @@ -1609,7 +1781,7 @@ static int render_lens_map(const Settings *s, StarCatalog *catalog) { if (prepare_movie_output_job(s, &map.frames[i].mesh, hdr, map.width, map.height, &output_paths, (size_t)map.frames[i].frame_id, images, - catalog->count, "; imported lens map", + catalog->count, replay_incomplete ? "; INCOMPLETE imported lens map" : "; imported lens map", &psf_stats, &frame_timing, &job)) { free(hdr); result = -1; @@ -1632,7 +1804,7 @@ static int render_lens_map(const Settings *s, StarCatalog *catalog) { const int write_result = write_frame_outputs( s, &map.frames[i].mesh, hdr, map.width, map.height, map.horizontal_fov_deg, &output_paths, images, catalog->count, - "; imported lens map", NULL); + replay_incomplete ? "; INCOMPLETE imported lens map" : "; imported lens map", NULL); free(hdr); report_psf_splat(s, &psf_stats); if (write_result) { result = -1; break; } @@ -1688,7 +1860,8 @@ int main(int argc, char **argv) { "Usage: %s [--catalog PATH | --all-sky-catalog DIR] [--output PATH] [--width N] [--height " "N] [--fov-deg D] [--look-ra-deg D] [--look-dec-deg D] " "[--lens-map-input FILE | --lens-map-output FILE] " - "[--exposure E] [--tone-map softclip|reinhard] [--tone-map-p P] " + "[--exposure E] [--dark-threshold T] " + "[--tone-map softclip|reinhard] [--tone-map-p P] " "[--sensor-bloom-limit E --sensor-bloom-transfer e] " "[--observer-radius R | --observer-position X Y Z] " "[--observer-velocity VX VY VZ] [--camera-roll-deg ANGLE] " @@ -1696,7 +1869,7 @@ int main(int argc, char **argv) { "[--max-magnification M] [--max-cache-psf-flux F] " "[--psf-relative-tail R] [--psf-min-y Y] " "[--psf-direct] [--fast-mode --fast-supersample N " - "--fast-deposit nearest|bilinear] [--verbose] " + "--fast-deposit nearest|bilinear] [--verbose] [--allow-incomplete] " #ifdef ENABLE_HDR_OUTPUT "[--hdr-output] " #endif diff --git a/src/ray.c b/src/ray.c index 68902ab..6241640 100644 --- a/src/ray.c +++ b/src/ray.c @@ -14,7 +14,9 @@ int ray_pool_init(RayPool *p, size_t capacity) { if (!(RAY_ALLOC(t) && RAY_ALLOC(x0) && RAY_ALLOC(x1) && RAY_ALLOC(x2) && RAY_ALLOC(p0) && RAY_ALLOC(p1) && RAY_ALLOC(p2) && RAY_ALLOC(observer) && RAY_ALLOC(direction0) && RAY_ALLOC(direction1) && RAY_ALLOC(direction2) && - RAY_ALLOC(log_alpha_p0) && RAY_ALLOC(activate_t) && RAY_ALLOC(steps) && + RAY_ALLOC(log_alpha_p0) && RAY_ALLOC(log_alpha_p0_0) && + RAY_ALLOC(activate_t) && RAY_ALLOC(steps) && + RAY_ALLOC(step_limit) && RAY_ALLOC(continuation) && RAY_ALLOC(frame_id) && RAY_ALLOC(vertex_id) && RAY_ALLOC(status) && RAY_ALLOC(endpoint))) { ray_pool_destroy(p); @@ -43,7 +45,55 @@ int ray_pool_append(RayPool *p, const ObserverState *observer, p->status[i] = RAY_POOL_PENDING; p->endpoint[i] = (RayEndpoint){.magnification = 1.0, .end_id = SPACETIME_END_NONE, - .status = RAY_ENDPOINT_INTEGRATION_FAILURE}; + .outcome = RAY_OUTCOME_INCOMPLETE, + .reason = RAY_REASON_NONE, + .stop_coordinate_time = NAN, + .accepted_steps = 0, + .final_x = {NAN, NAN, NAN}, + .final_Pi = {NAN, NAN, NAN}, + .final_log_alpha_p0 = NAN, + .final_log_alpha_p0_0 = NAN, + .threshold_value = NAN}; + p->step_limit[i] = 0; + p->continuation[i] = 0; + p->log_alpha_p0_0[i] = 0.0; + ++p->count; + return 0; +} + +int ray_pool_append_continuation(RayPool *p, size_t frame_id, + size_t vertex_id, double t, const double x[3], + const double Pi[3], double log_alpha_p0, + double log_alpha_p0_0, unsigned int steps, + unsigned int limit) { + if (p == NULL || p->count == p->capacity || x == NULL || Pi == NULL) + return -1; + const size_t i = p->count; + p->t[i] = t; + p->activate_t[i] = t; + p->observer[i] = NULL; + p->direction0[i] = p->direction1[i] = p->direction2[i] = 0.0; + p->x0[i] = x[0]; p->x1[i] = x[1]; p->x2[i] = x[2]; + p->p0[i] = Pi[0]; p->p1[i] = Pi[1]; p->p2[i] = Pi[2]; + p->log_alpha_p0[i] = log_alpha_p0; + p->log_alpha_p0_0[i] = log_alpha_p0_0; + p->steps[i] = steps; + p->step_limit[i] = limit; + p->continuation[i] = 1; + p->frame_id[i] = frame_id; + p->vertex_id[i] = vertex_id; + p->status[i] = RAY_POOL_PENDING; + p->endpoint[i] = (RayEndpoint){.magnification = 1.0, + .end_id = SPACETIME_END_NONE, + .outcome = RAY_OUTCOME_INCOMPLETE, + .reason = RAY_REASON_NONE, + .stop_coordinate_time = NAN, + .accepted_steps = steps, + .final_x = {NAN, NAN, NAN}, + .final_Pi = {NAN, NAN, NAN}, + .final_log_alpha_p0 = NAN, + .final_log_alpha_p0_0 = log_alpha_p0_0, + .threshold_value = NAN}; ++p->count; return 0; } @@ -53,7 +103,7 @@ void ray_pool_preroute(RayPool *p, const SpacetimeSource *source) { return; #pragma omp parallel for schedule(static) for (size_t i = 0; i < p->count; ++i) { - if (p->status[i] != RAY_POOL_PENDING) + if (p->status[i] != RAY_POOL_PENDING || p->continuation[i]) continue; AsymptoticRoute route; const AsymptoticStatus status = asymptotic_route_camera( @@ -61,13 +111,17 @@ void ray_pool_preroute(RayPool *p, const SpacetimeSource *source) { (double[]){p->direction0[i], p->direction1[i], p->direction2[i]}, &route); if (status == ASYMPTOTIC_UNSUPPORTED || status == ASYMPTOTIC_INVALID) { - p->endpoint[i].status = RAY_ENDPOINT_INVALID; + p->endpoint[i].outcome = RAY_OUTCOME_INCOMPLETE; + p->endpoint[i].reason = status == ASYMPTOTIC_UNSUPPORTED + ? RAY_REASON_UNSUPPORTED + : RAY_REASON_PROTOCOL_ERROR; p->endpoint[i].end_id = route.end_id; p->status[i] = RAY_POOL_FAILED; continue; } if (status == ASYMPTOTIC_TIME_RANGE_EXHAUSTED) { - p->endpoint[i].status = RAY_ENDPOINT_TIME_RANGE_EXHAUSTED; + p->endpoint[i].outcome = RAY_OUTCOME_INCOMPLETE; + p->endpoint[i].reason = RAY_REASON_TIME_RANGE_EXHAUSTED; p->endpoint[i].end_id = route.end_id; p->status[i] = RAY_POOL_TERMINATED; continue; @@ -77,19 +131,22 @@ void ray_pool_preroute(RayPool *p, const SpacetimeSource *source) { p->endpoint[i].n_infinity[axis] = route.n_infinity[axis]; p->endpoint[i].frequency_ratio = route.frequency_ratio; p->endpoint[i].end_id = route.end_id; - p->endpoint[i].status = RAY_ENDPOINT_ESCAPED; + p->endpoint[i].outcome = RAY_OUTCOME_ESCAPED; + p->endpoint[i].reason = RAY_REASON_NONE; p->status[i] = RAY_POOL_TERMINATED; continue; } if (route.kind == ASYMPTOTIC_ROUTE_TIME_RANGE_EXHAUSTED) { - p->endpoint[i].status = RAY_ENDPOINT_TIME_RANGE_EXHAUSTED; + p->endpoint[i].outcome = RAY_OUTCOME_INCOMPLETE; + p->endpoint[i].reason = RAY_REASON_TIME_RANGE_EXHAUSTED; p->endpoint[i].end_id = route.end_id; p->status[i] = RAY_POOL_TERMINATED; continue; } if (route.kind != ASYMPTOTIC_ROUTE_INSIDE && route.kind != ASYMPTOTIC_ROUTE_ENTRY) { - p->endpoint[i].status = RAY_ENDPOINT_INVALID; + p->endpoint[i].outcome = RAY_OUTCOME_INCOMPLETE; + p->endpoint[i].reason = RAY_REASON_PROTOCOL_ERROR; p->status[i] = RAY_POOL_FAILED; continue; } @@ -101,6 +158,8 @@ void ray_pool_preroute(RayPool *p, const SpacetimeSource *source) { p->p1[i] = route.Pi[1]; p->p2[i] = route.Pi[2]; p->log_alpha_p0[i] = route.log_alpha_p0; + /* Camera-event reference, distinct from the entry-state L. */ + p->log_alpha_p0_0[i] = route.log_alpha_p0_camera; } } @@ -110,7 +169,9 @@ void ray_pool_activate_in_time_range(RayPool *p, const MetricSlab *slab) { p->activate_t[i] <= slab->t_lo) continue; p->t[i] = p->activate_t[i]; - p->steps[i] = 0; + /* Continuation rays keep the accepted-step count they already consumed. */ + if (!p->continuation[i]) + p->steps[i] = 0; p->status[i] = RAY_POOL_ACTIVE; } } @@ -129,16 +190,22 @@ void ray_pool_advance_active(RayPool *p, const MetricSlab *slab, .x = {p->x0[i], p->x1[i], p->x2[i]}, .Pi = {p->p0[i], p->p1[i], p->p2[i]}, .log_alpha_p0 = p->log_alpha_p0[i], + .log_alpha_p0_0 = p->log_alpha_p0_0[i], .steps = p->steps[i]}; + GeodesicTraceConfig per_ray = *config; + if (p->step_limit[i] != 0) + per_ray.max_steps = p->step_limit[i]; const GeodesicAdvanceResult result = - geodesic_advance_past_ray(slab, &s, slab->t_lo, config, &p->endpoint[i]); + geodesic_advance_past_ray(slab, &s, slab->t_lo, &per_ray, &p->endpoint[i]); p->t[i] = s.coordinate_time; p->x0[i] = s.x[0]; p->x1[i] = s.x[1]; p->x2[i] = s.x[2]; p->p0[i] = s.Pi[0]; p->p1[i] = s.Pi[1]; p->p2[i] = s.Pi[2]; p->log_alpha_p0[i] = s.log_alpha_p0; p->steps[i] = s.steps; if (result == GEODESIC_ADVANCE_TERMINATED) - p->status[i] = RAY_POOL_TERMINATED; + p->status[i] = p->endpoint[i].outcome == RAY_OUTCOME_UNRESOLVED + ? RAY_POOL_UNRESOLVED + : RAY_POOL_TERMINATED; else if (result == GEODESIC_ADVANCE_FAILED) p->status[i] = RAY_POOL_FAILED; } @@ -159,7 +226,9 @@ void ray_pool_destroy(RayPool *p) { free(p->t); free(p->x0); free(p->x1); free(p->x2); free(p->observer); free(p->direction0); free(p->direction1); free(p->direction2); free(p->p0); free(p->p1); free(p->p2); free(p->log_alpha_p0); - free(p->activate_t); free(p->steps); free(p->frame_id); free(p->vertex_id); + free(p->log_alpha_p0_0); + free(p->activate_t); free(p->steps); free(p->step_limit); + free(p->continuation); free(p->frame_id); free(p->vertex_id); free(p->status); free(p->endpoint); *p = (RayPool){0}; diff --git a/src/ray.h b/src/ray.h index eb38a70..f7439f6 100644 --- a/src/ray.h +++ b/src/ray.h @@ -10,11 +10,12 @@ typedef enum { RAY_POOL_PENDING, RAY_POOL_ACTIVE, RAY_POOL_TERMINATED, + RAY_POOL_UNRESOLVED, /* trustworthy but budget-exhausted; retryable */ RAY_POOL_FAILED } RayPoolStatus; typedef struct { - double *t, *x0, *x1, *x2, *p0, *p1, *p2, *log_alpha_p0; + double *t, *x0, *x1, *x2, *p0, *p1, *p2, *log_alpha_p0, *log_alpha_p0_0; /* Coordinate time at which the pre-routed interior state becomes valid. * For a camera inside a worldtube this equals the camera time; for an * exterior hit it is the earlier entry time. */ @@ -22,6 +23,11 @@ typedef struct { const ObserverState **observer; double *direction0, *direction1, *direction2; unsigned int *steps; + /* Per-ray total accepted-step limit; zero means use the trace config. */ + unsigned int *step_limit; + /* Nonzero for a retry that resumes from a saved state instead of from the + * camera; such rays are not pre-routed. */ + uint8_t *continuation; size_t *frame_id, *vertex_id; uint8_t *status; RayEndpoint *endpoint; @@ -32,6 +38,13 @@ int ray_pool_init(RayPool *pool, size_t capacity); int ray_pool_append(RayPool *pool, const ObserverState *observer, const double direction[3], size_t frame_id, size_t vertex_id); +/* Append a retry that resumes an UNRESOLVED ray from its last accepted state. + * `limit` is the new total accepted-step budget for this ray. */ +int ray_pool_append_continuation(RayPool *pool, size_t frame_id, + size_t vertex_id, double t, const double x[3], + const double Pi[3], double log_alpha_p0, + double log_alpha_p0_0, unsigned int steps, + unsigned int limit); /* Pre-route every still-PENDING ray once, before the slab sweep. */ void ray_pool_preroute(RayPool *pool, const SpacetimeSource *source); void ray_pool_activate_in_time_range(RayPool *pool, const MetricSlab *slab); diff --git a/src/spacetime.h b/src/spacetime.h index c5e0228..34ad702 100644 --- a/src/spacetime.h +++ b/src/spacetime.h @@ -14,10 +14,23 @@ typedef struct { double d_gamma[3][3][3]; /* d_gamma[spatial derivative][j][k] */ } MetricData; +/* Result of evaluating the metric at one event. `OK` is zero so that legacy + * `if (eval(...))` call sites keep working. These codes describe data + * availability only; they never express a physical capture. */ +typedef enum { + SPACETIME_POINT_OK = 0, + SPACETIME_POINT_TIME_UNAVAILABLE, + SPACETIME_POINT_OUT_OF_DOMAIN, + SPACETIME_POINT_INVALID_METRIC, + SPACETIME_POINT_INTERNAL_ERROR +} SpacetimePointStatus; + +/* Optional legacy region query for backends that declare no asymptotic end. + * It can only report ACTIVE or ESCAPED; it can never report a physical + * capture, and it is not required by spacetime_source_finalize(). */ typedef enum { SPACETIME_RAY_ACTIVE, - SPACETIME_RAY_ESCAPED, - SPACETIME_RAY_CAPTURED + SPACETIME_RAY_ESCAPED } SpacetimeRayStatus; /* Stable identifier for one asymptotic end (infinity) of a backend. Backends @@ -66,15 +79,15 @@ struct MetricSlab { }; typedef struct { - int (*eval)(const SpacetimeSource *source, double t, const double x[3], - MetricData *metric); + SpacetimePointStatus (*eval)(const SpacetimeSource *source, double t, + const double x[3], MetricData *metric); SpacetimeRayStatus (*classify)(const SpacetimeSource *source, double t, const double x[3]); int (*load_slab)(const SpacetimeSource *source, double t_hi, double t_lo, MetricSlab **out); void (*free_slab)(MetricSlab *slab); - int (*eval_slab)(const MetricSlab *slab, double t, const double x[3], - MetricData *metric); + SpacetimePointStatus (*eval_slab)(const MetricSlab *slab, double t, + const double x[3], MetricData *metric); SpacetimeRayStatus (*classify_slab)(const MetricSlab *slab, double t, const double x[3]); /* Declared asymptotic ends and their moving escape worldtubes. Backends @@ -108,8 +121,7 @@ struct SpacetimeSource { int spacetime_create_default(SpacetimeSource *source); int spacetime_create_minkowski(SpacetimeSource *source, double escape_radius); int spacetime_create_schwarzschild_ks(SpacetimeSource *source, double mass, - double escape_radius, - double capture_radius); + double escape_radius); /* Moving Alcubierre bubble with x_s(t) = vs*t and x_s(0) = 0. Requires * |vs| < 1, R > 0, and sigma > 0. */ int spacetime_create_alcubierre(SpacetimeSource *source, double vs, @@ -118,15 +130,15 @@ int spacetime_create_alcubierre(SpacetimeSource *source, double vs, * callers size their integration step budget. */ double spacetime_alcubierre_escape_radius(double radius, double sigma); void spacetime_destroy(SpacetimeSource *source); -int spacetime_eval(const SpacetimeSource *source, double t, const double x[3], - MetricData *metric); +SpacetimePointStatus spacetime_eval(const SpacetimeSource *source, double t, + const double x[3], MetricData *metric); SpacetimeRayStatus spacetime_classify(const SpacetimeSource *source, double t, const double x[3]); int spacetime_load_slab(const SpacetimeSource *source, double t_hi, double t_lo, MetricSlab **out); void spacetime_free_slab(MetricSlab *slab); -int spacetime_slab_eval(const MetricSlab *slab, double t, const double x[3], - MetricData *metric); +SpacetimePointStatus spacetime_slab_eval(const MetricSlab *slab, double t, + const double x[3], MetricData *metric); SpacetimeRayStatus spacetime_slab_classify(const MetricSlab *slab, double t, const double x[3]); size_t spacetime_asymptotic_end_count(const SpacetimeSource *source); diff --git a/src/spacetime_alcubierre.c b/src/spacetime_alcubierre.c index 0d7843a..6a75be2 100644 --- a/src/spacetime_alcubierre.c +++ b/src/spacetime_alcubierre.c @@ -72,8 +72,9 @@ static double alcubierre_shape_derivative(double r, double radius, * = -v_s (delta_jx d_i f + delta_ix d_j f) / 2, * where d_i differentiates at fixed t (only the spatial argument of f moves * with t). K encodes the time dependence required by the 3+1 null-ray RHS. */ -static int alcubierre_eval(const SpacetimeSource *source, double t, - const double x[3], MetricData *metric) { +static SpacetimePointStatus alcubierre_eval(const SpacetimeSource *source, + double t, const double x[3], + MetricData *metric) { const AlcubierreContext *context = source->context; const double vs = context->vs; const double dx = x[0] - vs * t; @@ -81,7 +82,7 @@ static int alcubierre_eval(const SpacetimeSource *source, double t, double df[3] = {0.0, 0.0, 0.0}; double f; if (!isfinite(r2)) - return -1; + return SPACETIME_POINT_INVALID_METRIC; const double r = sqrt(r2); *metric = (MetricData){ .alpha = 1.0, @@ -103,7 +104,7 @@ static int alcubierre_eval(const SpacetimeSource *source, double t, metric->K[i][j] = -0.5 * vs * ((j == 0 ? df[i] : 0.0) + (i == 0 ? df[j] : 0.0)); } - return 0; + return SPACETIME_POINT_OK; } /* A warp bubble has no curvature singularity or horizon for |v_s| < 1, so diff --git a/src/spacetime_common.c b/src/spacetime_common.c index 364e38a..49e231d 100644 --- a/src/spacetime_common.c +++ b/src/spacetime_common.c @@ -9,17 +9,18 @@ void spacetime_destroy(SpacetimeSource *source) { source->ops->destroy(source); } -int spacetime_eval(const SpacetimeSource *source, double t, const double x[3], - MetricData *metric) { - return source == NULL || source->ops == NULL - ? -1 +SpacetimePointStatus spacetime_eval(const SpacetimeSource *source, double t, + const double x[3], MetricData *metric) { + return source == NULL || source->ops == NULL || source->ops->eval == NULL + ? SPACETIME_POINT_INTERNAL_ERROR : source->ops->eval(source, t, x, metric); } SpacetimeRayStatus spacetime_classify(const SpacetimeSource *source, double t, const double x[3]) { - return source == NULL || source->ops == NULL - ? SPACETIME_RAY_CAPTURED + /* A missing or incomplete source is never reported as escaped. */ + return source == NULL || source->ops == NULL || source->ops->classify == NULL + ? SPACETIME_RAY_ACTIVE : source->ops->classify(source, t, x); } @@ -47,10 +48,12 @@ void spacetime_free_slab(MetricSlab *slab) { free(slab); } -int spacetime_slab_eval(const MetricSlab *slab, double t, const double x[3], - MetricData *metric) { - if (slab == NULL || t < slab->t_lo || t > slab->t_hi) - return -1; +SpacetimePointStatus spacetime_slab_eval(const MetricSlab *slab, double t, + const double x[3], MetricData *metric) { + if (slab == NULL || slab->source == NULL || slab->source->ops == NULL) + return SPACETIME_POINT_INTERNAL_ERROR; + if (t < slab->t_lo || t > slab->t_hi) + return SPACETIME_POINT_TIME_UNAVAILABLE; if (slab->source->ops->eval_slab != NULL) return slab->source->ops->eval_slab(slab, t, x, metric); return spacetime_eval(slab->source, t, x, metric); @@ -58,8 +61,10 @@ int spacetime_slab_eval(const MetricSlab *slab, double t, const double x[3], SpacetimeRayStatus spacetime_slab_classify(const MetricSlab *slab, double t, const double x[3]) { - if (slab == NULL || t < slab->t_lo || t > slab->t_hi) - return SPACETIME_RAY_CAPTURED; + if (slab == NULL || slab->source == NULL) + return SPACETIME_RAY_ACTIVE; + if (t < slab->t_lo || t > slab->t_hi) + return SPACETIME_RAY_ACTIVE; if (slab->source->ops->classify_slab != NULL) return slab->source->ops->classify_slab(slab, t, x); return spacetime_classify(slab->source, t, x); @@ -102,7 +107,7 @@ int spacetime_source_finalize(SpacetimeSource *source) { if (source == NULL || source->ops == NULL || source->context == NULL) return -1; const SpacetimeOps *ops = source->ops; - if (ops->eval == NULL || ops->classify == NULL || ops->destroy == NULL) + if (ops->eval == NULL || ops->destroy == NULL) return -1; const size_t count = spacetime_asymptotic_end_count(source); if (count == 0) diff --git a/src/spacetime_minkowski.c b/src/spacetime_minkowski.c index 3b6007b..4fcc40c 100644 --- a/src/spacetime_minkowski.c +++ b/src/spacetime_minkowski.c @@ -6,15 +6,16 @@ typedef struct { double escape_radius; } MinkowskiContext; -static int minkowski_eval(const SpacetimeSource *source, double t, - const double x[3], MetricData *metric) { +static SpacetimePointStatus minkowski_eval(const SpacetimeSource *source, + double t, const double x[3], + MetricData *metric) { (void)source; (void)t; (void)x; *metric = (MetricData){ .alpha = 1.0, .gamma = {{1.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0.0, 0.0, 1.0}}}; - return 0; + return SPACETIME_POINT_OK; } static SpacetimeRayStatus minkowski_classify(const SpacetimeSource *source, diff --git a/src/spacetime_schwarzschild.c b/src/spacetime_schwarzschild.c index 9137c25..42a12a3 100644 --- a/src/spacetime_schwarzschild.c +++ b/src/spacetime_schwarzschild.c @@ -6,22 +6,26 @@ typedef struct { double mass; double escape_radius; - double capture_radius; } SchwarzschildKsContext; /* Schwarzschild in ingoing Cartesian Kerr--Schild coordinates: * g_mu_nu = eta_mu_nu + (2 M / r) l_mu l_nu, l_mu = (1, x_i / r). - * These slices are regular at r = 2 M; only the physical r = 0 singularity - * is excluded by the conservative capture cutoff. */ -static int schwarzschild_ks_eval(const SpacetimeSource *source, double t, - const double x[3], MetricData *metric) { + * These slices are regular at r = 2 M. Only r = 0 is a coordinate + * singularity; it is reported as a data/domain status, not as a physical + * capture. Normal dark endpoints come from the redshift threshold in the + * geodesic layer (see design section 18). */ +static SpacetimePointStatus schwarzschild_ks_eval(const SpacetimeSource *source, + double t, const double x[3], + MetricData *metric) { const SchwarzschildKsContext *context = source->context; double r2 = 0.0; (void)t; for (int i = 0; i < 3; ++i) r2 += x[i] * x[i]; - if (!isfinite(r2) || r2 <= 0.0) - return -1; + if (!isfinite(r2)) + return SPACETIME_POINT_INVALID_METRIC; + if (r2 <= 0.0) + return SPACETIME_POINT_OUT_OF_DOMAIN; /* r = 0 coordinate singularity */ const double r = sqrt(r2); const double m = context->mass; const double f = 2.0 * m / r; @@ -79,11 +83,14 @@ static int schwarzschild_ks_eval(const SpacetimeSource *source, double t, } metric->K[i][j] = (d_beta_cov_i_j - connection_term_ij + d_beta_cov_j_i - connection_term_ji) / - (2.0 * alpha); + (2.0 * alpha); } - return 0; + return SPACETIME_POINT_OK; } +/* Optional legacy region test: reports the escape sphere only. It never + * reports a physical capture; the normal dark terminal is the redshift + * threshold in the geodesic layer. */ static SpacetimeRayStatus schwarzschild_ks_classify( const SpacetimeSource *source, double t, const double x[3]) { const SchwarzschildKsContext *context = source->context; @@ -91,8 +98,8 @@ static SpacetimeRayStatus schwarzschild_ks_classify( (void)t; for (int i = 0; i < 3; ++i) r2 += x[i] * x[i]; - if (!isfinite(r2) || r2 <= context->capture_radius * context->capture_radius) - return SPACETIME_RAY_CAPTURED; + if (!isfinite(r2)) + return SPACETIME_RAY_ACTIVE; return r2 >= context->escape_radius * context->escape_radius ? SPACETIME_RAY_ESCAPED : SPACETIME_RAY_ACTIVE; @@ -150,16 +157,13 @@ static const SpacetimeOps schwarzschild_ks_ops = { }; int spacetime_create_schwarzschild_ks(SpacetimeSource *source, double mass, - double escape_radius, - double capture_radius) { - if (source == NULL || mass <= 0.0 || escape_radius <= 2.0 * mass || - capture_radius <= 0.0 || capture_radius >= 2.0 * mass || - capture_radius >= escape_radius) + double escape_radius) { + if (source == NULL || mass <= 0.0 || escape_radius <= 2.0 * mass) return -1; SchwarzschildKsContext *context = malloc(sizeof *context); if (context == NULL) return -1; - *context = (SchwarzschildKsContext){mass, escape_radius, capture_radius}; + *context = (SchwarzschildKsContext){mass, escape_radius}; source->ops = &schwarzschild_ks_ops; source->context = context; if (spacetime_source_finalize(source)) { @@ -170,5 +174,5 @@ int spacetime_create_schwarzschild_ks(SpacetimeSource *source, double mass, } int spacetime_create_default(SpacetimeSource *source) { - return spacetime_create_schwarzschild_ks(source, 1.0, 256.0, 1.5); + return spacetime_create_schwarzschild_ks(source, 1.0, 256.0); } diff --git a/tests/capture_psf.c b/tests/capture_psf.c index 498e752..9e86b05 100644 --- a/tests/capture_psf.c +++ b/tests/capture_psf.c @@ -98,7 +98,7 @@ int main(int argc, char **argv) { } printf("selection: raw_magnification=[%.17g,%.17g) exposure=%.17g max_cache_flux=%.17g min_y=%.17g\n",selected_min_mag,selected_max_mag,exposure,max_flux,min_y); LensMap map={0}; StarCatalog catalog={0}; - if (lens_map_read(argv[1],&map) || map.frame_count!=1 || + if (lens_map_read(argv[1],NULL,&map) || map.frame_count!=1 || last>map.frames[0].mesh.triangle_count || catalog_load_csv(&catalog,argv[2]) || catalog.count>32768) return 2; if (blackbody_backend_init(NULL, 0, NAN, NAN, NULL, stderr)) return 1; diff --git a/tests/test_alcubierre.c b/tests/test_alcubierre.c index 8bf1716..a63d4c6 100644 --- a/tests/test_alcubierre.c +++ b/tests/test_alcubierre.c @@ -214,8 +214,8 @@ int main(void) { n[k] /= norm; const RayEndpoint r0 = geodesic_trace_past(&source, &o0, n, &trace); const RayEndpoint r1 = geodesic_trace_past(&source, &o1, n, &trace); - CHECK(r0.status == RAY_ENDPOINT_ESCAPED); - CHECK(r1.status == RAY_ENDPOINT_ESCAPED); + CHECK(r0.outcome == RAY_OUTCOME_ESCAPED); + CHECK(r1.outcome == RAY_OUTCOME_ESCAPED); for (int k = 0; k < 3; ++k) CHECK(fabs(r0.n_infinity[k] - r1.n_infinity[k]) < 1e-6); CHECK(fabs(r0.frequency_ratio - r1.frequency_ratio) < 1e-6); @@ -236,7 +236,7 @@ int main(void) { const ObserverState observer = observer_fixed_at_origin(); const RayEndpoint ray = geodesic_trace_past( &flat, &observer, (double[]){1, 0, 0}, &trace); - CHECK(ray.status == RAY_ENDPOINT_ESCAPED); + CHECK(ray.outcome == RAY_OUTCOME_ESCAPED); CHECK(fabs(ray.n_infinity[0]) < 1e-12); CHECK(fabs(ray.n_infinity[1]) < 1e-12); CHECK(fabs(ray.n_infinity[2] + 1.0) < 1e-12); @@ -277,7 +277,7 @@ int main(void) { .max_steps = 2500000u}; const RayEndpoint ray = geodesic_trace_past( &fast, &observer, (double[]){-1, 0, 0}, &trace); - CHECK(ray.status == RAY_ENDPOINT_ESCAPED); + CHECK(ray.outcome == RAY_OUTCOME_ESCAPED); spacetime_destroy(&fast); } @@ -305,7 +305,7 @@ int main(void) { .coordinate_time_step = 0.08 / (1 << level), .max_steps = 1u << 20}; const RayEndpoint ray = geodesic_trace_past(&source, &observer, n, &trace); - CHECK(ray.status == RAY_ENDPOINT_ESCAPED); + CHECK(ray.outcome == RAY_OUTCOME_ESCAPED); if (level > 0) { double error = 0.0; for (int k = 0; k < 3; ++k) { diff --git a/tests/test_asymptotic.c b/tests/test_asymptotic.c index b4b786f..44d85a5 100644 --- a/tests/test_asymptotic.c +++ b/tests/test_asymptotic.c @@ -52,8 +52,9 @@ typedef struct { int sample_call_count; } SyntheticContext; -static int synthetic_eval(const SpacetimeSource *source, double t, - const double x[3], MetricData *metric) { +static SpacetimePointStatus synthetic_eval(const SpacetimeSource *source, + double t, const double x[3], + MetricData *metric) { (void)source; (void)t; (void)x; @@ -61,7 +62,7 @@ static int synthetic_eval(const SpacetimeSource *source, double t, .gamma = {{1.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0.0, 0.0, 1.0}}}; - return 0; + return SPACETIME_POINT_OK; } static SpacetimeRayStatus synthetic_classify(const SpacetimeSource *source, @@ -233,7 +234,7 @@ static void test_fixed_sphere(void) { CHECK(asymptotic_finish_escape(&source, 0, 0.0, (double[]){10.0, 0.0, 0.0}, (double[]){-1.0, 0.0, 0.0}, 0.0, &endpoint) == ASYMPTOTIC_OK && - endpoint.status == RAY_ENDPOINT_ESCAPED && endpoint.end_id == 0, + endpoint.outcome == RAY_OUTCOME_ESCAPED && endpoint.end_id == 0, "finish outward crossing"); CHECK(fabs(endpoint.n_infinity[0] - 1.0) < 1e-12 && fabs(endpoint.frequency_ratio - 1.0) < 1e-12, @@ -493,10 +494,11 @@ static void test_end_protocol_error(void) { RayEndpoint endpoint = {.frequency_ratio = 0.0, .magnification = 1.0, .end_id = SPACETIME_END_NONE, - .status = RAY_ENDPOINT_INVALID}; + .outcome = RAY_OUTCOME_INCOMPLETE}; CHECK(geodesic_advance_past_ray(slab, &state, -10.0, &config, &endpoint) == GEODESIC_ADVANCE_FAILED && - endpoint.status == RAY_ENDPOINT_INVALID, + endpoint.outcome == RAY_OUTCOME_INCOMPLETE && + endpoint.reason == RAY_REASON_PROTOCOL_ERROR, "advance rejects a declared-but-broken end without legacy"); spacetime_free_slab(slab); } @@ -517,7 +519,8 @@ static void test_interior_crossing_bisection_failure(void) { .max_steps = 2048}; const RayEndpoint endpoint = geodesic_trace_past( &source, &observer, (double[]){-1.0, 0.0, 0.0}, &config); - CHECK(endpoint.status == RAY_ENDPOINT_TIME_RANGE_EXHAUSTED && + CHECK(endpoint.outcome == RAY_OUTCOME_INCOMPLETE && + endpoint.reason == RAY_REASON_TIME_RANGE_EXHAUSTED && endpoint.end_id == 0, "interior crossing bisection propagates history exhaustion"); } @@ -550,7 +553,8 @@ static void test_interior_history_exhaustion(void) { .max_steps = 2048}; const RayEndpoint endpoint = geodesic_trace_past( &source, &observer, (double[]){-1.0, 0.0, 0.0}, &config); - CHECK(endpoint.status == RAY_ENDPOINT_TIME_RANGE_EXHAUSTED && + CHECK(endpoint.outcome == RAY_OUTCOME_INCOMPLETE && + endpoint.reason == RAY_REASON_TIME_RANGE_EXHAUSTED && endpoint.end_id == 0, "interior worldtube history exhaustion on a single trace"); @@ -567,7 +571,8 @@ static void test_interior_history_exhaustion(void) { ray_pool_activate_in_time_range(&pool, slab); CHECK(pool.status[0] == RAY_POOL_ACTIVE, "exhaustion ray activates"); ray_pool_advance_active(&pool, slab, &config); - CHECK(pool.endpoint[0].status == RAY_ENDPOINT_TIME_RANGE_EXHAUSTED && + CHECK(pool.endpoint[0].outcome == RAY_OUTCOME_INCOMPLETE && + pool.endpoint[0].reason == RAY_REASON_TIME_RANGE_EXHAUSTED && pool.endpoint[0].end_id == 0 && pool.status[0] == RAY_POOL_TERMINATED, "interior worldtube history exhaustion on a RayPool"); @@ -700,7 +705,7 @@ static void test_ray_pool_lifecycle(void) { pool.activate_t[0] < observer.coordinate_time - 1.0, "entry ray stays pending until entry time"); CHECK(pool.status[1] == RAY_POOL_TERMINATED && - pool.endpoint[1].status == RAY_ENDPOINT_ESCAPED, + pool.endpoint[1].outcome == RAY_OUTCOME_ESCAPED, "miss ray escapes during pre-route"); MetricSlab *early = NULL; diff --git a/tests/test_asymptotic_schwarzschild.c b/tests/test_asymptotic_schwarzschild.c index fd8bfb0..8754aa1 100644 --- a/tests/test_asymptotic_schwarzschild.c +++ b/tests/test_asymptotic_schwarzschild.c @@ -26,7 +26,7 @@ static double angle_between(const double a[3], const double b[3]) { static void test_round_trip(void) { SpacetimeSource source = {0}; - CHECK(spacetime_create_schwarzschild_ks(&source, 1.0, 256.0, 1.5) == 0, + CHECK(spacetime_create_schwarzschild_ks(&source, 1.0, 256.0) == 0, "create schwarzschild"); SpacetimeAsymptoticEnd end; CHECK(spacetime_asymptotic_end(&source, 0, &end) == 0, "end descriptor"); @@ -58,9 +58,9 @@ static void test_round_trip(void) { static void test_finish_matches_integration(void) { SpacetimeSource near = {0}, far = {0}; - CHECK(spacetime_create_schwarzschild_ks(&near, 1.0, 256.0, 1.5) == 0, + CHECK(spacetime_create_schwarzschild_ks(&near, 1.0, 256.0) == 0, "create near"); - CHECK(spacetime_create_schwarzschild_ks(&far, 1.0, 1.0e5, 1.5) == 0, + CHECK(spacetime_create_schwarzschild_ks(&far, 1.0, 1.0e5) == 0, "create far"); SpacetimeAsymptoticEnd end; CHECK(spacetime_asymptotic_end(&near, 0, &end) == 0, "near end"); @@ -97,12 +97,12 @@ static void test_finish_matches_integration(void) { CHECK(spacetime_load_slab(&far, 0.0, -1.0e6, &slab) == 0, "far slab"); RayEndpoint endpoint = {.frequency_ratio = 0, .magnification = 1.0, .end_id = SPACETIME_END_NONE, - .status = RAY_ENDPOINT_INVALID}; + .outcome = RAY_OUTCOME_INCOMPLETE}; const GeodesicAdvanceResult result = geodesic_advance_past_ray(slab, &state, -1.0e6, &config, &endpoint); spacetime_free_slab(slab); CHECK(result == GEODESIC_ADVANCE_TERMINATED && - endpoint.status == RAY_ENDPOINT_ESCAPED, + endpoint.outcome == RAY_OUTCOME_ESCAPED, "far integration escapes"); /* Pipeline check only: the far integration at step 5 and escape radius * 1e5 has its own O(1e-5..1e-3) error. Quantitative accuracy is checked @@ -152,7 +152,7 @@ static double transfer_dphi(double r, const void *context) { static void test_preroute_entry(void) { SpacetimeSource source = {0}; - CHECK(spacetime_create_schwarzschild_ks(&source, 1.0, 256.0, 1.5) == 0, + CHECK(spacetime_create_schwarzschild_ks(&source, 1.0, 256.0) == 0, "create schwarzschild"); const ObserverCamera camera = {.look_ra_deg = 0.0, .look_dec_deg = 0.0}; ObserverCamera positioned = camera; @@ -232,6 +232,51 @@ static void test_preroute_entry(void) { spacetime_destroy(&source); } +/* An external camera's dark-threshold reference must be the L at the camera + * event, not the worldtube entry energy. Changing only the worldtube radius + * must not change the reference but may change the entry L. */ +static void test_camera_reference_radius_independent(void) { + const double radii[2] = {128.0, 256.0}; + double reference[2], entry[2]; + for (int k = 0; k < 2; ++k) { + SpacetimeSource source = {0}; + CHECK(spacetime_create_schwarzschild_ks(&source, 1.0, radii[k]) == 0, + "create reference source"); + ObserverCamera positioned = {.look_ra_deg = 180.0, .look_dec_deg = 0.0}; + positioned.position[0] = 500.0; + MetricData metric; + CHECK(spacetime_eval(&source, 0.0, positioned.position, &metric) == 0, + "reference camera metric"); + ObserverState observer; + CHECK(observer_from_coordinate_camera(&metric, &positioned, &observer, + NULL) == OBSERVER_BUILD_OK, + "reference camera observer"); + const double direction[3] = {cos(0.05), sin(0.05), 0.0}; + MetricSlab *slab = NULL; + GeodesicRayState camera_state; + CHECK(spacetime_load_slab(&source, 0.0, -1.0, &slab) == 0, + "reference camera slab"); + CHECK(geodesic_initialize_past_ray(slab, &observer, direction, + &camera_state) == 0, + "reference camera state"); + spacetime_free_slab(slab); + AsymptoticRoute route; + CHECK(asymptotic_route_camera(&source, &observer, direction, &route) == + ASYMPTOTIC_OK && + route.kind == ASYMPTOTIC_ROUTE_ENTRY, + "reference ray enters the worldtube"); + CHECK(fabs(route.log_alpha_p0_camera - camera_state.log_alpha_p0) < 1e-12, + "route reference is the camera-event L"); + reference[k] = route.log_alpha_p0_camera; + entry[k] = route.log_alpha_p0; + spacetime_destroy(&source); + } + CHECK(fabs(reference[0] - reference[1]) < 1e-12, + "camera reference is worldtube-radius independent"); + CHECK(fabs(entry[0] - entry[1]) > 1e-6, + "entry energy depends on the worldtube radius"); +} + /* High-precision (mpmath, 60 digits) reference values fixed into the ordinary * C test: radial, complex-pair, three-real, grazing, and large-radius angle * cases. */ @@ -261,7 +306,7 @@ static void test_phi_reference_constants(void) { * measured ~1e-13. */ static void test_finish_reference_constants(void) { SpacetimeSource source = {0}; - CHECK(spacetime_create_schwarzschild_ks(&source, 1.0, 256.0, 1.5) == 0, + CHECK(spacetime_create_schwarzschild_ks(&source, 1.0, 256.0) == 0, "create schwarzschild"); SpacetimeAsymptoticEnd end; CHECK(spacetime_asymptotic_end(&source, 0, &end) == 0, "end"); @@ -313,7 +358,7 @@ static void test_turning_reference(void) { * both the KS time transfer and the entry direction construction. */ static void test_time_reference(void) { SpacetimeSource source = {0}; - CHECK(spacetime_create_schwarzschild_ks(&source, 1.0, 256.0, 1.5) == 0, + CHECK(spacetime_create_schwarzschild_ks(&source, 1.0, 256.0) == 0, "create schwarzschild"); SpacetimeAsymptoticEnd end; CHECK(spacetime_asymptotic_end(&source, 0, &end) == 0, "end"); @@ -365,7 +410,7 @@ static void test_time_reference(void) { * reachable from a camera outside R/M >= 64, so it is not tested here.) */ static void test_grazing_reference(void) { SpacetimeSource source = {0}; - CHECK(spacetime_create_schwarzschild_ks(&source, 1.0, 256.0, 1.5) == 0, + CHECK(spacetime_create_schwarzschild_ks(&source, 1.0, 256.0) == 0, "create schwarzschild"); SpacetimeAsymptoticEnd end; CHECK(spacetime_asymptotic_end(&source, 0, &end) == 0, "end"); @@ -411,7 +456,7 @@ static void test_grazing_reference(void) { * past-inward hit, and past-inward miss (turn before the worldtube). */ static void test_preroute_branches(void) { SpacetimeSource source = {0}; - CHECK(spacetime_create_schwarzschild_ks(&source, 1.0, 256.0, 1.5) == 0, + CHECK(spacetime_create_schwarzschild_ks(&source, 1.0, 256.0) == 0, "create schwarzschild"); SpacetimeAsymptoticEnd end; CHECK(spacetime_asymptotic_end(&source, 0, &end) == 0, "end"); @@ -480,6 +525,7 @@ int main(void) { test_round_trip(); test_finish_matches_integration(); test_preroute_entry(); + test_camera_reference_radius_independent(); test_phi_reference_constants(); test_finish_reference_constants(); test_turning_reference(); diff --git a/tests/test_camera_cli.py b/tests/test_camera_cli.py index 616a5a3..f5dde69 100644 --- a/tests/test_camera_cli.py +++ b/tests/test_camera_cli.py @@ -91,14 +91,16 @@ def image_payload(path, dimensions=(64, 48), allow_black=False): def map_vertices(path): data = path.read_bytes() assert data[:8] == b'GRLENS\x01\x00' - assert struct.unpack_from('= 0.01)) { + fprintf(stderr, + "last-step dark threshold regression failed (outcome=%d reason=%d " + "value=%.12g)\n", + endpoint.outcome, endpoint.reason, endpoint.threshold_value); + goto done; + } + } + /* A budget-exhausted ray is UNRESOLVED (retryable), keeps its last trusted + * state, and resolves when resumed from that state. */ + { + ObserverCamera camera = {.position = {30,0,0}, + .velocity = {-0.99999999,0,0}, .look_ra_deg = 0}; + ObserverState boosted; + MetricData m; + if (spacetime_eval(&spacetime, 0, camera.position, &m) || + observer_from_coordinate_camera(&m, &camera, &boosted, NULL)) goto done; + GeodesicRayState initial; + if (geodesic_initialize_past_ray_metric(&m, &boosted, + (double[]){1,0,0}, &initial) || + initial.log_alpha_p0 <= 8) goto done; + GeodesicTraceConfig disabled = trace; + disabled.threshold.kind = THRESHOLD_DISABLED; + RayEndpoint enabled = geodesic_trace_past(&spacetime, &boosted, + (double[]){1,0,0}, &trace); + RayEndpoint reference = geodesic_trace_past(&spacetime, &boosted, + (double[]){1,0,0}, &disabled); + if (enabled.outcome != RAY_OUTCOME_ESCAPED || + reference.outcome != RAY_OUTCOME_ESCAPED || + fabs(enabled.frequency_ratio/reference.frequency_ratio-1) > 1e-10) { + fputs("initial high-energy false-dark regression failed\n", stderr); goto done; + } + } + const GeodesicTraceConfig tiny = { + .coordinate_time_step = 0.1, + .max_steps = 30, + .threshold = {.kind = THRESHOLD_LOG_ENERGY_GROWTH, + .value = 8.0, + .policy_version = 3}}; + const RayEndpoint unresolved = geodesic_trace_past( + &spacetime, &observer, (double[]){cos(0.30), sin(0.30), 0.0}, &tiny); + if (unresolved.outcome != RAY_OUTCOME_UNRESOLVED || + unresolved.reason != RAY_REASON_BUDGET_EXHAUSTED || + unresolved.end_id != SPACETIME_END_NONE) { + fputs("budget-exhausted ray classification regression failed\n", stderr); + goto done; + } + const GeodesicRayState continuation = { + .coordinate_time = unresolved.stop_coordinate_time, + .x = {unresolved.final_x[0], unresolved.final_x[1], + unresolved.final_x[2]}, + .Pi = {unresolved.final_Pi[0], unresolved.final_Pi[1], + unresolved.final_Pi[2]}, + .log_alpha_p0 = unresolved.final_log_alpha_p0, + .log_alpha_p0_0 = unresolved.final_log_alpha_p0_0, + .steps = unresolved.accepted_steps}; + GeodesicTraceConfig more = tiny; + more.max_steps = 8192; + const RayEndpoint resumed = + geodesic_trace_past_from_state(&spacetime, &continuation, &more); + if (resumed.outcome != RAY_OUTCOME_ESCAPED) { + fprintf(stderr, + "resumed ray classification regression failed (outcome=%d)\n", + (int)resumed.outcome); goto done; } /* A coarse field covering the shadow must genuinely refine: its initial diff --git a/tests/test_termination_oracle.c b/tests/test_termination_oracle.c new file mode 100644 index 0000000..e23bedc --- /dev/null +++ b/tests/test_termination_oracle.c @@ -0,0 +1,293 @@ +/* + * Independent physics oracle for the ray-termination policy (plan P0). + * + * This test does not read production endpoints for its central assertions. + * It builds Schwarzschild-KS states independently and checks: + * 1. the two radial null branches dr/ds = 1 and dr/ds = (2M-r)/(2M+r); + * 2. the critical impact parameter b = 3 sqrt(3) M and photon sphere r = 3M; + * 3. camera energy normalization E_camera = 1 and tetrad orthonormality; + * 4. the threshold proxy identity ln(p^0) = L - ln(alpha) with + * L = ln(alpha p^0). + * + * The radial-branch assertions are integrated with the production RK4 RHS in + * src/geodesic.c so that an error in the 3+1 reduction is caught against a + * closed-form invariant rather than against a second copy of the same algebra. + */ +#include "asymptotic_schwarzschild.h" +#include "geodesic.h" +#include "observer.h" +#include "spacetime.h" + +#include +#include + +static int failures = 0; +#define CHECK(condition, message) \ + do { \ + if (!(condition)) { \ + fprintf(stderr, "FAIL %s:%d: %s\n", __FILE__, __LINE__, message); \ + ++failures; \ + } \ + } while (0) + +/* Static Eulerian orthonormal tetrad at x = (r0, 0, 0) for r0 > 0. At this + * point the KS spatial metric is diagonal, so the principal axes are already + * orthonormal (up to the radial scale sqrt(gamma_rr)). */ +static void radial_static_observer(const MetricData *metric, double r0, + ObserverState *out) { + *out = (ObserverState){.coordinate_time = 0.0, + .coordinate_position = {r0, 0.0, 0.0}}; + const double alpha = metric->alpha; + out->tetrad[0][0] = 1.0 / alpha; + for (int i = 0; i < 3; ++i) + out->tetrad[0][i + 1] = -metric->beta[i] / alpha; + const double radial_scale = sqrt(metric->gamma[0][0]); + out->tetrad[1][1] = 1.0 / radial_scale; + out->tetrad[2][2] = 1.0; + out->tetrad[3][3] = 1.0; +} + +/* Integrate a purely radial past ray with the production stepper and return + * its final state. Output endpoint is not inspected. */ +static int trace_radial(const SpacetimeSource *source, const ObserverState *o, + double direction, GeodesicRayState *state) { + MetricData metric; + if (spacetime_eval(source, o->coordinate_time, o->coordinate_position, + &metric) != SPACETIME_POINT_OK) + return -1; + const double n[3] = {direction, 0.0, 0.0}; + if (geodesic_initialize_past_ray_metric(&metric, o, n, state)) + return -1; + const GeodesicTraceConfig config = {.coordinate_time_step = 0.02, + .max_steps = 400, + .threshold = {.kind = THRESHOLD_DISABLED, .value = 0.0, .policy_version = 0}}; + MetricSlab *slab = NULL; + if (spacetime_load_slab(source, 0.0, -1000.0, &slab)) + return -1; + RayEndpoint endpoint = {.end_id = SPACETIME_END_NONE, + .outcome = RAY_OUTCOME_INCOMPLETE}; + const GeodesicAdvanceResult result = + geodesic_advance_past_ray(slab, state, -1000.0, &config, &endpoint); + spacetime_free_slab(slab); + return result == GEODESIC_ADVANCE_FAILED ? -1 : 0; +} + +static void test_radial_branches(void) { + SpacetimeSource source = {0}; + CHECK(spacetime_create_schwarzschild_ks(&source, 1.0, 256.0) == 0, + "create schwarzschild"); + const double r0 = 10.0; + MetricData metric; + CHECK(spacetime_eval(&source, 0.0, (double[]){r0, 0.0, 0.0}, &metric) == + SPACETIME_POINT_OK, + "metric at r0"); + ObserverState observer; + radial_static_observer(&metric, r0, &observer); + + /* Branch dr/ds = +1: the closed-form solution is r = r0 + s. */ + GeodesicRayState outward; + CHECK(trace_radial(&source, &observer, 1.0, &outward) == 0, "outward trace"); + const double s_out = -outward.coordinate_time; + const double invariant_out = outward.x[0] - r0 - s_out; + CHECK(fabs(invariant_out) < 1e-6, "outward branch r = r0 + s"); + + /* Branch dr/ds = (2M-r)/(2M+r): the closed-form invariant is + * (r-2M) + 4M ln(r-2M) + s = const. */ + GeodesicRayState inward; + CHECK(trace_radial(&source, &observer, -1.0, &inward) == 0, "inward trace"); + const double s_in = -inward.coordinate_time; + const double c0 = (r0 - 2.0) + 4.0 * log(r0 - 2.0); + const double c1 = (inward.x[0] - 2.0) + 4.0 * log(inward.x[0] - 2.0) + s_in; + CHECK(inward.x[0] > 2.0, "inward branch stays outside the horizon"); + CHECK(inward.x[0] < r0, "inward branch decreases r"); + CHECK(fabs(c1 - c0) < 1e-6, "inward branch closed-form invariant"); + + /* Both branches are time-reversal partners: the outward and inward states + * reach the same |dr/ds| magnitude in opposite senses at r0. */ + CHECK(outward.x[0] > r0, "outward branch increases r"); + spacetime_destroy(&source); +} + +/* The production dark policy is the camera-relative growth A_0 = L - L_0, + * independent of the backend. This oracle retains the stationary-KS + * conserved Killing energy A_K = L - ln|E_K| as an independent cross-check of + * the same ray: it verifies E_K conservation and the identity + * A_K - A_0 = -ln|alpha_0 - beta_0.Pi_0|. It is not the production + * criterion. */ +static void test_killing_energy_reference(void) { + SpacetimeSource source = {0}; + CHECK(spacetime_create_schwarzschild_ks(&source, 1.0, 256.0) == 0, + "create schwarzschild for killing reference"); + const double r0 = 10.0; + MetricData start_metric; + CHECK(spacetime_eval(&source, 0.0, (double[]){r0, 0.0, 0.0}, &start_metric) == + SPACETIME_POINT_OK, + "metric for killing reference"); + ObserverState observer; + radial_static_observer(&start_metric, r0, &observer); + GeodesicRayState start; + CHECK(geodesic_initialize_past_ray_metric(&start_metric, &observer, + (double[]){1.0, 0.0, 0.0}, + &start) == 0, + "initialize killing ray"); + double beta0 = 0.0; + for (int i = 0; i < 3; ++i) + beta0 += start_metric.beta[i] * start.Pi[i]; + const double ek0 = exp(start.log_alpha_p0) * (start_metric.alpha - beta0); + CHECK(isfinite(ek0) && fabs(ek0) > 0.0, "nonzero Killing energy"); + + GeodesicRayState end; + CHECK(trace_radial(&source, &observer, 1.0, &end) == 0, + "trace killing reference ray"); + MetricData end_metric; + CHECK(spacetime_eval(&source, end.coordinate_time, end.x, &end_metric) == + SPACETIME_POINT_OK, + "metric at killing reference end"); + double beta1 = 0.0; + for (int i = 0; i < 3; ++i) + beta1 += end_metric.beta[i] * end.Pi[i]; + const double ek1 = exp(end.log_alpha_p0) * (end_metric.alpha - beta1); + CHECK(fabs(ek1 / ek0 - 1.0) < 1e-6, + "Killing energy conserved along the geodesic"); + const double a0 = end.log_alpha_p0 - start.log_alpha_p0; + const double ak = end.log_alpha_p0 - log(fabs(ek1)); + const double predicted = -log(fabs(start_metric.alpha - beta0)); + CHECK(fabs((ak - a0) - predicted) < 1e-9, + "A_K - A_0 equals the initial boost factor"); + spacetime_destroy(&source); +} + +static void test_critical_parameters(void) { + const double b_crit = 3.0 * sqrt(3.0); + CHECK(!isfinite(asymptotic_schwarzschild_turning_rho(b_crit - 1e-6)), + "no turning point below b_crit"); + CHECK(!isfinite(asymptotic_schwarzschild_turning_rho(3.0)), + "no turning point for a deeply plunging ray"); + const double just_above = asymptotic_schwarzschild_turning_rho(b_crit + 1e-6); + CHECK(isfinite(just_above) && just_above > 3.0 && just_above < 3.01, + "turning radius approaches the photon sphere at b_crit"); + const double b6 = asymptotic_schwarzschild_turning_rho(6.0); + CHECK(isfinite(b6) && b6 > 3.0, "turning radius above the photon sphere"); + /* Verify the turning radius is an independent root of + * f(rho) = rho^3 - b^2 rho + 2 b^2. */ + const double residual = b6 * b6 * b6 - 36.0 * b6 + 72.0; + CHECK(fabs(residual) < 1e-9, "turning radius satisfies the radial equation"); +} + +static void test_observer_normalization(void) { + SpacetimeSource source = {0}; + CHECK(spacetime_create_schwarzschild_ks(&source, 1.0, 256.0) == 0, + "create schwarzschild for observer"); + ObserverCamera camera = {.coordinate_time = 0.0, + .position = {30.0, 0.0, 0.0}, + .velocity = {0.0, 0.0, 0.0}, + .look_ra_deg = 0.0, + .look_dec_deg = 0.0}; + MetricData metric; + ObserverState observer; + CHECK(spacetime_eval(&source, 0.0, camera.position, &metric) == + SPACETIME_POINT_OK, + "metric at camera"); + CHECK(observer_from_coordinate_camera(&metric, &camera, &observer, NULL) == + OBSERVER_BUILD_OK, + "build observer"); + /* Orthonormality of the production tetrad, independently of the geodesic + * layer: g(e_a, e_b) = diag(-1, 1, 1, 1). */ + for (int a = 0; a < 4; ++a) { + for (int b = 0; b < 4; ++b) { + const double *ea = observer.tetrad[a]; + const double *eb = observer.tetrad[b]; + double inner = -metric.alpha * metric.alpha * ea[0] * eb[0]; + for (int i = 0; i < 3; ++i) + for (int j = 0; j < 3; ++j) + inner += metric.gamma[i][j] * (ea[i + 1] + metric.beta[i] * ea[0]) * + (eb[j + 1] + metric.beta[j] * eb[0]); + const double expected = a == b ? (a == 0 ? -1.0 : 1.0) : 0.0; + CHECK(fabs(inner - expected) < 1e-10, "tetrad orthonormal"); + } + } + const double local[3] = {0.3, 0.5, 0.9}; + const double norm = sqrt(local[0] * local[0] + local[1] * local[1] + + local[2] * local[2]); + const double direction[3] = {local[0] / norm, local[1] / norm, + local[2] / norm}; + GeodesicRayState state; + CHECK(geodesic_initialize_past_ray_metric(&metric, &observer, direction, + &state) == 0, + "initialize past ray"); + /* gamma is diagonal at (30, 0, 0): gamma_xx = 1 + 2/r. */ + double gamma_inv[3][3]; + for (int i = 0; i < 3; ++i) + for (int j = 0; j < 3; ++j) + gamma_inv[i][j] = (i == j) ? 1.0 / metric.gamma[i][j] : 0.0; + /* Null constraint gamma^{ij} Pi_i Pi_j = 1. */ + double null_residual = 0.0; + for (int i = 0; i < 3; ++i) + for (int j = 0; j < 3; ++j) + null_residual += gamma_inv[i][j] * state.Pi[i] * state.Pi[j]; + CHECK(fabs(null_residual - 1.0) < 1e-10, "null constraint preserved"); + /* Observed energy -g(k, e0) = 1 for the unit observer four-velocity. The + * photon four-momentum is reconstructed from the stored state: + * p^0 = exp(L)/alpha and p^i = alpha p^0 gamma^{ij} Pi_j - beta^i p^0. */ + const double k0 = exp(state.log_alpha_p0) / metric.alpha; + double k[4] = {k0, 0.0, 0.0, 0.0}; + for (int i = 0; i < 3; ++i) { + double covariant = 0.0; + for (int j = 0; j < 3; ++j) + covariant += gamma_inv[i][j] * state.Pi[j]; + k[i + 1] = metric.alpha * k0 * covariant - metric.beta[i] * k0; + } + const double *e0 = observer.tetrad[0]; + double inner = -metric.alpha * metric.alpha * k[0] * e0[0]; + for (int i = 0; i < 3; ++i) + for (int j = 0; j < 3; ++j) + inner += metric.gamma[i][j] * (k[i + 1] + metric.beta[i] * k[0]) * + (e0[j + 1] + metric.beta[j] * e0[0]); + CHECK(fabs(inner + 1.0) < 1e-10, "camera energy normalized to one"); + spacetime_destroy(&source); +} + +static void test_threshold_proxies(void) { + SpacetimeSource source = {0}; + CHECK(spacetime_create_schwarzschild_ks(&source, 1.0, 256.0) == 0, + "create schwarzschild for proxies"); + const double r0 = 30.0; + MetricData metric; + CHECK(spacetime_eval(&source, 0.0, (double[]){r0, 0.0, 0.0}, &metric) == + SPACETIME_POINT_OK, + "metric for proxies"); + ObserverState observer; + radial_static_observer(&metric, r0, &observer); + const double direction[3] = {1.0, 0.0, 0.0}; + GeodesicRayState state; + CHECK(geodesic_initialize_past_ray_metric(&metric, &observer, direction, + &state) == 0, + "initialize proxy ray"); + /* L = ln(alpha p^0) is stored; ln(p^0) = L - ln(alpha). Recompute p^0 from + * the tetrad and direction independently. */ + const double k0 = observer.tetrad[0][0] - direction[0] * observer.tetrad[1][0] - + direction[1] * observer.tetrad[2][0] - + direction[2] * observer.tetrad[3][0]; + const double log_p0 = state.log_alpha_p0 - log(metric.alpha); + CHECK(fabs(log_p0 - log(k0)) < 1e-12, + "ln(p^0) = L - ln(alpha) with L = ln(alpha p^0)"); + /* For a static observer far outside, alpha -> 1 and the two proxies agree + * to O(M/r); this documents why a fixed L threshold is not a fixed p^0 + * threshold. */ + CHECK(fabs(state.log_alpha_p0 - log_p0) > 1e-3, + "L and ln(p^0) differ near the hole"); + spacetime_destroy(&source); +} + +int main(void) { + test_radial_branches(); + test_killing_energy_reference(); + test_critical_parameters(); + test_observer_normalization(); + test_threshold_proxies(); + if (failures != 0) { + fprintf(stderr, "termination oracle: %d failure(s)\n", failures); + return 1; + } + return 0; +} diff --git a/usage.md b/usage.md index 76ad529..77bec4e 100644 --- a/usage.md +++ b/usage.md @@ -54,12 +54,24 @@ Single-frame camera options cannot be combined with `--observer-track`, `--frames-dir`, or `--lens-map-input`. Schwarzschild uses Cartesian ingoing Kerr–Schild coordinates with `M=1`. -Cameras at and inside the horizon `r=2` are allowed with a valid timelike -coordinate velocity. The current backend excludes camera positions at or -inside its capture cutoff `r=1.5`; its finite escape radius is `256`. -These remain analytic demonstration settings, not criteria for future NR data. +Cameras at and inside the horizon `r=2` (and inside the old `r=1.5` guard) are +allowed with a valid timelike coordinate velocity; position never decides a ray +endpoint. Its finite escape radius is `256`. These remain analytic demonstration +settings, not criteria for future NR data. Zero coordinate velocity at or inside the horizon is not timelike and is rejected. +The normal dark terminal, for every backend, is the camera-relative local energy +growth `L - L0 >= T` (default `T = 8`, overridable with `--dark-threshold`), +where `L = ln(alpha p^0)` and `L0` is the photon's `L` at the **camera event** +(kept distinct from the escape-worldtube entry energy for an external camera). +A constant camera boost cancels, so a large initial `L` alone does not produce a +dark ray. Neither the photon energy nor the frequency ratio is reset; +budget-exhausted and data/integration failures are separate +unresolved/incomplete outcomes. +Failed and unresolved midpoint probes are retained as diagnostic samples, not +discarded after refinement. Lens-map replay uses its saved geometric policy and +the same incomplete-output check as live tracing. + The following complete examples use the bundled synthetic catalog: ```sh @@ -108,11 +120,12 @@ The bubble therefore propagates through the coordinates, and the metric is time-dependent: the renderer evaluates `f(r_s)` and its spatial derivatives at each coordinate time, while the extrinsic curvature supplies the required `d_t gamma` information to the 3+1 null-ray equations. The exotic matter that -would source the bubble is treated as optically transparent, so there is -**no capture**: rays are only active or escaped. This is why the backend -requires a sub-luminal `|v_s| < 1`; at or above `1` the metric develops an -ergoregion/event horizon and static observers cease to exist, which is outside -the current no-capture scope. +would source the bubble is treated as optically transparent and there is no +horizon, so in practice rays are active or escaped: the bubble has no causal +boundary at which `L - L0` can diverge, and the shared dark policy is not +expected to trigger. This is why the backend requires a sub-luminal +`|v_s| < 1`; at or above `1` the metric develops an ergoregion/event horizon and +static observers cease to exist, which is outside the current scope. | Option | Meaning / default | | --- | --- | @@ -257,9 +270,13 @@ source-sky/lens-map length. A locally escaped triangle is split only when `e / max(s, 1e-15) > --refine-angle-rel`. `P` and `A` prevent selecting a leaf already at or below the requested image-plane long-edge and area scales. -Triangles whose three vertices disagree between capture and escape are split -independently of the direction-error thresholds, allowing the mesh to follow a -shadow boundary. +Triangles whose three vertices straddle a dark/escape or unresolved/dark +boundary are split independently of the direction-error thresholds, allowing the +mesh to follow a shadow boundary. Unresolved vertices with an escape vertex (or +three unresolved vertices) are retried with more step budget before any split; +at the configured total cap the render is reported incomplete unless +`--allow-incomplete` is given. A UUD/UDD boundary triangle at the geometric stop +scale is approximately blackened and recorded with its image-plane area. Independently of the midpoint geometry test, an all-escaped triangle also computes the discrete lens Jacobian