# Traceback: What Did It Actually Align?

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

A score is a number; a biologist wants the alignment. That means walking backwards
from the best cell, asking at each step *which of the three candidates won here*, until
the score drops to `0` and the local alignment began.

The walk itself is hopelessly serial — one cell at a time, and the path is only as long as
the alignment. Putting it on a GPU would be silly. But the *question* it asks at every
step is embarrassingly parallel: "which predecessor explains this cell?" has the same answer
whenever you ask it, so compute the answer for all 90 cells in one launch and store it. A
**pointer matrix** — 1 for diagonal, 2 for up, 3 for left, 0 for "stop here" —
is what every real aligner keeps beside the scores for exactly this reason.

Knowing what *not* to move onto the GPU is part of the craft. The fill is
`|A| × |B|` cells of work; the walk is a few dozen steps of pointer chasing that
finishes before a kernel launch would have been scheduled.

## Goal

**Goal:** write the pointer kernel. Every cell reports which of its three
predecessors produced its score, or `0` if the score is `0` and the
alignment starts there.

## Requirements

- A cell whose score is `0` returns `0` — that is where the walk stops, and it also keeps the boundary safe
- Return `1` when `H[i][j]` equals `H[i - 1][j - 1] + s`
- Return `2` when it equals `H[i - 1][j] - gap`, and `3` otherwise
- Test the candidates in that order, so ties prefer the diagonal

## Hint 1 — the stop test comes first

Not just because the walk needs it: row 0 and column 0 hold zeros, and putting
the test first means those threads return before they ever look at
`a[i - 1]`.

```js
if (H[i][j] === 0) return 0;
```

## Hint 2 — comparing scores for equality is safe here

Every score in this matrix is a whole number — match, mismatch and gap are all
integers — so `===` between two of them is exact on every backend. That is a
property of this scoring scheme, not a general licence: with fractional weights you
would compare against a tolerance instead.

## Hint 3 — the whole body

```js
if (H[i][j] === 0) return 0;
let s = this.constants.mismatch;
if (a[i - 1] === b[j - 1]) s = this.constants.match;
if (H[i][j] === H[i - 1][j - 1] + s) return 1;
if (H[i][j] === H[i - 1][j] - this.constants.gap) return 2;
return 3;
```

## Same idea elsewhere

Storing a decision alongside a value so the reconstruction is cheap is the standard
move in every parallel dynamic program — the same reason a CUDA Viterbi decoder writes a
backpointer array and a GPU video encoder writes mode decisions. The serial tail stays on the
host, and the honest question on every platform is not "can this run on the GPU" but "is this
the part worth moving".

## Starter code

```js
// The fill from the last task, then one launch that says WHERE each score came from.
const gpu = new GPU({ mode });

const diagonal = gpu.createKernel(function (H, a, b, d, iStart) {
  const i = iStart + this.thread.x;
  const j = d - i;
  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)));
}, { dynamicOutput: true, constants: { match: 3, mismatch: -3, gap: 2 } });

// Anti-diagonal d covers rows i = iStart … iStart + len - 1, and j = d - i.
// Clipped at both ends: i never passes 8 (= |A|), and j never passes 9 (= |B|).
function diagonalRange(d) {
  const iStart = Math.max(1, d - 9);
  const iEnd = Math.min(8, d - 1);
  return { iStart, len: iEnd - iStart + 1 };
}

// The wavefront fill you built in the last task, folded into one helper.
async function fillMatrix() {
  const H = [];
  for (let i = 0; i <= 8; i++) H.push(new Array(10).fill(0));
  let best = 0;
  let bi = 0;
  let bj = 0;
  for (let d = 2; d <= 17; d++) {
    const range = diagonalRange(d);
    diagonal.setOutput([range.len]);
    const values = await diagonal(H, codesA, codesB, d, range.iStart);
    for (let t = 0; t < range.len; t++) {
      const i = range.iStart + t;
      H[i][d - i] = values[t];
      if (values[t] > best) {
        best = values[t];
        bi = i;
        bj = d - i;
      }
    }
  }
  return { H, best, bi, bj };
}

const pointers = gpu.createKernel(function (H, a, b) {
  const i = this.thread.y;
  const j = this.thread.x;
  // TODO 1: a cell scoring 0 is where a local alignment starts — return 0.
  // TODO 2: score the pair as before, then report which candidate won:
  //         1 = H[i - 1][j - 1] + s, 2 = H[i - 1][j] - gap, 3 = H[i][j - 1] - gap.
  return 0;
}, { output: [10, 9], constants: { match: 3, mismatch: -3, gap: 2 } });

const filled = await fillMatrix();
const dir = await pointers(filled.H, codesA, codesB);

// Walk the pointers back. Serial, tiny, and perfectly happy on the CPU.
let i = filled.bi;
let j = filled.bj;
let alignA = '';
let alignB = '';
while (dir[i][j] !== 0) {
  if (dir[i][j] === 1) {
    alignA = seqA[i - 1] + alignA;
    alignB = seqB[j - 1] + alignB;
    i--;
    j--;
  } else if (dir[i][j] === 2) {
    alignA = seqA[i - 1] + alignA;
    alignB = '-' + alignB;
    i--;
  } else {
    alignA = '-' + alignA;
    alignB = seqB[j - 1] + alignB;
    j--;
  }
}

console.log('score', filled.best);
console.log('A: ' + alignA);
console.log('B: ' + alignB);
```

---

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

[Previous task](https://gpu.rocks/learn/sequence-alignment-a85ca6d9/4.md) · [Next task](https://gpu.rocks/learn/sequence-alignment-a85ca6d9/6.md)