CRP

CRP#

brainmaze_eeg.crp.ccep_proj(V)#

Semi-normalized pairwise cross-projections between trials.

Each trial is L2-normalized and projected onto every other (raw) trial; the self-projections (diagonal) are discarded. For K trials this returns the K^2 - K off-diagonal magnitudes, flattened in row-major order so that element i * (K - 1) + (j or j - 1) is the projection of the normalized trial i onto the raw trial j.

Parameters:

V (numpy.ndarray) – Trial matrix of shape (n_samples, K) (each column is a trial).

Returns:

Off-diagonal cross-projection magnitudes, shape (K^2 - K,).

Return type:

numpy.ndarray

brainmaze_eeg.crp.crp_method(v, t_win, prune_the_data=False, rejection_coeff=3.5, min_trials=5)#

Canonical Response Parameterization (CRP) of single-pulse stimulation responses.

Python implementation of the method described in:

Miller, K. J., Müller, K.-R., Valencia, G. O., Huang, H., Gregg, N. M., Worrell, G. A., & Hermes, D. (2023). Canonical Response Parameterization: Quantifying the structure of responses to single-pulse intracranial electrical brain stimulation. PLoS Computational Biology, 19(5), e1011105. https://doi.org/10.1371/journal.pcbi.1011105

The method takes a set of K single-trial voltage traces of an evoked response (e.g. a cortico-cortical evoked potential, CCEP), and:

  1. Builds a temporal profile of the semi-normalized pairwise cross-projection magnitude S(t) and takes its peak as the response duration tau_R.

  2. Extracts a canonical response shape C(t) over [0, tau_R] via linear kernel PCA (the 1st principal component of the truncated trial matrix).

  3. Parameterizes each trial k as V_k(t) = alpha_k * C(t) + eps_k(t), yielding per-trial projection weight alpha_k, residual eps_k and derived quantities (SNR, explained variance).

  4. Optionally rejects artefactual trials as outliers of the per-trial sub-distributions of S(tau_R) (Miller et al. 2023, “extraction significance” section).

Parameters:

  • v (np.ndarray) – Voltage matrix of shape (n_samples, n_trials). Each column is one trial; trials must be baseline-corrected (the method is sensitive to baseline offsets).

  • t_win (np.ndarray) – Time vector of shape (n_samples,) in seconds, aligned with the rows of v.

  • prune_the_data (bool, optional) – If True, run one pass of artefact-trial rejection and re-parameterize on the surviving trials. Default False.

  • rejection_coeff (float, optional) – Robust (median/MAD) z-score threshold for the two-sided artefact-trial rejection. Larger is more permissive. Default 3.5.

  • min_trials (int, optional) – Minimum number of trials required to attempt rejection, and minimum number that must remain after rejection. Default 5.

Returns:

  • crp_parameters (dict) –

    • V_tR : truncated voltage matrix (T_R, K) over [0, tau_R].

    • C : canonical response shape C(t), unit-norm vector of length T_R. Sign-normalized so the response is positively represented across trials.

    • al : projection weights alpha_k (length K), in the same units as v times sqrt(#samples).

    • al_p : alpha_k normalized by sqrt(T_R) (paper’s alpha'_k, uV).

    • ep : residual eps_k(t) after removing alpha_k * C(t), shape (T_R, K).

    • erp : fitted evoked response alpha_k * C(t), shape (T_R, K) – the part of each trial captured by the canonical shape (over the response window).

    • erp_full : erp embedded in a full-length (n_samples, K) array, zero after tau_R. Subtract from v to remove the evoked response across the whole epoch.

    • Vsnr : per-trial signal-to-noise magnitude |alpha_k| / ||eps_k||.

    • expl_var : per-trial fraction of variance explained by C(t).

    • tR : response duration tau_R in seconds.

    • epep_root : sqrt(diag(ep.T @ ep)), per-trial residual norm.

    • avg_trace_tR / std_trace_tR : mean / std trace over [0, tau_R].

    • parms_times : time vector over [0, tau_R].

    • kept_trials : indices (into the original trial axis) of retained trials.

    • rejected_trials : indices (into the original trial axis) of rejected trials.

    • rejection_stat : per-original-trial anomaly statistic used for rejection (mean of the trial’s S(tau_R) sub-distribution; NaN if not computed).

  • crp_projections (dict) –

    • proj_tpts : profile time points (seconds).

    • s_all : all cross-projection magnitudes, shape (K^2 - K, n_tpts).

    • mean_proj_profile / var_proj_profile : mean / variance of S vs time.

    • tR_index : index into proj_tpts of the response duration.

    • avg_trace_input / std_trace_input : mean / std of the full input traces.

    • stat_indices : the non-overlapping half-selection of s_all rows used for significance and rejection.

    • t_value_tR / p_value_tR : extraction significance at tau_R.

    • t_value_full / p_value_full : extraction significance at the full window.

Notes

The sign of C (and hence al / al_p) is fixed so that the projection weights are predominantly positive; the product alpha_k * C(t) (erp) is sign-invariant regardless.

Raises:

ValueError – If v has 10 or fewer samples (the projection profile cannot be formed).

Example

Remove the evoked response from every trial across the full epoch, and inspect which trials were rejected as artefacts:

params, proj = crp_method(v, t_win, prune_the_data=True)
clean = v[:, params['kept_trials']] - params['erp_full']
print('rejected trials:', params['rejected_trials'])
brainmaze_eeg.crp.get_stat_indices(N)#

Non-overlapping half-selection of the cross-projection magnitudes for statistics.

ccep_proj returns both orientations of every trial pair (i normalized onto j, and j normalized onto i). Using both double-counts each interaction and inflates significance (Miller et al. 2023). This function selects one orientation per unordered pair – exactly N * (N - 1) / 2 of the N^2 - N projections – balancing which trial plays the normalized role so that each trial is the normalized one in as close to half of its pairs as possible (exactly half when N is odd).

Parameters:

N – Number of trials.

Returns:

Sorted flat indices into ccep_proj’s output, shape (N * (N - 1) / 2,).

brainmaze_eeg.crp.kt_pca(X)#

Linear kernel PCA (“kernel trick”) of X.

Implements the trick from Schölkopf et al., ICANN 1998, needed when the number of timepoints T greatly exceeds the number of trials N: the eigenvectors of the T x T matrix X @ X.T are recovered from the eigendecomposition of the much smaller N x N matrix X.T @ X.

X.T @ X is symmetric positive-semidefinite, so np.linalg.eigh is used: it returns real, ascending eigenvalues and orthonormal eigenvectors (no complex round-off, no manual re-orthogonalization).

Parameters:

X – Data matrix of shape (T, N).

Returns:

(E, S) – eigenvectors of X @ X.T (columns) and S, the singular values of X (square roots of the eigenvalues of X.T @ X), both in descending order.