GR-raytracing/tests/capture_psf.c

/* Standalone bounded diagnostic. Include the production producer so its
 * private query/mapping path stays identical; renderer builds omit the hook. */
#include "optics.h"
#include <math.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;
  if (blackbody_backend_init(NULL, 0, NAN, NAN, NULL, stderr)) return 1;
  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,NULL);
        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);blackbody_backend_destroy();
    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,NULL);
  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);blackbody_backend_destroy();
  return 0;
}