EEGLAB drop-in output parity (#92)

.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.

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).

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

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``.

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."""

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:

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):