Task: smith_waterman
Category: Spatial--temporal algorithms
Difficulty: 3
Tags: dynamic_programming,sequence_alignment,wavefront,batched,local_alignment
Source: CUDASW++4.0 reference: https://github.com/asbschmidt/CUDASW4
Description: Smith--Waterman local sequence-alignment dynamic program for sequence-pair scoring.

Prompt template (YAML):
# Smith-Waterman Local Sequence Alignment -- Prompt Template

task_description: |
  Given `N` independent pairs of DNA sequences, compute the Smith-Waterman local alignment
  score for each pair.

  For a query sequence `Q` of length `m` and a target sequence `T` of length `n`, define
  a dynamic-programming table `H` with shape `(m+1) x (n+1)` and initialize row 0 and
  column 0 to zero. For `i >= 1, j >= 1`:

  ```
  s = (Q[i-1] == T[j-1]) ? match_score : mismatch_penalty
  H[i][j] = max(
      0,
      H[i-1][j-1] + s,
      H[i-1][j]   + gap_penalty,
      H[i][j-1]   + gap_penalty
  )
  ```

  The output for one pair is `max(H[i][j])` over the whole table.
  All sequence pairs are independent, so the task is naturally batched.

interface: |
  ```c
  extern "C" void solution_compute(
      int N,
      int total_query_len,
      int total_target_len,
      int match_score,
      int mismatch_penalty,
      int gap_penalty,
      const int* query_seqs,      // host, [total_query_len], symbols in {0,1,2,3}
      const int* target_seqs,     // host, [total_target_len], symbols in {0,1,2,3}
      const int* query_offsets,   // host, [N+1], CSR offsets into query_seqs
      const int* target_offsets,  // host, [N+1], CSR offsets into target_seqs
      int* scores                 // host, [N], output: best local-alignment score per pair
  );

  ```

input_size_notes: |
  - `query_offsets[p] ... query_offsets[p+1]-1` is the p-th query sequence.
  - `target_offsets[p] ... target_offsets[p+1]-1` is the p-th target sequence.
  - All symbols are encoded as `0=A, 1=C, 2=G, 3=T`.
  - Sequence lengths vary across pairs.
  - Scoring parameters are small integers; expected outputs are exact integers.

output_contract: |
  - `scores[p]` must equal the exact optimal Smith-Waterman local-alignment score for pair `p`.
  - Output is one integer per pair.
  - `solution_compute` must be deterministic, idempotent, and synchronize before returning.
  - `scores` is a HOST pointer and must contain valid values when the function returns.

algorithm_background: |
  The Smith-Waterman algorithm is the standard dynamic-programming method for local sequence
  alignment. Unlike global alignment, it allows the optimal alignment to start and end at
  arbitrary positions inside the two sequences.

  Dependency structure:
  - Each DP cell depends only on its left, up, and upper-left neighbors.
  - Therefore, cells on the same anti-diagonal are independent and can be computed in parallel.
  - Across the batch, different sequence pairs are completely independent.

  Classical complexity is `O(m*n)` time and `O(n)` or `O(m*n)` memory per pair.

hints_l2: |
  - Exploit **inter-pair parallelism** first: many sequence pairs can be processed independently.
  - Within one pair, use **anti-diagonal wavefront parallelism** because cells on the same anti-diagonal are independent.
  - Avoid allocating and freeing GPU memory repeatedly inside inner loops.
  - Variable sequence lengths can cause load imbalance; consider grouping pairs with similar lengths.
  - You only need the maximum score, not the traceback path.

hints_l1: |
  ### 1. Batch-level parallelism
  A practical strategy is to assign one block / one warp / a small cooperative group to one
  sequence pair, then process many pairs concurrently.

  ### 2. Anti-diagonal wavefront
  For one pair, diagonal number `d` contains cells `(i, j)` such that `i + j = d`.
  All cells on the same anti-diagonal can be computed in parallel once the previous
  anti-diagonal has finished.

  ### 3. Memory reduction
  The CPU reference uses only two rows. On GPU, alternatives include:
  - two rolling diagonals,
  - two rolling rows,
  - tiled shared-memory wavefronts.

  ### 4. Workspace reuse
  Since `solution_compute` may be called repeatedly, persistent device buffers managed by
  deterministic.

  ### 5. Variable-length handling
  Pairs have different lengths. A simple correct baseline is fine, but for better utilization
  you may bucket pairs by length or launch separate kernels for small / medium / large pairs.

  ### 6. Synchronization and exactness
  Integer scores must match the CPU reference exactly. Do not approximate, band, prune, or
  change scoring semantics.


CPU reference implementation (C):
// cpu_reference.c -- smith_waterman CPU baseline (compute_only interface)
//
// No solution_init. All work happens in solution_compute.
// NO file I/O. All I/O handled by task_io_cpu.c.

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

static int max2(int a, int b) { return (a > b) ? a : b; }
static int max4(int a, int b, int c, int d) { return max2(max2(a, b), max2(c, d)); }

void solution_compute(
    int N,
    int total_query_len,
    int total_target_len,
    int match_score,
    int mismatch_penalty,
    int gap_penalty,
    const int* query_seqs,
    const int* target_seqs,
    const int* query_offsets,
    const int* target_offsets,
    int* scores
) {
    (void)total_query_len;
    (void)total_target_len;

    for (int p = 0; p < N; p++) {
        int q_start = query_offsets[p];
        int q_end   = query_offsets[p + 1];
        int t_start = target_offsets[p];
        int t_end   = target_offsets[p + 1];
        int m = q_end - q_start;
        int n = t_end - t_start;

        int* prev = (int*)calloc((size_t)(n + 1), sizeof(int));
        int* curr = (int*)calloc((size_t)(n + 1), sizeof(int));
        int best = 0;

        if (!prev || !curr) {
            free(prev);
            free(curr);
            return;
        }

        for (int i = 1; i <= m; i++) {
            int qi = query_seqs[q_start + i - 1];
            for (int j = 1; j <= n; j++) {
                int tj = target_seqs[t_start + j - 1];
                int subst = (qi == tj) ? match_score : mismatch_penalty;

                int val = max4(
                    0,
                    prev[j - 1] + subst,
                    prev[j] + gap_penalty,
                    curr[j - 1] + gap_penalty
                );
                curr[j] = val;
                if (val > best) best = val;
            }

            {
                int* tmp = prev;
                prev = curr;
                curr = tmp;
            }
            memset(curr, 0, (size_t)(n + 1) * sizeof(int));
        }

        scores[p] = best;
        free(prev);
        free(curr);
    }
}

void solution_free(void) {
    // CPU baseline has no persistent state.
}