EEGLAB interoperability¶

pamica is a drop-in replacement for EEGLAB's AMICA: a fit written to disk loads directly with the same reader EEGLAB uses (loadmodout15.m), with the components in the same order and orientation, so no manual re-sorting, sign-flipping, or reformatting is needed.

Writing EEGLAB-readable output¶

After a fit, call write_amica_output with a destination directory:

from pamica import AMICA

model = AMICA(n_models=1, n_mix=3)
model.fit(X)                      # X is (n_channels, n_samples)
model.write_amica_output("amicaout")

This writes the raw binary files EEGLAB's AMICA loader reads:

File	Contents
`gm`	model probabilities
`W`	unmixing weights (post-sphering)
`S`	sphering matrix
`mean`	data mean
`c`	per-model centers
`alpha`, `mu`, `sbeta`, `rho`	source mixture-density parameters
`A`	mixing matrix in sphered space, one column per component (the reference's layout, for any number of models)
`comp_list`	component ids (for component sharing)
`LL`	log-likelihood per iteration
`LLt`	per-timepoint, per-model log-likelihood, plus the per-timepoint total

For a single model the bytes are identical to the reference Fortran binary's amicaout files, so the directory is interchangeable with a native AMICA run. A directory written with do_approx_sphere=False (a full-rank, genuinely asymmetric S) by a pamica older than issue #336's fix has S transposed; re-run write_amica_output to regenerate it (the default symmetric sphere and a rank-reduced fit were unaffected). loadmodout15.m does not read A; pamica's own load_results does, and refuses an A that does not invert the W beside it. A multi-model directory written by a pamica older than issue #334 stored A in another layout and must be written again to load there (a single-model A was already in the reference layout).

The MLX backend writes the same directory: an Apple-Silicon fit exports to EEGLAB directly, with no torch round trip (epic #278). The export is validated against a torch twin in the test suite to float32 precision (max abs diff < 1e-5; comp_list is the only field compared exactly -- MLX computes in float32, so its files are not bit-identical to a float64 export the way a torch fit's are to Fortran's):

from pamica import AMICA

model = AMICA(n_mix=3, backend="mlx")   # requires the mlx extra (issue #313)
model.fit(X)                            # X is (n_channels, n_samples)
model.write_amica_output("amicaout")

LLt is what loadmodout15.m turns into Lht/Lt and the model-probability odds v; it is written after a fresh fit(), and omitted (with a warning) for a model restored from load(), which carries no E-step stash. Under do_reject, rejected samples are written as exactly 0.0 to match the reference: AMICA's own load_rej reconstructs the rejection mask from those zeros, so they are load-bearing rather than padding.

Like the reference, pamica writes the E-step that produced the last entry of LL, so LLt describes the parameters as they stood one M-step before the W written beside it (issue #157; see amica-differences). Use model.model_loglik(X) if you want the log-likelihood of the written parameters themselves.

Loading in EEGLAB / MATLAB¶

In MATLAB with the AMICA plugin on the path:

mod = loadmodout15('amicaout');
% mod.W   : unmixing weights (n x n x num_models)
% mod.A   : component scalp maps, columns ordered IC1..ICn by variance
% mod.S   : sphering matrix
% mod.svar: back-projected variance per component

loadmodout15 applies the EEGLAB conventions on load: it orders components by back-projected variance (IC1 has the highest), derives the sensor-space mixing A = pinv(W * S), and normalizes each map to unit norm. Because pamica writes the same format, the components you get in EEGLAB match a native AMICA run.

Variance ordering in Python¶

To get the EEGLAB display order without a disk round-trip, use variance_order, which ranks sources by the same back-projected variance (IC1 = highest):

order = model.variance_order()           # source indices, highest variance first
A = model.get_sensor_mixing_matrix()[:, order]  # scalp maps in EEGLAB order
W = model.get_unmixing_matrix()[order]   # unmixing rows in EEGLAB order

Pass return_svar=True to also get the per-component variances.

get_sensor_mixing_matrix() gives the maps in input-channel space; get_mixing_matrix() is the sphered-space matrix, which is not a scalp map.

Multi-model note¶

Every file is written in the reference's layout for any number of models: W with the model axis slowest (issue #159) and A with one column per component (issue #334), so loadmodout15 reads a multi-model directory the same way it reads a native one. The values themselves are another matter: multi-model AMICA is not partition-identifiable, so two runs, native or pamica, do not produce the same models; see the multi-model discussion in Validation & Parity.