diff --git a/Makefile b/Makefile index 6097a01..514efd1 100644 --- a/Makefile +++ b/Makefile @@ -29,7 +29,7 @@ else IMAGE_EXT := ppm endif -COMMON_SOURCES := $(filter-out src/main.c src/spacetime_minkowski.c src/spacetime_schwarzschild.c,$(wildcard src/*.c)) +COMMON_SOURCES := $(filter-out src/main.c src/dummy_psf.c src/spacetime_minkowski.c src/spacetime_schwarzschild.c,$(wildcard src/*.c)) PROVIDER_SOURCE := src/spacetime_$(SPACETIME).c BUILD_DIR := build/$(BUILD_TYPE) TARGET_BASENAME := $(SPACETIME)_sky @@ -56,8 +56,11 @@ RENDER_LINKER := $(CC) else ifeq ($(PSF_BACKEND),hip) RENDER_LINKER := $(HIPCC) BUILD_CPPFLAGS += -DPSF_BACKEND_HIP +else ifeq ($(PSF_BACKEND),dummy) +RENDER_LINKER := $(CC) +BUILD_CPPFLAGS += -DPSF_BACKEND_DUMMY else -$(error Unknown PSF_BACKEND '$(PSF_BACKEND)'; choose cpu or hip) +$(error Unknown PSF_BACKEND '$(PSF_BACKEND)'; choose cpu, hip, or dummy) endif ifeq ($(SPACETIME),minkowski) @@ -80,10 +83,15 @@ endif ifeq ($(PSF_BACKEND),hip) TARGET := $(BUILD_DIR)/$(TARGET_BASENAME)_hip +else ifeq ($(PSF_BACKEND),dummy) +TARGET := $(BUILD_DIR)/$(TARGET_BASENAME)_dummy else TARGET := $(BUILD_DIR)/$(TARGET_BASENAME) endif RENDER_SOURCES := $(COMMON_SOURCES) $(PROVIDER_SOURCE) src/main.c +ifeq ($(PSF_BACKEND),dummy) +RENDER_SOURCES += src/dummy_psf.c +endif RENDER_OBJECTS := $(patsubst %.c,$(OBJECT_DIR)/$(HDR_BUILD_TAG)/%.o,$(RENDER_SOURCES)) ifeq ($(PSF_BACKEND),hip) HIP_PSF_OBJECT := $(OBJECT_DIR)/$(HDR_BUILD_TAG)/src/hip_psf.o diff --git a/benchmarks/dummy_psf_chunks_2026-09-11/.gitignore b/benchmarks/dummy_psf_chunks_2026-09-11/.gitignore new file mode 100644 index 0000000..c18dd8d --- /dev/null +++ b/benchmarks/dummy_psf_chunks_2026-09-11/.gitignore @@ -0,0 +1 @@ +__pycache__/ diff --git a/benchmarks/dummy_psf_chunks_2026-09-11/run.py b/benchmarks/dummy_psf_chunks_2026-09-11/run.py new file mode 100644 index 0000000..f9a5e82 --- /dev/null +++ b/benchmarks/dummy_psf_chunks_2026-09-11/run.py @@ -0,0 +1,99 @@ +#!/usr/bin/env python3 +"""Run the fixed full-catalog dummy diagnostic once without touching images.""" + +import fcntl +import hashlib +import os +from pathlib import Path +import resource +import shlex +import signal +import subprocess +import time + +ROOT = Path(__file__).resolve().parents[2] +LOG = Path(__file__).with_name("full_catalog.log") +BINARY = ROOT / "build/Release/schwarzschild_sky_dummy" +LENS = ROOT / "output/lens/schwarzschild_galactic_center_R100_45deg_16_4_4k_.grlens" +PNG = ROOT / "output/imgs/2mass_galactic_center_schwarzschild_R100_4k_1e13_beta4.5_16_4-full-gpu.png" +FITS = ROOT / "output/imgs/2mass_galactic_center_schwarzschild_R100_4k_1e13_beta4.5_16_4-full-gpu_HDR.fits" + + +def digest(path): + return hashlib.sha256(path.read_bytes()).hexdigest() + + +def main(): + if LOG.exists(): + raise RuntimeError(f"refusing to overwrite {LOG}") + for path in (BINARY, LENS, PNG, FITS): + if not path.is_file(): + raise RuntimeError(f"required file is missing: {path}") + protected_before = { + path: (digest(path), path.stat().st_size, path.stat().st_mtime_ns) + for path in (PNG, FITS) + } + command = [ + str(BINARY), "--verbose", + "--psf-relative-tail", "1e-8", "--psf-min-y", "0", + "--max-cache-psf-flux", "1e8", + "--all-sky-catalog", "assets/2mass/processed/all_sky", + "--exposure", "1e13", "--psf-fwhm-pixels", "2.7", + "--psf-moffat-beta", "4.5", "--catalog-load-workers", "4", + "--hdr-output", "--lens-map-input", + "output/lens/schwarzschild_galactic_center_R100_45deg_16_4_4k_.grlens", + "--output", + "output/imgs/2mass_galactic_center_schwarzschild_R100_4k_1e13_beta4.5_16_4-full-gpu.png", + ] + guard = open("/tmp/gr_dummy_psf_chunks.lock", "w") + fcntl.flock(guard, fcntl.LOCK_EX | fcntl.LOCK_NB) + env = dict(os.environ, OMP_NUM_THREADS="16", OMP_DYNAMIC="FALSE") + with LOG.open("x") as out: + out.write("COMMAND OMP_NUM_THREADS=16 OMP_DYNAMIC=FALSE " + + shlex.join(command) + "\n") + out.write("GIT " + subprocess.check_output( + ["git", "rev-parse", "HEAD"], cwd=ROOT, text=True)) + out.write("BINARY_SHA256 " + digest(BINARY) + "\n") + out.write("LENS_SHA256 " + digest(LENS) + "\n") + for path, state in protected_before.items(): + out.write( + f"PROTECTED_BEFORE {path.relative_to(ROOT)} sha256={state[0]} " + f"size={state[1]} mtime_ns={state[2]}\n") + out.flush() + start = time.monotonic() + usage_start = resource.getrusage(resource.RUSAGE_CHILDREN) + process = subprocess.Popen(command, cwd=ROOT, env=env, stdout=out, + stderr=subprocess.STDOUT, + start_new_session=True) + try: + returncode = process.wait(timeout=600) + except subprocess.TimeoutExpired: + os.killpg(process.pid, signal.SIGKILL) + process.wait() + out.write("TIMEOUT killed and reaped process group\n") + raise + out.write( + f"PROCESS_WALL {time.monotonic() - start:.9f} EXIT {returncode}\n") + usage_end = resource.getrusage(resource.RUSAGE_CHILDREN) + out.write( + f"PROCESS_USER {usage_end.ru_utime - usage_start.ru_utime:.9f} " + f"PROCESS_SYSTEM {usage_end.ru_stime - usage_start.ru_stime:.9f} " + f"MAX_RSS_KIB {usage_end.ru_maxrss}\n") + protected_after = { + path: (digest(path), path.stat().st_size, path.stat().st_mtime_ns) + for path in (PNG, FITS) + } + if protected_after != protected_before: + raise RuntimeError("a protected render output changed") + with LOG.open("a") as out: + for path, state in protected_after.items(): + out.write( + f"PROTECTED_AFTER {path.relative_to(ROOT)} sha256={state[0]} " + f"size={state[1]} mtime_ns={state[2]} unchanged=yes\n") + if returncode: + raise RuntimeError(f"dummy diagnostic exited {returncode}") + print(LOG) + + +if __name__ == "__main__": + main() diff --git a/build.md b/build.md index dfbac17..dc55da9 100644 --- a/build.md +++ b/build.md @@ -104,6 +104,29 @@ times must not be added to GPU timings or interpreted as process CPU time. The two upload slots bound staging memory and allow CPU production ahead of GPU completion; copies and kernels still execute in order on one stream. +### Diagnostic dummy PSF backend + +To inspect production chunk geometry without starting HIP or allocating an HDR +framebuffer, build the diagnostic-only dummy backend: + +```sh +make -j1 PSF_BACKEND=dummy SPACETIME=schwarzschild ENABLE_HDR=1 backend +OMP_NUM_THREADS=16 OMP_DYNAMIC=FALSE \ + ./build/Release/schwarzschild_sky_dummy --lens-map-input MAP.grlens \ + --all-sky-catalog CATALOG_DIR --output PROTECTED_OR_UNUSED.png [normal PSF options] +``` + +The dummy executable requires an imported lens map and the PSF cache; it rejects +`--psf-direct`. It runs the real OpenMP triangle scheduling, catalog query, +inverse mapping, redshift/color/flux preparation, and cache/direct/min-Y +classification, but records chunk statistics instead of accumulating pixels. +`--output` and `--hdr-output` are accepted only for command parity: neither PNG +nor FITS is created or modified. Reports separate full 16,384-event chunks from +the final partial chunk of each worker and include event-producing triangle +counts, occupied 32-pixel center tiles, bounding-box diagonals, and successive +triangle-centroid jumps. This is a scheduling diagnostic, not a rendering or +backend-performance benchmark. + With both Minkowski binaries built using `ENABLE_HDR=1`, `python3 tests/test_hip_renderer.py build/Release/minkowski_sky build/Release/minkowski_sky_hip` runs small CPU/HIP fallback and min-Y comparisons sequentially. Each renderer diff --git a/src/dummy_psf.c b/src/dummy_psf.c new file mode 100644 index 0000000..25e58d5 --- /dev/null +++ b/src/dummy_psf.c @@ -0,0 +1,309 @@ +#include "dummy_psf.h" + +#include +#include +#include +#include +#include +#include + +enum { DUMMY_SELECTOR_TILE = 32 }; + +typedef struct { + size_t worker, events, triangles, occupied_tiles; + double event_bbox_diagonal, triangle_bbox_diagonal; + double maximum_triangle_jump; + int adaptive_tile16; +} DummyChunkSample; + +struct DummyPsfSink { + DummyChunkSample *samples; + size_t count, capacity, event_capacity; + size_t total_events; + int width, height, tiles_x, tiles_y; + int failed; + omp_lock_t lock; +}; + +struct DummyPsfChunk { + DummyPsfSink *sink; + uint32_t *tile_generation; + uint32_t generation; + size_t worker, events, triangles, occupied_tiles; + size_t last_triangle; + int have_event, have_triangle; + double event_min_x, event_max_x, event_min_y, event_max_y; + double triangle_min_x, triangle_max_x, triangle_min_y, triangle_max_y; + double last_triangle_x, last_triangle_y, maximum_triangle_jump; +}; + +static int append_sample(DummyPsfSink *sink, DummyChunkSample sample) { + int result = 0; + omp_set_lock(&sink->lock); + if (sink->failed) { + result = -1; + } else if (sink->count == sink->capacity) { + const size_t capacity = sink->capacity ? 2 * sink->capacity : 1024; + DummyChunkSample *samples = + capacity < sink->capacity || capacity > SIZE_MAX / sizeof *samples + ? NULL + : realloc(sink->samples, + capacity * sizeof *samples); + if (samples == NULL) { + sink->failed = 1; + result = -1; + } else { + sink->samples = samples; + sink->capacity = capacity; + } + } + if (!result) { + sink->samples[sink->count++] = sample; + sink->total_events += sample.events; + } + omp_unset_lock(&sink->lock); + return result; +} + +int dummy_psf_sink_create(DummyPsfSink **out, int width, int height, + size_t event_capacity) { + if (out == NULL || width <= 0 || height <= 0 || event_capacity == 0) + return -1; + DummyPsfSink *sink = calloc(1, sizeof *sink); + if (sink == NULL) + return -1; + sink->width = width; + sink->height = height; + sink->tiles_x = (width + DUMMY_SELECTOR_TILE - 1) / DUMMY_SELECTOR_TILE; + sink->tiles_y = (height + DUMMY_SELECTOR_TILE - 1) / DUMMY_SELECTOR_TILE; + sink->event_capacity = event_capacity; + omp_init_lock(&sink->lock); + *out = sink; + return 0; +} + +int dummy_psf_chunk_create(DummyPsfChunk **out, DummyPsfSink *sink, + size_t worker_id) { + if (out == NULL || sink == NULL) + return -1; + DummyPsfChunk *chunk = calloc(1, sizeof *chunk); + if (chunk == NULL) + return -1; + const size_t tile_count = (size_t)sink->tiles_x * sink->tiles_y; + chunk->tile_generation = calloc(tile_count, sizeof *chunk->tile_generation); + if (chunk->tile_generation == NULL) { + free(chunk); + return -1; + } + chunk->sink = sink; + chunk->worker = worker_id; + chunk->generation = 1; + *out = chunk; + return 0; +} + +static void reset_chunk(DummyPsfChunk *chunk) { + chunk->events = 0; + chunk->triangles = 0; + chunk->occupied_tiles = 0; + chunk->have_event = 0; + chunk->have_triangle = 0; + chunk->maximum_triangle_jump = 0.0; + if (++chunk->generation == 0) { + const size_t count = (size_t)chunk->sink->tiles_x * chunk->sink->tiles_y; + memset(chunk->tile_generation, 0, count * sizeof *chunk->tile_generation); + chunk->generation = 1; + } +} + +int dummy_psf_chunk_flush(DummyPsfChunk *chunk) { + if (chunk == NULL) + return -1; + if (chunk->events == 0) + return 0; + const double event_diagonal = + hypot(chunk->event_max_x - chunk->event_min_x, + chunk->event_max_y - chunk->event_min_y); + const double triangle_diagonal = + hypot(chunk->triangle_max_x - chunk->triangle_min_x, + chunk->triangle_max_y - chunk->triangle_min_y); + const DummyChunkSample sample = { + .worker = chunk->worker, + .events = chunk->events, + .triangles = chunk->triangles, + .occupied_tiles = chunk->occupied_tiles, + .event_bbox_diagonal = event_diagonal, + .triangle_bbox_diagonal = triangle_diagonal, + .maximum_triangle_jump = chunk->maximum_triangle_jump, + .adaptive_tile16 = chunk->events >= 8192 && chunk->occupied_tiles != 0 && + chunk->events / chunk->occupied_tiles >= 32}; + const int result = append_sample(chunk->sink, sample); + reset_chunk(chunk); + return result; +} + +int dummy_psf_chunk_emit(DummyPsfChunk *chunk, const PsfCachedEvent *event, + size_t triangle, double triangle_center_x, + double triangle_center_y) { + if (chunk == NULL || event == NULL) + return -1; + if (chunk->events == chunk->sink->event_capacity && + dummy_psf_chunk_flush(chunk)) + return -1; + + if (!chunk->have_event) { + chunk->event_min_x = chunk->event_max_x = event->x; + chunk->event_min_y = chunk->event_max_y = event->y; + chunk->have_event = 1; + } else { + chunk->event_min_x = fmin(chunk->event_min_x, event->x); + chunk->event_max_x = fmax(chunk->event_max_x, event->x); + chunk->event_min_y = fmin(chunk->event_min_y, event->y); + chunk->event_max_y = fmax(chunk->event_max_y, event->y); + } + + const int pixel_x = (int)floor(event->x); + const int pixel_y = (int)floor(event->y); + if (pixel_x >= 0 && pixel_x < chunk->sink->width && pixel_y >= 0 && + pixel_y < chunk->sink->height) { + const size_t tile = (size_t)(pixel_y / DUMMY_SELECTOR_TILE) * + chunk->sink->tiles_x + + pixel_x / DUMMY_SELECTOR_TILE; + if (chunk->tile_generation[tile] != chunk->generation) { + chunk->tile_generation[tile] = chunk->generation; + ++chunk->occupied_tiles; + } + } + + if (!chunk->have_triangle || triangle != chunk->last_triangle) { + if (!chunk->have_triangle) { + chunk->triangle_min_x = chunk->triangle_max_x = triangle_center_x; + chunk->triangle_min_y = chunk->triangle_max_y = triangle_center_y; + chunk->have_triangle = 1; + } else { + chunk->triangle_min_x = fmin(chunk->triangle_min_x, triangle_center_x); + chunk->triangle_max_x = fmax(chunk->triangle_max_x, triangle_center_x); + chunk->triangle_min_y = fmin(chunk->triangle_min_y, triangle_center_y); + chunk->triangle_max_y = fmax(chunk->triangle_max_y, triangle_center_y); + chunk->maximum_triangle_jump = + fmax(chunk->maximum_triangle_jump, + hypot(triangle_center_x - chunk->last_triangle_x, + triangle_center_y - chunk->last_triangle_y)); + } + chunk->last_triangle = triangle; + chunk->last_triangle_x = triangle_center_x; + chunk->last_triangle_y = triangle_center_y; + ++chunk->triangles; + } + ++chunk->events; + return 0; +} + +void dummy_psf_chunk_destroy(DummyPsfChunk *chunk) { + if (chunk == NULL) + return; + (void)dummy_psf_chunk_flush(chunk); + free(chunk->tile_generation); + free(chunk); +} + +static int compare_double(const void *left, const void *right) { + const double a = *(const double *)left; + const double b = *(const double *)right; + return (a > b) - (a < b); +} + +static double percentile(const double *values, size_t count, double fraction) { + if (count == 0) + return 0.0; + const double position = fraction * (count - 1); + const size_t low = (size_t)floor(position); + const size_t high = low + 1 < count ? low + 1 : low; + return values[low] + (values[high] - values[low]) * (position - low); +} + +static void report_group(const DummyPsfSink *sink, int full_only, + const char *label) { + size_t count = 0, adaptive = 0; + for (size_t i = 0; i < sink->count; ++i) + if (!full_only || sink->samples[i].events == sink->event_capacity) + ++count; + if (count == 0) { + fprintf(stderr, "Dummy PSF %s chunks: none\n", label); + return; + } + double *events = count > SIZE_MAX / (6 * sizeof *events) + ? NULL + : malloc(6 * count * sizeof *events); + if (events == NULL) { + fputs("Dummy PSF report allocation failed\n", stderr); + return; + } + double *triangles = events + count; + double *tiles = triangles + count; + double *event_distance = tiles + count; + double *triangle_distance = event_distance + count; + double *maximum_jump = triangle_distance + count; + size_t next = 0; + for (size_t i = 0; i < sink->count; ++i) { + const DummyChunkSample *sample = &sink->samples[i]; + if (full_only && sample->events != sink->event_capacity) + continue; + events[next] = sample->events; + triangles[next] = sample->triangles; + tiles[next] = sample->occupied_tiles; + event_distance[next] = sample->event_bbox_diagonal; + triangle_distance[next] = sample->triangle_bbox_diagonal; + maximum_jump[next] = sample->maximum_triangle_jump; + adaptive += sample->adaptive_tile16; + ++next; + } + double *series[] = {events, triangles, tiles, event_distance, + triangle_distance, maximum_jump}; + const char *names[] = {"events", "distinct_triangles", "occupied_32px_tiles", + "event_bbox_diagonal_px", + "triangle_centroid_bbox_diagonal_px", + "maximum_successive_triangle_jump_px"}; + for (size_t field = 0; field < 6; ++field) { + qsort(series[field], count, sizeof **series, compare_double); + fprintf(stderr, + "Dummy PSF %s %s: min=%.3f p10=%.3f p25=%.3f p50=%.3f " + "p75=%.3f p90=%.3f p99=%.3f max=%.3f\n", + label, names[field], series[field][0], + percentile(series[field], count, 0.10), + percentile(series[field], count, 0.25), + percentile(series[field], count, 0.50), + percentile(series[field], count, 0.75), + percentile(series[field], count, 0.90), + percentile(series[field], count, 0.99), series[field][count - 1]); + } + fprintf(stderr, + "Dummy PSF %s adaptive criterion: %zu/%zu chunks tile16-eligible " + "(events>=8192 and events/occupied_32px_tiles>=32)\n", + label, adaptive, count); + free(events); +} + +void dummy_psf_sink_report(const DummyPsfSink *sink) { + if (sink == NULL) + return; + size_t full = 0; + for (size_t i = 0; i < sink->count; ++i) + full += sink->samples[i].events == sink->event_capacity; + fprintf(stderr, + "Dummy PSF diagnostic only: no HDR or image was accumulated/written.\n" + "Dummy PSF chunks: total=%zu full=%zu partial=%zu events=%zu " + "capacity=%zu selector_tile=%dpx\n", + sink->count, full, sink->count - full, sink->total_events, + sink->event_capacity, DUMMY_SELECTOR_TILE); + report_group(sink, 0, "all"); + report_group(sink, 1, "full"); +} + +void dummy_psf_sink_destroy(DummyPsfSink *sink) { + if (sink == NULL) + return; + omp_destroy_lock(&sink->lock); + free(sink->samples); + free(sink); +} diff --git a/src/dummy_psf.h b/src/dummy_psf.h new file mode 100644 index 0000000..8c16729 --- /dev/null +++ b/src/dummy_psf.h @@ -0,0 +1,25 @@ +#ifndef DUMMY_PSF_H +#define DUMMY_PSF_H + +#include "optics.h" + +#include + +typedef struct DummyPsfSink DummyPsfSink; +typedef struct DummyPsfChunk DummyPsfChunk; + +/* Diagnostic-only consumer for the production PSF event stream. It preserves + * per-worker chunk boundaries but never accumulates or writes an HDR image. */ +int dummy_psf_sink_create(DummyPsfSink **out, int width, int height, + size_t event_capacity); +int dummy_psf_chunk_create(DummyPsfChunk **out, DummyPsfSink *sink, + size_t worker_id); +int dummy_psf_chunk_emit(DummyPsfChunk *chunk, const PsfCachedEvent *event, + size_t triangle, double triangle_center_x, + double triangle_center_y); +int dummy_psf_chunk_flush(DummyPsfChunk *chunk); +void dummy_psf_chunk_destroy(DummyPsfChunk *chunk); +void dummy_psf_sink_report(const DummyPsfSink *sink); +void dummy_psf_sink_destroy(DummyPsfSink *sink); + +#endif diff --git a/src/frame.c b/src/frame.c index 38d1857..9f33894 100644 --- a/src/frame.c +++ b/src/frame.c @@ -3,6 +3,9 @@ #ifdef PSF_BACKEND_HIP #include "hip_psf.h" #endif +#ifdef PSF_BACKEND_DUMMY +#include "dummy_psf.h" +#endif #include "optics.h" @@ -138,6 +141,11 @@ typedef struct { double submit_seconds; char hip_message[256]; HipPsfTiming hip_timing; +#endif +#ifdef PSF_BACKEND_DUMMY + DummyPsfSink *dummy; + DummyPsfChunk *dummy_chunk; + int borrowed_dummy; #endif int failed; } PsfEventSink; @@ -155,6 +163,16 @@ static int psf_event_sink_init(PsfEventSink *sink, double *hdr, int width, *sink = (PsfEventSink){.hdr = hdr, .width = width, .height = height, .cache = cache}; if (cache != NULL) { +#ifdef PSF_BACKEND_DUMMY + if (dummy_psf_sink_create(&sink->dummy, width, height, + PSF_EVENT_SINK_CAPACITY)) { + fputs("Dummy PSF backend initialization failed\n", stderr); + dummy_psf_sink_destroy(sink->dummy); + sink->dummy = NULL; + sink->failed = 1; + return -1; + } +#else sink->events = malloc(PSF_EVENT_SINK_CAPACITY * sizeof *sink->events); if (sink->events == NULL) { #ifdef PSF_BACKEND_HIP @@ -174,6 +192,7 @@ static int psf_event_sink_init(PsfEventSink *sink, double *hdr, int width, sink->events = NULL; return -1; } +#endif #endif } return 0; @@ -182,6 +201,10 @@ static int psf_event_sink_init(PsfEventSink *sink, double *hdr, int width, static int psf_event_sink_flush_unlocked(PsfEventSink *sink) { if (sink->failed) return -1; +#ifdef PSF_BACKEND_DUMMY + if (sink->dummy_chunk != NULL) + return dummy_psf_chunk_flush(sink->dummy_chunk); +#endif #ifdef PSF_BACKEND_HIP if (sink->hip != NULL) { if (hip_psf_sink_submit(sink->hip, sink->events, sink->count, @@ -233,6 +256,7 @@ static int psf_event_sink_finish_for_cpu(PsfEventSink *sink) { return 0; } +#ifndef PSF_BACKEND_DUMMY static int psf_event_sink_resume_gpu(PsfEventSink *sink) { #ifdef PSF_BACKEND_HIP if (sink->hip != NULL && hip_psf_sink_load_hdr(sink->hip, sink->hdr, @@ -246,8 +270,22 @@ static int psf_event_sink_resume_gpu(PsfEventSink *sink) { #endif return 0; } +#endif static int psf_event_sink_destroy(PsfEventSink *sink) { +#ifdef PSF_BACKEND_DUMMY + int dummy_result = sink->dummy_chunk != NULL + ? dummy_psf_chunk_flush(sink->dummy_chunk) + : 0; + dummy_psf_chunk_destroy(sink->dummy_chunk); + sink->dummy_chunk = NULL; + if (!sink->borrowed_dummy) { + dummy_psf_sink_report(sink->dummy); + dummy_psf_sink_destroy(sink->dummy); + sink->dummy = NULL; + } + return dummy_result; +#endif #ifdef PSF_BACKEND_HIP if (sink->borrowed_hip) { const int result = psf_event_sink_flush(sink); @@ -267,10 +305,22 @@ static int psf_event_sink_destroy(PsfEventSink *sink) { return result; } -#if FRAME_PSF_EVENT_SINK || defined(PSF_BACKEND_HIP) -static void psf_event_sink_emit(PsfEventSink *sink, const PsfCachedEvent *event) { +#if FRAME_PSF_EVENT_SINK || defined(PSF_BACKEND_HIP) || defined(PSF_BACKEND_DUMMY) +static void psf_event_sink_emit(PsfEventSink *sink, const PsfCachedEvent *event, + size_t triangle, double triangle_center_x, + double triangle_center_y) { if (sink->failed) return; +#ifdef PSF_BACKEND_DUMMY + if (dummy_psf_chunk_emit(sink->dummy_chunk, event, triangle, + triangle_center_x, triangle_center_y)) + sink->failed = 1; + return; +#else + (void)triangle; + (void)triangle_center_x; + (void)triangle_center_y; +#endif if (sink->events == NULL) { splat_prepared_cached_event(sink->hdr, sink->width, sink->height, event, sink->cache); return; @@ -1129,6 +1179,9 @@ static int splat_catalog_tile(const Star *stars, size_t count, #endif if (direct_fallback == 1) { /* Exclude every other submit and fallback until HDR is back on device. */ +#ifdef PSF_BACKEND_DUMMY + if (psf_event_sink_flush(context->event_sink)) return -1; +#else #ifdef PSF_BACKEND_HIP const double submit_start = omp_get_wtime(); if (context->event_sink->hip_lock) @@ -1147,9 +1200,17 @@ static int splat_catalog_tile(const Star *stars, size_t count, context->event_sink->submit_seconds += omp_get_wtime() - submit_start; #endif if (failed) return -1; +#endif } else if (direct_fallback != 3) { -#if FRAME_PSF_EVENT_SINK || defined(PSF_BACKEND_HIP) - psf_event_sink_emit(context->event_sink, &event); +#if FRAME_PSF_EVENT_SINK || defined(PSF_BACKEND_HIP) || defined(PSF_BACKEND_DUMMY) + const double triangle_center_x = + (context->vertex[0]->image_x + context->vertex[1]->image_x + + context->vertex[2]->image_x) / 3.0; + const double triangle_center_y = + (context->vertex[0]->image_y + context->vertex[1]->image_y + + context->vertex[2]->image_y) / 3.0; + psf_event_sink_emit(context->event_sink, &event, context->triangle_index, + triangle_center_x, triangle_center_y); #else splat_prepared_cached_event(context->hdr, context->width, context->height, &event, context->psf_cache); @@ -1319,6 +1380,83 @@ size_t frame_splat_catalog(const FrameLensMesh *mesh, progress->callback(progress->context, FRAME_SPLAT_PROGRESS_BEGIN, 0, mesh->triangle_count); +#ifdef PSF_BACKEND_DUMMY + PsfEventSink owner; + if (psf_event_sink_init(&owner, hdr, width, height, psf_cache)) return SIZE_MAX; + size_t dummy_images = 0, dummy_direct = 0, dummy_clipped = 0, + dummy_discarded = 0; + int dummy_failed = 0, dummy_workers = 0; + const double dummy_start = omp_get_wtime(); +#ifdef GR_DEBUG + double dummy_max_magnification = 0.0; + size_t dummy_clamped = 0; +#pragma omp parallel reduction(+ : dummy_images, dummy_direct, dummy_clipped, dummy_discarded, dummy_clamped) reduction(max : dummy_failed, dummy_max_magnification) +#else +#pragma omp parallel reduction(+ : dummy_images, dummy_direct, dummy_clipped, dummy_discarded) reduction(max : dummy_failed) +#endif + { +#pragma omp single + dummy_workers = omp_get_num_threads(); + const size_t worker = (size_t)omp_get_thread_num(); + PsfEventSink local = {.hdr = hdr, .width = width, .height = height, + .cache = psf_cache, .dummy = owner.dummy, .borrowed_dummy = 1}; + if (dummy_psf_chunk_create(&local.dummy_chunk, owner.dummy, worker)) { + fputs("Dummy PSF worker chunk initialization failed\n", stderr); + local.failed = 1; + } + size_t completed = 0, next_report = 8; + if (progress && progress->worker_callback) + progress->worker_callback(progress->context, worker, + (size_t)dummy_workers, 0, 0); +#pragma omp for schedule(dynamic, 1) nowait + for (size_t t = 0; t < mesh->triangle_count; ++t) { + if (local.failed) continue; + const CatalogSplatStats stats = splat_catalog_triangles( + mesh, catalog, hdr, width, height, exposure, psf, psf_cache, + max_magnification, max_cache_psf_flux, psf_relative_tail, psf_min_y, + t, t + 1, &local); + dummy_images += stats.images; + dummy_direct += stats.direct_fallbacks; + dummy_clipped += stats.cached_wing_clipped; + dummy_discarded += stats.discarded_below_min_y; + dummy_failed |= stats.failed; +#ifdef GR_DEBUG + dummy_max_magnification = fmax(dummy_max_magnification, + stats.max_raw_magnification); + dummy_clamped += stats.magnification_clamped_triangles; +#endif + ++completed; + if (progress && progress->worker_callback && completed == next_report) { + progress->worker_callback(progress->context, worker, + (size_t)dummy_workers, completed, 0); + if (next_report <= SIZE_MAX / 2) next_report *= 2; + } + } + if (psf_event_sink_destroy(&local)) dummy_failed = 1; + if (progress && progress->worker_callback) + progress->worker_callback(progress->context, worker, + (size_t)dummy_workers, completed, 1); + } + if (psf_event_sink_destroy(&owner)) dummy_failed = 1; + fprintf(stderr, + "Dummy PSF producers: %d workers; classification/chunk wall %.3f s\n", + dummy_workers, omp_get_wtime() - dummy_start); + copy_psf_splat_stats(psf_stats, (CatalogSplatStats){ + .images = dummy_images, .direct_fallbacks = dummy_direct, + .cached_wing_clipped = dummy_clipped, + .discarded_below_min_y = dummy_discarded, +#ifdef GR_DEBUG + .max_raw_magnification = dummy_max_magnification, + .magnification_clamped_triangles = dummy_clamped, +#endif + }); + if (dummy_failed) return SIZE_MAX; + if (progress != NULL && progress->callback != NULL) + progress->callback(progress->context, FRAME_SPLAT_PROGRESS_END, + mesh->triangle_count, mesh->triangle_count); + return dummy_images; +#endif + #ifdef PSF_BACKEND_HIP /* CPU workers own only bounded event chunks. A single device cache/HDR is * borrowed under a coarse submission lock, never duplicated per worker. */ diff --git a/src/main.c b/src/main.c index ae8e5b6..87711db 100644 --- a/src/main.c +++ b/src/main.c @@ -445,8 +445,13 @@ static void report_splat_progress(void *context, FrameSplatProgressStage stage, progress->frame_id, total); break; case FRAME_SPLAT_PROGRESS_END: +#ifdef PSF_BACKEND_DUMMY + fprintf(stderr, "Frame %zu: catalog classification finished in %.1f s; reporting chunk statistics...\n", + progress->frame_id, omp_get_wtime() - progress->splat_start); +#else fprintf(stderr, "Frame %zu: catalog splatting finished in %.1f s; writing image...\n", progress->frame_id, omp_get_wtime() - progress->splat_start); +#endif break; } } @@ -922,7 +927,15 @@ static int render_lens_map(const Settings *s, StarCatalog *catalog) { } output_path = movie_path; } +#ifdef PSF_BACKEND_DUMMY + (void)output_path; +#endif +#ifdef PSF_BACKEND_DUMMY + /* The dummy backend classifies real events but never touches pixels. */ + double *hdr = calloc(1, sizeof *hdr); +#else double *hdr = calloc((size_t)map.width * map.height * 3, sizeof *hdr); +#endif if (hdr == NULL) { result = -1; break; } CatalogPrefetchStats prefetch = {0}; PsfSplatStats psf_stats = {0}; @@ -945,7 +958,12 @@ static int render_lens_map(const Settings *s, StarCatalog *catalog) { break; } if (s->draw_mesh) +#ifdef PSF_BACKEND_DUMMY + fputs("Dummy PSF backend ignores --draw-mesh.\n", stderr); +#else frame_draw_mesh(&map.frames[i].mesh, hdr, map.width, map.height, 0.5, 0.5); +#endif +#ifndef PSF_BACKEND_DUMMY #ifdef ENABLE_HDR_OUTPUT if (s->write_hdr_output && write_hdr_fits(s->hdr_output_path, hdr, map.width, map.height, @@ -960,6 +978,13 @@ static int render_lens_map(const Settings *s, StarCatalog *catalog) { images, catalog->count, output_path, write_result == 0 ? "ok" : "write failed"); psf_kernel_cache_report(&s->psf_cache, &psf_stats, stderr); if (write_result) { result = -1; break; } +#else + free(hdr); + fprintf(stderr, + "Dummy PSF classified %zu images from %zu catalog stars; no HDR, PNG, or PPM was written.\n", + images, catalog->count); + psf_kernel_cache_report(&s->psf_cache, &psf_stats, stderr); +#endif } lens_map_destroy(&map); return result; @@ -1026,6 +1051,20 @@ int main(int argc, char **argv) { fputs("--lens-map-input and --lens-map-output are mutually exclusive.\n", stderr); return 2; } +#ifdef PSF_BACKEND_DUMMY + if (settings.lens_map_input_path == NULL) { + fputs("The dummy PSF backend requires --lens-map-input and never traces or writes images.\n", + stderr); + return 2; + } + if (settings.psf_direct) { + fputs("The dummy PSF backend measures cache-event chunks and does not support --psf-direct.\n", + stderr); + return 2; + } + fputs("Dummy PSF backend: --output and --hdr-output are accepted for command parity but no image files will be written.\n", + stderr); +#endif #ifdef ENABLE_HDR_OUTPUT if (settings.frames_dir != NULL && settings.write_hdr_output) { fputs("--hdr-output is available only for a single-frame render.\n", stderr); @@ -1073,8 +1112,17 @@ int main(int argc, char **argv) { fprintf(stderr, "Created test catalog: %s\n", settings.catalog_path); } if (!settings.psf_direct && psf_kernel_cache_init(&settings.psf_cache, &settings.psf, - settings.psf_relative_tail)) + settings.psf_relative_tail)) { +#ifdef PSF_BACKEND_DUMMY + fputs("PSF cache construction failed; dummy chunk statistics are unavailable.\n", + stderr); + catalog_destroy(&catalog); + spacetime_destroy(&spacetime); + return 1; +#else fputs("PSF cache construction failed; using direct evaluator.\n", stderr); +#endif + } psf_kernel_cache_report_ready(&settings.psf_cache, stderr); if (settings.lens_map_input_path != NULL) { const int result = render_lens_map(&settings, &catalog);