GR-raytracing/src/dummy_psf.c

#include "dummy_psf.h"

#include <math.h>
#include <omp.h>
#include <stdint.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

enum { DUMMY_SELECTOR_TILE = 32, DUMMY_SUPPORT_TILE = 16,
       DUMMY_TASK_EVENTS = 256 };

typedef struct {
  size_t worker, events, triangles, occupied_tiles;
  double event_bbox_diagonal, triangle_bbox_diagonal;
  double maximum_triangle_jump;
  int adaptive_tile16;
  size_t support_refs, support_tiles, support_tasks, maximum_tile_list;
  double support_pass_seconds;
} DummyChunkSample;

struct DummyPsfSink {
  DummyChunkSample *samples;
  size_t count, capacity, event_capacity;
  size_t total_events;
  int width, height, tiles_x, tiles_y;
  int support_tiles_x, support_tiles_y, support_stats;
  double max_radius_pixels;
  int failed;
  omp_lock_t lock;
};

struct DummyPsfChunk {
  DummyPsfSink *sink;
  uint32_t *tile_generation;
  PsfCachedEvent *support_events;
  uint32_t *support_counts, *support_touched;
  size_t support_touched_count;
  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, double max_radius_pixels) {
  if (out == NULL || width <= 0 || height <= 0 || event_capacity == 0 ||
      !isfinite(max_radius_pixels) || max_radius_pixels < 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;
  sink->support_tiles_x = (width + DUMMY_SUPPORT_TILE - 1) / DUMMY_SUPPORT_TILE;
  sink->support_tiles_y = (height + DUMMY_SUPPORT_TILE - 1) / DUMMY_SUPPORT_TILE;
  sink->max_radius_pixels = max_radius_pixels;
  const char *support = getenv("DUMMY_PSF_SUPPORT_STATS");
  if (support && strcmp(support, "0") && strcmp(support, "1")) {
    free(sink); return -1;
  }
  sink->support_stats = support && !strcmp(support, "1");
  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;
  chunk->sink = sink;
  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;
  }
  if (sink->support_stats) {
    const size_t count = (size_t)sink->support_tiles_x * sink->support_tiles_y;
    chunk->support_events = malloc(sink->event_capacity * sizeof *chunk->support_events);
    chunk->support_counts = calloc(count, sizeof *chunk->support_counts);
    chunk->support_touched = malloc(count * sizeof *chunk->support_touched);
    if (!chunk->support_events || !chunk->support_counts ||
        !chunk->support_touched) {
      dummy_psf_chunk_destroy(chunk);
      return -1;
    }
  }
  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);
  size_t support_refs = 0, support_tasks = 0, maximum_tile_list = 0;
  double support_pass_seconds = 0.0;
  if (chunk->sink->support_stats) {
    const double start = omp_get_wtime();
    for (size_t i = 0; i < chunk->support_touched_count; ++i)
      chunk->support_counts[chunk->support_touched[i]] = 0;
    chunk->support_touched_count = 0;
    const int tiles_x = chunk->sink->support_tiles_x;
    for (size_t i = 0; i < chunk->events; ++i) {
      const PsfCachedEvent *event = &chunk->support_events[i];
      const int support = (int)fmin(ceil(event->support_radius),
                                    chunk->sink->max_radius_pixels);
      const int bx = (int)floor(event->x), by = (int)floor(event->y);
      const int x0 = fmax(0, bx - support);
      const int x1 = fmin(chunk->sink->width - 1, bx + support);
      const int y0 = fmax(0, by - support);
      const int y1 = fmin(chunk->sink->height - 1, by + support);
      if (x0 > x1 || y0 > y1) continue;
      for (int y = y0 / DUMMY_SUPPORT_TILE; y <= y1 / DUMMY_SUPPORT_TILE; ++y)
        for (int x = x0 / DUMMY_SUPPORT_TILE; x <= x1 / DUMMY_SUPPORT_TILE; ++x) {
          const size_t tile = (size_t)y * tiles_x + x;
          if (!chunk->support_counts[tile])
            chunk->support_touched[chunk->support_touched_count++] = tile;
          ++chunk->support_counts[tile];
          ++support_refs;
        }
    }
    for (size_t i = 0; i < chunk->support_touched_count; ++i) {
      const size_t count = chunk->support_counts[chunk->support_touched[i]];
      if (count > maximum_tile_list) maximum_tile_list = count;
      support_tasks += (count + DUMMY_TASK_EVENTS - 1) / DUMMY_TASK_EVENTS;
    }
    support_pass_seconds = omp_get_wtime() - start;
  }
  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,
      .support_refs = support_refs,
      .support_tiles = chunk->support_touched_count,
      .support_tasks = support_tasks,
      .maximum_tile_list = maximum_tile_list,
      .support_pass_seconds = support_pass_seconds};
  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->sink->support_stats)
    chunk->support_events[chunk->events] = *event;

  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->support_events);
  free(chunk->support_counts);
  free(chunk->support_touched);
  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);
}

static void report_support_group(const DummyPsfSink *sink, int full_only,
                                 const char *label) {
  size_t count = 0, total_refs = 0, total_tasks = 0;
  double total_pass = 0.0;
  for (size_t i = 0; i < sink->count; ++i)
    if (!full_only || sink->samples[i].events == sink->event_capacity)
      ++count;
  if (!count || count > SIZE_MAX / (5 * sizeof(double))) return;
  double *values = malloc(5 * count * sizeof *values);
  if (!values) { fputs("Dummy PSF support report allocation failed\n", stderr); return; }
  const char *names[] = {"support_16px_references", "support_16px_tiles",
                         "support_16px_tasks", "maximum_tile_list",
                         "support_first_pass_ms"};
  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;
    values[next] = sample->support_refs;
    values[count + next] = sample->support_tiles;
    values[2 * count + next] = sample->support_tasks;
    values[3 * count + next] = sample->maximum_tile_list;
    values[4 * count + next] = 1000.0 * sample->support_pass_seconds;
    total_refs += sample->support_refs;
    total_tasks += sample->support_tasks;
    total_pass += sample->support_pass_seconds;
    ++next;
  }
  for (size_t field = 0; field < 5; ++field) {
    double *series = values + field * count;
    qsort(series, 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[0],
            percentile(series, count, 0.10), percentile(series, count, 0.25),
            percentile(series, count, 0.50), percentile(series, count, 0.75),
            percentile(series, count, 0.90), percentile(series, count, 0.99),
            series[count - 1]);
  }
  fprintf(stderr,
          "Dummy PSF %s support totals: references=%zu tasks=%zu summed_first_pass=%.3f s\n",
          label, total_refs, total_tasks, total_pass);
  free(values);
}

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");
  if (sink->support_stats) {
    fputs("Dummy PSF support diagnostic: 16px rectangular cache-limited first pass only; no references were materialized or uploaded.\n", stderr);
    report_support_group(sink, 0, "all");
    report_support_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);
}