diff --git a/.context/scratch_history.md b/.context/scratch_history.md index 22743b8..b994cb2 100644 --- a/.context/scratch_history.md +++ b/.context/scratch_history.md @@ -45,6 +45,30 @@ python -m pyAMICA.amica_cli pyAMICA/sample_data/sample_params.json --verbose --o gfortran -O3 -fopenmp amica17.f90 funmod2.f90 -o amica -llapack -lblas ``` +## Issue #92 (EEGLAB drop-in output): the column-major layout fix (KEY) +Fortran/EEGLAB store arrays **column-major**. The MATLAB round-trip exposed that +pyAMICA wrote the non-square mixture params (`alpha`/`mu`/`sbeta`/`rho`, shape +`(num_mix, num_comps)`) and `c`/`comp_list` in **C-order**, so real MATLAB +`loadmodout15.m` (column-major reads) got scrambled mixture params -- e.g. the +per-component mixture proportions did NOT sum to 1. FIXED (issue #92): the writer +`write_amicaout` and BOTH numpy readers (`loadmodout`, `data.py:load_results`) now +use `order="F"` for those arrays. Diagnostic that nails the layout: read genuine +`sample_data/amicaout/alpha` -- `reshape(3,32,order='F').sum(0)` is all 1.0 +(correct); C-order is garbage `[0.49..1.31]`. +- `W` (square) stays C-order in the writer AND is byte-identical to Fortran: the + internal-vs-true-unmixing transpose (#24) cancels against Fortran's column-major + storage (`self.W` C-order bytes == Fortran column-major `W_true` bytes). `S` is + symmetric (order-agnostic); `mean`/`gm`/`LL` are 1-D. +- Remaining, deliberately-out-of-scope quirk: `loadmodout`/`load_results` still + read the **W** file C-order, so the port's `mod.W` is the transpose of MATLAB's + and its derived `A`/`svar`/`origord` use `pinv(self.W_int @ S)`. This does NOT + corrupt values (unlike the mixture bug, which did) and nothing consumes those + fields for correctness; every parity test matches `W` via transpose-tolerant + Hungarian |corr|. Do NOT flip the W read without re-checking #24/#37 parity. + Consequence: `AMICATorchNG.variance_order()` (uses `W_fort=self.W.T`, matching + real MATLAB) is validated against the MATLAB-faithful column-major reader, NOT + the numpy port's `origord`. + ## Lessons / check first next time - [ ] Positive LL almost always means a wrong PDF normalization constant. - [ ] NaN at a phase transition (e.g. Newton start) points at unclipped gradients or zero denominators. diff --git a/docs/guides/eeglab.md b/docs/guides/eeglab.md new file mode 100644 index 0000000..d06ea40 --- /dev/null +++ b/docs/guides/eeglab.md @@ -0,0 +1,72 @@ +# EEGLAB interoperability + +pyAMICA 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: + +```python +from pyAMICA 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 | +| `comp_list` | component ids (for component sharing) | +| `LL` | log-likelihood per iteration | + +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. + +## Loading in EEGLAB / MATLAB + +In MATLAB with the AMICA plugin on the path: + +```matlab +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 pyAMICA 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): + +```python +order = model.variance_order() # source indices, highest variance first +A = model.get_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. + +## Multi-model note + +Single-model output is byte-identical to the Fortran reference. For +`n_models > 1` the per-model axis layout is self-consistent (it round-trips +through `loadmodout15` and pyAMICA's own reader) but is not byte-identical to a +native multi-model AMICA run; see the multi-model equivalence discussion in +[Validation & Parity](validation.md). diff --git a/mkdocs.yml b/mkdocs.yml index 0e83a9e..3013d20 100644 --- a/mkdocs.yml +++ b/mkdocs.yml @@ -93,6 +93,7 @@ nav: - User Guide: - Overview: guides/index.md - Backends & Devices: guides/backends.md + - EEGLAB interoperability: guides/eeglab.md - Validation & Parity: guides/validation.md - API Reference: - Overview: api/index.md diff --git a/pyAMICA/amica.py b/pyAMICA/amica.py index 2ad70f6..9866050 100644 --- a/pyAMICA/amica.py +++ b/pyAMICA/amica.py @@ -319,6 +319,54 @@ def get_unmixing_matrix(self, model_idx: int = 0) -> np.ndarray: return self.model_.get_unmixing_matrix(model_idx=model_idx) + def variance_order( + self, model_idx: int = 0, return_svar: bool = False + ) -> Union[np.ndarray, tuple]: + """ + Component order by EEGLAB back-projected variance (IC1 = highest). + + Reports the display order EEGLAB's ``loadmodout15.m`` applies on load, + without mutating the fitted parameters. Apply it to the columns of + :meth:`get_mixing_matrix` (or rows of :meth:`get_unmixing_matrix`) to get + EEGLAB-ordered components in Python. + + Parameters + ---------- + model_idx : int, default=0 + Which model's components to order. + return_svar : bool, default=False + If True, also return the per-source variance sorted to ``order``. + + Returns + ------- + order : np.ndarray of int + Source indices, highest back-projected variance first. + """ + self._check_usable("compute the variance order") + + return self.model_.variance_order(model_idx=model_idx, return_svar=return_svar) + + def write_amica_output(self, outdir: str) -> None: + """ + Write the fitted model as an EEGLAB-readable AMICA output directory. + + Emits the raw binary files that EEGLAB's ``loadmodout15.m`` reads (``W``, + ``S``, ``gm``, ``mean``, ``c``, ``alpha``, ``mu``, ``sbeta``, ``rho``, + ``comp_list``, ``LL``), so a pyAMICA fit drops directly into an EEGLAB + workflow (``mod = loadmodout15(outdir)``). ``loadmodout15`` applies the + variance-ordering and normalization on load, so no manual re-ordering or + sign-flipping is needed. Single-model output is byte-compatible with the + Fortran reference (issue #92). + + Parameters + ---------- + outdir : str + Destination directory (created if absent). + """ + self._check_usable("write EEGLAB output") + + self.model_.write_amica_output(outdir) + def save(self, filepath: str) -> None: """ Save the fitted model to ``filepath`` via ``torch.save``. diff --git a/pyAMICA/numpy_impl/core.py b/pyAMICA/numpy_impl/core.py index 294aec6..a28a004 100644 --- a/pyAMICA/numpy_impl/core.py +++ b/pyAMICA/numpy_impl/core.py @@ -1367,44 +1367,28 @@ def _write_results(self): but is not byte-identical to multi-model Fortran output. Multi-model Fortran interop is out of scope here (see #27). """ - if not self.outdir.exists(): - self.outdir.mkdir(parents=True) - - def _w(name, arr, dtype=np.float64): - np.ascontiguousarray(arr, dtype=dtype).tofile(self.outdir / name) - - # gm: (num_models,) - _w("gm", self.gm) - # A: mixing matrix (data_dim, num_comps). Not part of the Fortran output - # (loadmodout derives A from W and S), but written here so load_results - # can restore it directly for the viz helpers; loadmodout ignores it. - _w("A", self.A) - # W: internal (nw, nw, num_models). For a single model the C-order dump - # is byte-identical to Fortran's 'W' (the internal-vs-true-unmixing - # transpose of issue #24 cancels against Fortran's column-major storage); - # for num_models>1 only the model-axis nesting differs (see the docstring - # note), and loadmodout reads it back with the matching C-order. - _w("W", self.W) - # Sphering and mean. - _w("S", self.sphere) - _w("mean", self.mean) - # Per-model bias: (nw, num_models). - _w("c", self.c) - # Mixture params: (num_mix, num_comps). loadmodout maps columns to - # sources via comp_list, so column order is the component index. - _w("alpha", self.alpha) - _w("mu", self.mu) - _w("sbeta", self.beta) # Fortran's 'sbeta' is pyAMICA's beta (scale) - _w("rho", self.rho) - # comp_list is 1-based in the Fortran format (loadmodout subtracts 1 - # when indexing); pyAMICA stores it 0-based. - _w("comp_list", self.comp_list + 1, dtype=np.int32) - # Log-likelihood history (per iteration). - _w("LL", np.asarray(self.ll)) - # Note: the Fortran 'nd' file is a per-component weight-change history - # (max_iter, nw, num_models); pyAMICA's self.nd is a per-iteration - # gradient-norm scalar, a different quantity, so it is not emitted in - # the Fortran format (loadmodout treats 'nd' as optional). + # A is written (Fortran output omits it; loadmodout derives A from W and + # S) only so load_results can restore it directly for the viz helpers. + # The Fortran 'nd' file (per-component weight-change history) is a + # different quantity from pyAMICA's scalar self.nd, so it is not emitted + # (loadmodout treats 'nd' as optional). + from .load import write_amicaout + + write_amicaout( + self.outdir, + gm=self.gm, + W=self.W, + sphere=self.sphere, + mean=self.mean, + c=self.c, + alpha=self.alpha, + mu=self.mu, + sbeta=self.beta, # Fortran's 'sbeta' is pyAMICA's beta (scale) + rho=self.rho, + comp_list=self.comp_list, + ll=np.asarray(self.ll), + A=self.A, + ) def _write_history(self): """Write optimization history at current iteration.""" diff --git a/pyAMICA/numpy_impl/data.py b/pyAMICA/numpy_impl/data.py index e03b921..b9bb2fa 100644 --- a/pyAMICA/numpy_impl/data.py +++ b/pyAMICA/numpy_impl/data.py @@ -245,7 +245,9 @@ def _read(name, dtype: type = np.float64): len(A) // num_comps, num_comps ) # (data_dim, num_comps) - # Mixture params are stored (num_mix, num_comps); Fortran names 'sbeta'. + # Mixture params are stored (num_mix, num_comps) column-major (Fortran names + # 'sbeta'); reshape order="F" matches the write_amicaout writer and Fortran + # output (a C-order read would scramble the non-square layout, issue #92). for fname, key in ( ("alpha", "alpha"), ("mu", "mu"), @@ -254,16 +256,16 @@ def _read(name, dtype: type = np.float64): ): arr = _read(fname) if arr is not None: - results[key] = arr.reshape(-1, num_comps) + results[key] = arr.reshape(-1, num_comps, order="F") comp_list = _read("comp_list", dtype=np.int32) if comp_list is not None: # Stored 1-based (Fortran convention); restore AMICA's 0-based indices. - results["comp_list"] = comp_list.reshape(nw, num_models) - 1 + results["comp_list"] = comp_list.reshape(nw, num_models, order="F") - 1 c = _read("c") if c is not None: - results["c"] = c.reshape(nw, num_models) + results["c"] = c.reshape(nw, num_models, order="F") mean = _read("mean") if mean is not None: diff --git a/pyAMICA/numpy_impl/load.py b/pyAMICA/numpy_impl/load.py index 4ca3f1c..8fccf4f 100644 --- a/pyAMICA/numpy_impl/load.py +++ b/pyAMICA/numpy_impl/load.py @@ -48,6 +48,84 @@ def read_binary_file( return None +def write_amicaout( + outdir: Union[str, Path], + *, + gm, + W, + sphere, + mean, + c, + alpha, + mu, + sbeta, + rho, + comp_list, + ll, + A=None, +): + """Write a fitted AMICA model as the Fortran/EEGLAB binary output directory. + + Emits the raw little-endian files that :func:`loadmodout` and EEGLAB's + ``loadmodout15.m`` read: ``gm``, ``W``, ``S``, ``mean``, ``c``, ``alpha``, + ``mu``, ``sbeta``, ``rho``, ``comp_list`` (1-based ``int32``) and ``LL``. + This is the write counterpart of :func:`loadmodout`, so a pyAMICA fit (either + backend) drops into an EEGLAB workflow (issue #92). + + Both backends store these arrays in the same convention, so for a single + model the bytes are identical to the Fortran reference's ``amicaout`` files; + for ``num_models > 1`` the per-model axis nesting is self-consistent (it + round-trips through :func:`loadmodout`) but not byte-identical to genuine + multi-model Fortran output (issue #27). + + Parameters + ---------- + outdir : str or path-like + Destination directory (created if absent). + gm, W, sphere, mean, c, alpha, mu, sbeta, rho : array-like + Model weights, unmixing, sphere, data mean, per-model centers and the + mixture-density parameters (``sbeta`` is the scale, pyAMICA's ``beta``). + comp_list : array-like of int + 0-based component ids; written 1-based to match the Fortran format. + ll : array-like + Per-iteration log-likelihood history. + A : array-like, optional + Mixing matrix. ``loadmodout15`` derives ``A`` from ``W`` and ``S`` and + ignores this file; it is written (when given) only so pyAMICA's own + ``load_results`` can restore ``A`` directly for the viz helpers. + """ + outdir = Path(outdir) + outdir.mkdir(parents=True, exist_ok=True) + + def _w(name, arr, dtype=np.float64, order="C"): + # Fortran dumps arrays column-major; ``order="F"`` reproduces that byte + # layout so real EEGLAB ``loadmodout15.m`` reads the file correctly. + np.asarray(arr, dtype=dtype).ravel(order=order).tofile(outdir / name) + + _w("gm", gm) + if A is not None: + _w("A", A) + # W is byte-identical to Fortran in C order: the internal-vs-true-unmixing + # transpose (issue #24) cancels against Fortran's column-major storage, so a + # square W written C-order equals the Fortran/EEGLAB column-major bytes. The + # symmetric sphere S is order-agnostic; mean/gm/LL are 1-D. + _w("W", W) + _w("S", sphere) + _w("mean", mean) + # The (num_mix, num_comps) mixture params and (num_comps, num_models) c / + # comp_list are non-square, so their byte layout DOES depend on order: they + # must be column-major (Fortran) for loadmodout15 to read them correctly + # (e.g. mixture proportions per component sum to 1). Issue #92. + _w("c", c, order="F") + _w("alpha", alpha, order="F") + _w("mu", mu, order="F") + _w("sbeta", sbeta, order="F") + _w("rho", rho, order="F") + # comp_list is 1-based on disk (loadmodout subtracts 1 when indexing). + _w("comp_list", np.asarray(comp_list) + 1, dtype=np.int32, order="F") + _w("LL", np.asarray(ll)) + + def loadmodout(outdir: Union[str, Path]) -> AmicaOutput: """Load AMICA output files from directory. @@ -106,9 +184,9 @@ def loadmodout(outdir: Union[str, Path]) -> AmicaOutput: # else is a corrupt or unexpected file and should fail loudly # rather than silently truncating. if comp_list.size == expected: - comp_list = comp_list.reshape(nw, num_models) + comp_list = comp_list.reshape(nw, num_models, order="F") elif comp_list.size == 2 * expected and not np.any(comp_list[expected:]): - comp_list = comp_list[:expected].reshape(nw, num_models) + comp_list = comp_list[:expected].reshape(nw, num_models, order="F") else: raise ValueError( f"comp_list has {comp_list.size} elements; expected {expected} " @@ -136,13 +214,15 @@ def loadmodout(outdir: Union[str, Path]) -> AmicaOutput: if c is None: c = np.zeros((nw, num_models)) else: - c = c.reshape(nw, num_models) + c = c.reshape(nw, num_models, order="F") - # Read mixture parameters + # Read mixture parameters. Stored column-major (Fortran), so reshape order="F" + # (matches loadmodout15.m and the write_amicaout writer); a C-order read would + # scramble the (num_mix, num_comps) layout. Issue #92. alpha_tmp = read_binary_file(outdir / "alpha") if alpha_tmp is not None: num_mix = len(alpha_tmp) // (nw * num_models) - alpha_tmp = alpha_tmp.reshape(num_mix, nw * num_models) + alpha_tmp = alpha_tmp.reshape(num_mix, nw * num_models, order="F") alpha = np.zeros((num_mix, nw, num_models)) for h in range(num_models): for i in range(nw): @@ -158,7 +238,7 @@ def loadmodout(outdir: Union[str, Path]) -> AmicaOutput: # Read mu, sbeta, rho mu_tmp = read_binary_file(outdir / "mu") if mu_tmp is not None: - mu_tmp = mu_tmp.reshape(num_mix, nw * num_models) + mu_tmp = mu_tmp.reshape(num_mix, nw * num_models, order="F") mu = np.zeros((num_mix, nw, num_models)) for h in range(num_models): for i in range(nw): diff --git a/pyAMICA/tests/torch_tests/test_amica_ng_wrapper.py b/pyAMICA/tests/torch_tests/test_amica_ng_wrapper.py index b0066f7..bac7fba 100644 --- a/pyAMICA/tests/torch_tests/test_amica_ng_wrapper.py +++ b/pyAMICA/tests/torch_tests/test_amica_ng_wrapper.py @@ -139,6 +139,119 @@ def test_ng_mps_float32_escape_hatch(real_data): assert model.model_.dtype == torch.float32 +# --- EEGLAB drop-in output (issue #92) ------------------------------------- +# Files EEGLAB's loadmodout15.m / the numpy port loadmodout() read. +_AMICAOUT_FILES = ("gm", "W", "S", "mean", "c", "alpha", "mu", "sbeta", "rho", + "comp_list", "LL") # fmt: skip + + +def test_write_amica_output_requires_fit(tmp_path): + """An unfit model must refuse to write output, mirroring save() (#50).""" + model = AMICA(verbose=False) + with pytest.raises(ValueError, match="fitted"): + model.write_amica_output(str(tmp_path / "amicaout")) + + +def test_variance_order_requires_fit(): + """variance_order also refuses an unfit model (same _check_usable guard).""" + model = AMICA(verbose=False) + with pytest.raises(ValueError, match="fitted"): + model.variance_order() + + +def test_write_amica_output_bytes(fitted_ng, tmp_path): + """The written files are the model's exact float64 parameters: the on-disk + EEGLAB directory is a lossless serialization, not a lossy export (#92). W and + the symmetric sphere are byte-identical in C order; the non-square mixture + params and c/comp_list are column-major (Fortran layout), so read order="F". + """ + outdir = tmp_path / "amicaout" + fitted_ng.write_amica_output(str(outdir)) + ng = fitted_ng.model_ + + for name, attr in [("gm", ng.gm), ("W", ng.W), ("S", ng.sphere), + ("mean", ng.mean)]: # fmt: skip + got = np.fromfile(outdir / name).reshape(attr.shape) # C order + np.testing.assert_array_equal(got, attr.cpu().numpy(), err_msg=name) + for name, attr in [("c", ng.c), ("alpha", ng.alpha), ("mu", ng.mu), + ("sbeta", ng.beta), ("rho", ng.rho)]: # fmt: skip + got = np.fromfile(outdir / name).reshape(attr.shape, order="F") + np.testing.assert_array_equal(got, attr.cpu().numpy(), err_msg=name) + # comp_list is written 1-based int32, column-major. + np.testing.assert_array_equal( + np.fromfile(outdir / "comp_list", np.int32).reshape( + ng.comp_list.shape, order="F" + ), + ng.comp_list.cpu().numpy() + 1, + ) + # LL ends at the exported (kept) iterate; the fixture fit is monotone, so the + # full trajectory is written. + ll = np.fromfile(outdir / "LL") + assert ll[-1] == ng.final_ll_ + np.testing.assert_array_equal( + ll, np.asarray(ng.ll_history[: len(ll)], dtype=np.float64) + ) + + +def test_write_amica_output_loadmodout_readable(fitted_ng, tmp_path): + """A PyTorch NG fit written with write_amica_output() is a directory the + EEGLAB reader (loadmodout / loadmodout15) loads with the expected shapes and + correct Fortran (column-major) mixture-param layout (issue #92).""" + from pyAMICA.numpy_impl.load import loadmodout + + outdir = tmp_path / "amicaout" + fitted_ng.write_amica_output(str(outdir)) + for name in _AMICAOUT_FILES: + assert (outdir / name).exists(), f"missing {name}" + + mod = loadmodout(outdir) + assert mod.num_models == 1 + assert mod.W.shape == (NW, NW, 1) + assert mod.A.shape == (NW, NW, 1) + assert mod.S.shape == (NW, NW) + # Mixture proportions per component must sum to 1: a meaningful check that the + # (num_mix, n_comp) params were read back with the correct column-major layout + # (a C-order write would scramble them and break this). + np.testing.assert_allclose(mod.alpha[:, :, 0].sum(axis=0), 1.0, atol=1e-9) + + +def test_variance_order(fitted_ng): + """variance_order() returns a permutation of the sources ranked by descending + back-projected variance (EEGLAB IC1 = highest), for use in Python without a + disk round-trip (issue #92).""" + order, svar = fitted_ng.variance_order(return_svar=True) + assert sorted(order.tolist()) == list(range(NW)) # a permutation + assert np.all(np.diff(svar) <= 1e-9) # descending + + +def test_write_amica_output_ll_matches_kept_iterate(real_data, tmp_path): + """When keep_best (#51) restores an earlier iterate, the written LL trajectory + ends at that kept iterate (LL[-1] == final_ll_), not at a later discarded + overshoot -- so a user reading mod.LL(end) in EEGLAB sees the loaded model's + likelihood (review finding, #92).""" + model = AMICA(n_models=2, n_mix=3, device="cpu", verbose=False) + model.fit( + real_data[:, :4096], + max_iter=60, + do_newton=True, + newt_start=1, + lrate=0.5, + seed=0, + block_size=1024, + ) + if not model.is_fitted_: + pytest.skip("aggressive run ended degenerate; not the case under test") + ng = model.model_ + if np.isclose(ng.ll_history[-1], ng.final_ll_): + pytest.skip("run was monotone; keep_best restore did not fire") + + outdir = tmp_path / "amicaout" + model.write_amica_output(str(outdir)) + ll = np.fromfile(outdir / "LL") + assert np.isclose(ll[-1], ng.final_ll_) # ends at the kept iterate + assert len(ll) < len(ng.ll_history) # the later overshoot is dropped + + def test_ng_wrapper_fit_transform_real_data(fitted_ng, real_data): assert fitted_ng.is_fitted_ assert len(fitted_ng.ll_history_) >= 1 @@ -202,6 +315,8 @@ def test_degenerate_fit_refuses_output(real_data, tmp_path, caplog): lambda: model.get_mixing_matrix(), lambda: model.get_unmixing_matrix(), lambda: model.save(str(tmp_path / "degenerate.pt")), + lambda: model.write_amica_output(str(tmp_path / "degenerate_out")), + lambda: model.variance_order(), ): with pytest.raises(RuntimeError, match="degenerate.*nan_ll"): action() diff --git a/pyAMICA/torch_impl/core.py b/pyAMICA/torch_impl/core.py index 1c78ac7..3c1d796 100644 --- a/pyAMICA/torch_impl/core.py +++ b/pyAMICA/torch_impl/core.py @@ -1767,6 +1767,107 @@ def get_unmixing_matrix(self, model_idx: int = 0) -> np.ndarray: """True unmixing matrix ``W_fort`` = (stored W)^T (issue #24 convention).""" return self.W[:, :, model_idx].T.cpu().numpy() + # ------------------------------------------------------------------ + # EEGLAB drop-in output (issue #92) + # ------------------------------------------------------------------ + def variance_order( + self, model_idx: int = 0, return_svar: bool = False + ) -> Union[np.ndarray, tuple]: + """EEGLAB back-projected-variance component order (IC1 = highest variance). + + Returns the source indices sorted by descending back-projected variance, + the ordering EEGLAB's ``loadmodout15.m`` applies on load (so ``order[0]`` + is IC1). The de-sphered sensor-space mixing column ``a_i = pinv(W S)[:, i]`` + contributes ``||a_i||^2 * sum_k alpha_ki (mu_ki^2 + r_ki / sbeta_ki^2)`` + with ``r_ki = gamma(3/rho_ki)/gamma(1/rho_ki)`` (the source's mixture + variance), matching ``loadmodout15`` exactly. Non-mutating: the stored + parameters keep their fit order; this only reports the display order. + + Parameters + ---------- + model_idx : int, default=0 + Which model's components to order. + return_svar : bool, default=False + If True, also return the per-source variance sorted to ``order``. + + Returns + ------- + order : np.ndarray of int, shape (n_sources,) + Source indices, highest back-projected variance first. + svar : np.ndarray, optional + Present only when ``return_svar``; the sorted variances. + """ + from scipy.special import gamma + + cl = self.comp_list[:, model_idx].cpu().numpy() + alpha = self.alpha[:, cl].cpu().numpy() + mu = self.mu[:, cl].cpu().numpy() + sbeta = self.beta[:, cl].cpu().numpy() + rho = self.rho[:, cl].cpu().numpy() + # source mixture variance (sum over the mixture components); unused + # mixtures carry alpha == 0 and drop out, matching loadmodout15. + ratio = gamma(3.0 / rho) / gamma(1.0 / rho) + mix_var = (alpha * (mu**2 + ratio / sbeta**2)).sum(axis=0) + # de-sphered sensor-space mixing: A = pinv(W_fort @ S), columns = maps. + w_fort = self.W[:, :, model_idx].T.cpu().numpy() + sphere = self.sphere.cpu().numpy() + a_sensor = np.linalg.pinv(w_fort @ sphere) + svar = mix_var * (a_sensor**2).sum(axis=0) + order = np.argsort(-svar) + if return_svar: + return order, svar[order] + return order + + def write_amica_output(self, outdir) -> None: + """Write this fitted model as the Fortran/EEGLAB AMICA output directory. + + Produces the raw binary files that EEGLAB's ``loadmodout15.m`` (and the + Python port :func:`pyAMICA.numpy_impl.load.loadmodout`) read, so a + PyTorch NG fit drops directly into an EEGLAB workflow (issue #92). + ``loadmodout15`` performs the variance-ordering and unit-norm + normalization on load, so the on-disk parameters are written in fit + order. Single-model output is byte-compatible with the Fortran reference. + + Parameters + ---------- + outdir : str or path-like + Destination directory (created if absent). + """ + from ..numpy_impl.load import write_amicaout + + def _np(t): + return t.detach().cpu().numpy() + + # The exported parameters are the fit()-kept iterate (LL == final_ll_). + # Under the keep_best safeguard (#51) that can be an earlier iterate than + # the last, so end the written LL trajectory at that iterate rather than + # at a later, discarded overshoot -- otherwise LL[-1] would not match the + # model just written. Monotone runs keep the full trajectory unchanged. + ll = np.asarray(self.ll_history, dtype=np.float64) + if ( + self.final_ll_ is not None + and np.isfinite(self.final_ll_) + and ll.size + and not np.isclose(ll[-1], self.final_ll_) + ): + ll = ll[: int(np.argmax(ll)) + 1] + + write_amicaout( + outdir, + gm=_np(self.gm), + W=_np(self.W), + sphere=_np(self.sphere), + mean=_np(self.mean), + c=_np(self.c), + alpha=_np(self.alpha), + mu=_np(self.mu), + sbeta=_np(self.beta), # Fortran's 'sbeta' is pyAMICA's beta (scale) + rho=_np(self.rho), + comp_list=_np(self.comp_list), + ll=ll, + A=_np(self.A), + ) + # ------------------------------------------------------------------ # Persistence (issue #36) # ------------------------------------------------------------------