Support custom diagonal source covariance matrix for inverse operator

.git-blame-ignore-revs

@@ -14,3 +14,4 @@ e39995d9be6fc831c7a4a59f09b7a7c0a41ae315  # 12588, percent formatting
 b8b168088cb474f27833f5f9db9d60abe00dca83  # 12779, PR JSONs
 ee64eba6f345e895e3d5e7d2804fa6aa2dac2e6d  # 12781, Header unification
 362f9330925fb79a6adc19a42243672676dec63e  # 12799, UP038
+c166d8c191d07b0efa7d6772bb0fffcf4cdb9252  # 14316, Fully populate docstrings

doc/changes/dev/14315.newfeature.rst

@@ -0,0 +1 @@
+Add parameter ``source_cov`` for :func:`mne.minimum_norm.make_inverse_operator`, allowing use of a custom (diagonal) source covariance matrix for the inverse solution, by `Teemu Taivainen`_.

doc/sphinxext/credit_tools.py

@@ -521,9 +521,11 @@ def generate_credit_rst(
         logger.info(f"{pr.ljust(5)} @ {commits[(name, pr)]:5d} by {name}")
 
     logger.info("\nIgnored commits:")
-    for pr, files in ignores.items():  # should have found one of each
+    # should have found one of each, except PRs whose JSON isn't fetched yet
+    last_pr = max(int(pr) for _, pr in commits)
+    for pr, files in ignores.items():
         logger.info(f"ignored {len(files):3d} files for {pr}")
-        assert len(files) >= 1, (pr, files)
+        assert len(files) >= 1 or int(pr) > last_pr, (pr, files)
 
     mod_stats, link_overrides, mod_file_map = _aggregate_module_stats(stats)
     _write_credit_rst(mod_stats, link_overrides, mod_file_map, urls)

mne/minimum_norm/inverse.py

@@ -1846,6 +1846,7 @@ def _prepare_forward(
     combine_xyz,
     allow_fixed_depth,
     limit,
+    source_cov=None,
 ):
     """Prepare a gain matrix and noise covariance for localization."""
     # Steps (according to MNE-C, we change the order of various steps
@@ -1983,6 +1984,18 @@ def _prepare_forward(
 
     logger.info("Creating the source covariance matrix")
     source_std = np.ones(gain.shape[1], dtype=gain.dtype)
+    if source_cov is not None:
+        source_cov = np.asarray(source_cov, dtype=np.float64)
+        n_sources = forward["nsource"]
+        if source_cov.shape != (n_sources,):
+            raise ValueError(
+                f"source_cov must have shape ({n_sources},), got {source_cov.shape}"
+            )
+        if not (np.isfinite(source_cov).all() and (source_cov > 0).all()):
+            raise ValueError("source_cov must contain finite, positive variances")
+        logger.info("    Using user-specified source variances")
+        # one variance per source location, applied to all orientations
+        source_std *= np.repeat(source_cov, gain.shape[1] // n_sources)
     if depth_prior is not None:
         source_std *= depth_prior
     if orient_prior is not None:
@@ -2021,6 +2034,7 @@ def make_inverse_operator(
     fixed="auto",
     rank=None,
     use_cps=True,
+    source_cov=None,
     verbose=None,
 ):
     """Assemble inverse operator.
@@ -2113,6 +2127,14 @@ def make_inverse_operator(
     use_cps : bool
         Whether to use cortical patch statistics to define normal orientations for
         surfaces (default True).
+    source_cov : array-like, shape (n_sources,) | None
+        Prior variance of each source location, i.e., the diagonal of a custom
+        source covariance matrix. It is applied to all orientations of a source
+        and multiplied by the depth and orientation priors determined by
+        ``depth``, ``loose``, and ``fixed``. If None (default), all sources get
+        the same prior variance.
+
+        .. versionadded:: 1.14
     verbose : bool | str | int | None
         Control verbosity of the logging output. If ``None``, use the default
         verbosity level. See the :ref:`logging documentation <tut-logging>` and
@@ -2188,6 +2210,7 @@ def make_inverse_operator(
         rank,
         pca="white",
         use_cps=use_cps,
+        source_cov=source_cov,
         **depth,
     )
     # no need to copy any attributes of forward here because there is

mne/minimum_norm/tests/test_inverse.py

@@ -879,7 +879,29 @@ def test_make_inverse_operator_fixed(evoked, noise_cov):
     assert "EEG channels: 0" in repr(inv_op)
     assert "MEG channels: 305" in repr(inv_op)
     assert "Fixed" in repr(inv_op)
-    del fwd_fixed
+
+    # uniform source_cov should be equivalent to the default (None)
+    kwargs = dict(depth=0.0, fixed=True, use_cps=False)
+    inv_op_ones = make_inverse_operator(
+        evoked.info, fwd, noise_cov, source_cov=np.ones(fwd["nsource"]), **kwargs
+    )
+    _compare_inverses_approx(inv_op, inv_op_ones, evoked, rtol=1e-5, atol=1e-8)
+    # non-uniform source_cov should scale the final source covariance
+    source_cov = np.ones(fwd["nsource"])
+    source_cov[0] = 4.0
+    inv_op_scaled = make_inverse_operator(
+        evoked.info, fwd, noise_cov, source_cov=source_cov, **kwargs
+    )
+    got = inv_op_scaled["source_cov"]["data"]
+    assert_allclose(got / got[1], source_cov)
+    with pytest.raises(ValueError, match="must have shape"):
+        make_inverse_operator(evoked.info, fwd, noise_cov, source_cov=[1.0], **kwargs)
+    with pytest.raises(ValueError, match="finite, positive"):
+        make_inverse_operator(
+            evoked.info, fwd, noise_cov, source_cov=-source_cov, **kwargs
+        )
+    del fwd_fixed, inv_op_ones, inv_op_scaled, source_cov, kwargs
+
     inverse_operator_nodepth = read_inverse_operator(fname_inv_fixed_nodepth)
     # XXX We should have this but we don't (MNE-C doesn't restrict info):
     # assert 'EEG channels: 0' in repr(inverse_operator_nodepth)
@@ -932,6 +954,25 @@ def test_make_inverse_operator_free(evoked, noise_cov):
         stc_surf = apply_inverse(evoked, inv_surf, pick_ori=pick_ori)
         assert_allclose(stc_surf.data, stc.data, atol=1e-2)
 
+    # uniform source_cov should be equivalent to the default (None)
+    inv_ones = make_inverse_operator(
+        evoked.info,
+        fwd,
+        noise_cov,
+        depth=None,
+        loose=1.0,
+        source_cov=np.ones(fwd["nsource"]),
+    )
+    _compare_inverses_approx(inv, inv_ones, evoked, rtol=1e-5, atol=1e-8)
+    # non-uniform source_cov should scale all three orientations of a source
+    source_cov = np.ones(fwd["nsource"])
+    source_cov[-1] = 5.0
+    inv_scaled = make_inverse_operator(
+        evoked.info, fwd, noise_cov, depth=None, loose=1.0, source_cov=source_cov
+    )
+    got = inv_scaled["source_cov"]["data"]
+    assert_allclose(got / got[0], np.repeat(source_cov, 3))
+
 
 @pytest.mark.slowtest
 def test_make_inverse_operator_vector(evoked, noise_cov):

mne/viz/backends/_jupyterlite.py

@@ -31,7 +31,8 @@
 _orig = dict()  # the functions patched below that are imported lazily, by name
 
 
-def setup_notebook(data_path="/tmp/mne_data"):
+# /tmp here is the Pyodide in-browser virtual filesystem, not a shared host dir
+def setup_notebook(data_path="/tmp/mne_data"):  # nosec B108
     """Patch MNE for the browser kernel; the docs' first cell calls this.
 
     Parameters