Task: nash_flows_over_time
Category: Operations research
Difficulty: 3
Tags: simulation,graphs,shortest-path
Source: Olver, Sering \
Description: Vargas Koch 2026, "Continuity, Uniqueness, and Long-Term Behavior of Nash Flows over Time," Operations Research: https://doi.org/10.1287/opre.2024.1007 & Discrete-time dynamic traffic assignment on a time-expanded network.

Prompt template (YAML):
task_description: |
  You are tasked with implementing a discrete-time dynamic traffic assignment simulation on a GPU.
  The network is a directed graph where node 0 is the source and node `num_nodes - 1` is the destination.
  Flow enters the source at a constant `inflow_rate` per time step.
  
  At each time step `t` (from 0 to `num_steps - 1`), the simulation proceeds in 5 phases:
  
  1. **Pipeline Arrivals**: For each edge `e`, flow that has completed its transit time arrives at the destination node `v`.
  2. **External Inflow**: `inflow_rate` is added to the flow waiting at node 0.
  3. **Routing**: A shortest path tree to the destination is computed using Bellman-Ford.
     The weight of an edge is calculated as:
     `w = edge_transit_time[e] * 1000 + (int)((queue_volume[e] / edge_capacity[e]) * 1000.0f)`
     Distances are integers. Ties are broken by choosing the edge with the smaller index `e`.
     After computing the shortest paths, any flow waiting at node `u` (except the destination) is moved into the queue of the chosen outgoing edge `next_edge[u]`.
  4. **Queue Processing**: For each edge `e`, up to `edge_capacity[e]` flow leaves the queue and enters the transit pipeline. It will arrive at the destination node at step `t + edge_transit_time[e]`.
  5. **Destination Accumulation**: Any flow that has reached `num_nodes - 1` is added to `total_arrived` and removed from the node.
  
  The maximum transit time for any edge is 31, so you can use a circular buffer of size 32 for the transit pipeline of each edge.

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 num_nodes,
      int num_edges,
      int num_steps,
      const int* edge_u,
      const int* edge_v,
      const float* edge_capacity,
      const int* edge_transit_time,
      int num_requests,                  // number of inflow-rate scenarios
      const float* inflow_rates,         // HOST [num_requests]
      float* results                     // HOST [num_requests] output
  );
  ```

input_size_notes: |
  - `num_nodes` is up to 5,000.
  - `num_edges` is up to 20,000.
  - `num_steps` is up to 10,000.
  - `edge_transit_time` is between 1 and 31.

output_contract: |
  - `solution_compute` setup phase should allocate any necessary GPU memory and copy the graph structure.
  - `solution_compute` must perform the simulation and write the final `total_arrived` to `out_total_arrived`.
  - `solution_compute` must be idempotent and reset all simulation state (queues, pipelines, node flows) to zero at the beginning of each call.
  - The Bellman-Ford routing must exactly match the integer arithmetic and tie-breaking logic of the reference implementation to ensure deterministic results.

algorithm_background: |
  The simulation models the Vickrey bottleneck model (fluid queueing) in discrete time.
  Flow traverses links in a fixed transit time, but if the inflow exceeds capacity, a point queue forms at the entrance of the link.
  Users dynamically route themselves along the instantaneous shortest path to the destination, considering both transit time and current queue delays.

hints_l2: |
  - Use a 1D array of size `num_edges * 32` for the transit pipelines. Index it as `e * 32 + (t % 32)`.
  - For Bellman-Ford on the GPU, you can use a boolean flag in device memory to track if any distance changed, allowing early termination.
  - To handle the tie-breaking deterministically in parallel, you can pack the distance and edge index into a single 64-bit integer (`(uint64_t)dist << 32 | (uint32_t)edge_index`) and use `atomicMin` on `unsigned long long`. This perfectly matches the CPU's sequential tie-breaking logic.

hints_l1: |
  - In `solution_compute` setup phase, allocate device arrays for the graph structure, node flows, edge queues, pipelines, distances, and next edges.
  - In `solution_compute`, launch a kernel to clear the state arrays (node flows, queues, pipelines).
  - Implement the 5 phases as separate CUDA kernels launched sequentially in a loop for `num_steps`.
  - For Bellman-Ford, initialize distances to 1,000,000,000 and the destination to 0. Launch a kernel over edges that reads `dist[v]`, computes `new_dist`, and uses `atomicMin` on a 64-bit packed value at `dist_packed[u]` to update both distance and `next_edge`.


CPU reference implementation (C):
#include <stdlib.h>
#include <string.h>
#include <math.h>

static int g_num_nodes;
static int g_num_edges;
static int g_num_steps;
static int* g_edge_u;
static int* g_edge_v;
static float* g_edge_capacity;
static int* g_edge_transit_time;

static void _orbench_old_init(
    int num_nodes,
    int num_edges,
    int num_steps,
    const int* edge_u,
    const int* edge_v,
    const float* edge_capacity,
    const int* edge_transit_time
) {
    g_num_nodes = num_nodes;
    g_num_edges = num_edges;
    g_num_steps = num_steps;
    
    g_edge_u = (int*)malloc(num_edges * sizeof(int));
    g_edge_v = (int*)malloc(num_edges * sizeof(int));
    g_edge_capacity = (float*)malloc(num_edges * sizeof(float));
    g_edge_transit_time = (int*)malloc(num_edges * sizeof(int));
    
    memcpy(g_edge_u, edge_u, num_edges * sizeof(int));
    memcpy(g_edge_v, edge_v, num_edges * sizeof(int));
    memcpy(g_edge_capacity, edge_capacity, num_edges * sizeof(float));
    memcpy(g_edge_transit_time, edge_transit_time, num_edges * sizeof(int));
}

static void _orbench_old_compute(float inflow_rate, float* out_total_arrived) {
    float* node_flow = (float*)calloc(g_num_nodes, sizeof(float));
    float* queue_volume = (float*)calloc(g_num_edges, sizeof(float));
    float* pipeline = (float*)calloc(g_num_edges * 32, sizeof(float));
    
    int* dist = (int*)malloc(g_num_nodes * sizeof(int));
    int* dist_next = (int*)malloc(g_num_nodes * sizeof(int));
    int* next_edge = (int*)malloc(g_num_nodes * sizeof(int));
    int* next_edge_next = (int*)malloc(g_num_nodes * sizeof(int));
    
    float total_arrived = 0.0f;
    
    for (int t = 0; t < g_num_steps; t++) {
        // 1. Process arrivals from pipelines
        for (int e = 0; e < g_num_edges; e++) {
            int p_idx = e * 32 + (t % 32);
            float arrive = pipeline[p_idx];
            pipeline[p_idx] = 0.0f;
            node_flow[g_edge_v[e]] += arrive;
        }
        
        // 2. Add external inflow
        node_flow[0] += inflow_rate;
        
        // 3. Route node flows
        for (int i = 0; i < g_num_nodes; i++) {
            dist[i] = 1000000000;
            next_edge[i] = -1;
        }
        dist[g_num_nodes - 1] = 0;
        
        for (int iter = 0; iter < g_num_nodes; iter++) {
            int changed = 0;
            for (int i = 0; i < g_num_nodes; i++) {
                dist_next[i] = dist[i];
                next_edge_next[i] = next_edge[i];
            }
            
            for (int e = 0; e < g_num_edges; e++) {
                int u = g_edge_u[e];
                int v = g_edge_v[e];
                int w = g_edge_transit_time[e] * 1000 + (int)((queue_volume[e] / g_edge_capacity[e]) * 1000.0f);
                
                if (dist[v] < 500000000) {
                    int new_dist = dist[v] + w;
                    if (new_dist < dist_next[u]) {
                        dist_next[u] = new_dist;
                        next_edge_next[u] = e;
                        changed = 1;
                    } else if (new_dist == dist_next[u]) {
                        if (next_edge_next[u] == -1 || e < next_edge_next[u]) {
                            if (next_edge_next[u] != e) {
                                next_edge_next[u] = e;
                                changed = 1;
                            }
                        }
                    }
                }
            }
            
            for (int i = 0; i < g_num_nodes; i++) {
                dist[i] = dist_next[i];
                next_edge[i] = next_edge_next[i];
            }
            
            if (!changed) break;
        }
        
        // Move flow to queues
        for (int u = 0; u < g_num_nodes - 1; u++) {
            if (node_flow[u] > 0.0f && next_edge[u] != -1) {
                int e = next_edge[u];
                queue_volume[e] += node_flow[u];
                node_flow[u] = 0.0f;
            }
        }
        
        // 4. Process edge queues into pipelines
        for (int e = 0; e < g_num_edges; e++) {
            float leave = queue_volume[e];
            if (leave > g_edge_capacity[e]) {
                leave = g_edge_capacity[e];
            }
            queue_volume[e] -= leave;
            
            int arrival_time = t + g_edge_transit_time[e];
            int p_idx = e * 32 + (arrival_time % 32);
            pipeline[p_idx] += leave;
        }
        
        // 5. Accumulate destination flow
        total_arrived += node_flow[g_num_nodes - 1];
        node_flow[g_num_nodes - 1] = 0.0f;
    }
    
    *out_total_arrived = total_arrived;
    
    free(node_flow);
    free(queue_volume);
    free(pipeline);
    free(dist);
    free(dist_next);
    free(next_edge);
    free(next_edge_next);
}

// ── Unified compute_only wrapper: batched over num_requests inflow rates ──
void solution_compute(int num_nodes, int num_edges, int num_steps,
                      const int* edge_u, const int* edge_v,
                      const float* edge_capacity, const int* edge_transit_time,
                      int num_requests, const float* inflow_rates, float* results) {
    _orbench_old_init(num_nodes, num_edges, num_steps,
                      edge_u, edge_v, edge_capacity, edge_transit_time);
    for (int i = 0; i < num_requests; i++) {
        _orbench_old_compute(inflow_rates[i], &results[i]);
    }
}