Catalog: parallelize all-sky tile prefetch

src/catalog.c

@@ -6,6 +6,8 @@
 #include <stdlib.h>
 #include <string.h>
 
+#include <omp.h>
+
 #define PI 3.14159265358979323846
 
 enum {
@@ -108,26 +110,40 @@ static size_t tile_index(int ra_index, int dec_index)
     return (size_t)dec_index * CATALOG_ALL_SKY_RA_TILES + ra_index;
 }
 
+static int load_tile_file(const StarCatalog *catalog, int ra_index,
+                          int dec_index, Star **stars, size_t *count)
+{
+    char path[PATH_MAX];
+    int written;
+    if (catalog == NULL || stars == NULL || count == NULL)
+        return -1;
+    *stars = NULL;
+    *count = 0;
+    written = snprintf(path, sizeof path, "%s/tile_ra%03d_dec%03d.csv",
+                       catalog->all_sky_root, ra_index, dec_index);
+    if (written < 0 || (size_t)written >= sizeof path)
+        return -1;
+    StarCatalog temporary = {0};
+    if (catalog_load_csv(&temporary, path))
+        return -1;
+    *stars = temporary.stars;
+    *count = temporary.count;
+    return 0;
+}
+
 static int load_tile(StarCatalog *catalog, int ra_index, int dec_index)
 {
     CatalogTile *tile = &catalog->tiles[tile_index(ra_index, dec_index)];
-    char path[PATH_MAX];
-    int written;
+    Star *stars;
+    size_t count;
     if (tile->state != 0)
         return tile->state == 1 ? 0 : -1;
-    written = snprintf(path, sizeof path, "%s/tile_ra%03d_dec%03d.csv",
-                       catalog->all_sky_root, ra_index, dec_index);
-    if (written < 0 || (size_t)written >= sizeof path) {
+    if (load_tile_file(catalog, ra_index, dec_index, &stars, &count)) {
         tile->state = -1;
         return -1;
     }
-    StarCatalog temporary = {0};
-    if (catalog_load_csv(&temporary, path)) {
-        tile->state = -1;
-        return -1;
-    }
-    tile->stars = temporary.stars;
-    tile->count = temporary.count;
+    tile->stars = stars;
+    tile->count = count;
     tile->state = 1;
     catalog->count += tile->count;
     return 0;
@@ -255,15 +271,15 @@ static int tile_is_fully_contained(const double direction[3][3], int ra_index,
     return 1;
 }
 
-int catalog_visit_source_triangle(StarCatalog *catalog,
-                                  const double direction[3][3],
-                                  int load_missing, CatalogTileVisitor visitor,
-                                  void *context)
+typedef int (*CatalogTileIndexVisitor)(int ra_index, int dec_index,
+                                       void *context);
+
+static int visit_source_triangle_tile_indices(
+    const double direction[3][3], CatalogTileIndexVisitor visitor,
+    void *context)
 {
-    if (catalog == NULL || direction == NULL || visitor == NULL)
+    if (direction == NULL || visitor == NULL)
         return -1;
-    if (catalog->kind == STAR_CATALOG_MEMORY)
-        return visitor(catalog->stars, catalog->count, 0, context);
     double longitude[3], latitude[3], unwrapped[3];
     for (int i = 0; i < 3; ++i) {
         lon_lat_from_direction(direction[i], &longitude[i], &latitude[i]);
@@ -311,17 +327,149 @@ int catalog_visit_source_triangle(StarCatalog *catalog,
         for (int raw_ra = all_ra ? 0 : ra_first;
              raw_ra <= (all_ra ? 359 : ra_last); ++raw_ra) {
             const int ra = (raw_ra % 360 + 360) % 360;
-            CatalogTile *tile = &catalog->tiles[tile_index(ra, dec)];
-            if (load_missing && load_tile(catalog, ra, dec))
-                continue; /* Downloader has not finished this tile yet. */
-            if (tile->state != 1 ||
-                visitor(tile->stars, tile->count,
-                        tile_is_fully_contained(direction, ra, dec), context))
+            if (visitor(ra, dec, context))
                 return -1;
         }
     return 0;
 }
 
+typedef struct {
+    StarCatalog *catalog;
+    const double (*direction)[3];
+    int load_missing;
+    CatalogTileVisitor visitor;
+    void *context;
+} CatalogVisitContext;
+
+static int visit_catalog_tile(int ra_index, int dec_index, void *opaque)
+{
+    CatalogVisitContext *context = opaque;
+    CatalogTile *tile =
+        &context->catalog->tiles[tile_index(ra_index, dec_index)];
+    if (context->load_missing && load_tile(context->catalog, ra_index, dec_index))
+        return 0; /* Downloader has not finished this tile yet. */
+    if (tile->state != 1)
+        return 0;
+    return context->visitor(
+        tile->stars, tile->count,
+        tile_is_fully_contained(context->direction, ra_index, dec_index),
+        context->context);
+}
+
+typedef struct {
+    unsigned char *requested;
+} CatalogMarkContext;
+
+static int mark_catalog_tile(int ra_index, int dec_index, void *opaque)
+{
+    CatalogMarkContext *context = opaque;
+    context->requested[tile_index(ra_index, dec_index)] = 1;
+    return 0;
+}
+
+int catalog_mark_source_triangle_tiles(
+    const double direction[3][3],
+    unsigned char requested[CATALOG_ALL_SKY_TILE_COUNT])
+{
+    CatalogMarkContext context = {.requested = requested};
+    if (requested == NULL)
+        return -1;
+    return visit_source_triangle_tile_indices(direction, mark_catalog_tile,
+                                              &context);
+}
+
+typedef struct {
+    size_t tile_id;
+    Star *stars;
+    size_t count;
+    int loaded;
+} CatalogPendingTile;
+
+int catalog_prefetch_marked_tiles(
+    StarCatalog *catalog,
+    const unsigned char requested[CATALOG_ALL_SKY_TILE_COUNT],
+    int worker_count, CatalogPrefetchStats *stats)
+{
+    CatalogPendingTile *pending = NULL;
+    size_t pending_count = 0;
+    double load_start;
+    int result = -1;
+    if (stats != NULL)
+        *stats = (CatalogPrefetchStats){0};
+    if (catalog == NULL || requested == NULL || worker_count <= 0)
+        return -1;
+    if (catalog->kind != STAR_CATALOG_ALL_SKY)
+        return 0;
+    for (size_t tile_id = 0; tile_id < CATALOG_ALL_SKY_TILE_COUNT; ++tile_id) {
+        if (!requested[tile_id])
+            continue;
+        if (stats != NULL)
+            ++stats->requested_tiles;
+        if (catalog->tiles[tile_id].state == 0)
+            ++pending_count;
+    }
+    if (pending_count == 0)
+        return 0;
+    pending = calloc(pending_count, sizeof *pending);
+    if (pending == NULL)
+        return -1;
+    size_t pending_index = 0;
+    for (size_t tile_id = 0; tile_id < CATALOG_ALL_SKY_TILE_COUNT; ++tile_id)
+        if (requested[tile_id] && catalog->tiles[tile_id].state == 0)
+            pending[pending_index++].tile_id = tile_id;
+    if (worker_count > (int)pending_count)
+        worker_count = (int)pending_count;
+    load_start = omp_get_wtime();
+#pragma omp parallel for num_threads(worker_count) schedule(static)
+    for (size_t i = 0; i < pending_count; ++i) {
+        const int ra_index =
+            (int)(pending[i].tile_id % CATALOG_ALL_SKY_RA_TILES);
+        const int dec_index =
+            (int)(pending[i].tile_id / CATALOG_ALL_SKY_RA_TILES);
+        pending[i].loaded = !load_tile_file(
+            catalog, ra_index, dec_index, &pending[i].stars, &pending[i].count);
+    }
+    if (stats != NULL)
+        stats->load_seconds = omp_get_wtime() - load_start;
+    /* Only this serial commit mutates the shared catalog cache. */
+    for (size_t i = 0; i < pending_count; ++i) {
+        CatalogTile *tile = &catalog->tiles[pending[i].tile_id];
+        if (pending[i].loaded) {
+            tile->stars = pending[i].stars;
+            tile->count = pending[i].count;
+            tile->state = 1;
+            catalog->count += tile->count;
+            if (stats != NULL) {
+                ++stats->newly_loaded_tiles;
+                stats->newly_loaded_stars += tile->count;
+            }
+        } else {
+            tile->state = -1;
+            if (stats != NULL)
+                ++stats->unavailable_tiles;
+        }
+    }
+    result = 0;
+    free(pending);
+    return result;
+}
+
+int catalog_visit_source_triangle(StarCatalog *catalog,
+                                  const double direction[3][3],
+                                  int load_missing, CatalogTileVisitor visitor,
+                                  void *context)
+{
+    if (catalog == NULL || direction == NULL || visitor == NULL)
+        return -1;
+    if (catalog->kind == STAR_CATALOG_MEMORY)
+        return visitor(catalog->stars, catalog->count, 0, context);
+    CatalogVisitContext visit_context = {
+        .catalog = catalog, .direction = direction, .load_missing = load_missing,
+        .visitor = visitor, .context = context};
+    return visit_source_triangle_tile_indices(direction, visit_catalog_tile,
+                                              &visit_context);
+}
+
 void catalog_destroy(StarCatalog *catalog)
 {
     if (catalog->tiles != NULL)

src/catalog.h

@@ -9,7 +9,12 @@ typedef struct {
     double amplitude;
 } Star;
 
-enum { CATALOG_ALL_SKY_RA_TILES = 360, CATALOG_ALL_SKY_DEC_TILES = 180 };
+enum {
+    CATALOG_ALL_SKY_RA_TILES = 360,
+    CATALOG_ALL_SKY_DEC_TILES = 180,
+    CATALOG_ALL_SKY_TILE_COUNT =
+        CATALOG_ALL_SKY_RA_TILES * CATALOG_ALL_SKY_DEC_TILES
+};
 
 typedef struct {
     Star *stars;
@@ -33,6 +38,14 @@ typedef struct {
 typedef int (*CatalogTileVisitor)(const Star *stars, size_t count,
                                   int fully_contained, void *context);
 
+typedef struct {
+    size_t requested_tiles;
+    size_t newly_loaded_tiles;
+    size_t unavailable_tiles;
+    size_t newly_loaded_stars;
+    double load_seconds;
+} CatalogPrefetchStats;
+
 /*
  * Synthetic lensing fixture: stars lie on the union of 10-degree longitude
  * and latitude lines, sampled every 2 degrees.  The eight longitude/hemisphere
@@ -44,6 +57,18 @@ int catalog_load_csv(StarCatalog *catalog, const char *path);
 /* The directory contains tile_raRRR_decDDD.csv plus optional .done markers.
  * Tiles are loaded only after a source triangle intersects them. */
 int catalog_load_all_sky(StarCatalog *catalog, const char *directory);
+/* Add every one-degree tile which can intersect a source triangle to the
+ * caller-owned CATALOG_ALL_SKY_TILE_COUNT-byte request bitmap. */
+int catalog_mark_source_triangle_tiles(
+    const double direction[3][3],
+    unsigned char requested[CATALOG_ALL_SKY_TILE_COUNT]);
+/* Load each requested, previously unseen tile once. File reads use at most
+ * worker_count OpenMP workers, but cache mutation happens serially afterward,
+ * leaving the cache immutable for later splatting. */
+int catalog_prefetch_marked_tiles(
+    StarCatalog *catalog,
+    const unsigned char requested[CATALOG_ALL_SKY_TILE_COUNT],
+    int worker_count, CatalogPrefetchStats *stats);
 /* Visit only 1-degree all-sky tiles that can intersect this source triangle.
  * Call with load_missing=1 before parallel rendering, then 0 inside workers. */
 int catalog_visit_source_triangle(StarCatalog *catalog,

src/frame.c

@@ -259,19 +259,13 @@ static size_t splat_catalog_triangles(const FrameLensMesh *mesh,
   return images;
 }
 
-static int prefetch_catalog_tile(const Star *stars, size_t count,
-                                 int fully_contained, void *context) {
-  (void)stars;
-  (void)count;
-  (void)fully_contained;
-  (void)context;
-  return 0;
-}
-
 static void prefetch_catalog_for_mesh(const FrameLensMesh *mesh,
-                                      StarCatalog *catalog) {
+                                      StarCatalog *catalog,
+                                      int worker_count,
+                                      CatalogPrefetchStats *stats) {
   if (catalog->kind != STAR_CATALOG_ALL_SKY)
     return;
+  unsigned char requested[CATALOG_ALL_SKY_TILE_COUNT] = {0};
   for (size_t t = 0; t < mesh->triangle_count; ++t) {
     const LensVertex *vertex[3];
     if (!usable_triangle(mesh, &mesh->triangles[t], vertex))
@@ -280,22 +274,24 @@ static void prefetch_catalog_for_mesh(const FrameLensMesh *mesh,
         {vertex[0]->n_infinity[0], vertex[0]->n_infinity[1], vertex[0]->n_infinity[2]},
         {vertex[1]->n_infinity[0], vertex[1]->n_infinity[1], vertex[1]->n_infinity[2]},
         {vertex[2]->n_infinity[0], vertex[2]->n_infinity[1], vertex[2]->n_infinity[2]}};
-    (void)catalog_visit_source_triangle(catalog, direction, 1,
-                                        prefetch_catalog_tile, NULL);
+    (void)catalog_mark_source_triangle_tiles(direction, requested);
   }
+  (void)catalog_prefetch_marked_tiles(catalog, requested, worker_count, stats);
 }
 
 size_t frame_splat_catalog(const FrameLensMesh *mesh,
                            StarCatalog *catalog, double *hdr, int width,
                            int height, double exposure,
-                           const PointSpreadFunction *psf) {
+                           const PointSpreadFunction *psf,
+                           int catalog_load_workers,
+                           CatalogPrefetchStats *prefetch_stats) {
   if (mesh == NULL || catalog == NULL || hdr == NULL || exposure <= 0.0 ||
-      psf == NULL || width <= 0 || height <= 0)
+      psf == NULL || width <= 0 || height <= 0 || catalog_load_workers <= 0)
     return 0;
 
-  /* Tile I/O is deliberately serial and complete before OpenMP workers start.
-   * The parallel splat pass then reads an immutable tile cache. */
-  prefetch_catalog_for_mesh(mesh, catalog);
+  /* A bounded parallel read phase completes before splatting. Its serial cache
+   * commit leaves immutable tile data for the OpenMP splat workers. */
+  prefetch_catalog_for_mesh(mesh, catalog, catalog_load_workers, prefetch_stats);
 
   const size_t pixel_count = (size_t)width * height * 3;
   if (pixel_count > SIZE_MAX / sizeof(double) ||

src/frame.h

@@ -38,7 +38,9 @@ int frame_lens_mesh_trace(FrameLensMesh *mesh, const SpacetimeSource *spacetime,
 size_t frame_splat_catalog(const FrameLensMesh *mesh,
                            StarCatalog *catalog, double *hdr, int width,
                            int height, double exposure,
-                           const PointSpreadFunction *psf);
+                           const PointSpreadFunction *psf,
+                           int catalog_load_workers,
+                           CatalogPrefetchStats *prefetch_stats);
 void frame_draw_mesh(const FrameLensMesh *mesh, double *hdr, int width,
                      int height, double gray, double opacity);
 void frame_lens_mesh_destroy(FrameLensMesh *mesh);

src/main.c

@@ -31,6 +31,7 @@ typedef struct {
   double movie_start_time, movie_duration, movie_fps;
   double slab_duration;
   double minkowski_proper_acceleration;
+  int catalog_load_workers;
 } Settings;
 
 static int parse_int(const char *text, int *value) {
@@ -108,7 +109,8 @@ static int parse_args(int argc, char **argv, Settings *s,
                   .movie_duration = 2.0,
                   .movie_fps = 30.0,
                   .slab_duration = 64.0,
-                  .minkowski_proper_acceleration = 1.52};
+                  .minkowski_proper_acceleration = 1.52,
+                  .catalog_load_workers = 4};
   *write_path = NULL;
   for (int i = 1; i < argc; ++i) {
     if (!strcmp(argv[i], "--catalog") && i + 1 < argc)
@@ -158,6 +160,8 @@ static int parse_args(int argc, char **argv, Settings *s,
     } else if (!strcmp(argv[i], "--proper-acceleration") && i + 1 < argc &&
                !parse_nonnegative(argv[++i],
                                   &s->minkowski_proper_acceleration)) {
+    } else if (!strcmp(argv[i], "--catalog-load-workers") && i + 1 < argc &&
+               !parse_int(argv[++i], &s->catalog_load_workers)) {
     } else if (!strcmp(argv[i], "--write-minkowski-accel-track") &&
                i + 1 < argc)
       s->write_minkowski_accel_track_path = argv[++i];
@@ -204,14 +208,24 @@ static int render_observer_frame(const Settings *s, StarCatalog *catalog,
     free(hdr);
     return -1;
   }
-  size_t images = frame_splat_catalog(&mesh, catalog, hdr, s->width, s->height,
-                                      s->exposure, &s->psf);
+  CatalogPrefetchStats prefetch = {0};
+  size_t images = frame_splat_catalog(
+      &mesh, catalog, hdr, s->width, s->height, s->exposure, &s->psf,
+      s->catalog_load_workers, &prefetch);
   if (s->draw_mesh)
     frame_draw_mesh(&mesh, hdr, s->width, s->height, 0.5, 0.5);
   int result = write_tonemapped_image(output_path, hdr, s->width, s->height);
   fprintf(stderr, "Rendered %zu images from %zu catalog stars to %s (%s)\n",
           images, catalog->count, output_path,
           result == 0 ? "ok" : "write failed");
+  if (catalog->kind == STAR_CATALOG_ALL_SKY)
+    fprintf(stderr,
+            "Catalog prefetch: %zu requested, %zu newly loaded (%zu stars), "
+            "%zu unavailable in %.3f s; %d loader workers\n",
+            prefetch.requested_tiles, prefetch.newly_loaded_tiles,
+            prefetch.newly_loaded_stars, prefetch.unavailable_tiles,
+            prefetch.load_seconds,
+            s->catalog_load_workers);
   frame_lens_mesh_destroy(&mesh);
   free(hdr);
   return result;
@@ -289,14 +303,24 @@ static int render_movie(const Settings *s, StarCatalog *catalog,
       free(hdr);
       goto done;
     }
-    const size_t images = frame_splat_catalog(&movie.frames[i].mesh, catalog, hdr,
-                                              s->width, s->height, s->exposure, &s->psf);
+    CatalogPrefetchStats prefetch = {0};
+    const size_t images = frame_splat_catalog(
+        &movie.frames[i].mesh, catalog, hdr, s->width, s->height, s->exposure,
+        &s->psf, s->catalog_load_workers, &prefetch);
     if (s->draw_mesh)
       frame_draw_mesh(&movie.frames[i].mesh, hdr, s->width, s->height, 0.5, 0.5);
     const int write_result = write_tonemapped_image(output_path, hdr, s->width, s->height);
     free(hdr);
     fprintf(stderr, "Rendered %zu images from %zu catalog stars to %s (%s)\n",
             images, catalog->count, output_path, write_result == 0 ? "ok" : "write failed");
+    if (catalog->kind == STAR_CATALOG_ALL_SKY)
+      fprintf(stderr,
+              "Catalog prefetch: %zu requested, %zu newly loaded (%zu stars), "
+              "%zu unavailable in %.3f s; %d loader workers\n",
+              prefetch.requested_tiles, prefetch.newly_loaded_tiles,
+              prefetch.newly_loaded_stars, prefetch.unavailable_tiles,
+              prefetch.load_seconds,
+              s->catalog_load_workers);
     if (write_result)
       goto done;
   }
@@ -331,6 +355,7 @@ int main(int argc, char **argv) {
             "[--exposure E] [--observer-inward-speed V] "
             "[--psf-fwhm-pixels N] [--psf-moffat-beta N] "
             "[--coarse-cell-pixels N] [--draw-mesh] [--write-catalog PATH] "
+            "[--catalog-load-workers N] "
             "[--observer-track PATH --frames-dir DIR --frames-prefix NAME "
             "--start-time T --duration T --fps N] "
             "[--proper-acceleration A --write-minkowski-accel-track PATH]\n",