# Watch the Wavefront, Then Read the Bill

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

Thirty-two bases against thirty-six, sharing a mutated stretch in the middle.
Sixty-seven launches, and this time each one renders — so the console gives you a scrubber
and you can watch the wavefront march from the top-left corner to the bottom-right, lighting
the ridge where the two sequences agree.

Then the bill. Diagonal 2 has one cell. Diagonal 34 has thirty-two. Diagonal 68 has one
again. Average occupancy across the sweep is around `|A| / 2` threads, and the
launches at either end are launches for a handful of work — which is why nobody ships this.
Real aligners cut the matrix into **tiles** and run a wavefront over the tiles
instead, so each launch has thousands of cells to chew on; or they abandon the wavefront
entirely and align a whole database at once, one sequence per thread, which is
embarrassingly parallel and needs no cleverness at all.

Both of those are engineering. The reframing — that a dynamic program nobody could
parallelise turns out to be parallel along a diagonal — is the idea, and it survives every
one of those refinements intact.

## Goal

**Goal:** paint each frame of the sweep, and plot how much work each
launch actually had against how many threads it started.

## Requirements

- The cell on the live diagonal — `x + y === d` — is painted `this.color(1, 0.6, 0.1, 1)`
- Every other cell is painted from its score: `const v = H[y][x] / this.constants.top;` then `this.color(v * 0.35, v * 0.95, v * 0.75, 1)`
- Fill `lengths` with the number of cells on each anti-diagonal, `d = 2 … 68`

## Hint 1 — the paint kernel is a two-way branch

Read the score once, then decide which of the two colours this pixel gets:

```js
const v = H[this.thread.y][this.thread.x] / this.constants.top;
if (this.thread.x + this.thread.y === d) {
  this.color(1, 0.6, 0.1, 1);
} else {
  this.color(v * 0.35, v * 0.95, v * 0.75, 1);
}
```

## Hint 2 — how long is a diagonal?

The same clipping as before, one sequence length apiece:

```js
const iStart = Math.max(1, d - 36);
const iEnd = Math.min(32, d - 1);
lengths.push(iEnd - iStart + 1);
```

Note the `d - 1` and the `1`: row 0 and column 0 are never
computed, so they are not cells this algorithm has work for.

## Hint 3 — reading the plot

The two series differ by two orders of magnitude at the ends, which is why the
plot asks for a log axis. The flat line is what every launch costs; the triangle under
it is what every launch is worth.

## Same idea elsewhere

Occupancy — how much of the machine a launch actually keeps busy — is the number
that decides whether a GPU port was worth it, and a triangular work profile is a classic way
to lose. The standard answers are the ones real aligners use: tile the domain so every
launch is full (CUDA's blocked wavefront, the same idea as cache blocking), or find a
coarser axis of parallelism and use that instead. Modern tools like WFA go further and
change the algorithm so the wavefront is short in the first place.

## Starter code

```js
// 32 × 36 bases, 67 launches, one rendered frame each.
const gpu = new GPU({ mode });

const sweep = gpu.createKernel(function (H, a, b, d) {
  const i = this.thread.y;
  const j = this.thread.x;
  if (i === 0 || j === 0) return 0;
  if (i + j !== d) return H[i][j];
  let s = this.constants.mismatch;
  if (a[i - 1] === b[j - 1]) s = this.constants.match;
  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)));
}, { output: [37, 33], constants: { match: 3, mismatch: -3, gap: 2 } });

const paint = gpu.createKernel(function (H, d) {
  // TODO 1: cells on the live diagonal (x + y === d) get this.color(1, 0.6, 0.1, 1).
  //         Everything else is shaded by its score: v = H[y][x] / this.constants.top,
  //         then this.color(v * 0.35, v * 0.95, v * 0.75, 1).
  this.color(1, 0, 1, 1);
}, { output: [37, 33], graphical: true, constants: { top: 40 } });

let H = [];
for (let i = 0; i <= 32; i++) H.push(new Array(37).fill(0));

for (let d = 2; d <= 68; d++) {
  H = await sweep(H, codesA, codesB, d);
  await paint(H, d);
  render(paint.canvas);
}

let best = 0;
for (let i = 0; i <= 32; i++) {
  for (let j = 0; j <= 36; j++) if (H[i][j] > best) best = H[i][j];
}
console.log('best score', best);

// TODO 2: how many cells does each anti-diagonal actually have?
const lengths = [];
for (let d = 2; d <= 68; d++) {
  lengths.push(0);
}

const launched = lengths.map(() => 33 * 37);
plot(
  { 'cells with work': lengths, 'threads launched': launched },
  { title: 'work vs. threads, per launch', log: true }
);
```

---

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

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