MATLAB: full expanded feature matrix, color-coded viz, and EFA

The existing readers kept scalar channels only. Add tooling to load ALL
~2,768 feature variables (every vector channel expanded per component),
visualize them color-coded by category, and run exploratory factor analysis.

New matlab/ functions:
- featureInfo            label/category table for every expanded variable
                         (class, subclass, level, model, per-class color).
- featuresToTable        one clip -> wide timetable, vectors expanded.
- readAnnotationCorpusFull  concatenate the whole corpus -> C.X + C.info + C.T.
- plotFeatureMatrix      heatmap / class-mean view, color-coded by category.
- factorAnalysisCorpus   EFA (factoran, varimax) on interpretable features;
                         pairwise-correlation fit handles audio-only NaNs;
                         returns loadings + per-timepoint factor scores + plots.

Docs: CONTENTS.md MATLAB table + recipe; walkthrough.m section 8. Verified
end-to-end on the 83-clip corpus (28,237 timepoints) via the MATLAB engine.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>

Author: torwager

docs/CONTENTS.md

@@ -145,10 +145,15 @@ Details: [`../analysis/web/README.md`](../analysis/web/README.md).
 | `readAnnotations(path)` | load one `.h5`/folder/JSON → struct (stimulus, time_sec, features) |
 | `getFeature(ann, "audio/low_level/mfcc")` | one channel + metadata by hierarchical path |
 | `featuresToTimetable(ann)` | scalar channels → timetable on the common grid |
-| `readAnnotationCorpus(folder)` | stack the whole corpus → `C.X [timepoints × channels]` |
+| `readAnnotationCorpus(folder)` | stack the whole corpus → `C.X [timepoints × channels]` (**scalar** channels only) |
 | `annotationMovieViewer(movie, ann)` | play movie with synced annotation time series + marker |
 | `analyzeCorpus(C)` | correlation heatmap, PCA, channel + class network graphs |
 | `selectStimulusSet(C)` | D-optimal high-variance / low-redundancy segment selection |
+| `featureInfo()` | label/category table for **all** expanded variables (class, subclass, level, model, color) |
+| `featuresToTable(ann)` | one clip → wide table with **every vector expanded** into per-component columns (~2.7k vars) |
+| `readAnnotationCorpusFull(folder)` | stack the whole corpus with vectors expanded → `C.X [timepoints × ~2.7k vars]` + `C.info` |
+| `plotFeatureMatrix(C)` | heatmap of the full feature time series, color-coded by category |
+| `factorAnalysisCorpus(C)` | exploratory factor analysis (EFA / `factoran`) + color-coded loadings plot |
 
 ---
 
@@ -167,14 +172,28 @@ mf  = getFeature(ann, "audio/low_level/mfcc");   % a vector channel [n × 13]
 m = "data/movies/spacetop/videos/ses-01/ses-01_run-01_order-04_content-parkour.mp4";
 annotationMovieViewer(m, "annotations/corpus/ses-01_run-01_order-04_content-parkour")
 
-% 3) Analyze the whole CORPUS
+% 3) Analyze the whole CORPUS (scalar channels — correlation / PCA / design tool)
 C   = readAnnotationCorpus("annotations/corpus");
 res = analyzeCorpus(C);                   % figures → analysis/figures/
 sel = selectStimulusSet(C, "K", 20);      % design a stimulus set; sel.table
 
+% 4) Load the FULL expanded feature set (~2.7k variables: every vector expanded),
+%    visualize it color-coded by category, and extract factors with EFA
+F   = readAnnotationCorpusFull("annotations/corpus");   % F.X [timepoints × ~2768], F.info labels
+plotFeatureMatrix(F, "Clip", "BigBuckBunny");           % heatmap, color-coded by class
+fa  = factorAnalysisCorpus(F, "NumFactors", 10);        % EFA + loadings plot; fa.scores = per-timepoint factors
+
 % Python inspection (alternative): h5py / pandas on <id>.h5 and corpus_index.csv
 ```
 
+> **Two corpus readers.** `readAnnotationCorpus` returns only the **scalar** channels
+> (what `analyzeCorpus`/`selectStimulusSet` expect). `readAnnotationCorpusFull` expands
+> every **vector** channel (SigLIP/DINOv2/CLAP embeddings, AudioSet/action posteriors,
+> EmoNet, MFCC, …) into one column per component — the full ~2,768-variable matrix with a
+> companion `featureInfo` label table (class / subclass / level / model / color) for
+> color-coded plotting and factor analysis. See the
+> [feature map](FEATURE_MAP.md) for how the variables are organized.
+
 ```bash
 # 4) SEARCH segments by feature in the browser (serve from project root)
 python3 tools/serve.py           # then open http://localhost:8000/analysis/web/index.html

docs/walkthrough.m

@@ -78,7 +78,28 @@
 disp(sel.table)                          % rank, stimulus, t_start, t_end, dur_s
 % sel.objTrace vs sel.randTrace quantifies the gain over random selection.
 
-%% 8. Where to go next
+%% 8. FULL feature set: all ~2,768 variables, color-coded, + factor analysis
+% readAnnotationCorpus (above) keeps only scalar channels. To get EVERY variable —
+% each vector channel (SigLIP/DINOv2/CLAP embeddings, AudioSet/action posteriors,
+% EmoNet, MFCC, ...) expanded into one column per component — use the FULL reader.
+F = readAnnotationCorpusFull("annotations/corpus");    % F.X [timepoints x ~2768]
+% F.info is the label table: one row per column, with class/subclass/level/model/color.
+fprintf("full set: %d variables (%d embedding, %d interpretable) over %d timepoints\n", ...
+        size(F.X,2), sum(F.info.IsEmbedding), sum(~F.info.IsEmbedding), size(F.X,1));
+
+% Visualize the whole feature matrix for one clip, color-coded by category:
+plotFeatureMatrix(F, "Clip", "BigBuckBunny");          % heatmap; class color strip on the left
+% plotFeatureMatrix(F, "Mode","classmean");            % or six class-mean trajectories
+
+% Exploratory factor analysis (interpretable features by default; excludes the opaque
+% embeddings and the big AudioSet/Kinetics taxonomies). Returns rotated loadings +
+% a per-timepoint factor-score time series you can plot like any feature.
+fa = factorAnalysisCorpus(F, "NumFactors", 10);        % draws loadings heatmap + scree
+r  = F.stim == "BigBuckBunny";
+figure; plot(F.time_sec(r), fa.scores(r, 1:3)); legend("F1","F2","F3");
+xlabel("time (s)"); ylabel("factor score"); title("Factor time series — BigBuckBunny");
+
+%% 9. Where to go next
 % - Inspect any other stimulus: change `stimId` in section 0.
 % - Annotate a NEW movie (Python):
 %     PYTHONPATH=src .venv/bin/python -m nfe.run <movie> --vision --audio-hl --events \

matlab/factorAnalysisCorpus.m

@@ -0,0 +1,180 @@
+function F = factorAnalysisCorpus(C, opts)
+%FACTORANALYSISCORPUS Exploratory factor analysis of the feature corpus.
+%   F = FACTORANALYSISCORPUS(C) runs EFA (MATLAB's FACTORAN, maximum-likelihood with
+%   varimax rotation) on the expanded feature matrix from READANNOTATIONCORPUSFULL
+%   and returns the factor solution plus a color-coded loadings plot.
+%
+%   By default it analyzes the INTERPRETABLE feature set — it excludes the opaque
+%   SigLIP/DINOv2/CLAP embeddings and the two large fixed-taxonomy posteriors
+%   (AudioSet 527 tags, Kinetics 400 actions), which otherwise dominate the solution
+%   with hundreds of near-collinear columns. Set the options below to include them.
+%
+%   Missing data: features are z-scored, then the factor MODEL is fit from the
+%   pairwise correlation matrix (using all available data per pair, so audio-only
+%   clips still contribute their audio/language features) after a nearest
+%   positive-definite adjustment. Per-timepoint FACTOR SCORES are then computed by
+%   the regression method with mean-imputation (standardized NaN -> 0), so every
+%   timepoint gets a score and the factors can be plotted as time series like the
+%   raw features.
+%
+%   Name-value options:
+%       "NumFactors"        (8)     number of factors to extract
+%       "ExcludeEmbeddings" (true)  drop SigLIP/DINOv2/CLAP dimensions
+%       "ExcludeTaxonomies" (true)  drop audioset_tags (527) + action_posteriors (400)
+%       "Rotate"            ("varimax")  any FACTORAN rotation ("none","promax",...)
+%       "MaxNaNFrac"        (0.8)   drop variables missing in > this fraction of rows
+%                                   (0.8 keeps visual features, which are absent on the
+%                                   ~67% of timepoints that come from audio-only clips)
+%       "Plot"              (true)  draw the loadings heatmap + scree
+%       "SavePng"           ("")    export the loadings figure to this path
+%
+%   Returns struct F with:
+%       F.loadings   [p x m] rotated factor loadings
+%       F.vars       [p x k] FEATUREINFO subset for the retained variables (row i <-> loadings row i)
+%       F.scores     [T x m] per-timepoint factor scores (regression method)
+%       F.specificVar[p x 1] uniquenesses; F.stats FACTORAN stats struct
+%       F.eigs       eigenvalues of the correlation matrix (for the scree plot)
+%       F.varExplained [1 x m] % variance each rotated factor accounts for
+%       F.m, F.rotation, F.nEff
+%
+%   Example:
+%       C = readAnnotationCorpusFull("annotations/corpus");
+%       F = factorAnalysisCorpus(C, "NumFactors", 8);
+%       % plot factor-3 time series for one clip:
+%       r = C.stim=="BigBuckBunny"; plot(C.time_sec(r), F.scores(r,3));
+%
+%   See also READANNOTATIONCORPUSFULL, PLOTFEATUREMATRIX, FACTORAN, PCA.
+
+arguments
+    C (1,1) struct
+    opts.NumFactors (1,1) double {mustBePositive, mustBeInteger} = 8
+    opts.ExcludeEmbeddings (1,1) logical = true
+    opts.ExcludeTaxonomies (1,1) logical = true
+    opts.Rotate (1,1) string = "varimax"
+    opts.MaxNaNFrac (1,1) double = 0.8
+    opts.Plot (1,1) logical = true
+    opts.SavePng (1,1) string = ""
+end
+
+info = C.info; X = C.X;
+
+% ---- column selection ----
+keep = info.Numeric;
+if opts.ExcludeEmbeddings, keep = keep & ~info.IsEmbedding; end
+if opts.ExcludeTaxonomies
+    keep = keep & ~ismember(info.Leaf, ["audioset_tags", "action_posteriors"]);
+end
+keep = keep & (mean(isnan(X),1)' <= opts.MaxNaNFrac);      % not-too-missing
+X = X(:, keep); info = info(keep, :);
+
+% z-score; drop zero-variance columns
+mu = mean(X,1,"omitnan"); sd = std(X,0,1,"omitnan");
+good = sd > 0 & isfinite(sd);
+X = X(:, good); info = info(good,:); mu = mu(good); sd = sd(good);
+Z = (X - mu) ./ sd;
+
+% ---- correlation matrix (pairwise) + nearest-PD repair ----
+R = corr(Z, "rows", "pairwise");
+R(~isfinite(R)) = 0; R(1:size(R,1)+1:end) = 1;
+% drop columns still fully collinear (|r|=1 with an earlier column)
+[R, info, Z] = dropCollinear(R, info, Z);
+R = nearestPD(R);
+p = size(R,1);
+nEff = size(Z,1);
+
+m = min(opts.NumFactors, p-1);
+if m < opts.NumFactors
+    warning("factorAnalysisCorpus:m", "Reduced NumFactors to %d (only %d variables).", m, p);
+end
+
+% ---- fit factor model ----
+% factoran needs complete data; audio-only clips make visual columns NaN, so when
+% any NaNs are present we fit from the (repaired) correlation matrix instead.
+if any(isnan(Z(:)))
+    [L, psi, ~, stats] = factoran(R, m, "Xtype","covariance", "Nobs",nEff, ...
+                                  "Rotate",char(opts.Rotate));
+else
+    [L, psi, ~, stats] = factoran(Z, m, "Xtype","data", "Rotate",char(opts.Rotate), ...
+                                  "Maxit",1000);
+end
+
+% ---- per-timepoint factor scores (regression method, mean-imputed) ----
+Zi = Z; Zi(isnan(Zi)) = 0;                 % standardized mean-imputation
+W = R \ L;                                 % regression weights  (p x m)
+scores = Zi * W;
+
+% variance explained by each rotated factor (sum of squared loadings / p)
+ve = 100 * sum(L.^2, 1) / p;
+[ve, si] = sort(ve, "descend"); L = L(:, si); scores = scores(:, si);
+
+F.loadings = L; F.vars = info; F.scores = scores; F.specificVar = psi;
+F.stats = stats; F.varExplained = ve; F.m = m; F.rotation = string(opts.Rotate); F.nEff = nEff;
+F.eigs = sort(eig(R), "descend");
+
+fprintf("EFA: %d variables, %d factors (%s), n_eff=%d. Variance explained: %s%%\n", ...
+    p, m, opts.Rotate, nEff, join(compose("%.1f", ve), "/"));
+
+if opts.Plot, plotLoadings(F, opts.SavePng); end
+end
+
+% -------------------------------------------------------------------------
+function [R, info, Z] = dropCollinear(R, info, Z)
+p = size(R,1); keep = true(p,1);
+for i = 2:p
+    if any(abs(R(i, 1:i-1)) > 0.999 & keep(1:i-1)'), keep(i) = false; end
+end
+R = R(keep,keep); info = info(keep,:); Z = Z(:,keep);
+end
+
+function A = nearestPD(A)
+A = (A + A')/2;
+[V,D] = eig(A); d = diag(D);
+d(d < 1e-6) = 1e-6;
+A = V*diag(d)*V';
+A = (A + A')/2;
+s = sqrt(diag(A)); A = A ./ (s*s');       % rescale to unit diagonal (correlation)
+A = (A + A')/2;
+A(1:size(A,1)+1:end) = 1;                 % exact unit diagonal
+end
+
+function plotLoadings(F, savePng)
+fig = figure("Color","w","Position",[100 100 980 760]);
+% sort variables by class then by dominant factor for a readable block structure
+classes = ["visual","audio","language","social","situation","affect"];
+key = double(arrayfun(@(c) find(classes==c,1), F.vars.Class));
+[~, dom] = max(abs(F.loadings), [], 2);
+[~, ord] = sortrows([key, dom]);
+L = F.loadings(ord,:); vinfo = F.vars(ord,:);
+
+ax = axes(fig, "Position",[0.30 0.08 0.60 0.84]);
+imagesc(ax, L); set(ax, "CLim",[-1 1]);
+colormap(ax, localDiv()); cb = colorbar(ax); cb.Label.String = "loading";
+ax.XTick = 1:F.m; ax.XTickLabel = compose("F%d", 1:F.m);
+xlabel(ax, sprintf("factors (%s rotation)", F.rotation));
+title(ax, sprintf("EFA loadings — %d variables x %d factors", size(L,1), F.m));
+% label a subset of rows (too many to show all)
+step = max(1, round(size(L,1)/60));
+ax.YTick = 1:step:size(L,1);
+lab = vinfo.Leaf + string(compose(" [%s]", extractBefore(vinfo.Level+"   ",4)));
+comp = vinfo.Component; hasC = comp~="";
+lab(hasC) = vinfo.Leaf(hasC) + ":" + comp(hasC);
+ax.YTickLabel = lab(1:step:end); ax.FontSize = 7;
+
+% class color strip
+strip = axes(fig, "Position",[0.27 0.08 0.02 0.84]);
+image(strip, reshape(vinfo.Color,[height(vinfo) 1 3]));
+strip.XTick=[]; strip.YTick=[]; box(strip,"on"); linkaxes([ax strip],"y");
+
+% scree inset
+sc = axes(fig, "Position",[0.06 0.62 0.16 0.30]);
+plot(sc, F.eigs(1:min(20,numel(F.eigs))), "-o","MarkerSize",3); grid(sc,"on");
+title(sc,"scree"); xlabel(sc,"factor"); ylabel(sc,"eigenvalue"); sc.FontSize=7;
+
+if savePng ~= ""
+    exportgraphics(fig, savePng, "Resolution",150); fprintf("wrote %s\n", savePng);
+end
+end
+
+function m = localDiv()
+n=256; x=linspace(0,1,n)'; m=[min(1,1.8*x), 1-abs(2*x-1), min(1,1.8*(1-x))];
+end

matlab/featureInfo.m

@@ -0,0 +1,139 @@
+function info = featureInfo(templatePath)
+%FEATUREINFO Metadata/label table for every expanded feature variable.
+%   INFO = FEATUREINFO() reads schema/channel_template.json and returns one row
+%   per *variable* in the fully-expanded feature set: every scalar channel is one
+%   row, and every vector channel is expanded into one row per component (so the
+%   768-D SigLIP embedding contributes 768 rows, EmoNet 20 rows, etc.). This is the
+%   companion label table for the wide feature matrix built by FEATURESTOTABLE /
+%   READANNOTATIONCORPUSFULL — INFO row i describes column i of that matrix.
+%
+%   INFO = FEATUREINFO(TEMPLATEPATH) uses a specific channel template.
+%
+%   Columns of INFO:
+%       VarName     valid, unique MATLAB name used as the table column name
+%       Path        hierarchical channel path (e.g. "affect/depicted/emonet")
+%       Class       top-level class (visual/audio/language/social/situation/affect)
+%       Subclass    second level (e.g. "high_level_static"), or "(direct)"
+%       Leaf        channel leaf name (e.g. "emonet")
+%       Component   component label for a vector element (e.g. "Fear"), else ""
+%       CompIndex   1-based index within the channel (1 for scalars)
+%       Model       model/tool that produced it
+%       Dtype       scalar | bool | event | vector | categorical | label
+%       Level       "low" | "mid" | "high"  (perceptual-abstraction level)
+%       Numeric     true if it enters the numeric matrix (scalar/bool/event/vector)
+%       IsEmbedding true for the opaque SigLIP/DINOv2/CLAP embedding dimensions
+%       Color       1x3 RGB for the class (for color-coded plots)
+%
+%   Only Numeric==true rows appear in the numeric matrix X; categorical "*_top"
+%   codes and free-text label channels (asr_text, scene_description, ...) are listed
+%   here for reference but excluded from X (their codes are not magnitudes).
+%
+%   See also FEATURESTOTABLE, READANNOTATIONCORPUSFULL, PLOTFEATUREMATRIX.
+
+arguments
+    templatePath (1,1) string = fullfile(fileparts(fileparts(mfilename("fullpath"))), ...
+                                          "schema", "channel_template.json")
+end
+
+t = jsondecode(fileread(templatePath));
+chans = t.channels;
+if ~iscell(chans), chans = num2cell(chans); end       % normalize to cell of structs
+
+Path=strings(0,1); Class=strings(0,1); Subclass=strings(0,1); Leaf=strings(0,1);
+Component=strings(0,1); CompIndex=zeros(0,1); Model=strings(0,1); Dtype=strings(0,1);
+Level=strings(0,1); Numeric=false(0,1); IsEmbedding=false(0,1);
+
+for i = 1:numel(chans)
+    c = chans{i};
+    p = string(c.path);
+    parts = split(p, "/");
+    cls  = parts(1);
+    sub  = "(direct)"; if numel(parts) > 2, sub = parts(2); end
+    leaf = parts(end);
+    dt   = string(c.dtype);
+    isNum = ismember(dt, ["scalar", "bool", "event", "vector"]);
+    isEmb = ismember(leaf, ["siglip_embedding", "dino_embedding", "clap_embedding"]);
+    lev   = localLevel(cls, sub, leaf);
+    mdl   = ""; if isfield(c, "model") && ~isempty(c.model), mdl = string(c.model); end
+
+    if dt == "vector"
+        comps = c.components;
+        D = double(c.dim);
+        if isempty(comps), comps = arrayfun(@(k) "d"+string(k-1), (1:D)'); end
+        for k = 1:D
+            comp = string(comps{k});
+            Path(end+1,1)=p; Class(end+1,1)=cls; Subclass(end+1,1)=sub; Leaf(end+1,1)=leaf; %#ok<*AGROW>
+            Component(end+1,1)=comp; CompIndex(end+1,1)=k; Model(end+1,1)=mdl;
+            Dtype(end+1,1)=dt; Level(end+1,1)=lev; Numeric(end+1,1)=isNum; IsEmbedding(end+1,1)=isEmb;
+        end
+    else
+        Path(end+1,1)=p; Class(end+1,1)=cls; Subclass(end+1,1)=sub; Leaf(end+1,1)=leaf;
+        Component(end+1,1)=""; CompIndex(end+1,1)=1; Model(end+1,1)=mdl;
+        Dtype(end+1,1)=dt; Level(end+1,1)=lev; Numeric(end+1,1)=isNum; IsEmbedding(end+1,1)=isEmb;
+    end
+end
+
+info = table(Path, Class, Subclass, Leaf, Component, CompIndex, Model, Dtype, ...
+             Level, Numeric, IsEmbedding);
+
+% Stable, unique, valid column names: base = path with "/"->"__"; components get a
+% "__<comp>" (or "__d<k-1>" for embeddings) suffix.
+base = replace(info.Path, "/", "__");
+suffix = strings(height(info), 1);
+isVec = info.Dtype == "vector";
+suffix(isVec & info.IsEmbedding) = "__d" + string(info.CompIndex(isVec & info.IsEmbedding) - 1);
+mask = isVec & ~info.IsEmbedding;
+suffix(mask) = "__" + info.Component(mask);
+raw = base + suffix;
+vn = matlab.lang.makeValidName(raw);
+vn = matlab.lang.makeUniqueStrings(vn, {}, namelengthmax);
+info.VarName = vn;
+
+% Per-class color (matches analysis/figures/feature_map.svg).
+info.Color = localColor(info.Class);
+
+info = movevars(info, "VarName", "Before", "Path");
+end
+
+% -------------------------------------------------------------------------
+function lev = localLevel(cls, sub, leaf)
+% Perceptual-abstraction level per channel; mirrors the FEATURE_MAP summary table.
+lev = "high";                                            % default (high-level semantics)
+switch cls
+    case "visual"
+        if ismember(sub, ["low_level_static", "dynamic_motion"])
+            lev = "low";
+        elseif ismember(sub, ["saliency_aesthetics_depth", "faces_bodies_gaze"])
+            lev = "mid";
+        end                                              % else high (high_level_static, action)
+    case "audio"
+        if sub == "low_level"
+            lev = "low";
+        elseif leaf == "speech_present" || leaf == "word_rate"
+            lev = "mid";
+        end                                              % else high (high_level, voice_*)
+    case "language"
+        if ismember(leaf, ["freq_zipf", "word_length"])
+            lev = "low";
+        elseif sub == "syntax" || ismember(leaf, ["valence","arousal","dominance","concreteness","aoa"])
+            lev = "mid";
+        end                                              % else high (surprisal, entropy)
+    case "social"
+        if ismember(leaf, ["n_agents", "min_pair_distance"]), lev = "mid"; end
+    case "affect"
+        if ismember(leaf, ["face_emotion", "face_valence", "face_arousal"]), lev = "mid"; end
+end
+end
+
+% -------------------------------------------------------------------------
+function C = localColor(classCol)
+map = ["visual", "#6366f1"; "audio", "#06b6d4"; "language", "#f59e0b"; ...
+       "social", "#ec4899"; "situation", "#10b981"; "affect", "#ef4444"];
+C = zeros(numel(classCol), 3);
+for i = 1:numel(classCol)
+    j = find(map(:,1) == classCol(i), 1);
+    hex = "#808080"; if ~isempty(j), hex = map(j,2); end
+    C(i,:) = double([hex2dec(extractBetween(hex,2,3)), hex2dec(extractBetween(hex,4,5)), ...
+                     hex2dec(extractBetween(hex,6,7))]) / 255;
+end
+end