Task 1 of 6

Three Ways Into a Cell

Two DNA sequences, and the question every genome browser asks: where do they say the same thing? Smith-Waterman answers it by scoring a matrix. H[i][j] is the best score of any alignment that ends at base i of A and base j of B, and there are exactly four candidates for it:

H[i][j] = max(
  0,                               // start fresh here
  H[i - 1][j - 1] + s(A[i], B[j]), // pair the bases
  H[i - 1][j]     - gap,           // skip a base of A
  H[i][j - 1]     - gap            // skip a base of B
)

s is +3 when the bases match and −3 when they do not; a gap costs 2. The 0 is the whole point of local alignment: an alignment that has gone badly is abandoned rather than carried, so a good match buried inside two otherwise unrelated sequences still surfaces. Row 0 and column 0 are the empty prefixes, and stay 0.

Now read the three predecessors again, because the trouble is right there. Every cell wants its neighbour up, its neighbour left, and its neighbour up-left — so no row can start before the row above it has finished, and the natural loop is about as serial as code gets. Nobody parallelises this. (Nobody parallelises it row by row, at least. The next task is the escape.)

four candidates, one winner — and the 0 is the escape hatch
Goal: write the kernel for one full pass of the recurrence — every interior cell recomputed from the matrix it was handed, with row 0 and column 0 pinned to 0.

Requirements

Hint 1 — the substitution score, without a boolean

gpu.js cannot keep a boolean in a kernel variable on the WebGL backend, so pick the score with a plain if:

let s = this.constants.mismatch;
if (a[i - 1] === b[j - 1]) s = this.constants.match;

The − 1 is because row i of the matrix belongs to base i − 1 of the sequence: the matrix has one extra row for the empty prefix.

Hint 2 — three candidates, then the floor

Math.max takes two arguments inside a kernel, so nest it:

const diag = H[i - 1][j - 1] + s;
const up = H[i - 1][j] - this.constants.gap;
const left = H[i][j - 1] - this.constants.gap;
return Math.max(0, Math.max(diag, Math.max(up, left)));

The 0 goes on the outside — it is a floor under the whole thing, not a fourth predecessor.

Same idea elsewhere

A dynamic program is a dependency graph wearing a grid costume, and this shape — each cell reading its three earlier neighbours — is shared by edit distance, dynamic time warping, the Viterbi decoder and pairwise HMMs. Whether any of them can go on a GPU is decided entirely by that graph, not by the arithmetic in the cell: CUDA, WGSL and Metal all give you thousands of threads and no way whatsoever to make one wait for another.

All tasks in Wavefronts: Aligning DNA on the Diagonal

  1. Three Ways Into a Cell
  2. The Anti-Diagonal Goes All At Once
  3. Which Cell Am I?
  4. Ride the Wavefront Down
  5. Traceback: What Did It Actually Align?
  6. Watch the Wavefront, Then Read the Bill

This page is an interactive exercise — the editor, the GPU runner and your saved progress need JavaScript. The text above is the full brief.