# Three Ways Into a Cell

*Task 1 of 6 · [Wavefronts: Aligning DNA on the Diagonal](https://gpu.rocks/learn/sequence-alignment-a85ca6d9.md) · GPU.js Learn*

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:

```js
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.)

## Figures

- **four candidates, one winner — and the 0 is the escape hatch**

## Goal

**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

- Return `0` for any cell with `i === 0` or `j === 0`
- Score the pair with `this.constants.match` / `this.constants.mismatch`, comparing `a[i - 1]` against `b[j - 1]`
- Take the largest of the three predecessors, each with its own adjustment — and *then* the max against `0`

## 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`:

```js
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:

```js
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.

## Starter code

```js
// One pass of the Smith-Waterman recurrence, over the whole matrix.
// The matrix is 9 rows (|A| + 1) by 10 columns (|B| + 1).
const gpu = new GPU({ mode });

const score = gpu.createKernel(function (H, a, b) {
  const i = this.thread.y;
  const j = this.thread.x;
  // TODO 1: row 0 and column 0 are the empty prefixes — return 0.
  // TODO 2: score the pair. s is this.constants.match when a[i - 1]
  //         equals b[j - 1], and this.constants.mismatch otherwise.
  // TODO 3: return the largest of 0, diag + s, up - gap and left - gap.
  return H[i][j];
}, { output: [10, 9], constants: { match: 3, mismatch: -3, gap: 2 } });

// An empty matrix: every alignment score starts life unknown.
const empty = [];
for (let i = 0; i <= 8; i++) empty.push(new Array(10).fill(0));

const once = await score(empty, codesA, codesB);
console.log('A =', seqA, '  B =', seqB);
for (let i = 1; i <= 8; i++) console.log('row', i, Array.from(once[i]).join(' '));
```

---

Interactive version: https://gpu.rocks/learn/sequence-alignment-a85ca6d9/1

[Next task](https://gpu.rocks/learn/sequence-alignment-a85ca6d9/2.md)