PCIST-1.0 annex

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.