Optics: add bounded HIP PSF replay and tile reduction experiments

Capture prepared events through a test-only producer consumer and measure serial and indexed parallel generation independently. Compare the production atomic kernel with bounded pixel-owned tile reductions while preserving double HDR, cache support, and direct fallback ordering.

Validate captured inputs, edge and multichunk fixtures, CPU/HIP renderer regressions, and 45 isolated final replays. Keep the production HIP algorithm unchanged because gains depend on event distribution.

Makefile

@@ -182,3 +182,17 @@ clean:
 
 -include $(RENDER_DEPS)
 include mk/reference_images.mk
+
+# Test-only producer consumer: never linked into a renderer.
+.PHONY: psf-capture
+psf-capture: $(BUILD_DIR)/capture_psf
+$(BUILD_DIR)/capture_psf: tests/capture_psf.c $(CORE_MINKOWSKI_SOURCES) | $(BUILD_DIR)
+	$(CC) $(CPPFLAGS) $(CFLAGS) $(BUILD_CFLAGS) $(OPENMP_FLAGS) -Isrc $< $(filter-out src/frame.c,$(CORE_MINKOWSKI_SOURCES)) $(LDLIBS) -o $@
+
+.PHONY: hip-psf-replay
+hip-psf-replay: $(BUILD_DIR)/replay_psf
+$(BUILD_DIR)/replay_psf: tests/replay_psf.hip src/hip_psf.hip src/hip_psf.h src/optics.h $(OBJECT_DIR)/$(HDR_BUILD_TAG)/src/optics.o | $(BUILD_DIR)
+	$(HIPCC) $(HIP_CXXFLAGS) $(OPENMP_FLAGS) -Isrc -x hip $< -x none $(filter %.o,$^) $(LDLIBS) $(HDR_LDLIBS) -o $@
+
+$(BUILD_DIR)/make_psf_fixture: tests/make_psf_fixture.c src/optics.c src/optics.h | $(BUILD_DIR)
+	$(CC) $(CPPFLAGS) $(CFLAGS) $(BUILD_CFLAGS) $(OPENMP_FLAGS) -Isrc $< src/optics.c $(LDLIBS) -o $@

src/frame.c

@@ -1114,6 +1114,19 @@ static int splat_catalog_tile(const Star *stars, size_t count,
     const int direct_fallback = psf_prepare_cached_event(
         &event, image_x, image_y, color, flux, context->psf, context->psf_cache,
         context->max_cache_psf_flux, context->psf_relative_tail, context->psf_min_y);
+#ifdef FRAME_PSF_DIAGNOSTIC
+    /* Test-only consumer: normal query/mapping/colour/classification above.
+     * Never produces an HDR image; not compiled into renderer binaries. */
+    if (frame_psf_diagnostic_visit(&event, direct_fallback,
+                                  context->triangle_index, image_x, image_y,
+                                  color, flux))
+      return -1;
+    context->direct_fallbacks += direct_fallback == 1;
+    context->cached_wing_clipped += direct_fallback == 2;
+    context->discarded_below_min_y += direct_fallback == 3;
+    ++context->images;
+    continue;
+#endif
     if (direct_fallback == 1) {
       /* Exclude every other submit and fallback until HDR is back on device. */
 #ifdef PSF_BACKEND_HIP
@@ -1180,6 +1193,9 @@ static CatalogSplatStats splat_catalog_triangles(
     event_sink = &owned_sink;
   }
   for (size_t t = first_triangle; t < last_triangle; ++t) {
+#ifdef FRAME_PSF_DIAGNOSTIC
+    if (frame_psf_diagnostic_stopped()) break;
+#endif
     const LensVertex *vertex[3];
     if (!usable_triangle(mesh, &mesh->triangles[t], vertex))
       continue;
@@ -1189,6 +1205,9 @@ static CatalogSplatStats splat_catalog_triangles(
         spherical_area(vertex[0]->camera_direction, vertex[1]->camera_direction,
                        vertex[2]->camera_direction);
     const double raw_magnification = image_area / source_area;
+#ifdef FRAME_PSF_DIAGNOSTIC
+    if (!frame_psf_diagnostic_select(raw_magnification)) continue;
+#endif
     double magnification = raw_magnification;
 #ifdef GR_DEBUG
     stats.max_raw_magnification = fmax(stats.max_raw_magnification,

tests/capture_psf.c

@@ -0,0 +1,151 @@
+/* Standalone bounded diagnostic. Include the production producer so its
+ * private query/mapping path stays identical; renderer builds omit the hook. */
+#include "optics.h"
+#include <stdint.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <sys/stat.h>
+static PsfCachedEvent captured[65536];
+static size_t limit;
+static int diagnosis;
+static _Thread_local size_t captured_count, classified[4], triangles_seen, last_triangle;
+static _Thread_local uint64_t checksum;
+static double selected_min_mag=0, selected_max_mag=1e300;
+static int frame_psf_diagnostic_select(double mag) {
+  return mag>=selected_min_mag && mag<selected_max_mag;
+}
+static _Thread_local uint64_t digest = UINT64_C(14695981039346656037);
+static int frame_psf_diagnostic_stopped(void) {
+  return captured_count >= limit || classified[0]+classified[1]+classified[2]+classified[3] >= 131072;
+}
+static int frame_psf_diagnostic_visit(const PsfCachedEvent *event, int state,
+    size_t triangle, double x, double y, LinearRgb color, double flux) {
+  if (frame_psf_diagnostic_stopped()) return -1;
+  if (state < 0 || state > 3) abort();
+  ++classified[state];
+  if (!triangles_seen || triangle != last_triangle) ++triangles_seen;
+  last_triangle = triangle;
+  /* Hash explicit scalar fields, including rejected/direct events; no padding. */
+  const double fields[] = {x,y,color.r,color.g,color.b,flux,(double)state};
+  const unsigned char *bytes = (const unsigned char *)fields;
+  for (size_t i=0; i<sizeof fields; ++i) digest=(digest ^ bytes[i])*UINT64_C(1099511628211);
+  uint64_t item=UINT64_C(14695981039346656037);
+  for(size_t i=0;i<sizeof fields;++i)item=(item ^ bytes[i])*UINT64_C(1099511628211);
+  checksum+=item;
+  if (state == 0 || state == 2) {
+    if(!diagnosis)captured[captured_count]=*event;
+    ++captured_count;
+  }
+  return 0;
+}
+#define FRAME_PSF_DIAGNOSTIC 1
+#include "../src/frame.c"
+#include "lens_map.h"
+
+static size_t number(const char *s, size_t max) {
+  char *end; unsigned long n = strtoul(s,&end,10);
+  if (!*s || *end || *s=='-' || n>max) { fputs("invalid bound\n",stderr); exit(2); }
+  return n;
+}
+/* Build the same one-degree ownership index used by production, from only
+ * the already bounded CSV. All tile arrays become immutable before workers. */
+static size_t star_tile(const Star *star) {
+  double ra=atan2(star->direction[1],star->direction[0])*180/pi;
+  if(ra<0)ra+=360;
+  const double dec=asin(fmax(-1,fmin(1,star->direction[2])))*180/pi;
+  return (size_t)fmin(179,floor(dec+90))*360+(size_t)fmin(359,floor(ra));
+}
+static int index_subset(StarCatalog *catalog) {
+  catalog->tiles=calloc(CATALOG_ALL_SKY_TILE_COUNT,sizeof *catalog->tiles);
+  if(!catalog->tiles)return -1;
+  for(size_t i=0;i<catalog->count;++i)++catalog->tiles[star_tile(&catalog->stars[i])].count;
+  for(size_t t=0;t<CATALOG_ALL_SKY_TILE_COUNT;++t) {
+    CatalogTile *tile=&catalog->tiles[t];
+    if(tile->count) {
+      tile->stars=malloc(tile->count*sizeof *tile->stars);
+      if(!tile->stars)return -1;
+      tile->state=1;tile->count=0;
+    }
+  }
+  for(size_t i=0;i<catalog->count;++i) {
+    CatalogTile *tile=&catalog->tiles[star_tile(&catalog->stars[i])];
+    tile->stars[tile->count++]=catalog->stars[i];
+  }
+  catalog->kind=STAR_CATALOG_ALL_SKY;
+  return 0;
+}
+int main(int argc, char **argv) {
+  if (argc != 8 && argc != 10 && argc != 12) {
+    fprintf(stderr,"Usage: capture_psf MAP SUBSET.csv FIRST LAST MAX_EVENTS OUT.events|--diagnose|--diagnose-parallel EXPOSURE [MAX_CACHE_FLUX MIN_Y [MIN_MAG MAX_MAG]]\n"
+      "Indexed diagnostics: --diagnose-indexed (parallel), --diagnose-indexed-serial.\nBounds: CSV <=8 MiB / 32768 stars; MAP <=64 MiB / one frame; triangle range <=131072; events <=65536.\n");
+    return 2;
+  }
+  struct stat st;
+  if (stat(argv[1],&st) || st.st_size>64*1024*1024 ||
+      stat(argv[2],&st) || st.st_size>8*1024*1024) return 2;
+  const size_t first=number(argv[3],131072), last=number(argv[4],262144);
+  limit=number(argv[5],65536);
+  if (!limit || last<=first || last-first>131072) return 2;
+  const double exposure=strtod(argv[7],NULL);
+  const double max_flux=argc>=10 ? strtod(argv[8],NULL) : 1e8;
+  const double min_y=argc>=10 ? strtod(argv[9],NULL) : 0;
+  if (!isfinite(exposure) || exposure<=0 || !isfinite(max_flux) || max_flux<1 || !isfinite(min_y) || min_y<0) return 2;
+  if(argc==12) {
+    selected_min_mag=strtod(argv[10],NULL);selected_max_mag=strtod(argv[11],NULL);
+    if(!isfinite(selected_min_mag) || !isfinite(selected_max_mag) || selected_min_mag<0 || selected_max_mag<=selected_min_mag)return 2;
+  }
+  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 ||
+      last>map.frames[0].mesh.triangle_count || catalog_load_csv(&catalog,argv[2]) ||
+      catalog.count>32768) return 2;
+  const int indexed=!strcmp(argv[6],"--diagnose-indexed") || !strcmp(argv[6],"--diagnose-indexed-serial");
+  if(indexed && index_subset(&catalog))return 1;
+  const char *kind=indexed ? "indexed-subset" : "memory";
+  const PointSpreadFunction psf={2.7,4.5}; PsfKernelCache cache={0};
+  if (psf_kernel_cache_init(&cache,&psf,1e-8)) return 1;
+  diagnosis=indexed || !strcmp(argv[6],"--diagnose") || !strcmp(argv[6],"--diagnose-parallel");
+  if(!strcmp(argv[6],"--diagnose-parallel") || !strcmp(argv[6],"--diagnose-indexed")) {
+    size_t totals[4]={0},total_events=0; uint64_t total_checksum=0;int failed=0,workers=0;
+    const double start=omp_get_wtime();
+#pragma omp parallel num_threads(omp_get_max_threads()<16 ? omp_get_max_threads() : 16) reduction(+:totals[:4],total_events,total_checksum) reduction(|:failed)
+    {
+#pragma omp single
+      workers=omp_get_num_threads();
+      PsfEventSink local={.cache=&cache};
+#pragma omp for schedule(dynamic,1)
+      for(size_t t=first;t<last;++t) {
+        CatalogSplatStats result=splat_catalog_triangles(&map.frames[0].mesh,&catalog,NULL,
+          map.width,map.height,exposure,&psf,&cache,1e6,max_flux,1e-8,min_y,t,t+1,&local);
+        failed|=result.failed;
+      }
+      for(int k=0;k<4;++k)totals[k]+=classified[k];
+      total_events+=captured_count;total_checksum+=checksum;
+      failed|=frame_psf_diagnostic_stopped();
+    }
+    printf("DIAGNOSTIC ONLY, no HDR: workers=%d catalog=%s stars=%zu triangles=[%zu,%zu) events=%zu cached=%zu wing=%zu direct=%zu discarded=%zu checksum=%016llx producer_wall=%.9f capped_or_failed=%d\n",
+      workers,kind,catalog.count,first,last,total_events,totals[0]+totals[2],totals[2],totals[1],totals[3],(unsigned long long)total_checksum,omp_get_wtime()-start,failed);
+    psf_kernel_cache_destroy(&cache);catalog_destroy(&catalog);lens_map_destroy(&map);
+    return failed ? 1 : 0;
+  }
+  PsfEventSink sink={.cache=&cache};
+  double start=omp_get_wtime();
+  CatalogSplatStats stats=splat_catalog_triangles(&map.frames[0].mesh,&catalog,NULL,
+      map.width,map.height,exposure,&psf,&cache,1e6,max_flux,1e-8,min_y,first,last,&sink);
+  printf("DIAGNOSTIC ONLY, no HDR: workers=1 catalog=%s stars=%zu triangles=[%zu,%zu) emitting_triangles=%zu last_emitting_triangle=%zu cached=%zu wing=%zu direct=%zu discarded=%zu events=%zu capped=%d producer_wall=%.9f digest=%016llx\n",
+    kind,catalog.count,first,last,triangles_seen,last_triangle,classified[0]+classified[2],classified[2],classified[1],classified[3],captured_count,frame_psf_diagnostic_stopped(),omp_get_wtime()-start,(unsigned long long)digest);
+  printf("checksum=%016llx\n",(unsigned long long)checksum);
+  if (stats.failed) return 1;
+  if (!diagnosis) {
+    FILE *f=fopen(argv[6],"wx"); if (!f) {perror(argv[6]);return 1;}
+    fprintf(f,"PSFEVENTS1 %d %d %zu %.17g %.17g %.17g\n",map.width,map.height,captured_count,psf.fwhm_pixels,psf.moffat_beta,1e-8);
+    for(size_t i=0;i<captured_count;++i) {
+      PsfCachedEvent e=captured[i];
+      fprintf(f,"%.17g %.17g %.17g %.17g %.17g %.17g %.17g\n",e.x,e.y,e.color.r,e.color.g,e.color.b,e.flux,e.support_radius);
+    }
+    if (ferror(f) || fclose(f)) return 1;
+  }
+  psf_kernel_cache_destroy(&cache);catalog_destroy(&catalog);lens_map_destroy(&map);
+  return 0;
+}

tests/make_psf_fixture.c

@@ -0,0 +1,26 @@
+/* Deliberately synthetic edge/phase/hotspot/wing fixture, not a benchmark
+ * substitute for captured catalog events. Output round-trips all doubles. */
+#include "optics.h"
+#include <stdio.h>
+int main(int argc,char **argv) {
+  if(argc!=2)return 2;
+  PsfKernelCache cache={0};const PointSpreadFunction psf={2.7,4.5};
+  if(psf_kernel_cache_init(&cache,&psf,1e-8))return 1;
+  FILE *out=fopen(argv[1],"wx");if(!out)return 2;
+  fprintf(out,"PSFEVENTS1 128 96 32771 2.7 4.5 1e-8\n");
+  int clipped=0;
+  for(int i=0;i<32771;++i) {
+    const double phases[]={0,0.25,0.5,0.999999999};
+    double x=-12+(i*7)%152+phases[i%4], y=-10+(i*11)%116+phases[(i/4)%4];
+    if(i%5==0){x=63.25;y=47.75;}
+    PsfCachedEvent e;LinearRgb color={0.7,0.2,0.5};
+    const int state=psf_prepare_cached_event(&e,x,y,color,i%1024==0?1000:0.01,&psf,&cache,1e8,1e-8,0);
+    if(state!=0 && state!=2)return 1;
+    clipped+=state==2;
+    fprintf(out,"%.17g %.17g %.17g %.17g %.17g %.17g %.17g\n",e.x,e.y,e.color.r,e.color.g,e.color.b,e.flux,e.support_radius);
+  }
+  fprintf(stderr,"synthetic fixture: events=32771 wing_clipped=%d; three chunks, exact/in-between phases, edges, hotspots\n",clipped);
+  const int failed=ferror(out);const int close_result=fclose(out);
+  psf_kernel_cache_destroy(&cache);
+  return failed || close_result ? 1:0;
+}

tests/replay_psf.hip

@@ -0,0 +1,246 @@
+/* Experimental pixel-owned reduction; intentionally not a production backend.
+ * Sharing this TU reuses the exact production cache indexing and row bounds. */
+#include "../src/hip_psf.hip"
+#include <omp.h>
+#include <vector>
+#include <algorithm>
+#include <cstdlib>
+#include <climits>
+#include <sys/resource.h>
+
+struct TileTask { unsigned tile, first, count, partial; };
+struct Prepared {
+  int bx, by, support;
+  size_t phase_base;
+  double tx, ty, r, g, b;
+};
+__global__ static void prepare_tiles(const PsfCachedEvent *events, size_t count,
+    int phases, int radius, double max_radius, Prepared *prepared, int *bounds) {
+  const size_t thread=(size_t)blockIdx.x*blockDim.x+threadIdx.x;
+  const size_t i=thread/32;const int lane=thread%32;
+  if(i>=count)return;
+  const PsfCachedEvent e=events[i];
+  const int bx=(int)floor(e.x),by=(int)floor(e.y);
+  const double fx=e.x-bx,fy=e.y-by,ax=fx*phases,ay=fy*phases;
+  const int x0=(int)floor(ax),y0=(int)floor(ay);
+  const int support=(int)fmin(ceil(e.support_radius),max_radius);
+  if(!lane)prepared[i]={bx,by,support,weight_index(phases,radius,x0,y0,0,0),
+                         ax-x0,ay-y0,e.color.r*e.flux,e.color.g*e.flux,e.color.b*e.flux};
+  const int side=2*radius+1;
+  for(int row=lane;row<side;row+=32) {
+    int first=1,last=0;
+    const int dy=row-radius;
+    if(dy>=-support && dy<=support)row_range(e.support_radius,fx,fy,support,dy,&first,&last);
+    bounds[2*(i*side+row)]=first;bounds[2*(i*side+row)+1]=last;
+  }
+}
+__global__ static void tile_partial(const Prepared *events, const int *bounds, const unsigned *refs,
+    const TileTask *tasks, int tile_size, int tiles_x, int width, int height,
+    const float *weights, int phases, int radius, double *partial, double *hdr) {
+  const TileTask task=tasks[blockIdx.x];
+  const int local=threadIdx.x;
+  const int px=(task.tile%tiles_x)*tile_size+local%tile_size;
+  const int py=(task.tile/tiles_x)*tile_size+local/tile_size;
+  const int side=2*radius+1;
+  const size_t plane=(size_t)side*side;
+  double r=0,g=0,b=0;
+  if(px<width && py<height) for(unsigned k=0;k<task.count;++k) {
+    const unsigned id=refs[task.first+k];
+    const Prepared e=events[id];
+    const int dx=px-e.bx,dy=py-e.by;
+    if(dy < -e.support || dy > e.support) continue;
+    const size_t row=2*((size_t)id*side+dy+radius);
+    if(dx<bounds[row] || dx>bounds[row+1])continue;
+    const size_t p=e.phase_base+(ptrdiff_t)dy*side+dx;
+    const double w00=weights[p],w10=weights[p+plane];
+    const double w01=weights[p+(phases+1)*plane],w11=weights[p+(phases+2)*plane];
+    const double w=(1-e.ty)*((1-e.tx)*w00+e.tx*w10)+e.ty*((1-e.tx)*w01+e.tx*w11);
+    r+=e.r*w;g+=e.g*w;b+=e.b*w;
+  }
+  if(task.partial==UINT_MAX) {
+    if(px<width && py<height) {
+      const size_t p=3*((size_t)py*width+px);
+      hdr[p]+=r;hdr[p+1]+=g;hdr[p+2]+=b;
+    }
+  } else {
+    const size_t p=3*((size_t)task.partial*tile_size*tile_size+local);
+    partial[p]=r;partial[p+1]=g;partial[p+2]=b;
+  }
+}
+__global__ static void tile_merge(const unsigned *starts, const TileTask *tasks, int tile_size,int tiles_x,
+    int width,int height,const double *partial,double *hdr) {
+  const unsigned tile=blockIdx.x; const int local=threadIdx.x;
+  const int px=(tile%tiles_x)*tile_size+local%tile_size;
+  const int py=(tile/tiles_x)*tile_size+local/tile_size;
+  if(px>=width || py>=height || starts[tile+1]-starts[tile]<=1) return;
+  double r=0,g=0,b=0;
+  for(unsigned k=starts[tile];k<starts[tile+1];++k) {
+    const size_t p=3*((size_t)tasks[k].partial*tile_size*tile_size+local);
+    r+=partial[p];g+=partial[p+1];b+=partial[p+2];
+  }
+  const size_t p=3*((size_t)py*width+px);
+  hdr[p]+=r;hdr[p+1]+=g;hdr[p+2]+=b;
+}
+static void check(hipError_t e) {if(e!=hipSuccess){fprintf(stderr,"HIP: %s\n",hipGetErrorString(e));exit(1);}}
+static double sync_time(HipPsfSink *s,double start) {check(hipStreamSynchronize(s->stream));return omp_get_wtime()-start;}
+static size_t contributions(const std::vector<PsfCachedEvent>& events,int width,int height,const PsfKernelCache& cache) {
+  size_t n=0;
+  for(const auto&e:events) {
+    const int bx=(int)floor(e.x),by=(int)floor(e.y);
+    const double fx=e.x-bx,fy=e.y-by;
+    const int support=(int)fmin(ceil(e.support_radius),cache.max_radius_pixels);
+    for(int dy=-support;dy<=support;++dy) {
+      if(by+dy<0 || by+dy>=height) continue;
+      const double left=e.support_radius*e.support_radius-(dy+0.5-fy)*(dy+0.5-fy);
+      if(left<0) continue;
+      int lo=std::max(-support,(int)ceil(fx-0.5-sqrt(left)));
+      int hi=std::min(support,(int)floor(fx-0.5+sqrt(left)));
+      lo=std::max(lo,-bx);hi=std::min(hi,width-bx-1);
+      if(hi>=lo)n+=(size_t)(hi-lo+1);
+    }
+  }
+  return n;
+}
+int main(int argc,char **argv) {
+  if(argc<4 || argc>5) {puts("Usage: replay_psf INPUT.events EVENTS(1..65536) atomic|16|32 [disperse|mixed]\nExternal timeout <=45s; sequential processes only.");return 2;}
+  char *end=nullptr;const long requested=strtol(argv[2],&end,10);
+  if(!*argv[2] || *end || requested<1 || requested>65536) return 2;
+  const int tile=!strcmp(argv[3],"atomic") ? 0 : !strcmp(argv[3],"16") ? 16 : !strcmp(argv[3],"32") ? 32 : -1;
+  if(tile<0 || (argc==5 && strcmp(argv[4],"disperse") && strcmp(argv[4],"mixed")))return 2;
+  FILE *f=fopen(argv[1],"r");if(!f){perror(argv[1]);return 2;}
+  int width,height;size_t count;PointSpreadFunction psf;double tail;
+  if(fscanf(f,"PSFEVENTS1 %d %d %zu %lf %lf %lf",&width,&height,&count,&psf.fwhm_pixels,&psf.moffat_beta,&tail)!=6 ||
+      width<1 || height<1 || width>3840 || height>2160 || count>65536 || count<(size_t)requested ||
+      !std::isfinite(psf.fwhm_pixels) || psf.fwhm_pixels<=0 || psf.fwhm_pixels>2.7 ||
+      !std::isfinite(psf.moffat_beta) || psf.moffat_beta<4.5 || !std::isfinite(tail) || tail<1e-8 || tail>=1) return 2;
+  PsfKernelCache cache={};if(psf_kernel_cache_init(&cache,&psf,tail))return 1;
+  psf_kernel_cache_report_ready(&cache,stdout);
+  std::vector<PsfCachedEvent> events(requested);
+  for(auto &e:events) {
+    if(fscanf(f,"%lf %lf %lf %lf %lf %lf %lf",&e.x,&e.y,&e.color.r,&e.color.g,&e.color.b,&e.flux,&e.support_radius)!=7) return 2;
+    const double fields[]={e.x,e.y,e.color.r,e.color.g,e.color.b,e.flux,e.support_radius};
+    for(double v:fields)if(!std::isfinite(v))return 2;
+    if(fabs(e.x)>width+cache.max_radius_pixels || fabs(e.y)>height+cache.max_radius_pixels || e.support_radius<0 || e.support_radius>1e6)return 2;
+  }
+  fclose(f);
+  // Wing-clipped events retain their physical radius; traversal alone is cache-limited.
+  const size_t original_contributions=contributions(events,width,height,cache);
+  const bool mixed=argc==5 && !strcmp(argv[4],"mixed");
+  if(argc==5 && !mixed) {
+    unsigned state=12345;
+    for(auto&e:events) {
+      // Preserve fractional phase and clipping: move only fully interior events.
+      const int s=(int)ceil(e.support_radius)+1;
+      if(e.x<s || e.y<s || e.x>=width-s || e.y>=height-s)continue;
+      state=1664525u*state+1013904223u; e.x=s+state%(width-2*s)+(e.x-floor(e.x));
+      state=1664525u*state+1013904223u; e.y=s+state%(height-2*s)+(e.y-floor(e.y));
+    }
+  }
+  const size_t work=contributions(events,width,height,cache);
+  if(work!=original_contributions)return 1;
+  const size_t values=(size_t)width*height*3;
+  std::vector<double> cpu(values,0),gpu(values,0);
+  const LinearRgb direct_color={0.7,0.2,0.5};
+  auto direct=[&](double *hdr){splat_moffat_direct(hdr,width,height,12.25,20.75,direct_color,0.2,&psf,tail,0);};
+  double start=omp_get_wtime();
+  for(size_t i=0;i<events.size();++i) {
+    splat_prepared_cached_event(cpu.data(),width,height,&events[i],&cache);
+    if(mixed && i+1==std::min(events.size(),(size_t)16384))direct(cpu.data());
+  }
+  const double cpu_time=omp_get_wtime()-start;
+  char message[256]={};HipPsfSink*s=nullptr;
+  start=omp_get_wtime();if(hip_psf_sink_create(&s,width,height,&cache,16384,message,sizeof message)){fprintf(stderr,"%s\n",message);return 1;}
+  const double create=omp_get_wtime()-start;
+  hipDeviceProp_t prop;check(hipGetDeviceProperties(&prop,0));
+  printf("device=%s arch=%s events=%zu frame=%dx%d mode=%s distribution=%s contributions=%zu CPU_reference=%.9f create=%.9f\n",prop.name,prop.gcnArchName,events.size(),width,height,argv[3],mixed?"fixture-mixed-direct":argc==5?"synthetic-phase-preserving-dispersal":"captured",work,cpu_time,create);fflush(stdout);
+  double bin=0,upload=0,prepare=0,kernel=0,merge=0,clear=0;size_t peak_scratch=0,ref_total=0,task_total=0;
+  Prepared *prepared=nullptr;int *bounds=nullptr;
+  unsigned *drefs=nullptr,*dstarts=nullptr;TileTask *dtasks=nullptr;double *partial=nullptr;
+  auto direct_boundary=[&]() {
+    if(hip_psf_sink_finish(s,gpu.data(),message,sizeof message)) {fprintf(stderr,"%s\n",message);exit(1);}
+    direct(gpu.data());
+    if(hip_psf_sink_load_hdr(s,gpu.data(),message,sizeof message)) {fprintf(stderr,"%s\n",message);exit(1);}
+  };
+  const double replay_start=omp_get_wtime();
+  if(tile) {
+    // Fixed hard bounds; no full-frame event list or per-chunk allocation on GPU.
+    check(hipMalloc(&prepared,16384*sizeof(Prepared)));
+    check(hipMalloc(&bounds,16384*(size_t)(2*cache.radius_pixels+1)*2*sizeof(int)));
+    check(hipMalloc(&drefs,32*1024*1024));check(hipMalloc(&dtasks,2*1024*1024));
+    check(hipMalloc(&dstarts,131072));check(hipMalloc(&partial,128*1024*1024));
+  }
+  for(size_t offset=0;offset<events.size();offset+=16384) {
+    const size_t n=std::min((size_t)16384,events.size()-offset);
+    if(!tile) {if(hip_psf_sink_submit(s,events.data()+offset,n,message,sizeof message)) {fprintf(stderr,"%s\n",message);return 1;}if(mixed && offset==0)direct_boundary();continue;}
+    start=omp_get_wtime();
+    const int nx=(width+tile-1)/tile,ny=(height+tile-1)/tile;
+    std::vector<std::vector<unsigned>> lists(nx*ny);
+    size_t refs_count=0;
+    for(size_t i=0;i<n;++i) {
+      const auto&e=events[offset+i];const int support=(int)fmin(ceil(e.support_radius),cache.max_radius_pixels);
+      const int bx=(int)floor(e.x),by=(int)floor(e.y);
+      const int x0=std::max(0,bx-support),x1=std::min(width-1,bx+support);
+      const int y0=std::max(0,by-support),y1=std::min(height-1,by+support);
+      if(x0>x1 || y0>y1)continue;
+      for(int y=y0/tile;y<=y1/tile;++y)for(int x=x0/tile;x<=x1/tile;++x) {
+        if(++refs_count>32*1024*1024/sizeof(unsigned)){fputs("reference limit exceeded\n",stderr);return 2;}
+        lists[y*nx+x].push_back((unsigned)i);
+      }
+    }
+    std::vector<unsigned> refs,starts(nx*ny+1);refs.reserve(refs_count);
+    std::vector<TileTask> tasks;size_t partial_count=0;
+    for(unsigned t=0;t<lists.size();++t) {
+      starts[t]=tasks.size();
+      for(size_t k=0;k<lists[t].size();k+=256) {
+        const unsigned length=std::min((size_t)256,lists[t].size()-k);
+        const unsigned part=lists[t].size()<=256 ? UINT_MAX : (unsigned)partial_count++;
+        if(partial_count*(size_t)tile*tile*3*sizeof(double)>128*1024*1024 ||
+           (tasks.size()+1)*sizeof(TileTask)>2*1024*1024) {
+          fputs("partial/task limit exceeded\n",stderr);return 2;
+        }
+        tasks.push_back({t,(unsigned)refs.size(),length,part});
+        refs.insert(refs.end(),lists[t].begin()+k,lists[t].begin()+k+length);
+      }
+    }
+    starts.back()=tasks.size();
+    const size_t bytes=partial_count*(size_t)tile*tile*3*sizeof(double);
+    if(bytes>128*1024*1024 || tasks.size()*sizeof(TileTask)>2*1024*1024 || starts.size()*sizeof(unsigned)>131072){fputs("partial/task limit exceeded\n",stderr);return 2;}
+    peak_scratch=std::max(peak_scratch,bytes+refs.size()*sizeof(unsigned)+tasks.size()*sizeof(TileTask)+starts.size()*sizeof(unsigned)+n*(sizeof(Prepared)+(size_t)(2*cache.radius_pixels+1)*2*sizeof(int)));
+    ref_total+=refs.size();task_total+=tasks.size();bin+=omp_get_wtime()-start;
+    start=omp_get_wtime();
+    check(hipMemcpyAsync(s->slots[0].device,events.data()+offset,n*sizeof(PsfCachedEvent),hipMemcpyHostToDevice,s->stream));
+    if(!refs.empty())check(hipMemcpyAsync(drefs,refs.data(),refs.size()*sizeof(unsigned),hipMemcpyHostToDevice,s->stream));
+    if(!tasks.empty())check(hipMemcpyAsync(dtasks,tasks.data(),tasks.size()*sizeof(TileTask),hipMemcpyHostToDevice,s->stream));
+    check(hipMemcpyAsync(dstarts,starts.data(),starts.size()*sizeof(unsigned),hipMemcpyHostToDevice,s->stream));
+    upload+=sync_time(s,start);
+    start=omp_get_wtime();
+    hipLaunchKernelGGL(prepare_tiles,dim3((n*32+127)/128),dim3(128),0,s->stream,s->slots[0].device,n,s->phase_resolution,s->radius_pixels,s->max_radius_pixels,prepared,bounds);
+    check(hipGetLastError());prepare+=sync_time(s,start);
+    start=omp_get_wtime();
+    if(!tasks.empty())hipLaunchKernelGGL(tile_partial,dim3(tasks.size()),dim3(tile*tile),0,s->stream,prepared,bounds,drefs,dtasks,tile,nx,width,height,s->weights,s->phase_resolution,s->radius_pixels,partial,s->hdr);
+    check(hipGetLastError());kernel+=sync_time(s,start);
+    start=omp_get_wtime();hipLaunchKernelGGL(tile_merge,dim3(nx*ny),dim3(tile*tile),0,s->stream,dstarts,dtasks,tile,nx,width,height,partial,s->hdr);
+    check(hipGetLastError());merge+=sync_time(s,start);
+    if(mixed && offset==0)direct_boundary();
+    // Every partial pixel is written by exactly one thread: no memset required.
+  }
+  if(hip_psf_sink_finish(s,gpu.data(),message,sizeof message)){fprintf(stderr,"%s\n",message);return 1;}
+  if(tile) {start=omp_get_wtime();check(hipFree(prepared));check(hipFree(bounds));check(hipFree(drefs));check(hipFree(dtasks));check(hipFree(dstarts));check(hipFree(partial));clear=omp_get_wtime()-start;}
+  const double wall=omp_get_wtime()-replay_start;
+  HipPsfTiming timing={};hip_psf_sink_get_timing(s,&timing);
+  if(!tile){upload=timing.upload_seconds;kernel=timing.kernel_seconds;}
+  double max_abs=0,max_rel=0;long double sums[3]={},diffs[3]={};bool pass=true;
+  for(size_t i=0;i<values;++i) {
+    if(!std::isfinite(cpu[i]) || !std::isfinite(gpu[i]))pass=false;
+    const double d=fabs(gpu[i]-cpu[i]);max_abs=fmax(max_abs,d);
+    if(cpu[i]!=0)max_rel=fmax(max_rel,d/fabs(cpu[i]));
+    sums[i%3]+=cpu[i];diffs[i%3]+=(long double)gpu[i]-cpu[i];
+  }
+  long double ysum=.2126L*sums[0]+.7152L*sums[1]+.0722L*sums[2],ydiff=.2126L*diffs[0]+.7152L*diffs[1]+.0722L*diffs[2];
+  struct rusage usage;getrusage(RUSAGE_SELF,&usage);
+  printf("complete_replay=%.9f replay_wall=%.9f bin=%.9f upload=%.9f prepare=%.9f accumulation=%.9f merge=%.9f download=%.9f cleanup=%.9f events_s=%.3f contributions_s=%.3f refs=%zu tasks=%zu scratch_used_peak=%zu scratch_device_reserved=%zu RSS_KiB=%ld max_abs=%.17g max_rel=%.17g flux_R=%.17Lg flux_G=%.17Lg flux_B=%.17Lg flux_Y=%.17Lg\n",create+wall,wall,bin,upload,prepare,kernel,merge,timing.download_seconds,clear,events.size()/wall,work/wall,ref_total,task_total,peak_scratch,tile?162*1024*1024+131072+16384*(sizeof(Prepared)+(size_t)(2*cache.radius_pixels+1)*2*sizeof(int)):0,usage.ru_maxrss,max_abs,max_rel,sums[0]?diffs[0]/sums[0]:0,sums[1]?diffs[1]/sums[1]:0,sums[2]?diffs[2]/sums[2]:0,ysum?ydiff/ysum:0);
+  for(int c=0;c<3;++c)if(fabsl(diffs[c])>1e-10L*fmaxl(fabsl(sums[c]),1e-30L))pass=false;
+  if(fabsl(ydiff)>1e-10L*fmaxl(fabsl(ysum),1e-30L))pass=false;
+  hip_psf_sink_destroy(s);psf_kernel_cache_destroy(&cache);
+  return pass && max_abs<1e-10 && max_rel<1e-10 ? 0 : 1;
+}

tests/test_psf_capture.py

@@ -0,0 +1,45 @@
+#!/usr/bin/env python3
+"""Bounded diagnostic accounting, output refusal, and parallel checksum checks."""
+from pathlib import Path
+import os
+import re
+import subprocess
+import sys
+import tempfile
+
+renderer=Path(sys.argv[1]).resolve()
+capture=Path(sys.argv[2]).resolve()
+repo=Path(__file__).resolve().parents[1]
+env=dict(os.environ,OMP_NUM_THREADS='4')
+def call(command, success=True):
+    print('RUN',command,flush=True)
+    p=subprocess.run(list(map(str,command)),env=env,text=True,stdout=subprocess.PIPE,stderr=subprocess.STDOUT,timeout=30)
+    print(p.stdout,flush=True)
+    assert (p.returncode==0)==success,(p.returncode,p.stdout)
+    return p.stdout
+with tempfile.TemporaryDirectory(prefix='gr-capture-test-') as temp:
+    temp=Path(temp);lens=temp/'mixed.grlens';events=temp/'mixed.events'
+    catalog=repo/'tests/data/hip_psf_mixed_catalog.csv'
+    call([renderer,'--catalog',catalog,'--width','96','--height','64','--fov-deg','30',
+          '--look-ra-deg','2','--look-dec-deg','1','--coarse-cell-pixels','16','--refine-max-level','0',
+          '--exposure','1000','--lens-map-output',lens,'--output',temp/'reference.png'])
+    base=[capture,lens,catalog,'0','48','65536']
+    for min_y,expected in [('0',(2,1,0)),('1e30',(0,0,3))]:
+        summaries=[]
+        for mode in ['--diagnose','--diagnose-parallel','--diagnose-indexed-serial','--diagnose-indexed']:
+            text=call(base+[mode,'1000','1',min_y])
+            counts=tuple(int(re.search(key+r'=(\d+)',text)[1]) for key in ['cached','direct','discarded'])
+            assert counts==expected,(counts,expected)
+            assert 'no HDR' in text
+            summaries.append(re.search(r'checksum=([0-9a-f]+)',text)[1])
+        assert len(set(summaries))==1,summaries
+    call(base+[events,'1000','1','0'])
+    lines=events.read_text().splitlines()
+    assert lines[0].split()[3]=='2' and len(lines)==3
+    call(base+[events,'1000','1','0'],False) # exclusive creation
+    call([capture,lens,catalog,'0','49','1','--diagnose','1000'],False)
+    call([capture,lens,catalog,'0','48','65537','--diagnose','1000'],False)
+    capped=temp/'capped.events'
+    call([capture,lens,catalog,'0','48','1',capped,'1000'])
+    assert len(capped.read_text().splitlines())==2
+print('PASS: capture bounds, classifications, exclusive output, serial/parallel checksum')