Task: dbscan
Category: Spatial--temporal algorithms
Difficulty: 3
Tags: DBSCAN,neighbor_search,spatial_indexing,cluster_expansion,atomics
Source: scikit-learn DBSCAN reference: https://github.com/scikit-learn/scikit-learn
Description: DBSCAN density-based clustering from spatial neighborhood queries and dense-region connectivity.

Prompt template (YAML):
# DBSCAN Density-Based Spatial Clustering — Prompt Template

task_description: |
  Implement GPU-accelerated DBSCAN (Density-Based Spatial Clustering of Applications
  with Noise) for 2D point data.

  Given N 2D points and parameters epsilon (distance threshold) and minPts (minimum
  neighbor count), classify each point as belonging to a cluster (positive integer ID)
  or as noise (-2).

  DBSCAN algorithm:
  1. For each unprocessed point, find all neighbors within distance epsilon
  2. If neighbors >= minPts: it's a core point — start a new cluster, expand from seeds
  3. If neighbors < minPts: mark as noise (may later be absorbed into a cluster)
  4. Expansion: for each seed neighbor, repeat neighbor search; if also core, add its
     neighbors to the seed list

  The main computational bottleneck is the neighbor search: O(N) per point brute-force,
  O(N^2) total. GPU acceleration via spatial indexing (grid/tree) can reduce this
  dramatically.

  Based on CUDA-DClust+ (Poudel & Gowanlock, "CUDA-DClust+: Revisiting Early
  GPU-Accelerated DBSCAN Clustering Designs", IEEE HiPC 2021).

interface: |
  ```c
  // Single end-to-end call. Do H2D + kernel + D2H and synchronize before returning.
  // Called multiple times (warmup + timed); must be idempotent.
  extern "C" void solution_compute(
      int N,
      const float* xs,
      const float* ys,
      float eps,
      int minPts,
      int* labels
  );
  ```

input_size_notes: |
  | Size   | N (points)  | epsilon | minPts |
  |--------|-------------|---------|--------|
  | small  | 10,000      | 0.015   | 4      |
  | medium | 100,000     | 0.008   | 8      |
  | large  | 500,000     | 0.004   | 8      |

  - Points are 2D (x, y) in the range [0, 1]
  - Points form Gaussian clusters with ~5% uniform noise
  - Cluster labels must be exact match (not just same clustering structure)
  - Labels: positive integers for clusters, -2 for noise

output_contract: |
  - `labels[i]` = integer cluster ID (>0) or NOISE (-2)
  - Every non-noise point must be assigned to exactly one cluster
  - Cluster IDs must be consistent (same point → same cluster across calls)
  - `solution_compute` must be idempotent and synchronize before returning

algorithm_background: |
  **DBSCAN Algorithm (Ester et al., 1996)**

  ```
  For each point p in dataset:
    if p is already processed: skip

    neighbors = regionQuery(p, eps)  // find all points within eps distance

    if |neighbors| < minPts:
      mark p as NOISE
    else:
      create new cluster C
      add p to C
      seeds = neighbors
      while seeds is not empty:
        q = seeds.pop()
        if q is NOISE: add q to C  // border point
        if q is already in a cluster: skip
        add q to C
        q_neighbors = regionQuery(q, eps)
        if |q_neighbors| >= minPts:
          seeds.addAll(q_neighbors)  // q is also a core point
  ```

  **regionQuery** is the bottleneck: brute-force is O(N) per query.
  Total: O(N^2) for brute-force, O(N log N) with spatial index.

hints_l2: |
  - **Neighbor search parallelization**: The most impactful optimization
  - Build a uniform grid (cell size = epsilon) to limit distance calculations
  - Each point only checks the 9 neighboring cells (2D) for potential neighbors
  - Grid construction: compute cell index for each point, sort by cell, build offsets
  - Cluster expansion is inherently sequential — but neighbor search within expansion
    can be parallelized
  - Use atomics for cluster label assignment (multiple threads may try to assign)
  - Consider a two-phase approach: (1) parallel neighbor count, (2) sequential expansion

hints_l1: |
  ### 1. Grid-Based Neighbor Search
  ```cuda
  // Precompute grid cell for each point
  int cell_x = (int)(x / eps);
  int cell_y = (int)(y / eps);
  int cell_id = cell_x * grid_dim_y + cell_y;

  // Sort points by cell_id (thrust::sort_by_key)
  // Build cell offsets (scan to find start/end of each cell)
  ```

  ### 2. Parallel Neighbor Counting Kernel
  ```cuda
  __global__ void count_neighbors(int N, float* xs, float* ys,
                                   int* cell_offsets, int* sorted_indices,
                                   float eps2, int* neighbor_counts) {
      int i = blockIdx.x * blockDim.x + threadIdx.x;
      if (i >= N) return;
      // Check 9 neighboring cells
      // Count points within eps distance
      neighbor_counts[i] = count;
  }
  ```

  ### 3. Cluster Expansion Strategy
  Two approaches:
  - **CPU-driven expansion + GPU neighbor search**: CPU maintains seed queue,
    launches GPU kernel per expansion step (simpler but kernel launch overhead)
  - **Fully GPU expansion**: Multiple thread blocks each expand one cluster
    simultaneously, merge conflicts on CPU (matches CUDA-DClust+ approach)

  ### 4. Atomic Cluster Assignment
  ```cuda
  atomicCAS(&labels[j], UNCLASSIFIED, cluster_id);
  ```
  Use compare-and-swap to avoid race conditions when multiple clusters
  try to claim the same point.


CPU reference implementation (C):
// cpu_reference.c — DBSCAN density-based spatial clustering (CPU baseline)
//
// Faithfully ported from fast-cuda-gpu-dbscan/multicore-cpu/dbscan.cpp
// Preserves original variable names, struct definitions, and algorithm flow.
//
// Key simplification: single-threaded (no OpenMP), brute-force neighbor search.
// This matches the serial DBSCAN algorithm complexity O(N^2).
//
// NO file I/O, NO main(). All I/O handled by task_io_cpu.c.

#include <math.h>
#include <stdlib.h>
#include <string.h>

// ===== Constants from original (common.h / dbscan.cpp) =====
#define UNCLASSIFIED -1
#define NOISE        -2
#define MARKED        0

// ===== Structs from original (multicore-cpu/dbscan.cpp) =====
typedef float real;

typedef struct {
    real x;
    real y;
} point;

// ===== Module-level state =====
static int    g_N;
static real   g_eps;       // epsilon radius
static int    g_minPts;    // minimum points for core
static point* g_points;    // point array (owned, copied from input)

static void _orbench_old_init(int N, const float* xs, const float* ys,
                   float eps, int minPts)
{
    g_N = N;
    g_eps = eps;
    g_minPts = minPts;

    g_points = (point*)malloc(N * sizeof(point));
    for (int i = 0; i < N; i++) {
        g_points[i].x = xs[i];
        g_points[i].y = ys[i];
    }
}

// ===== Neighbor search — matches processPointCPU from original =====
// Brute-force: check all N points for distance <= eps
static int calcRow(int thisPoint, int numPoints, real radius,
                   char* oneRow, point* point_Data)
{
    real myX = point_Data[thisPoint].x;
    real myY = point_Data[thisPoint].y;
    int nOfNeighbors = 0;

    for (int x = 0; x < numPoints; x++) {
        real otherX = point_Data[x].x;
        real otherY = point_Data[x].y;

        real distX = myX - otherX;
        real distY = myY - otherY;

        real Distance = sqrtf(distX * distX + distY * distY);

        if (Distance <= radius) {
            oneRow[x] = 1;
            if (x != thisPoint) {
                nOfNeighbors++;
            }
        } else {
            oneRow[x] = 0;
        }
    }
    return nOfNeighbors;
}

// ===== Cluster expansion — matches expandCluster from original =====
static int expandCluster(int thisPoint, int* clusterId, int nextClusterId,
                         int numPoints, real Eps, int MinPts,
                         int* numPointsInCluster, char* oneRow, point* point_Data)
{
    int nOfNeighbors = calcRow(thisPoint, numPoints, Eps, oneRow, point_Data);

    if (nOfNeighbors < MinPts) {
        clusterId[thisPoint] = NOISE;
        return 0;  // false
    } else {
        (*numPointsInCluster)++;
        clusterId[thisPoint] = nextClusterId;

        // Build seed queue (simple array-based queue)
        int* seeds = (int*)malloc(numPoints * sizeof(int));
        int seedHead = 0, seedTail = 0;

        for (int i = 0; i < numPoints; i++) {
            if (oneRow[i] == 1) {
                seeds[seedTail++] = i;
                clusterId[i] = nextClusterId;
            }
        }

        while (seedHead < seedTail) {
            int currentPoint = seeds[seedHead++];

            nOfNeighbors = calcRow(currentPoint, numPoints, Eps, oneRow, point_Data);

            if (nOfNeighbors >= MinPts) {
                // Core point: expand its neighbors
                for (int i = 0; i < numPoints; i++) {
                    if (oneRow[i] == 1) {
                        if (clusterId[i] == UNCLASSIFIED || clusterId[i] == NOISE) {
                            if (clusterId[i] == UNCLASSIFIED) {
                                seeds[seedTail++] = i;
                            }
                            clusterId[i] = nextClusterId;
                        }
                    }
                }
            }
            // Non-core seed points keep their cluster assignment (border points)
        }

        free(seeds);
        return 1;  // true
    }
}

// ===== Main DBSCAN — matches clusterThread from original =====
static void clusterThread(int numPoints, point* point_Data,
                          int* clusterId, real radius, int MinPts)
{
    char* oneRow = (char*)malloc(numPoints * sizeof(char));

    for (int i = 0; i < numPoints; i++) {
        clusterId[i] = UNCLASSIFIED;
    }

    int nextClusterId = 1;

    for (int i = 0; i < numPoints; i++) {
        if (clusterId[i] == UNCLASSIFIED) {
            int numPointsInCluster = 0;
            if (expandCluster(i, clusterId, nextClusterId, numPoints,
                              radius, MinPts, &numPointsInCluster,
                              oneRow, point_Data)) {
                nextClusterId++;
            }
        }
    }

    // Finalize: MARKED points that are in a cluster but not core get cluster ID
    // (In original, MARKED points get the cluster ID from their expansion)
    // Fix up: MARKED points should have been assigned during expansion
    // Any remaining MARKED points → assign to cluster 0 (shouldn't happen)
    for (int i = 0; i < numPoints; i++) {
        if (clusterId[i] == MARKED) {
            clusterId[i] = NOISE;  // fallback
        }
    }

    free(oneRow);
}

static void _orbench_old_compute(int N, int* labels)
{
    clusterThread(N, g_points, labels, g_eps, g_minPts);
}

// ── Unified compute_only wrapper (auto-migrated) ──
void solution_compute(int N, const float* xs, const float* ys, float eps, int minPts, int* labels) {
    _orbench_old_init(N, xs, ys, eps, minPts);
    _orbench_old_compute(N, labels);
}