Determine AMICA API

[snesmaeili]

see #13819 (comment)

could be related: https://github.com/sccn/eegprep

[snesmaeili]

I have a working prototype of the models-as-views design now, so there is something concrete to argue about rather than a sketch: snesmaeili/jamica#33

I want to be explicit up front that this is deliberately not an MNE-side PR, and it does not touch _KNOWN_ICA_METHODS or add any dependency. Which AMICA implementation MNE depends on — if any — is a maintainer decision, and there are several (sccn/pAMICA, scott-huberty/amica-python, DerAndereJohannes/pyamica, and mine). What I wanted to find out was whether the API we discussed actually survives contact with mne.preprocessing.ICA, and where it breaks.

A bug the prototype turned up

Worth reporting on its own. My M-step gated its per-model centre update on the same flag that controls global mean removal in preprocessing. Those are different questions, and the MNE path centres the data externally — so every model centre was frozen at its zero initialisation and a mixture could only distinguish its models by their unmixing matrices. Anyone else wiring a multi-model AMICA behind MNE's PCA is likely to hit the same thing, since it only shows up when the caller does its own centring.

What the API looks like

amica = AmicaICA(n_models=3, n_components=20, random_state=42).fit(raw)

amica.models_[0].plot_components()      # ordinary mne.preprocessing.ICA
amica.models_[0].exclude = [2, 5]
clean = amica.models_[0].apply(raw.copy())

amica.model_weights_          # (M,) priors
amica.model_posteriors_       # p(h | x_t) on the input grid
amica.fit_sample_mask_        # which samples drove the fit
amica.get_model_probabilities(other_raw)

Things that turned out to matter in practice:

• Children are built by assignment, not ICA.fit(). fit() sorts components by explained variance afterwards, which breaks the correspondence between a component and its density parameters.
• The per-model centre folds into pca_mean_. AMICA's sources are W_h @ (x_white - c_h) and MNE computes U @ P @ (x - pca_mean_); since P has orthonormal rows the offset is P.T @ (comp_stds * c_h). So an individual model really is representable as a plain ICA, with no new class and no changes to MNE.
• Shared vs per-model state. pca_components_ / pre_whitener_ / pca_explained_variance_ describe the one preprocessing pipeline and are shared; pca_mean_, unmixing_matrix_, mixing_matrix_, exclude, labels_ are per-child. Children are cached, because interactive review mutates exclude and it would be bad to hand back a fresh object each access.
• apply() raises for M > 1, as agreed above, and names model_idx= and models_[h].apply() in the error.

Tests are against real RawArray / EpochsArray / ICA objects rather than mocks, including per-model source equivalence with non-zero centres, residual-PCA preservation, and child state isolation.

What I would like input on

1. Saving these objects. This is the one @larsoner flagged and it is a genuine blocker. FIF cannot round-trip M models today: read_ica does ica_data = dir_tree_find(tree, FIFF.FIFFB_MNE_ICA) and then takes ica_data[0], so writing M blocks would load back as model 0 with no error at all. Silent truncation seems worse than refusing to save. The closest in-tree precedent I found is EOGRegression.save / read_eog_regression in mne/preprocessing/_regress.py, which uses h5io; that would work here, with a decision needed about the posterior array since it is M x n_times. Would an HDF5 container be acceptable, or would you rather extend the FIF reader?

2. Where this should live. It may be that the mixture container belongs in the backend package and MNE only ever sees single-model ICA objects. I genuinely do not know which you prefer, and it changes what a future MNE-side PR would look like.

3. Naming. AmicaICA ties it to one algorithm. If other mixture methods are plausible, something neutral might age better.

4. Decimation semantics. I report posteriors on the original input grid regardless of decimation, and keep a separate mask for which samples drove the fit. One wrinkle: scipy.signal.decimate anti-alias filters before downsampling, so every input sample contributes to the retained ones and "these samples were excluded" is not really true. I left the mask all-true under decimation rather than claim otherwise, but I am not certain that is the most useful convention for users.

Happy to reshape any of this — the point of building it was to find the parts that do not work, and persistence is the one I could not resolve on my own.

[larsoner]

For FIF it does [0] today because there is only one model. This makes it maybe easy to extend to multiple: if only 1 available, it's traditional ICA, if multiple, it's AMICA. If we want to be safer, we encode also a FIFF_DESCRIPTION or similar that says what type it is somewhere (backward compatible: if missing, regular ICA). So FIFF storage doesn't seem insurmountable.

Backing with HDF5 also seems okay as well, though.

Okay with me to have this in MNE-Python or jamica. I'd actually lean slightly toward jamica and in MNE we can advertise / show the use in our ICA tutorial. So it doesn't become an optional dependency but a sibling package. Will allow faster iteration and releases etc. But not a strong feeling from me

[snesmaeili]

Thanks @larsoner that clears up the FIF point.
For now I kept the mixture container in jamica and used HDF5, mainly so it doesn’t require changes to MNE’s reader. AmicaICA.save() uses the same h5io helper as EOGRegression, and export_model_fifs() writes each model as a normal -ica.fif. I also checked that those load with plain mne.preprocessing.read_ica without jamica installed.
The only thing I’m still unsure about is the single-model case from #13819: when you say sibling package, do you also mean no ICA(method="jamica") entry in MNE? I had imagined a Picard-style optional backend for M=1, while keeping the multi-model container entirely in jamica.
If that is worth considering, I’d prefer method="jamica" rather than "amica", so MNE isn’t implicitly choosing one implementation as the AMICA.
Probably easiest to settle that at the sprint on Thursday. I won’t open an MNE-side PR unless there’s appetite for it.
Implementation for reference: snesmaeili/jamica#33

[larsoner]

Ahh yes I think that seems fairly trivial for MNE-Python to support (right?). If so then yes let's add it here with a note for that method that if people want multi-model support they should use jamica directly. That will also help draw attention to the other package because it'll show up explicitly in the method options