feat: manifold diagnostic toolkit — dimension selection, CV reconstruction/decoding, Procrustes alignment, within-session stability

+diagnostics/+compute/align_procrustes.m

@@ -0,0 +1,88 @@
+function results = align_procrustes(Z1, Z2, allowScale)
+%ALIGN_PROCRUSTES Orthogonal Procrustes alignment of two latent spaces.
+%
+%   RESULTS = diagnostics.compute.align_procrustes(Z1, Z2) finds the
+%   orthogonal matrix R that minimises ||Z1 - Z2*R||_F, and returns the
+%   aligned Z2 and related metrics.
+%
+%   RESULTS = diagnostics.compute.align_procrustes(Z1, Z2, ALLOWSCALE)
+%   also optimises an isotropic scaling factor when ALLOWSCALE is true.
+%
+%   Both matrices are centred (mean-subtracted) before alignment.
+%
+%   Inputs
+%   ------
+%   Z1         : T x d latent space (reference session, e.g. trials x dims).
+%   Z2         : T x d latent space (session to be aligned). Must have the
+%                same number of columns d as Z1.
+%   allowScale : logical (default false). Allow isotropic scaling.
+%
+%   Outputs
+%   -------
+%   results : struct with fields
+%       .Z2_aligned  - (T x d) Z2 after applying R (and scale if enabled).
+%       .R           - (d x d) orthogonal rotation matrix.
+%       .scale       - isotropic scale factor (1.0 when allowScale=false).
+%       .disparity   - Procrustes disparity = ||Z1 - Z2_aligned||_F^2 /
+%                      ||Z1||_F^2 (normalised; 0 = perfect alignment).
+%
+%   Notes
+%   -----
+%   * Implements the classic SVD solution: [U,~,V] = svd(Z1'*Z2);
+%     R = V * U'.
+%   * No toolboxes required.
+%   * For multi-session alignment see NeuralEmbedding.alignSessions.
+%
+%   Example
+%   -------
+%   res = diagnostics.compute.align_procrustes(Z1, Z2);
+%   fprintf('Disparity = %.4f\n', res.disparity);
+%
+%   See also diagnostics.compute.alignment_metrics,
+%            diagnostics.pars.ProcrustesAlignment
+
+if nargin < 3
+    allowScale = false;
+end
+
+% --- Validation ---
+if ~ismatrix(Z1) || ~isnumeric(Z1) || ~ismatrix(Z2) || ~isnumeric(Z2)
+    error('diagnostics:align_procrustes:badInput', ...
+        'Z1 and Z2 must be 2-D numeric matrices.');
+end
+if size(Z1, 2) ~= size(Z2, 2)
+    error('diagnostics:align_procrustes:dimMismatch', ...
+        'Z1 and Z2 must have the same number of columns (latent dims).');
+end
+
+% --- Centre ---
+Z1c = Z1 - mean(Z1, 1);
+Z2c = Z2 - mean(Z2, 1);
+
+% --- SVD solution ---
+M       = Z1c' * Z2c;             % d x d
+[U, S, V] = svd(M);
+R       = V * U';                  % orthogonal rotation
+
+% --- Optional scale ---
+if allowScale
+    scale = trace(S) / (norm(Z2c, 'fro')^2);
+else
+    scale = 1.0;
+end
+
+Z2_aligned = scale * Z2c * R;
+
+% --- Disparity ---
+norm_Z1 = norm(Z1c, 'fro');
+if norm_Z1 == 0
+    disparity = 0;
+else
+    disparity = norm(Z1c - Z2_aligned, 'fro')^2 / norm_Z1^2;
+end
+
+results.Z2_aligned = Z2_aligned;
+results.R          = R;
+results.scale      = scale;
+results.disparity  = disparity;
+end

+diagnostics/+compute/alignment_metrics.m

@@ -0,0 +1,114 @@
+function results = alignment_metrics(Z1, Z2_aligned)
+%ALIGNMENT_METRICS Quantitative metrics for latent-space alignment quality.
+%
+%   RESULTS = diagnostics.compute.alignment_metrics(Z1, Z2_ALIGNED)
+%   computes several metrics that quantify how well two latent spaces are
+%   aligned, assuming Z2_ALIGNED has already been transformed by
+%   diagnostics.compute.align_procrustes.
+%
+%   Inputs
+%   ------
+%   Z1         : T x d reference latent space (centred).
+%   Z2_aligned : T x d aligned latent space.
+%
+%   Outputs
+%   -------
+%   results : struct with fields
+%       .disparity       - Procrustes disparity ||Z1-Z2||_F^2 / ||Z1||_F^2.
+%       .principalAngles - (1 x d) principal angles (radians) between the
+%                          column spaces of Z1 and Z2_aligned, in ascending
+%                          order.
+%       .meanPrincipalAngle - mean of principalAngles (radians).
+%       .distCorr        - Pearson correlation between vectorised pairwise
+%                          Euclidean distance matrices of Z1 and Z2_aligned
+%                          (Mantel-style; 1 = identical structure).
+%
+%   Notes
+%   -----
+%   * Principal angles: let Q1, Q2 be column-space orthonormal bases (from
+%     QR decomposition). Angles = acos(svd(Q1'*Q2)).
+%   * Distance correlation is computed on the upper-triangular part of the
+%     pairwise distance matrix to avoid double-counting.
+%   * For large T, distance-correlation computation is O(T^2); set
+%     distCorr to NaN if T > 2000 to avoid excessive runtime.
+%   * No toolboxes required.
+%
+%   Example
+%   -------
+%   proc = diagnostics.compute.align_procrustes(Z1, Z2);
+%   met  = diagnostics.compute.alignment_metrics(Z1, proc.Z2_aligned);
+%   fprintf('Mean principal angle = %.2f deg\n', ...
+%       rad2deg(met.meanPrincipalAngle));
+%
+%   See also diagnostics.compute.align_procrustes
+
+% --- Validation ---
+if size(Z1, 2) ~= size(Z2_aligned, 2)
+    error('diagnostics:alignment_metrics:dimMismatch', ...
+        'Z1 and Z2_aligned must have the same number of columns.');
+end
+
+% Centre
+Z1 = Z1 - mean(Z1, 1);
+Z2 = Z2_aligned - mean(Z2_aligned, 1);
+
+d = size(Z1, 2);
+T = size(Z1, 1);
+
+% --- Disparity ---
+nZ1 = norm(Z1, 'fro');
+if nZ1 == 0
+    disparity = 0;
+else
+    disparity = norm(Z1 - Z2, 'fro')^2 / nZ1^2;
+end
+
+% --- Principal angles ---
+[Q1, ~] = qr(Z1, 0);   % economy QR
+[Q2, ~] = qr(Z2, 0);
+svs     = svd(Q1' * Q2);
+svs     = min(max(svs, -1), 1);   % clamp for numerical safety
+principalAngles = acos(svs(:)');
+
+% --- Distance correlation (Mantel-style) ---
+if T <= 2000
+    D1 = i_pdist_upper(Z1);
+    D2 = i_pdist_upper(Z2);
+    distCorr = i_pearson(D1, D2);
+else
+    distCorr = NaN;  % too expensive for large T
+end
+
+results.disparity           = disparity;
+results.principalAngles     = principalAngles;
+results.meanPrincipalAngle  = mean(principalAngles);
+results.distCorr            = distCorr;
+end
+
+% =========================================================================
+%  Local helpers
+% =========================================================================
+
+function v = i_pdist_upper(Z)
+% Upper-triangular pairwise Euclidean distances as a vector.
+T = size(Z, 1);
+v = zeros(1, T*(T-1)/2);
+k = 0;
+for ii = 1:T-1
+    diff = Z(ii+1:end, :) - Z(ii, :);
+    n    = T - ii;
+    v(k+1 : k+n) = sqrt(sum(diff.^2, 2));
+    k = k + n;
+end
+end
+
+function r = i_pearson(a, b)
+a = a(:); b = b(:);
+a = a - mean(a); b = b - mean(b);
+na = norm(a); nb = norm(b);
+if na == 0 || nb == 0
+    r = 0;
+    return;
+end
+r = (a' * b) / (na * nb);
+end

+diagnostics/+compute/circular_shift.m

@@ -0,0 +1,36 @@
+function Xout = circular_shift(X)
+%CIRCULAR_SHIFT Circular time-shift null for neural time-series data.
+%
+%   XOUT = diagnostics.compute.circular_shift(X) returns a copy of X where
+%   each column (neuron) has been independently shifted by a uniformly
+%   random time offset in 1 .. T-1 (wrapping at the boundaries).
+%
+%   This shuffle:
+%     - Preserves each neuron's autocorrelation structure.
+%     - Preserves each neuron's marginal distribution.
+%     - Destroys the alignment between the neural data and any external
+%       labels (behaviour, task variables).
+%     - Partially preserves instantaneous cross-neuron correlations at
+%       non-zero lags (unlike neuronwise shuffle).
+%
+%   Use this as a temporal null for decoding permutation tests when the
+%   data are sampled continuously in time (not independent trials).
+%
+%   Input
+%   -----
+%   X : T x N matrix, where rows are time-ordered samples.
+%
+%   Output
+%   ------
+%   Xout : T x N matrix with independently circularly shifted columns.
+%
+%   See also diagnostics.compute.shuffle_neuronwise,
+%            diagnostics.shufflers.circular_shift
+
+[T, N] = size(X);
+Xout   = zeros(T, N);
+shifts = randi([1, T-1], 1, N);   % independent shift per neuron
+for nn = 1:N
+    Xout(:, nn) = circshift(X(:, nn), shifts(nn));
+end
+end