# Algorithmic specification of the reference computation

**PCIST-1.0** · normative annex to `METHOD.md` §2.1 · in force from 2026-09-22

## What this document is and why it exists

`METHOD.md` §2.1 defined the reference computation by naming PCIst and pointing at a
library. `gate/RESULT-DUAL.md` showed what that costs: a second implementation, written
from everything the standard and the published paper actually say, produced values for
the same evoked responses differing from the reference by 3.76× to 35.67×, with ICC(2,1)
= 0.043. Four of the five points the implementer had to guess were guessed wrongly, and
none of the four could have been guessed correctly from the published sources.

This annex removes the guessing. It specifies the computation in full, so that a third
party can reimplement it and obtain the same number to numerical tolerance, which is what
§12 has been inviting them to do since 0.3.6 without supplying the means.

**Incorporation by reference into a versioned annex in this repository is not the defect
this replaces.** The defect was that the operative content lived in a third party's
source tree at an unpinned revision. Here the content is present, versioned with the
standard, and changed only by amendment under `GOVERNANCE.md`.

## Notation

Input is a trial-averaged evoked response **Y**, a real matrix of `n_ch` channels by
`n_t` samples, and a vector **t** of sample times in milliseconds with the perturbation
at `t = 0`. Baseline times are negative. Parameters are those named in §2.1's parameter
sets: `baseline_window = (b₀, b₁)`, `response_window = (r₀, r₁)`, `k`, `min_snr`,
`max_var`, `n_steps`, `resample`, `embed`, and `avgref`.

Define the index function

> **idx(t, x)** = the number of entries of **t** strictly less than `x`.

Every window below is taken as the **half-open** index range `[idx(t, w₀), idx(t, w₁))`.
The sample at the window's upper bound is **excluded**. This is a choice, it changes the
sample count, and it is stated because it cannot be inferred.

## Step 1 — Preprocessing

1. If `avgref` is true, subtract from every channel, at each sample, the mean across
   channels at that sample. The default is false; §2.1 records why the earlier
   prescription was withdrawn.
2. Baseline correction is **not** applied by this computation unless explicitly
   requested. When requested, subtract from each channel its own mean over all samples
   with `t < −50` ms. Readings in this register baseline-correct upstream instead, and
   either route must be declared.
3. Crop **Y** and **t** to the index range `[idx(t, b₀), idx(t, r₁))`. All later windows
   are computed against the cropped **t**.
4. If `resample` is set, resample to `n = ⌊(t_last − t_first) / 1000 × resample⌋`
   samples using **Fourier-method resampling** — the signal is transformed, the spectrum
   truncated or zero-padded, and transformed back — applied along the time axis jointly
   with the time vector. Polynomial or linear interpolation is **not** equivalent and
   will not reproduce the reference. `resample` is required by §2.1 and may not be
   omitted.

## Step 2 — Dimensionality reduction

1. Take the response-window block of the preprocessed signal, `Y_r`, of shape
   `n_ch × n_response`.
2. Obtain a channel-space basis **V** whose columns are the left-singular vectors of
   `Y_r`, ordered by descending singular value `s₁ ≥ s₂ ≥ …`. Decomposing `Y_rᵀ` and
   taking its right-singular matrix gives the same basis; the reference does it that
   way and the two routes are interchangeable.

   > An earlier draft of this annex warned that decomposing `Y_r` rather than its
   > transpose "yields a basis in the wrong space." **That warning was wrong** and the
   > conformance vector set found it: for `A = U S Vᵀ`, the right-singular matrix of
   > `Aᵀ` is `U`, so the two routes agree exactly. A per-component sign flip is
   > possible and is immaterial, because every later step uses `|x_i − x_j|`.
   > `spec/vectors/VECTORS-1.0.json` asserts the invariance on all eight vectors
   > rather than leaving it as an argument.
3. Project the **whole cropped signal** — baseline and response together — onto the
   basis: `S = Vᵀ Y`. Row `c` of **S** is component `c` over all cropped samples. The
   basis comes from the response alone; the projection covers everything.
4. Explained variance of component `c` is `100 · s_c² / Σ_j s_j²`.
5. Retain the leading `m` components, where `m` is the smallest count whose cumulative
   explained variance reaches `max_var` percent, computed as
   `m = n_ch − #{c : cumvar_c ≥ max_var} + 1`, and `m = n_ch` when `max_var = 100`.
6. For each retained component compute

   > **snr_c** = √( mean over the response window of `S_c²` ÷ mean over the baseline
   > window of `S_c²` )

   and keep the components with `snr_c` **strictly greater than** `min_snr`.
   **The parameter `k` does not appear in this ratio.** It enters at step 4.
7. If no component survives, PCIst is **0**.

## Step 3 — Distance matrices

For each surviving component `c`, and separately for its baseline block (length `n_b`)
and its response block (length `n_r`):

> **D[i, j] = | x_i − x_j |**

the absolute difference between the component's amplitude at samples `i` and `j`. The
diagonal is zero and both matrices are symmetric. Where `embed` is true the component is
first time-delay embedded with `L` dimensions and delay `tau`, samples `0 … (L−1)·tau`
are dropped from the front, and `D` becomes the Euclidean norm over the embedding
dimensions; `embed` is false in every parameter set §2.1 currently defines.

## Step 4 — State transitions

This is the step no published source specifies, and the step where an implementer will
otherwise go wrong. **A state transition is not a thresholded distance.** It is a change
in the recurrence relation between one sample and the next.

For a threshold `ε`:

1. Form the **recurrence mask** `R = [D ≤ ε]`, taking the value 1 where the two samples
   are within `ε` of each other and 0 otherwise, inclusive of the diagonal, which is
   therefore always 1.

   > The direction of the inequality does not matter, and an earlier draft of this
   > annex implied it did. If `R' = 1 − R` then `diff(R') = −diff(R)`, and step 2 takes
   > the absolute value of that difference, so complementing the mask cannot change the
   > transition count. `≤` is specified for definiteness, not because `≥` would be
   > wrong. The vector set asserts this invariance too.
2. Form the **transition matrix** by differencing the mask along its second index and
   taking absolute values:

   > **T[i, j] = | R[i, j+1] − R[i, j] |**   for `j = 0 … n−2`,
   > with a final column of zeros appended so that **T** is `n × n`.

   `T[i, j]` is 1 exactly where sample `i` changes between being inside and outside
   `ε` of sample `j` as `j` advances by one. Differencing is **one-sided**; the matrix
   is not symmetrized.
3. The number of state transitions is the total, normalized by the squared window
   length:

   > **NST(ε) = ( Σ_i Σ_j T[i, j] ) ÷ n²**

   The sum runs over the full `n × n` matrix, including the diagonal and the appended
   zero column. Normalizing by `n²` rather than by the count of off-diagonal pairs is a
   choice and it changes the value.

The threshold is swept over `n_steps` values spaced **linearly** from

> **ε_min = the median of every entry of the baseline distance matrix**, including its
> zero diagonal, to
> **ε_max = the maximum entry of the response distance matrix**, multiplied by
> `max_thr_p` (default 1.0).

The range is set by the baseline's median and the response's maximum. Using the pooled
distances, or the minimum rather than the median, will not reproduce the reference.

## Step 5 — Per-component complexity and the index

For each surviving component, over the swept thresholds:

> **ΔNST_c = n_r · max_ε [ NST_r(ε) − k · NST_b(ε) ]**,  and **0** where that is negative.

`k` is the **noise-control multiplier on the baseline count**, not a term in the SNR and
not a weight on the response. The maximum is taken over thresholds before multiplying by
the response length, and `n_r` is the only length factor applied — the baseline count is
not separately weighted by `n_b`.

Finally

> **PCIst = Σ_c ΔNST_c**

summed over the surviving components.

## Conformance

An implementation conforms to PCIST-1.0 where, on the eight vectors of
`spec/vectors/VECTORS-1.0.json`, it reproduces every expected value to a relative
tolerance of **1 × 10⁻⁶**. The vectors ship as generators — a seed and a few shape
parameters — so the input is reproducible in any language with a standard normal
generator and an FFT, with nothing to download and nothing to drift.

**The set is diagnostic, not merely pass or fail.** Each vector records what eight
specific mis-implementations return on it, one broken clause at a time, four of them
errors an implementation written from the pre-specification standard actually made. An
implementation that fails is matched against that table and told which clause it most
likely got wrong, by section. Every mutation is distinguished by at least one vector,
and the build refuses to publish a set with a blind spot — which is how the two
immaterial choices noted above were found. A reading declares the specification version it was computed under in
`computation.spec_version`, alongside the library, version and commit that §2.1 already
requires.

## What is still not fixed here

This annex specifies the computation. It does not settle the parameter sets, which are
per-paradigm and versioned separately in §2.1, nor the `resample` figure for
`mouse_spes_v1`, which §2.1 still marks provisional. It also does not make the computed
value comparable across implementations of *different* specifications: §2.2's finding
stands, and a value remains meaningful by comparison rather than in itself.