From 238dc2f053f5fabb2abf97c6b2e3b3010ba8512e Mon Sep 17 00:00:00 2001
From: olafhauk <olaf.hauk@mrc-cbu.cam.ac.uk>
Date: Wed, 24 Aug 2022 16:55:12 +0100
Subject: [PATCH 01/27] WIP: added covariance to resolution matrix

---
 mne/minimum_norm/resolution_matrix.py | 30 +++++++++++++++++++++++++--
 1 file changed, 28 insertions(+), 2 deletions(-)

diff --git a/mne/minimum_norm/resolution_matrix.py b/mne/minimum_norm/resolution_matrix.py
index a027867091f..1d21d0538d2 100644
--- a/mne/minimum_norm/resolution_matrix.py
+++ b/mne/minimum_norm/resolution_matrix.py
@@ -11,14 +11,15 @@
 from ..io.constants import FIFF
 from ..utils import logger, verbose
 from ..forward.forward import convert_forward_solution
-from ..minimum_norm import apply_inverse
+from ..minimum_norm import apply_inverse, apply_inverse_cov
 from ..source_estimate import _prepare_label_extraction
 from ..label import Label
 
 
 @verbose
 def make_inverse_resolution_matrix(forward, inverse_operator, method='dSPM',
-                                   lambda2=1. / 9., verbose=None):
+                                   lambda2=1. / 9., cov=None, snr=None,
+                                   verbose=None):
     """Compute resolution matrix for linear inverse operator.
 
     Parameters
@@ -31,6 +32,10 @@ def make_inverse_resolution_matrix(forward, inverse_operator, method='dSPM',
         Inverse method to use (MNE, dSPM, sLORETA).
     lambda2 : float
         The regularisation parameter.
+    cov : None | instance of Covariance
+        Noise covariance matrix to compute noise power in source space.
+    snr : None | float
+        Signal-to-noise ratio for source signal vs noise power.
     %(verbose)s
 
     Returns
@@ -40,9 +45,16 @@ def make_inverse_resolution_matrix(forward, inverse_operator, method='dSPM',
         The result of applying the inverse operator to the forward operator.
         If source orientations are not fixed, all source components will be
         computed (i.e. for n_orient_inv > 1 or n_orient_fwd > 1).
+        If 'cov' and 'snr' are not None, then 'resmat' will be of shape
+        (n_dipoles, n_orient_fwd * n_dipoles) with values pooled across
+        n_orient_inv source orientations per location.
         The columns of the resolution matrix are the point-spread functions
         (PSFs) and the rows are the cross-talk functions (CTFs).
     """
+    if (cov is None and snr is not None) or (cov is not None and snr is None):
+        msg = 'cov and snr must both be None or not be None.'
+        raise ValueError(msg)
+
     # make sure forward and inverse operator match
     inv = inverse_operator
     fwd = _convert_forward_match_inv(forward, inv)
@@ -60,6 +72,20 @@ def make_inverse_resolution_matrix(forward, inverse_operator, method='dSPM',
                                                method=method, lambda2=lambda2)
     resmat = invmat.dot(leadfield)
     logger.info('Dimensions of resolution matrix: %d by %d.' % resmat.shape)
+
+    if cov is not None:
+        print('\n###\nAdding noise to resolution matrix.\n###')
+        info_inv = _prepare_info(inverse_operator)
+        stc = apply_inverse_cov(
+            cov, info_inv, inverse_operator, nave=1, lambda2=lambda2,
+            method=method, pick_ori=None, prepared=False, label=None,
+            method_params=None, use_cps=True, verbose=None)
+
+        # Add square root of source power per column
+        print(stc.data.shape)
+        resmat = np.zeros(resmat.shape)
+        resmat += np.sqrt(stc.data)
+
     return resmat
 
 

From c40451d9de099ed98e7474d9bf620cd5eab8e54c Mon Sep 17 00:00:00 2001
From: olafhauk <olaf.hauk@mrc-cbu.cam.ac.uk>
Date: Thu, 25 Aug 2022 15:21:12 +0100
Subject: [PATCH 02/27] WIP: first reasonable results

---
 doc/changes/latest.inc                |  2 +
 mne/minimum_norm/resolution_matrix.py | 85 ++++++++++++++++++++-------
 2 files changed, 67 insertions(+), 20 deletions(-)

diff --git a/doc/changes/latest.inc b/doc/changes/latest.inc
index d095c0e571a..d540c5d96c7 100644
--- a/doc/changes/latest.inc
+++ b/doc/changes/latest.inc
@@ -27,6 +27,8 @@ Enhancements
 
 - Add ``starting_affine`` keyword argument to :func:`mne.transforms.compute_volume_registration` to initialize an alignment with an affine (:gh:`11020` by `Alex Rockhill`_)
 
+- Added noisy point-spread functions to mne.minimum_norm.make_inverse_resolution_matrix (:gh:`11058` by `Olaf Hauk`)
+
 Bugs
 ~~~~
 - Fix bug in :func:`mne.decoding.TimeFrequency` that prevented cloning if constructor arguments were modified (:gh:`11004` by :newcontrib:`Daniel Carlström Schad`)
diff --git a/mne/minimum_norm/resolution_matrix.py b/mne/minimum_norm/resolution_matrix.py
index 1d21d0538d2..79ed1d71586 100644
--- a/mne/minimum_norm/resolution_matrix.py
+++ b/mne/minimum_norm/resolution_matrix.py
@@ -7,18 +7,20 @@
 
 import numpy as np
 
-from .. import pick_channels_forward, EvokedArray, SourceEstimate
+from .. import pick_types, pick_channels_forward, EvokedArray, SourceEstimate
 from ..io.constants import FIFF
 from ..utils import logger, verbose
 from ..forward.forward import convert_forward_solution
-from ..minimum_norm import apply_inverse, apply_inverse_cov
+from ..minimum_norm import (apply_inverse, apply_inverse_cov,
+                            prepare_inverse_operator)
+from ..minimum_norm.spatial_resolution import _rectify_resolution_matrix
 from ..source_estimate import _prepare_label_extraction
 from ..label import Label
 
 
 @verbose
 def make_inverse_resolution_matrix(forward, inverse_operator, method='dSPM',
-                                   lambda2=1. / 9., cov=None, snr=None,
+                                   lambda2=1. / 9., noise_cov=None, snr=None,
                                    verbose=None):
     """Compute resolution matrix for linear inverse operator.
 
@@ -32,7 +34,7 @@ def make_inverse_resolution_matrix(forward, inverse_operator, method='dSPM',
         Inverse method to use (MNE, dSPM, sLORETA).
     lambda2 : float
         The regularisation parameter.
-    cov : None | instance of Covariance
+    noise_cov : None | instance of Covariance
         Noise covariance matrix to compute noise power in source space.
     snr : None | float
         Signal-to-noise ratio for source signal vs noise power.
@@ -41,18 +43,25 @@ def make_inverse_resolution_matrix(forward, inverse_operator, method='dSPM',
     Returns
     -------
     resmat: array, shape (n_orient_inv * n_dipoles, n_orient_fwd * n_dipoles)
-        Resolution matrix (inverse operator times forward operator).
+        Resolution matrix (inverse operator times forward operator). The
+        columns of the resolution matrix are the point-spread functions (PSFs)
+        and the rows are the cross-talk functions (CTFs).
         The result of applying the inverse operator to the forward operator.
         If source orientations are not fixed, all source components will be
         computed (i.e. for n_orient_inv > 1 or n_orient_fwd > 1).
-        If 'cov' and 'snr' are not None, then 'resmat' will be of shape
+        If 'noise_cov' and 'snr' are not None, then source noise estimated from
+        noise covariance matrix will be added to columns of the resolution
+        matrix (PSFs). In this case, 'resmat' will be of shape
         (n_dipoles, n_orient_fwd * n_dipoles) with values pooled across
         n_orient_inv source orientations per location.
-        The columns of the resolution matrix are the point-spread functions
-        (PSFs) and the rows are the cross-talk functions (CTFs).
+        Note: If the resolution matrix is computed with a noise covariance
+        matrix then only is columns, i.e. PSFs, can meaningfully be
+        interpreted. It must not be used to compute CTFs or resolution metrics
+        for CTFs!
     """
-    if (cov is None and snr is not None) or (cov is not None and snr is None):
-        msg = 'cov and snr must both be None or not be None.'
+    if ((noise_cov is None and snr is not None) or
+            (noise_cov is not None and snr is None)):
+        msg = 'noise_cov and snr must both be None or not be None.'
         raise ValueError(msg)
 
     # make sure forward and inverse operator match
@@ -71,20 +80,56 @@ def make_inverse_resolution_matrix(forward, inverse_operator, method='dSPM',
     invmat = _get_matrix_from_inverse_operator(inv, fwd,
                                                method=method, lambda2=lambda2)
     resmat = invmat.dot(leadfield)
+    resmat_ori = resmat.copy()
     logger.info('Dimensions of resolution matrix: %d by %d.' % resmat.shape)
 
-    if cov is not None:
-        print('\n###\nAdding noise to resolution matrix.\n###')
-        info_inv = _prepare_info(inverse_operator)
+    # add source noise power to columns of resolution matrix
+    if noise_cov is not None:
+        inv = prepare_inverse_operator(inv, 1, lambda2, method, copy='non-src')
+        info = _prepare_info(inv)
+        # compute source noise poweralph
         stc = apply_inverse_cov(
-            cov, info_inv, inverse_operator, nave=1, lambda2=lambda2,
-            method=method, pick_ori=None, prepared=False, label=None,
+            noise_cov, info, inv, nave=1, lambda2=lambda2,
+            method=method, pick_ori=None, prepared=True, label=None,
             method_params=None, use_cps=True, verbose=None)
 
-        # Add square root of source power per column
-        print(stc.data.shape)
-        resmat = np.zeros(resmat.shape)
-        resmat += np.sqrt(stc.data)
+        # output will be square-root of power, so take intensities before
+        # adding noise
+        shape = resmat.shape
+        if not shape[0] == shape[1]:
+            # pool loose/free orientations
+            resmat = _rectify_resolution_matrix(resmat)
+        else:
+            resmat = np.abs(resmat)
+
+        # scaling factor
+
+        idx_ch_types = {}
+        idx_ch_types['eeg'] = pick_types(info, eeg=True, meg=False)
+        idx_ch_types['mag'] = pick_types(info, eeg=False, meg='mag')
+        idx_ch_types['gra'] = pick_types(info, eeg=False, meg='grad')
+
+        # scaling similar to Samuelsson et al., Neuroimage 2021 (p. 4)
+        alphas = {}
+        for cht in idx_ch_types:
+            # signal intensity of all sources across channels of this type
+            alphas[cht] = (
+                np.abs(leadfield[idx_ch_types[cht], :]).mean())
+            # divided by noise intensity across channels of this type
+            alphas[cht] /= (
+                np.sqrt(np.diag(noise_cov.data)[idx_ch_types[cht]].mean()))
+
+        alpha = (alphas['eeg'] + alphas['mag'] + alphas['gra']) / (3 * snr)
+
+        # whitener = inv['whitener']
+        # gain = np.dot(whitener, leadfield)
+        # gram_gain = np.dot(gain, gain.T)
+        # # alpha = np.sqrt(np.trace(gram_gain) / np.trace(noise_cov.data))
+        # alpha = np.sqrt(np.trace(gram_gain))
+
+        # Add square root of source noise power to every column, scale
+        # depending on SNR
+        resmat += alpha * np.sqrt(stc.data)
 
     return resmat
 
@@ -368,7 +413,7 @@ def _get_matrix_from_inverse_operator(inverse_operator, forward, method='dSPM',
         The inverse operator.
     forward : instance of Forward
         The forward operator.
-    method : 'MNE' | 'dSPM' | 'sLORETA'
+    method : 'MNE' | 'dSPM' | 'sLORETA' | 'eLORETA'
         Inverse methods (for apply_inverse).
     lambda2 : float
         The regularization parameter (for apply_inverse).

From 5ca0c930eb43dede51413f90f2a1d355ee4cc8e3 Mon Sep 17 00:00:00 2001
From: olafhauk <olaf.hauk@mrc-cbu.cam.ac.uk>
Date: Thu, 25 Aug 2022 16:24:48 +0100
Subject: [PATCH 03/27] WIP: added tests for noisy resmat

---
 .../tests/test_resolution_matrix.py           | 55 +++++++++++++++++++
 1 file changed, 55 insertions(+)

diff --git a/mne/minimum_norm/tests/test_resolution_matrix.py b/mne/minimum_norm/tests/test_resolution_matrix.py
index 2ed97fd454f..1667ff1ab58 100644
--- a/mne/minimum_norm/tests/test_resolution_matrix.py
+++ b/mne/minimum_norm/tests/test_resolution_matrix.py
@@ -167,3 +167,58 @@ def test_resolution_matrix():
         stc_psf_label.data, stc_psf_label2[1].data)
     assert_array_equal(
         stc_psf_label.data, stc_psf_idx.data)
+
+    # test noisy resolution matrix
+    # resmats without noise
+    rm_mne1 = make_inverse_resolution_matrix(
+        forward_fxd, inverse_operator_fxd, method='MNE', lambda2=lambda2,
+        noise_cov=None, snr=None)
+    rm_mne2 = make_inverse_resolution_matrix(
+        forward_fxd, inverse_operator_fxd, method='dSPM', lambda2=lambda2,
+        noise_cov=None, snr=None)
+    rm_mne3 = make_inverse_resolution_matrix(
+        forward_fxd, inverse_operator_fxd, method='sLORETA', lambda2=lambda2,
+        noise_cov=None, snr=None)
+    rm_mne4 = make_inverse_resolution_matrix(
+        forward_fxd, inverse_operator_fxd, method='eLORETA', lambda2=lambda2,
+        noise_cov=None, snr=None)
+
+    # same resmats with noise but very high SNR
+    rm_mne1a = make_inverse_resolution_matrix(
+        forward_fxd, inverse_operator_fxd, method='MNE', lambda2=lambda2,
+        noise_cov=noise_cov, snr=1e30)
+    rm_mne2a = make_inverse_resolution_matrix(
+        forward_fxd, inverse_operator_fxd, method='dSPM', lambda2=lambda2,
+        noise_cov=noise_cov, snr=1e30)
+    rm_mne3a = make_inverse_resolution_matrix(
+        forward_fxd, inverse_operator_fxd, method='sLORETA', lambda2=lambda2,
+        noise_cov=noise_cov, snr=1e30)
+    rm_mne4a = make_inverse_resolution_matrix(
+        forward_fxd, inverse_operator_fxd, method='eLORETA', lambda2=lambda2,
+        noise_cov=noise_cov, snr=1e30)
+
+    # resmats should be very close without noise and at high SNR
+    assert_array_almost_equal(np.abs(rm_mne1), rm_mne1a)
+    assert_array_almost_equal(np.abs(rm_mne2), rm_mne2a)
+    assert_array_almost_equal(np.abs(rm_mne3), rm_mne3a)
+    assert_array_almost_equal(np.abs(rm_mne4), rm_mne4a)
+
+    # rows of resmats (CTFs) should be similar between methods
+    row_corrs = [np.corrcoef(rm_mne1[i, :], rm_mne2[i, :])[0, 1]
+                 for i in np.arange(rm_mne1.shape[0])]
+    assert_allclose(row_corrs, 1)
+    row_corrs = [np.corrcoef(rm_mne1[i, :], rm_mne3[i, :])[0, 1]
+                 for i in np.arange(rm_mne1.shape[0])]
+    assert_allclose(row_corrs, 1)
+    row_corrs = [np.corrcoef(rm_mne1[i, :], rm_mne4[i, :])[0, 1]
+                 for i in np.arange(rm_mne1.shape[0])]
+    assert_allclose(row_corrs, 1, 0.2)  # eLORETA is different
+    row_corrs = [np.corrcoef(rm_mne1a[i, :], rm_mne2a[i, :])[0, 1]
+                 for i in np.arange(rm_mne1.shape[0])]
+    assert_allclose(row_corrs, 1)
+    row_corrs = [np.corrcoef(rm_mne1a[i, :], rm_mne3a[i, :])[0, 1]
+                 for i in np.arange(rm_mne1.shape[0])]
+    assert_allclose(row_corrs, 1)
+    row_corrs = [np.corrcoef(rm_mne1a[i, :], rm_mne4a[i, :])[0, 1]
+                 for i in np.arange(rm_mne1.shape[0])]
+    assert_allclose(row_corrs, 1, 0.2)  # eLORETA is different

From 89ea2763bc3161572c5d82a850e1dea6a9867b91 Mon Sep 17 00:00:00 2001
From: olafhauk <olaf.hauk@mrc-cbu.cam.ac.uk>
Date: Thu, 25 Aug 2022 16:31:09 +0100
Subject: [PATCH 04/27] WIP: removed unnecessary line

---
 mne/minimum_norm/resolution_matrix.py | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/mne/minimum_norm/resolution_matrix.py b/mne/minimum_norm/resolution_matrix.py
index 79ed1d71586..f1bdaa7cb86 100644
--- a/mne/minimum_norm/resolution_matrix.py
+++ b/mne/minimum_norm/resolution_matrix.py
@@ -80,7 +80,6 @@ def make_inverse_resolution_matrix(forward, inverse_operator, method='dSPM',
     invmat = _get_matrix_from_inverse_operator(inv, fwd,
                                                method=method, lambda2=lambda2)
     resmat = invmat.dot(leadfield)
-    resmat_ori = resmat.copy()
     logger.info('Dimensions of resolution matrix: %d by %d.' % resmat.shape)
 
     # add source noise power to columns of resolution matrix
@@ -121,10 +120,10 @@ def make_inverse_resolution_matrix(forward, inverse_operator, method='dSPM',
 
         alpha = (alphas['eeg'] + alphas['mag'] + alphas['gra']) / (3 * snr)
 
+        # # previous attempt to compute scaling via whitened gain matrix
         # whitener = inv['whitener']
         # gain = np.dot(whitener, leadfield)
         # gram_gain = np.dot(gain, gain.T)
-        # # alpha = np.sqrt(np.trace(gram_gain) / np.trace(noise_cov.data))
         # alpha = np.sqrt(np.trace(gram_gain))
 
         # Add square root of source noise power to every column, scale

From fa0c24407069bbb51ba81a93463b3e3abed4e4f1 Mon Sep 17 00:00:00 2001
From: olafhauk <olaf.hauk@mrc-cbu.cam.ac.uk>
Date: Fri, 26 Aug 2022 16:09:42 +0100
Subject: [PATCH 05/27] WIP: changed scaling of noisy resmat

---
 mne/minimum_norm/resolution_matrix.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/mne/minimum_norm/resolution_matrix.py b/mne/minimum_norm/resolution_matrix.py
index f1bdaa7cb86..82af0f58648 100644
--- a/mne/minimum_norm/resolution_matrix.py
+++ b/mne/minimum_norm/resolution_matrix.py
@@ -118,17 +118,17 @@ def make_inverse_resolution_matrix(forward, inverse_operator, method='dSPM',
             alphas[cht] /= (
                 np.sqrt(np.diag(noise_cov.data)[idx_ch_types[cht]].mean()))
 
-        alpha = (alphas['eeg'] + alphas['mag'] + alphas['gra']) / (3 * snr)
+        alpha = (3 * snr) / (alphas['eeg'] + alphas['mag'] + alphas['gra'])
 
         # # previous attempt to compute scaling via whitened gain matrix
         # whitener = inv['whitener']
         # gain = np.dot(whitener, leadfield)
         # gram_gain = np.dot(gain, gain.T)
-        # alpha = np.sqrt(np.trace(gram_gain))
+        # alpha = np.sqrt(np.trace(gram_gain) / gram_gain.shape[0])
 
         # Add square root of source noise power to every column, scale
         # depending on SNR
-        resmat += alpha * np.sqrt(stc.data)
+        resmat = alpha * resmat + np.sqrt(stc.data)
 
     return resmat
 

From 09ed30f60692dab1a7a02009e3c99940b9a47e6c Mon Sep 17 00:00:00 2001
From: olafhauk <olaf.hauk@mrc-cbu.cam.ac.uk>
Date: Fri, 26 Aug 2022 16:24:01 +0100
Subject: [PATCH 06/27] WIP: linked snr to lambda2 in API

---
 mne/minimum_norm/resolution_matrix.py | 11 ++++++++---
 1 file changed, 8 insertions(+), 3 deletions(-)

diff --git a/mne/minimum_norm/resolution_matrix.py b/mne/minimum_norm/resolution_matrix.py
index 82af0f58648..abdd1b577c0 100644
--- a/mne/minimum_norm/resolution_matrix.py
+++ b/mne/minimum_norm/resolution_matrix.py
@@ -38,6 +38,7 @@ def make_inverse_resolution_matrix(forward, inverse_operator, method='dSPM',
         Noise covariance matrix to compute noise power in source space.
     snr : None | float
         Signal-to-noise ratio for source signal vs noise power.
+        If None, SNR is inferred from lambda2 (as 1/sqrt(lambda2)).
     %(verbose)s
 
     Returns
@@ -59,11 +60,15 @@ def make_inverse_resolution_matrix(forward, inverse_operator, method='dSPM',
         interpreted. It must not be used to compute CTFs or resolution metrics
         for CTFs!
     """
-    if ((noise_cov is None and snr is not None) or
-            (noise_cov is not None and snr is None)):
-        msg = 'noise_cov and snr must both be None or not be None.'
+    if (noise_cov is None and snr is not None):
+        msg = 'snr should be None if noise_cov is None.'
         raise ValueError(msg)
 
+    if (noise_cov is not None and snr is None):
+        snr = np.sqrt(1. / lambda2)
+        logger.info('Inferring snr from lamdba2 as: %f.' % snr)
+
+
     # make sure forward and inverse operator match
     inv = inverse_operator
     fwd = _convert_forward_match_inv(forward, inv)

From 77b34704c5f0540987cb85ad5abb1279341c0da5 Mon Sep 17 00:00:00 2001
From: olafhauk <olaf.hauk@mrc-cbu.cam.ac.uk>
Date: Fri, 26 Aug 2022 19:20:31 +0100
Subject: [PATCH 07/27] WIP: added example, plus minor changes

---
 mne/minimum_norm/resolution_matrix.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/mne/minimum_norm/resolution_matrix.py b/mne/minimum_norm/resolution_matrix.py
index abdd1b577c0..cb9961acfbf 100644
--- a/mne/minimum_norm/resolution_matrix.py
+++ b/mne/minimum_norm/resolution_matrix.py
@@ -118,7 +118,7 @@ def make_inverse_resolution_matrix(forward, inverse_operator, method='dSPM',
         for cht in idx_ch_types:
             # signal intensity of all sources across channels of this type
             alphas[cht] = (
-                np.abs(leadfield[idx_ch_types[cht], :]).mean())
+                np.sqrt((leadfield[idx_ch_types[cht], :]**2).mean()))
             # divided by noise intensity across channels of this type
             alphas[cht] /= (
                 np.sqrt(np.diag(noise_cov.data)[idx_ch_types[cht]].mean()))

From 0f30564b6d3987a27a20d3e8d3569052cca4ca0a Mon Sep 17 00:00:00 2001
From: olafhauk <olaf.hauk@mrc-cbu.cam.ac.uk>
Date: Mon, 13 Nov 2023 11:16:55 +0000
Subject: [PATCH 08/27] WIP: corrected typo

---
 mne/minimum_norm/resolution_matrix.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/mne/minimum_norm/resolution_matrix.py b/mne/minimum_norm/resolution_matrix.py
index f77a479d710..f02c24fcf96 100755
--- a/mne/minimum_norm/resolution_matrix.py
+++ b/mne/minimum_norm/resolution_matrix.py
@@ -61,7 +61,7 @@ def make_inverse_resolution_matrix(forward, inverse_operator, method='dSPM',
         (n_dipoles, n_orient_fwd * n_dipoles) with values pooled across
         n_orient_inv source orientations per location.
         Note: If the resolution matrix is computed with a noise covariance
-        matrix then only is columns, i.e. PSFs, can meaningfully be
+        matrix then only its columns, i.e. PSFs, can meaningfully be
         interpreted. It must not be used to compute CTFs or resolution metrics
         for CTFs!
     """

From 6f82bfdbce128a11a9ec303ad9fe3c355b7ed1cc Mon Sep 17 00:00:00 2001
From: olafhauk <olaf.hauk@mrc-cbu.cam.ac.uk>
Date: Mon, 13 Nov 2023 16:44:29 +0000
Subject: [PATCH 09/27] WIP: more typos, doc style

---
 mne/minimum_norm/resolution_matrix.py | 9 +++------
 1 file changed, 3 insertions(+), 6 deletions(-)

diff --git a/mne/minimum_norm/resolution_matrix.py b/mne/minimum_norm/resolution_matrix.py
index f02c24fcf96..eac9e73710d 100755
--- a/mne/minimum_norm/resolution_matrix.py
+++ b/mne/minimum_norm/resolution_matrix.py
@@ -7,12 +7,10 @@
 
 import numpy as np
 
-<<<<<<< HEAD
 from .. import pick_types, pick_channels_forward, EvokedArray, SourceEstimate
-=======
 from .. import (pick_channels_forward, EvokedArray, SourceEstimate,
                 VectorSourceEstimate)
->>>>>>> master
+
 from ..io.constants import FIFF
 from ..utils import logger, verbose
 from ..forward.forward import convert_forward_solution
@@ -52,11 +50,10 @@ def make_inverse_resolution_matrix(forward, inverse_operator, method='dSPM',
         Resolution matrix (inverse operator times forward operator). The
         columns of the resolution matrix are the point-spread functions (PSFs)
         and the rows are the cross-talk functions (CTFs).
-        The result of applying the inverse operator to the forward operator.
         If source orientations are not fixed, all source components will be
         computed (i.e. for n_orient_inv > 1 or n_orient_fwd > 1).
         If 'noise_cov' and 'snr' are not None, then source noise estimated from
-        noise covariance matrix will be added to columns of the resolution
+        the noise covariance matrix will be added to columns of the resolution
         matrix (PSFs). In this case, 'resmat' will be of shape
         (n_dipoles, n_orient_fwd * n_dipoles) with values pooled across
         n_orient_inv source orientations per location.
@@ -96,7 +93,7 @@ def make_inverse_resolution_matrix(forward, inverse_operator, method='dSPM',
     if noise_cov is not None:
         inv = prepare_inverse_operator(inv, 1, lambda2, method, copy='non-src')
         info = _prepare_info(inv)
-        # compute source noise poweralph
+        # compute source noise power
         stc = apply_inverse_cov(
             noise_cov, info, inv, nave=1, lambda2=lambda2,
             method=method, pick_ori=None, prepared=True, label=None,

From 840de7f921cb1b9158bd008582dab0d0f7339276 Mon Sep 17 00:00:00 2001
From: Marijn van Vliet <w.m.vanvliet@gmail.com>
Date: Wed, 15 Nov 2023 12:03:37 +0200
Subject: [PATCH 10/27] typo

---
 mne/minimum_norm/resolution_matrix.py | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/mne/minimum_norm/resolution_matrix.py b/mne/minimum_norm/resolution_matrix.py
index 74c4f45acce..b9160258645 100755
--- a/mne/minimum_norm/resolution_matrix.py
+++ b/mne/minimum_norm/resolution_matrix.py
@@ -7,8 +7,7 @@
 import numpy as np
 
 from .._fiff.constants import FIFF
-from .._fiff.pick import pick_channels_forward
-from ..channels import pick_types
+from .._fiff.pick import pick_channels_forward, pick_types
 from ..evoked import EvokedArray
 from ..forward.forward import Forward, convert_forward_solution
 from ..label import Label

From 415a4b9830271542dfafda25233fcf98c4f33bab Mon Sep 17 00:00:00 2001
From: olafhauk <olaf.hauk@mrc-cbu.cam.ac.uk>
Date: Wed, 15 Nov 2023 16:38:27 +0000
Subject: [PATCH 11/27] WIP: added tests

---
 mne/minimum_norm/resolution_matrix.py         | 354 ++++++--------
 .../tests/test_resolution_matrix.py           | 461 ++++++++----------
 2 files changed, 333 insertions(+), 482 deletions(-)

diff --git a/mne/minimum_norm/resolution_matrix.py b/mne/minimum_norm/resolution_matrix.py
index b9160258645..755d637b58c 100755
--- a/mne/minimum_norm/resolution_matrix.py
+++ b/mne/minimum_norm/resolution_matrix.py
@@ -1,3 +1,4 @@
+# -*- coding: utf-8 -*-
 """Compute resolution matrix for linear estimators."""
 # Authors: olaf.hauk@mrc-cbu.cam.ac.uk
 #
@@ -6,33 +7,24 @@
 
 import numpy as np
 
-from .._fiff.constants import FIFF
-from .._fiff.pick import pick_channels_forward, pick_types
-from ..evoked import EvokedArray
-from ..forward.forward import Forward, convert_forward_solution
-from ..label import Label
-from ..minimum_norm.inverse import (
-    InverseOperator,
-    apply_inverse_cov,
-    prepare_inverse_operator,
-)
+from .. import pick_types, pick_channels_forward, EvokedArray, SourceEstimate
+from .. import (pick_channels_forward, EvokedArray, SourceEstimate,
+                VectorSourceEstimate)
+
+from ..io.constants import FIFF
+from ..utils import logger, verbose
+from ..forward.forward import convert_forward_solution
+from ..minimum_norm import (apply_inverse, apply_inverse_cov,
+                            prepare_inverse_operator)
 from ..minimum_norm.spatial_resolution import _rectify_resolution_matrix
-from ..source_estimate import _get_src_type, _make_stc, _prepare_label_extraction
-from ..source_space._source_space import SourceSpaces, _get_vertno
-from ..utils import _validate_type, logger, verbose
-from .inverse import apply_inverse
+from ..source_estimate import _prepare_label_extraction
+from ..label import Label
 
 
 @verbose
-def make_inverse_resolution_matrix(
-    forward,
-    inverse_operator,
-    method="dSPM",
-    lambda2=1.0 / 9.0,
-    noise_cov=None,
-    snr=None,
-    verbose=None,
-):
+def make_inverse_resolution_matrix(forward, inverse_operator, method='dSPM',
+                                   lambda2=1. / 9., noise_cov=None, snr=None,
+                                   verbose=None):
     """Compute resolution matrix for linear inverse operator.
 
     Parameters
@@ -70,13 +62,14 @@ def make_inverse_resolution_matrix(
         interpreted. It must not be used to compute CTFs or resolution metrics
         for CTFs!
     """
-    if noise_cov is None and snr is not None:
-        msg = "snr should be None if noise_cov is None."
+    if (noise_cov is None and snr is not None):
+        msg = 'snr should be None if noise_cov is None.'
         raise ValueError(msg)
 
-    if noise_cov is not None and snr is None:
-        snr = np.sqrt(1.0 / lambda2)
-        logger.info("Inferring snr from lamdba2 as: %f." % snr)
+    if (noise_cov is not None and snr is None):
+        snr = np.sqrt(1. / lambda2)
+        logger.info('Inferring snr from lamdba2 as: %f.' % snr)
+
 
     # make sure forward and inverse operator match
     inv = inverse_operator
@@ -84,36 +77,27 @@ def make_inverse_resolution_matrix(
 
     # don't include bad channels
     # only use good channels from inverse operator
-    bads_inv = inv["info"]["bads"]
+    bads_inv = inv['info']['bads']
     # good channels
-    ch_names = [c for c in inv["info"]["ch_names"] if (c not in bads_inv)]
+    ch_names = [c for c in inv['info']['ch_names'] if (c not in bads_inv)]
     fwd = pick_channels_forward(fwd, ch_names, ordered=True)
 
     # get leadfield matrix from forward solution
-    leadfield = fwd["sol"]["data"]
-    invmat = _get_matrix_from_inverse_operator(inv, fwd, method=method, lambda2=lambda2)
+    leadfield = fwd['sol']['data']
+    invmat = _get_matrix_from_inverse_operator(inv, fwd,
+                                               method=method, lambda2=lambda2)
     resmat = invmat.dot(leadfield)
-    logger.info("Dimensions of resolution matrix: %d by %d." % resmat.shape)
+    logger.info('Dimensions of resolution matrix: %d by %d.' % resmat.shape)
 
     # add source noise power to columns of resolution matrix
     if noise_cov is not None:
-        inv = prepare_inverse_operator(inv, 1, lambda2, method, copy="non-src")
+        inv = prepare_inverse_operator(inv, 1, lambda2, method, copy='non-src')
         info = _prepare_info(inv)
         # compute source noise power
         stc = apply_inverse_cov(
-            noise_cov,
-            info,
-            inv,
-            nave=1,
-            lambda2=lambda2,
-            method=method,
-            pick_ori=None,
-            prepared=True,
-            label=None,
-            method_params=None,
-            use_cps=True,
-            verbose=None,
-        )
+            noise_cov, info, inv, nave=1, lambda2=lambda2,
+            method=method, pick_ori=None, prepared=True, label=None,
+            method_params=None, use_cps=True, verbose=None)
 
         # output will be square-root of power, so take intensities before
         # adding noise
@@ -127,19 +111,21 @@ def make_inverse_resolution_matrix(
         # scaling factor
 
         idx_ch_types = {}
-        idx_ch_types["eeg"] = pick_types(info, eeg=True, meg=False)
-        idx_ch_types["mag"] = pick_types(info, eeg=False, meg="mag")
-        idx_ch_types["gra"] = pick_types(info, eeg=False, meg="grad")
+        idx_ch_types['eeg'] = pick_types(info, eeg=True, meg=False)
+        idx_ch_types['mag'] = pick_types(info, eeg=False, meg='mag')
+        idx_ch_types['gra'] = pick_types(info, eeg=False, meg='grad')
 
         # scaling similar to Samuelsson et al., Neuroimage 2021 (p. 4)
         alphas = {}
         for cht in idx_ch_types:
             # signal intensity of all sources across channels of this type
-            alphas[cht] = np.sqrt((leadfield[idx_ch_types[cht], :] ** 2).mean())
+            alphas[cht] = (
+                np.sqrt((leadfield[idx_ch_types[cht], :]**2).mean()))
             # divided by noise intensity across channels of this type
-            alphas[cht] /= np.sqrt(np.diag(noise_cov.data)[idx_ch_types[cht]].mean())
+            alphas[cht] /= (
+                np.sqrt(np.diag(noise_cov.data)[idx_ch_types[cht]].mean()))
 
-        alpha = (3 * snr) / (alphas["eeg"] + alphas["mag"] + alphas["gra"])
+        alpha = (3 * snr) / (alphas['eeg'] + alphas['mag'] + alphas['gra'])
 
         # # previous attempt to compute scaling via whitened gain matrix
         # whitener = inv['whitener']
@@ -149,66 +135,63 @@ def make_inverse_resolution_matrix(
 
         # Add square root of source noise power to every column, scale
         # depending on SNR
-        resmat = alpha * resmat + np.sqrt(stc.data)
+        # EDIT
+        resmat = resmat + (1 / alpha) * np.sqrt(stc.data)
 
     return resmat
 
 
 @verbose
-def _get_psf_ctf(
-    resmat,
-    src,
-    idx,
-    *,
-    func,
-    mode,
-    n_comp,
-    norm,
-    return_pca_vars,
-    vector=False,
-    verbose=None,
-):
-    """Get point-spread (PSFs) or cross-talk (CTFs) functions."""
+def _get_psf_ctf(resmat, src, idx, *, func, mode, n_comp, norm,
+                 return_pca_vars, vector=False, verbose=None):
+    """Get point-spread (PSFs) or cross-talk (CTFs) functions.
+
+    Parameters
+    ----------
+    resmat : array, shape (n_orient_inv * n_dipoles, n_orient_fwd * n_dipoles)
+        Forward Operator.
+    src : Source Space
+        Source space used to compute resolution matrix.
+    %(idx_pctf)s
+    func : str ('psf' | 'ctf')
+        Whether to produce PSFs or CTFs. Defaults to psf.
+    %(mode_pctf)s
+    %(n_comp_pctf_n)s
+    %(norm_pctf)s
+    %(return_pca_vars_pctf)s
+    %(vector_pctf)s
+    %(verbose)s
+
+    Returns
+    -------
+    %(stcs_pctf)s
+    %(pca_vars_pctf)s
+    """
     # check for consistencies in input parameters
     _check_get_psf_ctf_params(mode, n_comp, return_pca_vars)
 
     # backward compatibility
     if norm is True:
-        norm = "max"
+        norm = 'max'
 
     # get relevant vertices in source space
-    src_orig = src
-    _validate_type(src_orig, (InverseOperator, Forward, SourceSpaces), "src")
-    if not isinstance(src, SourceSpaces):
-        src = src["src"]
     verts_all = _vertices_for_get_psf_ctf(idx, src)
-    vertno = _get_vertno(src)
-    n_verts = sum(len(v) for v in vertno)
-    src_type = _get_src_type(src, vertno)
-    subject = src._subject
-    if vector and src_type == "surface":
-        _validate_type(
-            src_orig,
-            (Forward, InverseOperator),
-            "src",
-            extra="when creating a vector surface source estimate",
-        )
-        nn = src_orig["source_nn"]
-    else:
-        nn = np.repeat(np.eye(3, 3)[np.newaxis], n_verts, 0)
+    vertno_lh = src[0]['vertno']
+    vertno_rh = src[1]['vertno']
+    vertno = [vertno_lh, vertno_rh]
+
+    n_verts = len(vertno[0]) + len(vertno[1])
 
     n_r, n_c = resmat.shape
-    if ((n_verts != n_r) and (n_r / 3 != n_verts)) or (
-        (n_verts != n_c) and (n_c / 3 != n_verts)
-    ):
-        msg = (
-            "Number of vertices (%d) and corresponding dimension of"
-            "resolution matrix (%d, %d) do not match" % (n_verts, n_r, n_c)
-        )
+    if (((n_verts != n_r) and (n_r / 3 != n_verts))
+            or ((n_verts != n_c) and (n_c / 3 != n_verts))):
+        msg = ('Number of vertices (%d) and corresponding dimension of'
+               'resolution matrix ((%d, %d) do not match' %
+               (n_verts, n_r, n_c))
         raise ValueError(msg)
 
     # the following will operate on columns of funcs
-    if func == "ctf":
+    if func == 'ctf':
         resmat = resmat.T
         n_r, n_c = n_c, n_r
 
@@ -223,12 +206,13 @@ def _get_psf_ctf(
         for [i, verts] in enumerate(verts_all):
             verts_vec = np.empty(3 * len(verts), dtype=int)
             for [j, v] in enumerate(verts):
-                verts_vec[3 * j : 3 * j + 3] = 3 * verts[j] + np.array([0, 1, 2])
+                verts_vec[3 * j: 3 * j + 3] = \
+                    3 * verts[j] + np.array([0, 1, 2])
             verts_all[i] = verts_vec  # use these as indices
 
     for verts in verts_all:
         # get relevant PSFs or CTFs for specified vertices
-        if isinstance(verts, int):
+        if type(verts) is int:
             verts = [verts]  # to keep array dimensions
         funcs = resmat[:, verts]
 
@@ -239,28 +223,24 @@ def _get_psf_ctf(
         # summarise PSFs/CTFs across vertices if requested
         pca_var = None  # variances computed only if return_pca_vars=True
         if mode is not None:
-            funcs, pca_var = _summarise_psf_ctf(
-                funcs, mode, n_comp, return_pca_vars, nn
-            )
+            funcs, pca_var = _summarise_psf_ctf(funcs, mode, n_comp,
+                                                return_pca_vars)
 
         if not vector:  # if one value per vertex requested
             if n_verts != n_r:  # if 3 orientations per vertex, combine
                 funcs_int = np.empty([int(n_r / 3), funcs.shape[1]])
                 for i in np.arange(0, n_verts):
-                    funcs_vert = funcs[3 * i : 3 * i + 3, :]
-                    funcs_int[i, :] = np.sqrt((funcs_vert**2).sum(axis=0))
-                funcs = funcs_int
-
-        stc = _make_stc(
-            funcs,
-            vertno,
-            src_type,
-            tmin=0.0,
-            tstep=1.0,
-            subject=subject,
-            vector=vector,
-            source_nn=nn,
-        )
+                    funcs_vert = funcs[3 * i:3 * i + 3, :]
+                    funcs_int[i, :] = np.sqrt((funcs_vert ** 2).sum(axis=0))
+                stc = SourceEstimate(funcs_int, vertno, tmin=0., tstep=1.)
+            else:  # use as is
+                stc = SourceEstimate(funcs, vertno, tmin=0., tstep=1.)
+        else:  # STC with orientations
+            # convert to vector source estimate
+            m, n = int(funcs.shape[0] / 3), int(funcs.shape[1])
+            data = funcs.reshape(m, 3, n)
+            stc = VectorSourceEstimate(data, vertno, tmin=0., tstep=1.)
+
         stcs.append(stc)
         pca_vars.append(pca_var)
 
@@ -277,11 +257,11 @@ def _get_psf_ctf(
 
 def _check_get_psf_ctf_params(mode, n_comp, return_pca_vars):
     """Check input parameters of _get_psf_ctf() for consistency."""
-    if mode in [None, "sum", "mean"] and n_comp > 1:
-        msg = "n_comp must be 1 for mode=%s." % mode
+    if mode in [None, 'sum', 'mean'] and n_comp > 1:
+        msg = 'n_comp must be 1 for mode=%s.' % mode
         raise ValueError(msg)
-    if mode != "pca" and return_pca_vars:
-        msg = "SVD variances can only be returned if mode=" "pca" "."
+    if mode != 'pca' and return_pca_vars:
+        msg = 'SVD variances can only be returned if mode=''pca''.'
         raise ValueError(msg)
 
 
@@ -292,13 +272,8 @@ def _vertices_for_get_psf_ctf(idx, src):
     if type(idx[0]) is Label:
         # specify without source time courses, gets indices per label
         verts_labs, _ = _prepare_label_extraction(
-            stc=None,
-            labels=idx,
-            src=src,
-            mode="mean",
-            allow_empty=False,
-            use_sparse=False,
-        )
+            stc=None, labels=idx, src=src, mode='mean', allow_empty=False,
+            use_sparse=False)
         # verts_labs can be list of lists
         # concatenate indices per label across hemispheres
         # one list item per label
@@ -306,7 +281,7 @@ def _vertices_for_get_psf_ctf(idx, src):
 
         for v in verts_labs:
             # if two hemispheres present
-            if isinstance(v, list):
+            if type(v) is list:
                 # indices for both hemispheres in one list
                 this_verts = np.concatenate((v[0], v[1]))
             else:
@@ -314,7 +289,7 @@ def _vertices_for_get_psf_ctf(idx, src):
             verts.append(this_verts)
     # check if list of list or just list
     else:
-        if isinstance(idx[0], list):  # if list of list of integers
+        if type(idx[0]) is list:  # if list of list of integers
             verts = idx
         else:  # if list of integers
             verts = [idx]
@@ -325,21 +300,21 @@ def _vertices_for_get_psf_ctf(idx, src):
 def _normalise_psf_ctf(funcs, norm):
     """Normalise PSFs/CTFs in _get_psf_ctf()."""
     # normalise PSFs/CTFs if specified
-    if norm == "max":
+    if norm == 'max':
         maxval = max(-funcs.min(), funcs.max())
         funcs = funcs / maxval
-    elif norm == "norm":  # normalise to maximum norm across columns
+    elif norm == 'norm':  # normalise to maximum norm across columns
         norms = np.linalg.norm(funcs, axis=0)
         funcs = funcs / norms.max()
 
     return funcs
 
 
-def _summarise_psf_ctf(funcs, mode, n_comp, return_pca_vars, nn):
+def _summarise_psf_ctf(funcs, mode, n_comp, return_pca_vars):
     """Summarise PSFs/CTFs across vertices."""
     s_var = None  # only computed for return_pca_vars=True
 
-    if mode == "maxval":  # pick PSF/CTF with maximum absolute value
+    if mode == 'maxval':  # pick PSF/CTF with maximum absolute value
         absvals = np.maximum(-np.min(funcs, axis=0), np.max(funcs, axis=0))
         if n_comp > 1:  # only keep requested number of sorted PSFs/CTFs
             sortidx = np.argsort(absvals)
@@ -348,7 +323,7 @@ def _summarise_psf_ctf(funcs, mode, n_comp, return_pca_vars, nn):
             maxidx = [absvals.argmax()]
         funcs = funcs[:, maxidx]
 
-    elif mode == "maxnorm":  # pick PSF/CTF with maximum norm
+    elif mode == 'maxnorm':  # pick PSF/CTF with maximum norm
         norms = np.linalg.norm(funcs, axis=0)
         if n_comp > 1:  # only keep requested number of sorted PSFs/CTFs
             sortidx = np.argsort(norms)
@@ -357,13 +332,13 @@ def _summarise_psf_ctf(funcs, mode, n_comp, return_pca_vars, nn):
             maxidx = [norms.argmax()]
         funcs = funcs[:, maxidx]
 
-    elif mode == "sum":  # sum across PSFs/CTFs
+    elif mode == 'sum':  # sum across PSFs/CTFs
         funcs = np.sum(funcs, axis=1, keepdims=True)
 
-    elif mode == "mean":  # mean of PSFs/CTFs
+    elif mode == 'mean':  # mean of PSFs/CTFs
         funcs = np.mean(funcs, axis=1, keepdims=True)
 
-    elif mode == "pca":  # SVD across PSFs/CTFs
+    elif mode == 'pca':  # SVD across PSFs/CTFs
         # compute SVD of PSFs/CTFs across vertices
         u, s, _ = np.linalg.svd(funcs, full_matrices=False, compute_uv=True)
         if n_comp > 1:
@@ -380,28 +355,16 @@ def _summarise_psf_ctf(funcs, mode, n_comp, return_pca_vars, nn):
 
 
 @verbose
-def get_point_spread(
-    resmat,
-    src,
-    idx,
-    mode=None,
-    *,
-    n_comp=1,
-    norm=False,
-    return_pca_vars=False,
-    vector=False,
-    verbose=None,
-):
+def get_point_spread(resmat, src, idx, mode=None, *, n_comp=1, norm=False,
+                     return_pca_vars=False, vector=False, verbose=None):
     """Get point-spread (PSFs) functions for vertices.
 
     Parameters
     ----------
     resmat : array, shape (n_dipoles, n_dipoles)
         Forward Operator.
-    src : instance of SourceSpaces | instance of InverseOperator | instance of Forward
+    src : instance of SourceSpaces
         Source space used to compute resolution matrix.
-        Must be an InverseOperator if ``vector=True`` and a surface
-        source space is used.
     %(idx_pctf)s
     %(mode_pctf)s
     %(n_comp_pctf_n)s
@@ -414,43 +377,23 @@ def get_point_spread(
     -------
     %(stcs_pctf)s
     %(pca_vars_pctf)s
-    """  # noqa: E501
-    return _get_psf_ctf(
-        resmat,
-        src,
-        idx,
-        func="psf",
-        mode=mode,
-        n_comp=n_comp,
-        norm=norm,
-        return_pca_vars=return_pca_vars,
-        vector=vector,
-    )
+    """
+    return _get_psf_ctf(resmat, src, idx, func='psf', mode=mode, n_comp=n_comp,
+                        norm=norm, return_pca_vars=return_pca_vars,
+                        vector=vector)
 
 
 @verbose
-def get_cross_talk(
-    resmat,
-    src,
-    idx,
-    mode=None,
-    *,
-    n_comp=1,
-    norm=False,
-    return_pca_vars=False,
-    vector=False,
-    verbose=None,
-):
+def get_cross_talk(resmat, src, idx, mode=None, *, n_comp=1, norm=False,
+                   return_pca_vars=False, vector=False, verbose=None):
     """Get cross-talk (CTFs) function for vertices.
 
     Parameters
     ----------
     resmat : array, shape (n_dipoles, n_dipoles)
         Forward Operator.
-    src : instance of SourceSpaces | instance of InverseOperator | instance of Forward
+    src : instance of SourceSpaces
         Source space used to compute resolution matrix.
-        Must be an InverseOperator if ``vector=True`` and a surface
-        source space is used.
     %(idx_pctf)s
     %(mode_pctf)s
     %(n_comp_pctf_n)s
@@ -463,18 +406,10 @@ def get_cross_talk(
     -------
     %(stcs_pctf)s
     %(pca_vars_pctf)s
-    """  # noqa: E501
-    return _get_psf_ctf(
-        resmat,
-        src,
-        idx,
-        func="ctf",
-        mode=mode,
-        n_comp=n_comp,
-        norm=norm,
-        return_pca_vars=return_pca_vars,
-        vector=vector,
-    )
+    """
+    return _get_psf_ctf(resmat, src, idx, func='ctf', mode=mode, n_comp=n_comp,
+                        norm=norm, return_pca_vars=return_pca_vars,
+                        vector=vector)
 
 
 def _convert_forward_match_inv(fwd, inv):
@@ -483,8 +418,6 @@ def _convert_forward_match_inv(fwd, inv):
     Inverse operator and forward operator must have same surface orientations,
     but can have different source orientation constraints.
     """
-    _validate_type(fwd, Forward, "fwd")
-    _validate_type(inv, InverseOperator, "inverse_operator")
     # did inverse operator use fixed orientation?
     is_fixed_inv = _check_fixed_ori(inv)
     # did forward operator use fixed orientation?
@@ -494,11 +427,13 @@ def _convert_forward_match_inv(fwd, inv):
     # if inv loose: surf_ori must be True
     # if inv free: surf_ori must be False
     if not is_fixed_inv and not is_fixed_fwd:
-        inv_surf_ori = inv._is_surf_ori
-        if inv_surf_ori != fwd["surf_ori"]:
-            fwd = convert_forward_solution(
-                fwd, surf_ori=inv_surf_ori, force_fixed=False
-            )
+        is_loose_inv = not (inv['orient_prior']['data'] == 1.).all()
+
+        if is_loose_inv:
+            if not fwd['surf_ori']:
+                fwd = convert_forward_solution(fwd, surf_ori=True)
+        elif fwd['surf_ori']:  # free orientation, change fwd
+            fwd = convert_forward_solution(fwd, surf_ori=False)
 
     return fwd
 
@@ -508,17 +443,15 @@ def _prepare_info(inverse_operator):
     # in order to convert sub-leadfield matrix to evoked data type (pretending
     # it's an epoch, see in loop below), uses 'info' from inverse solution
     # because this has all the correct projector information
-    info = deepcopy(inverse_operator["info"])
+    info = deepcopy(inverse_operator['info'])
     with info._unlock():
-        info["sfreq"] = 1000.0  # necessary
-        info["projs"] = inverse_operator["projs"]
-        info["custom_ref_applied"] = False
+        info['sfreq'] = 1000.  # necessary
+        info['projs'] = inverse_operator['projs']
     return info
 
 
-def _get_matrix_from_inverse_operator(
-    inverse_operator, forward, method="dSPM", lambda2=1.0 / 9.0
-):
+def _get_matrix_from_inverse_operator(inverse_operator, forward, method='dSPM',
+                                      lambda2=1. / 9.):
     """Get inverse matrix from an inverse operator.
 
     Currently works only for fixed/loose orientation constraints
@@ -549,9 +482,9 @@ def _get_matrix_from_inverse_operator(
     info_inv = _prepare_info(inverse_operator)
 
     # only use channels that are good for inverse operator and forward sol
-    ch_names_inv = info_inv["ch_names"]
+    ch_names_inv = info_inv['ch_names']
     n_chs_inv = len(ch_names_inv)
-    bads_inv = inverse_operator["info"]["bads"]
+    bads_inv = inverse_operator['info']['bads']
 
     # indices of bad channels
     ch_idx_bads = [ch_names_inv.index(ch) for ch in bads_inv]
@@ -561,7 +494,7 @@ def _get_matrix_from_inverse_operator(
     id_mat = np.eye(n_chs_inv)
 
     # convert identity matrix to evoked data type (pretending it's an epoch)
-    ev_id = EvokedArray(id_mat, info=info_inv, tmin=0.0)
+    ev_id = EvokedArray(id_mat, info=info_inv, tmin=0.)
 
     # apply inverse operator to identity matrix in order to get inverse matrix
     # free orientation constraint not possible because apply_inverse would
@@ -574,12 +507,11 @@ def _get_matrix_from_inverse_operator(
     if is_fixed_inv:
         pick_ori = None
     else:
-        pick_ori = "vector"
+        pick_ori = 'vector'
 
     # columns for bad channels will be zero
-    invmat_op = apply_inverse(
-        ev_id, inverse_operator, lambda2=lambda2, method=method, pick_ori=pick_ori
-    )
+    invmat_op = apply_inverse(ev_id, inverse_operator, lambda2=lambda2,
+                              method=method, pick_ori=pick_ori)
 
     # turn source estimate into numpy array
     invmat = invmat_op.data
@@ -606,5 +538,5 @@ def _get_matrix_from_inverse_operator(
 
 def _check_fixed_ori(inst):
     """Check if inverse or forward was computed for fixed orientations."""
-    is_fixed = inst["source_ori"] != FIFF.FIFFV_MNE_FREE_ORI
+    is_fixed = inst['source_ori'] != FIFF.FIFFV_MNE_FREE_ORI
     return is_fixed
diff --git a/mne/minimum_norm/tests/test_resolution_matrix.py b/mne/minimum_norm/tests/test_resolution_matrix.py
index 3e09bab83df..95bb7fccc8f 100755
--- a/mne/minimum_norm/tests/test_resolution_matrix.py
+++ b/mne/minimum_norm/tests/test_resolution_matrix.py
@@ -1,172 +1,96 @@
+# -*- coding: utf-8 -*-
 # Author: Olaf Hauk <olaf.hauk@mrc-cbu.cam.ac.uk>
 #         Alexandre Gramfort <alexandre.gramfort@inria.fr>
 #         Eric Larson <larson.eric.d@gmail.com>
 #
 # License: BSD-3-Clause
 
-from contextlib import nullcontext
-
+import os.path as op
 import numpy as np
-import pytest
-from numpy.testing import assert_allclose, assert_array_almost_equal, assert_array_equal
+from numpy.testing import (assert_equal, assert_array_almost_equal,
+                           assert_array_equal, assert_allclose)
 
 import mne
 from mne.datasets import testing
-from mne.minimum_norm.resolution_matrix import (
-    _vertices_for_get_psf_ctf,
-    get_cross_talk,
-    get_point_spread,
-    make_inverse_resolution_matrix,
-)
+from mne.minimum_norm import apply_inverse_cov, prepare_inverse_operator
+from mne.minimum_norm.resolution_matrix_tmp import (make_inverse_resolution_matrix,
+                                                    get_cross_talk,
+                                                    get_point_spread,
+                                                    _vertices_for_get_psf_ctf)
 
 data_path = testing.data_path(download=False)
-subjects_dir = data_path / "subjects"
-fname_evoked = data_path / "MEG" / "sample" / "sample_audvis_trunc-ave.fif"
-# Intentional mismatch here!
-fname_fwd = data_path / "MEG" / "sample" / "sample_audvis_trunc-meg-eeg-oct-4-fwd.fif"
-fname_cov = data_path / "MEG" / "sample" / "sample_audvis_trunc-cov.fif"
-fname_label = data_path / "subjects" / "sample" / "label" / "lh.V1.label"
+subjects_dir = op.join(data_path, 'subjects')
+fname_inv = op.join(data_path, 'MEG', 'sample',
+                    'sample_audvis_trunc-meg-eeg-oct-6-meg-inv.fif')
+fname_evoked = op.join(data_path, 'MEG', 'sample',
+                       'sample_audvis_trunc-ave.fif')
+fname_raw = op.join(data_path, 'MEG', 'sample', 'sample_audvis_trunc_raw.fif')
+fname_t1 = op.join(data_path, 'subjects', 'sample', 'mri', 'T1.mgz')
+fname_src = op.join(data_path, 'MEG', 'sample',
+                    'sample_audvis_trunc-meg-eeg-oct-6-fwd.fif')
+fname_src_fs = op.join(data_path, 'subjects', 'fsaverage', 'bem',
+                       'fsaverage-ico-5-src.fif')
+fname_src_3 = op.join(data_path, 'subjects', 'sample', 'bem',
+                      'sample-oct-4-src.fif')
+fname_stc = op.join(data_path, 'MEG', 'sample', 'sample_audvis_trunc-meg')
+fname_vol = op.join(data_path, 'MEG', 'sample',
+                    'sample_audvis_trunc-grad-vol-7-fwd-sensmap-vol.w')
+fname_vsrc = op.join(data_path, 'MEG', 'sample',
+                     'sample_audvis_trunc-meg-vol-7-fwd.fif')
+fname_inv_vol = op.join(data_path, 'MEG', 'sample',
+                        'sample_audvis_trunc-meg-vol-7-meg-inv.fif')
+rng = np.random.RandomState(0)
+
+fname_fwd = op.join(data_path, 'MEG', 'sample',
+                    'sample_audvis_trunc-meg-eeg-oct-4-fwd.fif')
+fname_cov = op.join(data_path, 'MEG', 'sample', 'sample_audvis_trunc-cov.fif')
+
+fname_label = op.join(data_path, 'subjects', 'sample', 'label', 'lh.V1.label')
 
 
 @testing.requires_testing_data
-@pytest.mark.parametrize("src_type", ("surface", "volume"))
-def test_resolution_matrix_free(src_type, fwd_volume_small):
-    """Test make_inverse_resolution_matrix on surfaces."""
+def test_resolution_matrix():
+    """Test make_inverse_resolution_matrix() function."""
     # read forward solution
-    if src_type == "surface":
-        forward = mne.read_forward_solution(fname_fwd)
-        verbose = False
-        # Some arbitrary vertex numbers
-        idx = [1, 100, 400]
-    else:
-        assert src_type == "volume"
-        forward = fwd_volume_small
-        verbose = "error"  # ignore missing chs in vol
-        idx = [1, 3, 8]
-    assert forward["src"].kind == src_type
+    forward = mne.read_forward_solution(fname_fwd)
+    # forward operator with fixed source orientations
+    forward_fxd = mne.convert_forward_solution(forward, surf_ori=True,
+                                               force_fixed=True)
+
+    # noise covariance matrix
     noise_cov = mne.read_cov(fname_cov)
+    # evoked data for info
     evoked = mne.read_evokeds(fname_evoked, 0)
 
     # make inverse operator from forward solution
     # free source orientation
     inverse_operator = mne.minimum_norm.make_inverse_operator(
-        info=evoked.info,
-        forward=forward,
-        noise_cov=noise_cov,
-        loose=1.0,
-        depth=None,
-        verbose=verbose,
-    )
-    assert_allclose(inverse_operator["source_nn"], forward["source_nn"])
+        info=evoked.info, forward=forward, noise_cov=noise_cov, loose=1.,
+        depth=None)
+    # fixed source orientation
+    inverse_operator_fxd = mne.minimum_norm.make_inverse_operator(
+        info=evoked.info, forward=forward, noise_cov=noise_cov, loose=0.,
+        depth=None, fixed=True)
 
     # regularisation parameter based on SNR
     snr = 3.0
-    lambda2 = 1.0 / snr**2
+    lambda2 = 1.0 / snr ** 2
     # resolution matrices for free source orientation
     # compute resolution matrix for MNE with free source orientations
-    rm_mne_free = make_inverse_resolution_matrix(
-        forward, inverse_operator, method="MNE", lambda2=lambda2, verbose=verbose
-    )
-    assert_array_almost_equal(rm_mne_free, rm_mne_free.T)
-    # check various summary and normalisation options
-    for mode in [None, "sum", "mean", "maxval", "maxnorm", "pca"]:
-        n_comps = [1, 3]
-        if mode in [None, "sum", "mean"]:
-            n_comps = [1]
-        for n_comp in n_comps:
-            for norm in [None, "max", "norm", True]:
-                # with free orientations and vector source estimates
-                if src_type == "surface":
-                    ctx = pytest.raises(TypeError, match="vector surface")
-                else:
-                    ctx = nullcontext()
-                with ctx:
-                    stc_psf_free = get_point_spread(
-                        rm_mne_free,
-                        forward["src"],
-                        idx,
-                        mode=mode,
-                        n_comp=n_comp,
-                        norm=norm,
-                        return_pca_vars=False,
-                        vector=True,
-                    )
-                stc_psf_free = get_point_spread(
-                    rm_mne_free,
-                    inverse_operator,
-                    idx,
-                    mode=mode,
-                    n_comp=n_comp,
-                    norm=norm,
-                    return_pca_vars=False,
-                    vector=True,
-                )
-                stc_ctf_free = get_cross_talk(
-                    rm_mne_free,
-                    forward,
-                    idx,
-                    mode=mode,
-                    n_comp=n_comp,
-                    norm=norm,
-                    return_pca_vars=False,
-                    vector=True,
-                )
-                err_msg = f"mode={mode}, n_comp={n_comp}, norm={norm}"
-                # There is an ambiguity in the sign flip from the PCA here.
-                # Ideally we would use the normals to fix it, but it's not
-                # trivial.
-                if mode == "pca" and n_comp == 3:
-                    stc_psf_free = abs(stc_psf_free)
-                    stc_ctf_free = abs(stc_psf_free)
-                assert_array_almost_equal(
-                    stc_psf_free.data, stc_ctf_free.data, err_msg=err_msg
-                )
-    # For MNE, PSF and CTF for same vertices should be the same
-    label = mne.read_label(fname_label)
-    label = [label]
-    stc_psf_label_free = get_point_spread(
-        rm_mne_free, inverse_operator, label, norm="norm", vector=True
-    )
-    stc_ctf_label_free = get_cross_talk(
-        rm_mne_free, inverse_operator, label, norm="norm", vector=True
-    )
-    assert_array_almost_equal(stc_psf_label_free.data, stc_ctf_label_free.data)
-
-
-@testing.requires_testing_data
-def test_resolution_matrix_fixed():
-    """Test resolution matrices with fixed orientations."""
+    rm_mne_free = make_inverse_resolution_matrix(forward, inverse_operator,
+                                                 method='MNE', lambda2=lambda2)
     # compute resolution matrix for MNE, fwd fixed and inv free
-    forward = mne.read_forward_solution(fname_fwd)
-    forward_fxd = mne.convert_forward_solution(forward, surf_ori=True, force_fixed=True)
-    noise_cov = mne.read_cov(fname_cov)
-    evoked = mne.read_evokeds(fname_evoked, 0)
-    inverse_operator = mne.minimum_norm.make_inverse_operator(
-        info=evoked.info, forward=forward, noise_cov=noise_cov, loose=1.0, depth=None
-    )
-    inverse_operator_fxd = mne.minimum_norm.make_inverse_operator(
-        info=evoked.info,
-        forward=forward,
-        noise_cov=noise_cov,
-        loose=0.0,
-        depth=None,
-        fixed=True,
-    )
-    snr = 3.0
-    lambda2 = 1.0 / snr**2
-    rm_mne_fxdfree = make_inverse_resolution_matrix(
-        forward_fxd, inverse_operator, method="MNE", lambda2=lambda2
-    )
-
+    rm_mne_fxdfree = make_inverse_resolution_matrix(forward_fxd,
+                                                    inverse_operator,
+                                                    method='MNE',
+                                                    lambda2=lambda2)
     # resolution matrices for fixed source orientation
     # compute resolution matrix for MNE
-    rm_mne = make_inverse_resolution_matrix(
-        forward_fxd, inverse_operator_fxd, method="MNE", lambda2=lambda2
-    )
+    rm_mne = make_inverse_resolution_matrix(forward_fxd, inverse_operator_fxd,
+                                            method='MNE', lambda2=lambda2)
     # compute resolution matrix for sLORETA
-    rm_lor = make_inverse_resolution_matrix(
-        forward_fxd, inverse_operator_fxd, method="sLORETA", lambda2=lambda2
-    )
+    rm_lor = make_inverse_resolution_matrix(forward_fxd, inverse_operator_fxd,
+                                            method='sLORETA', lambda2=lambda2)
     # rectify resolution matrix for sLORETA before determining maxima
     rm_lor_abs = np.abs(rm_lor)
 
@@ -180,65 +104,53 @@ def test_resolution_matrix_fixed():
     assert_array_equal(maxidxs, goodidxs)
     # MNE resolution matrices symmetric?
     assert_array_almost_equal(rm_mne, rm_mne.T)
+    assert_array_almost_equal(rm_mne_free, rm_mne_free.T)
 
     # Some arbitrary vertex numbers
     idx = [1, 100, 400]
     # check various summary and normalisation options
-    for mode in [None, "sum", "mean", "maxval", "maxnorm", "pca"]:
+    for mode in [None, 'sum', 'mean', 'maxval', 'maxnorm', 'pca']:
         n_comps = [1, 3]
-        if mode in [None, "sum", "mean"]:
+        if mode in [None, 'sum', 'mean']:
             n_comps = [1]
         for n_comp in n_comps:
-            for norm in [None, "max", "norm", True]:
+            for norm in [None, 'max', 'norm', True]:
                 stc_psf = get_point_spread(
-                    rm_mne,
-                    forward_fxd["src"],
-                    idx,
-                    mode=mode,
-                    n_comp=n_comp,
-                    norm=norm,
-                    return_pca_vars=False,
-                )
+                    rm_mne, forward_fxd['src'], idx, mode=mode, n_comp=n_comp,
+                    norm=norm, return_pca_vars=False)
                 stc_ctf = get_cross_talk(
-                    rm_mne,
-                    forward_fxd["src"],
-                    idx,
-                    mode=mode,
-                    n_comp=n_comp,
-                    norm=norm,
-                    return_pca_vars=False,
-                )
+                    rm_mne, forward_fxd['src'], idx, mode=mode, n_comp=n_comp,
+                    norm=norm, return_pca_vars=False)
                 # for MNE, PSF/CTFs for same vertices should be the same
                 assert_array_almost_equal(stc_psf.data, stc_ctf.data)
+                # with free orientations and vector source estimates
+                stc_psf_free = get_point_spread(
+                    rm_mne_free, forward_fxd['src'], idx, mode=mode,
+                    n_comp=n_comp, norm=norm, return_pca_vars=False,
+                    vector=True)
+                stc_ctf_free = get_cross_talk(
+                    rm_mne_free, forward_fxd['src'], idx, mode=mode,
+                    n_comp=n_comp, norm=norm, return_pca_vars=False,
+                    vector=True)
+                assert_array_almost_equal(stc_psf_free.data, stc_ctf_free.data)
 
     # check SVD variances
     n_comp = 3
     stc_psf, s_vars_psf = get_point_spread(
-        rm_mne,
-        forward_fxd["src"],
-        idx,
-        mode=mode,
-        n_comp=n_comp,
-        norm="norm",
-        return_pca_vars=True,
-    )
+        rm_mne, forward_fxd['src'], idx, mode=mode, n_comp=n_comp,
+        norm='norm', return_pca_vars=True)
     stc_ctf, s_vars_ctf = get_cross_talk(
-        rm_mne,
-        forward_fxd["src"],
-        idx,
-        mode=mode,
-        n_comp=n_comp,
-        norm="norm",
-        return_pca_vars=True,
-    )
+        rm_mne, forward_fxd['src'], idx, mode=mode, n_comp=n_comp,
+        norm='norm', return_pca_vars=True)
     assert_array_almost_equal(s_vars_psf, s_vars_ctf)
     # variances for SVD components should be ordered
     assert s_vars_psf[0] > s_vars_psf[1] > s_vars_psf[2]
     # all variances should sum up to 100
-    assert_allclose(s_vars_psf.sum(), 100.0)
+    assert_allclose(s_vars_psf.sum(), 100.)
 
     # Test application of free inv to fixed fwd
-    assert rm_mne_fxdfree.shape == (3 * rm_mne.shape[0], rm_mne.shape[0])
+    assert_equal(rm_mne_fxdfree.shape, (3 * rm_mne.shape[0],
+                 rm_mne.shape[0]))
 
     # Test PSF/CTF for labels
     label = mne.read_label(fname_label)
@@ -246,125 +158,132 @@ def test_resolution_matrix_fixed():
     label = [label]
     label2 = 2 * label
     # get relevant vertices in source space
-    verts = _vertices_for_get_psf_ctf(label, forward_fxd["src"])[0]
+    verts = _vertices_for_get_psf_ctf(label, forward_fxd['src'])[0]
 
-    stc_psf_label = get_point_spread(rm_mne, forward_fxd["src"], label, norm="max")
+    stc_psf_label = get_point_spread(rm_mne, forward_fxd['src'], label,
+                                     norm='max')
     # for list of indices
-    stc_psf_idx = get_point_spread(rm_mne, forward_fxd["src"], verts, norm="max")
-    stc_ctf_label = get_cross_talk(rm_mne, forward_fxd["src"], label, norm="max")
+    stc_psf_idx = get_point_spread(rm_mne, forward_fxd['src'], verts,
+                                   norm='max')
+    stc_ctf_label = get_cross_talk(rm_mne, forward_fxd['src'], label,
+                                   norm='max')
     # For MNE, PSF and CTF for same vertices should be the same
     assert_array_almost_equal(stc_psf_label.data, stc_ctf_label.data)
+    # same with free orientations and vector source estimates
+    stc_psf_label_free = get_point_spread(
+        rm_mne_free, forward_fxd['src'], label, norm='norm', vector=True)
+    stc_ctf_label_free = get_cross_talk(
+        rm_mne_free, forward_fxd['src'], label, norm='norm', vector=True)
+    assert_array_almost_equal(stc_psf_label_free.data, stc_ctf_label_free.data)
 
     # test multiple labels
-    stc_psf_label2 = get_point_spread(rm_mne, forward_fxd["src"], label2, norm="max")
+    stc_psf_label2 = get_point_spread(rm_mne, forward_fxd['src'], label2,
+                                      norm='max')
     m, n = stc_psf_label.data.shape
-    assert_array_equal(stc_psf_label.data, stc_psf_label2[0].data)
-    assert_array_equal(stc_psf_label.data, stc_psf_label2[1].data)
-    assert_array_equal(stc_psf_label.data, stc_psf_idx.data)
+    assert_array_equal(
+        stc_psf_label.data, stc_psf_label2[0].data)
+    assert_array_equal(
+        stc_psf_label.data, stc_psf_label2[1].data)
+    assert_array_equal(
+        stc_psf_label.data, stc_psf_idx.data)
 
     # test noisy resolution matrix
     # resmats without noise
     rm_mne1 = make_inverse_resolution_matrix(
-        forward_fxd,
-        inverse_operator_fxd,
-        method="MNE",
-        lambda2=lambda2,
-        noise_cov=None,
-        snr=None,
-    )
+        forward_fxd, inverse_operator_fxd, method='MNE', lambda2=lambda2,
+        noise_cov=None, snr=None)
     rm_mne2 = make_inverse_resolution_matrix(
-        forward_fxd,
-        inverse_operator_fxd,
-        method="dSPM",
-        lambda2=lambda2,
-        noise_cov=None,
-        snr=None,
-    )
+        forward_fxd, inverse_operator_fxd, method='dSPM', lambda2=lambda2,
+        noise_cov=None, snr=None)
     rm_mne3 = make_inverse_resolution_matrix(
-        forward_fxd,
-        inverse_operator_fxd,
-        method="sLORETA",
-        lambda2=lambda2,
-        noise_cov=None,
-        snr=None,
-    )
+        forward_fxd, inverse_operator_fxd, method='sLORETA', lambda2=lambda2,
+        noise_cov=None, snr=None)
     rm_mne4 = make_inverse_resolution_matrix(
-        forward_fxd,
-        inverse_operator_fxd,
-        method="eLORETA",
-        lambda2=lambda2,
-        noise_cov=None,
-        snr=None,
-    )
+        forward_fxd, inverse_operator_fxd, method='eLORETA', lambda2=lambda2,
+        noise_cov=None, snr=None)
 
     # same resmats with noise but very high SNR
-    rm_mne1a = make_inverse_resolution_matrix(
-        forward_fxd,
-        inverse_operator_fxd,
-        method="MNE",
-        lambda2=lambda2,
-        noise_cov=noise_cov,
-        snr=1e30,
-    )
-    rm_mne2a = make_inverse_resolution_matrix(
-        forward_fxd,
-        inverse_operator_fxd,
-        method="dSPM",
-        lambda2=lambda2,
-        noise_cov=noise_cov,
-        snr=1e30,
-    )
-    rm_mne3a = make_inverse_resolution_matrix(
-        forward_fxd,
-        inverse_operator_fxd,
-        method="sLORETA",
-        lambda2=lambda2,
-        noise_cov=noise_cov,
-        snr=1e30,
-    )
-    rm_mne4a = make_inverse_resolution_matrix(
-        forward_fxd,
-        inverse_operator_fxd,
-        method="eLORETA",
-        lambda2=lambda2,
-        noise_cov=noise_cov,
-        snr=1e30,
-    )
+    rm_mne1h = make_inverse_resolution_matrix(
+        forward_fxd, inverse_operator_fxd, method='MNE', lambda2=lambda2,
+        noise_cov=noise_cov, snr=1e20)
+    rm_mne2h = make_inverse_resolution_matrix(
+        forward_fxd, inverse_operator_fxd, method='dSPM', lambda2=lambda2,
+        noise_cov=noise_cov, snr=1e20)
+    rm_mne3h = make_inverse_resolution_matrix(
+        forward_fxd, inverse_operator_fxd, method='sLORETA', lambda2=lambda2,
+        noise_cov=noise_cov, snr=1e20)
+    rm_mne4h = make_inverse_resolution_matrix(
+        forward_fxd, inverse_operator_fxd, method='eLORETA', lambda2=lambda2,
+        noise_cov=noise_cov, snr=1e20)
+
+    # same resmats with noise and realistic/medium SNR
+    rm_mne1m = make_inverse_resolution_matrix(
+        forward_fxd, inverse_operator_fxd, method='MNE', lambda2=lambda2,
+        noise_cov=noise_cov, snr=3)
+    rm_mne2m = make_inverse_resolution_matrix(
+        forward_fxd, inverse_operator_fxd, method='dSPM', lambda2=lambda2,
+        noise_cov=noise_cov, snr=3)
+    rm_mne3m = make_inverse_resolution_matrix(
+        forward_fxd, inverse_operator_fxd, method='sLORETA', lambda2=lambda2,
+        noise_cov=noise_cov, snr=3)
+    rm_mne4m = make_inverse_resolution_matrix(
+        forward_fxd, inverse_operator_fxd, method='eLORETA', lambda2=lambda2,
+        noise_cov=noise_cov, snr=3)
+
+    # same resmats with noise and very low SNR
+    rm_mne1l = make_inverse_resolution_matrix(
+        forward_fxd, inverse_operator_fxd, method='MNE', lambda2=lambda2,
+        noise_cov=noise_cov, snr=1e-20)
 
     # resmats should be very close without noise and at high SNR
-    assert_array_almost_equal(np.abs(rm_mne1), rm_mne1a)
-    assert_array_almost_equal(np.abs(rm_mne2), rm_mne2a)
-    assert_array_almost_equal(np.abs(rm_mne3), rm_mne3a)
-    assert_array_almost_equal(np.abs(rm_mne4), rm_mne4a)
+    assert_array_almost_equal(np.abs(rm_mne1), rm_mne1h)
+    assert_array_almost_equal(np.abs(rm_mne2), rm_mne2h)
+    assert_array_almost_equal(np.abs(rm_mne3), rm_mne3h)
+    assert_array_almost_equal(np.abs(rm_mne4), rm_mne4h)
 
     # rows of resmats (CTFs) should be similar between methods
-    row_corrs = [
-        np.corrcoef(rm_mne1[i, :], rm_mne2[i, :])[0, 1]
-        for i in np.arange(rm_mne1.shape[0])
-    ]
+    row_corrs = [np.corrcoef(rm_mne1[i, :], rm_mne2[i, :])[0, 1]
+                 for i in np.arange(rm_mne1.shape[0])]
     assert_allclose(row_corrs, 1)
-    row_corrs = [
-        np.corrcoef(rm_mne1[i, :], rm_mne3[i, :])[0, 1]
-        for i in np.arange(rm_mne1.shape[0])
-    ]
+    row_corrs = [np.corrcoef(rm_mne1[i, :], rm_mne3[i, :])[0, 1]
+                 for i in np.arange(rm_mne1.shape[0])]
     assert_allclose(row_corrs, 1)
-    row_corrs = [
-        np.corrcoef(rm_mne1[i, :], rm_mne4[i, :])[0, 1]
-        for i in np.arange(rm_mne1.shape[0])
-    ]
+    row_corrs = [np.corrcoef(rm_mne1[i, :], rm_mne4[i, :])[0, 1]
+                 for i in np.arange(rm_mne1.shape[0])]
     assert_allclose(row_corrs, 1, 0.2)  # eLORETA is different
-    row_corrs = [
-        np.corrcoef(rm_mne1a[i, :], rm_mne2a[i, :])[0, 1]
-        for i in np.arange(rm_mne1.shape[0])
-    ]
+
+    # noisy:
+    row_corrs = [np.corrcoef(rm_mne1h[i, :], rm_mne2h[i, :])[0, 1]
+                 for i in np.arange(rm_mne1.shape[0])]
     assert_allclose(row_corrs, 1)
-    row_corrs = [
-        np.corrcoef(rm_mne1a[i, :], rm_mne3a[i, :])[0, 1]
-        for i in np.arange(rm_mne1.shape[0])
-    ]
+    row_corrs = [np.corrcoef(rm_mne1h[i, :], rm_mne3h[i, :])[0, 1]
+                 for i in np.arange(rm_mne1.shape[0])]
     assert_allclose(row_corrs, 1)
-    row_corrs = [
-        np.corrcoef(rm_mne1a[i, :], rm_mne4a[i, :])[0, 1]
-        for i in np.arange(rm_mne1.shape[0])
-    ]
+    row_corrs = [np.corrcoef(rm_mne1h[i, :], rm_mne4h[i, :])[0, 1]
+                 for i in np.arange(rm_mne1.shape[0])]
     assert_allclose(row_corrs, 1, 0.2)  # eLORETA is different
+
+    row_corrs = [np.corrcoef(rm_mne1m[i, :], rm_mne2m[i, :])[0, 1]
+                 for i in np.arange(rm_mne1.shape[0])]
+    assert_allclose(row_corrs, 1)
+    row_corrs = [np.corrcoef(rm_mne1m[i, :], rm_mne3m[i, :])[0, 1]
+                 for i in np.arange(rm_mne1.shape[0])]
+    assert_allclose(row_corrs, 1)
+    row_corrs = [np.corrcoef(rm_mne1m[i, :], rm_mne4m[i, :])[0, 1]
+                 for i in np.arange(rm_mne1.shape[0])]
+    assert_allclose(row_corrs, 1, 0.2)  # eLORETA is different
+
+    # for very low SNR, resmat should be equal to noise power (except for scaling)
+    inv = prepare_inverse_operator(inverse_operator_fxd, 1, lambda2, 'MNE', copy='non-src')
+    stc = apply_inverse_cov(
+        noise_cov, evoked.info, inv, nave=1, lambda2=lambda2,
+        method='MNE', pick_ori=None, prepared=True, label=None,
+        method_params=None, use_cps=True, verbose=None)
+
+    x = np.sqrt(stc.data)
+    y = rm_mne1l
+    x = x / x.max()
+    x = np.tile(x, x.shape[0])
+    y = y / y.max()
+
+    assert_array_almost_equal(x, y)

From be49e808df6bdb0c6aca5e8f4d7332e013e13fe1 Mon Sep 17 00:00:00 2001
From: "pre-commit-ci[bot]"
 <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Date: Wed, 15 Nov 2023 16:44:38 +0000
Subject: [PATCH 12/27] [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci
---
 mne/minimum_norm/resolution_matrix.py         | 256 +++++++----
 .../tests/test_resolution_matrix.py           | 407 ++++++++++++------
 2 files changed, 437 insertions(+), 226 deletions(-)

diff --git a/mne/minimum_norm/resolution_matrix.py b/mne/minimum_norm/resolution_matrix.py
index 755d637b58c..754684f2f8c 100755
--- a/mne/minimum_norm/resolution_matrix.py
+++ b/mne/minimum_norm/resolution_matrix.py
@@ -7,24 +7,32 @@
 
 import numpy as np
 
-from .. import pick_types, pick_channels_forward, EvokedArray, SourceEstimate
-from .. import (pick_channels_forward, EvokedArray, SourceEstimate,
-                VectorSourceEstimate)
-
-from ..io.constants import FIFF
-from ..utils import logger, verbose
+from .. import (
+    EvokedArray,
+    SourceEstimate,
+    VectorSourceEstimate,
+    pick_channels_forward,
+    pick_types,
+)
 from ..forward.forward import convert_forward_solution
-from ..minimum_norm import (apply_inverse, apply_inverse_cov,
-                            prepare_inverse_operator)
+from ..io.constants import FIFF
+from ..label import Label
+from ..minimum_norm import apply_inverse, apply_inverse_cov, prepare_inverse_operator
 from ..minimum_norm.spatial_resolution import _rectify_resolution_matrix
 from ..source_estimate import _prepare_label_extraction
-from ..label import Label
+from ..utils import logger, verbose
 
 
 @verbose
-def make_inverse_resolution_matrix(forward, inverse_operator, method='dSPM',
-                                   lambda2=1. / 9., noise_cov=None, snr=None,
-                                   verbose=None):
+def make_inverse_resolution_matrix(
+    forward,
+    inverse_operator,
+    method="dSPM",
+    lambda2=1.0 / 9.0,
+    noise_cov=None,
+    snr=None,
+    verbose=None,
+):
     """Compute resolution matrix for linear inverse operator.
 
     Parameters
@@ -62,14 +70,13 @@ def make_inverse_resolution_matrix(forward, inverse_operator, method='dSPM',
         interpreted. It must not be used to compute CTFs or resolution metrics
         for CTFs!
     """
-    if (noise_cov is None and snr is not None):
-        msg = 'snr should be None if noise_cov is None.'
+    if noise_cov is None and snr is not None:
+        msg = "snr should be None if noise_cov is None."
         raise ValueError(msg)
 
-    if (noise_cov is not None and snr is None):
-        snr = np.sqrt(1. / lambda2)
-        logger.info('Inferring snr from lamdba2 as: %f.' % snr)
-
+    if noise_cov is not None and snr is None:
+        snr = np.sqrt(1.0 / lambda2)
+        logger.info("Inferring snr from lamdba2 as: %f." % snr)
 
     # make sure forward and inverse operator match
     inv = inverse_operator
@@ -77,27 +84,36 @@ def make_inverse_resolution_matrix(forward, inverse_operator, method='dSPM',
 
     # don't include bad channels
     # only use good channels from inverse operator
-    bads_inv = inv['info']['bads']
+    bads_inv = inv["info"]["bads"]
     # good channels
-    ch_names = [c for c in inv['info']['ch_names'] if (c not in bads_inv)]
+    ch_names = [c for c in inv["info"]["ch_names"] if (c not in bads_inv)]
     fwd = pick_channels_forward(fwd, ch_names, ordered=True)
 
     # get leadfield matrix from forward solution
-    leadfield = fwd['sol']['data']
-    invmat = _get_matrix_from_inverse_operator(inv, fwd,
-                                               method=method, lambda2=lambda2)
+    leadfield = fwd["sol"]["data"]
+    invmat = _get_matrix_from_inverse_operator(inv, fwd, method=method, lambda2=lambda2)
     resmat = invmat.dot(leadfield)
-    logger.info('Dimensions of resolution matrix: %d by %d.' % resmat.shape)
+    logger.info("Dimensions of resolution matrix: %d by %d." % resmat.shape)
 
     # add source noise power to columns of resolution matrix
     if noise_cov is not None:
-        inv = prepare_inverse_operator(inv, 1, lambda2, method, copy='non-src')
+        inv = prepare_inverse_operator(inv, 1, lambda2, method, copy="non-src")
         info = _prepare_info(inv)
         # compute source noise power
         stc = apply_inverse_cov(
-            noise_cov, info, inv, nave=1, lambda2=lambda2,
-            method=method, pick_ori=None, prepared=True, label=None,
-            method_params=None, use_cps=True, verbose=None)
+            noise_cov,
+            info,
+            inv,
+            nave=1,
+            lambda2=lambda2,
+            method=method,
+            pick_ori=None,
+            prepared=True,
+            label=None,
+            method_params=None,
+            use_cps=True,
+            verbose=None,
+        )
 
         # output will be square-root of power, so take intensities before
         # adding noise
@@ -111,21 +127,19 @@ def make_inverse_resolution_matrix(forward, inverse_operator, method='dSPM',
         # scaling factor
 
         idx_ch_types = {}
-        idx_ch_types['eeg'] = pick_types(info, eeg=True, meg=False)
-        idx_ch_types['mag'] = pick_types(info, eeg=False, meg='mag')
-        idx_ch_types['gra'] = pick_types(info, eeg=False, meg='grad')
+        idx_ch_types["eeg"] = pick_types(info, eeg=True, meg=False)
+        idx_ch_types["mag"] = pick_types(info, eeg=False, meg="mag")
+        idx_ch_types["gra"] = pick_types(info, eeg=False, meg="grad")
 
         # scaling similar to Samuelsson et al., Neuroimage 2021 (p. 4)
         alphas = {}
         for cht in idx_ch_types:
             # signal intensity of all sources across channels of this type
-            alphas[cht] = (
-                np.sqrt((leadfield[idx_ch_types[cht], :]**2).mean()))
+            alphas[cht] = np.sqrt((leadfield[idx_ch_types[cht], :] ** 2).mean())
             # divided by noise intensity across channels of this type
-            alphas[cht] /= (
-                np.sqrt(np.diag(noise_cov.data)[idx_ch_types[cht]].mean()))
+            alphas[cht] /= np.sqrt(np.diag(noise_cov.data)[idx_ch_types[cht]].mean())
 
-        alpha = (3 * snr) / (alphas['eeg'] + alphas['mag'] + alphas['gra'])
+        alpha = (3 * snr) / (alphas["eeg"] + alphas["mag"] + alphas["gra"])
 
         # # previous attempt to compute scaling via whitened gain matrix
         # whitener = inv['whitener']
@@ -142,8 +156,19 @@ def make_inverse_resolution_matrix(forward, inverse_operator, method='dSPM',
 
 
 @verbose
-def _get_psf_ctf(resmat, src, idx, *, func, mode, n_comp, norm,
-                 return_pca_vars, vector=False, verbose=None):
+def _get_psf_ctf(
+    resmat,
+    src,
+    idx,
+    *,
+    func,
+    mode,
+    n_comp,
+    norm,
+    return_pca_vars,
+    vector=False,
+    verbose=None,
+):
     """Get point-spread (PSFs) or cross-talk (CTFs) functions.
 
     Parameters
@@ -172,26 +197,28 @@ def _get_psf_ctf(resmat, src, idx, *, func, mode, n_comp, norm,
 
     # backward compatibility
     if norm is True:
-        norm = 'max'
+        norm = "max"
 
     # get relevant vertices in source space
     verts_all = _vertices_for_get_psf_ctf(idx, src)
-    vertno_lh = src[0]['vertno']
-    vertno_rh = src[1]['vertno']
+    vertno_lh = src[0]["vertno"]
+    vertno_rh = src[1]["vertno"]
     vertno = [vertno_lh, vertno_rh]
 
     n_verts = len(vertno[0]) + len(vertno[1])
 
     n_r, n_c = resmat.shape
-    if (((n_verts != n_r) and (n_r / 3 != n_verts))
-            or ((n_verts != n_c) and (n_c / 3 != n_verts))):
-        msg = ('Number of vertices (%d) and corresponding dimension of'
-               'resolution matrix ((%d, %d) do not match' %
-               (n_verts, n_r, n_c))
+    if ((n_verts != n_r) and (n_r / 3 != n_verts)) or (
+        (n_verts != n_c) and (n_c / 3 != n_verts)
+    ):
+        msg = (
+            "Number of vertices (%d) and corresponding dimension of"
+            "resolution matrix ((%d, %d) do not match" % (n_verts, n_r, n_c)
+        )
         raise ValueError(msg)
 
     # the following will operate on columns of funcs
-    if func == 'ctf':
+    if func == "ctf":
         resmat = resmat.T
         n_r, n_c = n_c, n_r
 
@@ -206,8 +233,7 @@ def _get_psf_ctf(resmat, src, idx, *, func, mode, n_comp, norm,
         for [i, verts] in enumerate(verts_all):
             verts_vec = np.empty(3 * len(verts), dtype=int)
             for [j, v] in enumerate(verts):
-                verts_vec[3 * j: 3 * j + 3] = \
-                    3 * verts[j] + np.array([0, 1, 2])
+                verts_vec[3 * j : 3 * j + 3] = 3 * verts[j] + np.array([0, 1, 2])
             verts_all[i] = verts_vec  # use these as indices
 
     for verts in verts_all:
@@ -223,23 +249,22 @@ def _get_psf_ctf(resmat, src, idx, *, func, mode, n_comp, norm,
         # summarise PSFs/CTFs across vertices if requested
         pca_var = None  # variances computed only if return_pca_vars=True
         if mode is not None:
-            funcs, pca_var = _summarise_psf_ctf(funcs, mode, n_comp,
-                                                return_pca_vars)
+            funcs, pca_var = _summarise_psf_ctf(funcs, mode, n_comp, return_pca_vars)
 
         if not vector:  # if one value per vertex requested
             if n_verts != n_r:  # if 3 orientations per vertex, combine
                 funcs_int = np.empty([int(n_r / 3), funcs.shape[1]])
                 for i in np.arange(0, n_verts):
-                    funcs_vert = funcs[3 * i:3 * i + 3, :]
-                    funcs_int[i, :] = np.sqrt((funcs_vert ** 2).sum(axis=0))
-                stc = SourceEstimate(funcs_int, vertno, tmin=0., tstep=1.)
+                    funcs_vert = funcs[3 * i : 3 * i + 3, :]
+                    funcs_int[i, :] = np.sqrt((funcs_vert**2).sum(axis=0))
+                stc = SourceEstimate(funcs_int, vertno, tmin=0.0, tstep=1.0)
             else:  # use as is
-                stc = SourceEstimate(funcs, vertno, tmin=0., tstep=1.)
+                stc = SourceEstimate(funcs, vertno, tmin=0.0, tstep=1.0)
         else:  # STC with orientations
             # convert to vector source estimate
             m, n = int(funcs.shape[0] / 3), int(funcs.shape[1])
             data = funcs.reshape(m, 3, n)
-            stc = VectorSourceEstimate(data, vertno, tmin=0., tstep=1.)
+            stc = VectorSourceEstimate(data, vertno, tmin=0.0, tstep=1.0)
 
         stcs.append(stc)
         pca_vars.append(pca_var)
@@ -257,11 +282,11 @@ def _get_psf_ctf(resmat, src, idx, *, func, mode, n_comp, norm,
 
 def _check_get_psf_ctf_params(mode, n_comp, return_pca_vars):
     """Check input parameters of _get_psf_ctf() for consistency."""
-    if mode in [None, 'sum', 'mean'] and n_comp > 1:
-        msg = 'n_comp must be 1 for mode=%s.' % mode
+    if mode in [None, "sum", "mean"] and n_comp > 1:
+        msg = "n_comp must be 1 for mode=%s." % mode
         raise ValueError(msg)
-    if mode != 'pca' and return_pca_vars:
-        msg = 'SVD variances can only be returned if mode=''pca''.'
+    if mode != "pca" and return_pca_vars:
+        msg = "SVD variances can only be returned if mode=" "pca" "."
         raise ValueError(msg)
 
 
@@ -272,8 +297,13 @@ def _vertices_for_get_psf_ctf(idx, src):
     if type(idx[0]) is Label:
         # specify without source time courses, gets indices per label
         verts_labs, _ = _prepare_label_extraction(
-            stc=None, labels=idx, src=src, mode='mean', allow_empty=False,
-            use_sparse=False)
+            stc=None,
+            labels=idx,
+            src=src,
+            mode="mean",
+            allow_empty=False,
+            use_sparse=False,
+        )
         # verts_labs can be list of lists
         # concatenate indices per label across hemispheres
         # one list item per label
@@ -300,10 +330,10 @@ def _vertices_for_get_psf_ctf(idx, src):
 def _normalise_psf_ctf(funcs, norm):
     """Normalise PSFs/CTFs in _get_psf_ctf()."""
     # normalise PSFs/CTFs if specified
-    if norm == 'max':
+    if norm == "max":
         maxval = max(-funcs.min(), funcs.max())
         funcs = funcs / maxval
-    elif norm == 'norm':  # normalise to maximum norm across columns
+    elif norm == "norm":  # normalise to maximum norm across columns
         norms = np.linalg.norm(funcs, axis=0)
         funcs = funcs / norms.max()
 
@@ -314,7 +344,7 @@ def _summarise_psf_ctf(funcs, mode, n_comp, return_pca_vars):
     """Summarise PSFs/CTFs across vertices."""
     s_var = None  # only computed for return_pca_vars=True
 
-    if mode == 'maxval':  # pick PSF/CTF with maximum absolute value
+    if mode == "maxval":  # pick PSF/CTF with maximum absolute value
         absvals = np.maximum(-np.min(funcs, axis=0), np.max(funcs, axis=0))
         if n_comp > 1:  # only keep requested number of sorted PSFs/CTFs
             sortidx = np.argsort(absvals)
@@ -323,7 +353,7 @@ def _summarise_psf_ctf(funcs, mode, n_comp, return_pca_vars):
             maxidx = [absvals.argmax()]
         funcs = funcs[:, maxidx]
 
-    elif mode == 'maxnorm':  # pick PSF/CTF with maximum norm
+    elif mode == "maxnorm":  # pick PSF/CTF with maximum norm
         norms = np.linalg.norm(funcs, axis=0)
         if n_comp > 1:  # only keep requested number of sorted PSFs/CTFs
             sortidx = np.argsort(norms)
@@ -332,13 +362,13 @@ def _summarise_psf_ctf(funcs, mode, n_comp, return_pca_vars):
             maxidx = [norms.argmax()]
         funcs = funcs[:, maxidx]
 
-    elif mode == 'sum':  # sum across PSFs/CTFs
+    elif mode == "sum":  # sum across PSFs/CTFs
         funcs = np.sum(funcs, axis=1, keepdims=True)
 
-    elif mode == 'mean':  # mean of PSFs/CTFs
+    elif mode == "mean":  # mean of PSFs/CTFs
         funcs = np.mean(funcs, axis=1, keepdims=True)
 
-    elif mode == 'pca':  # SVD across PSFs/CTFs
+    elif mode == "pca":  # SVD across PSFs/CTFs
         # compute SVD of PSFs/CTFs across vertices
         u, s, _ = np.linalg.svd(funcs, full_matrices=False, compute_uv=True)
         if n_comp > 1:
@@ -355,8 +385,18 @@ def _summarise_psf_ctf(funcs, mode, n_comp, return_pca_vars):
 
 
 @verbose
-def get_point_spread(resmat, src, idx, mode=None, *, n_comp=1, norm=False,
-                     return_pca_vars=False, vector=False, verbose=None):
+def get_point_spread(
+    resmat,
+    src,
+    idx,
+    mode=None,
+    *,
+    n_comp=1,
+    norm=False,
+    return_pca_vars=False,
+    vector=False,
+    verbose=None,
+):
     """Get point-spread (PSFs) functions for vertices.
 
     Parameters
@@ -378,14 +418,32 @@ def get_point_spread(resmat, src, idx, mode=None, *, n_comp=1, norm=False,
     %(stcs_pctf)s
     %(pca_vars_pctf)s
     """
-    return _get_psf_ctf(resmat, src, idx, func='psf', mode=mode, n_comp=n_comp,
-                        norm=norm, return_pca_vars=return_pca_vars,
-                        vector=vector)
+    return _get_psf_ctf(
+        resmat,
+        src,
+        idx,
+        func="psf",
+        mode=mode,
+        n_comp=n_comp,
+        norm=norm,
+        return_pca_vars=return_pca_vars,
+        vector=vector,
+    )
 
 
 @verbose
-def get_cross_talk(resmat, src, idx, mode=None, *, n_comp=1, norm=False,
-                   return_pca_vars=False, vector=False, verbose=None):
+def get_cross_talk(
+    resmat,
+    src,
+    idx,
+    mode=None,
+    *,
+    n_comp=1,
+    norm=False,
+    return_pca_vars=False,
+    vector=False,
+    verbose=None,
+):
     """Get cross-talk (CTFs) function for vertices.
 
     Parameters
@@ -407,9 +465,17 @@ def get_cross_talk(resmat, src, idx, mode=None, *, n_comp=1, norm=False,
     %(stcs_pctf)s
     %(pca_vars_pctf)s
     """
-    return _get_psf_ctf(resmat, src, idx, func='ctf', mode=mode, n_comp=n_comp,
-                        norm=norm, return_pca_vars=return_pca_vars,
-                        vector=vector)
+    return _get_psf_ctf(
+        resmat,
+        src,
+        idx,
+        func="ctf",
+        mode=mode,
+        n_comp=n_comp,
+        norm=norm,
+        return_pca_vars=return_pca_vars,
+        vector=vector,
+    )
 
 
 def _convert_forward_match_inv(fwd, inv):
@@ -427,12 +493,12 @@ def _convert_forward_match_inv(fwd, inv):
     # if inv loose: surf_ori must be True
     # if inv free: surf_ori must be False
     if not is_fixed_inv and not is_fixed_fwd:
-        is_loose_inv = not (inv['orient_prior']['data'] == 1.).all()
+        is_loose_inv = not (inv["orient_prior"]["data"] == 1.0).all()
 
         if is_loose_inv:
-            if not fwd['surf_ori']:
+            if not fwd["surf_ori"]:
                 fwd = convert_forward_solution(fwd, surf_ori=True)
-        elif fwd['surf_ori']:  # free orientation, change fwd
+        elif fwd["surf_ori"]:  # free orientation, change fwd
             fwd = convert_forward_solution(fwd, surf_ori=False)
 
     return fwd
@@ -443,15 +509,16 @@ def _prepare_info(inverse_operator):
     # in order to convert sub-leadfield matrix to evoked data type (pretending
     # it's an epoch, see in loop below), uses 'info' from inverse solution
     # because this has all the correct projector information
-    info = deepcopy(inverse_operator['info'])
+    info = deepcopy(inverse_operator["info"])
     with info._unlock():
-        info['sfreq'] = 1000.  # necessary
-        info['projs'] = inverse_operator['projs']
+        info["sfreq"] = 1000.0  # necessary
+        info["projs"] = inverse_operator["projs"]
     return info
 
 
-def _get_matrix_from_inverse_operator(inverse_operator, forward, method='dSPM',
-                                      lambda2=1. / 9.):
+def _get_matrix_from_inverse_operator(
+    inverse_operator, forward, method="dSPM", lambda2=1.0 / 9.0
+):
     """Get inverse matrix from an inverse operator.
 
     Currently works only for fixed/loose orientation constraints
@@ -482,9 +549,9 @@ def _get_matrix_from_inverse_operator(inverse_operator, forward, method='dSPM',
     info_inv = _prepare_info(inverse_operator)
 
     # only use channels that are good for inverse operator and forward sol
-    ch_names_inv = info_inv['ch_names']
+    ch_names_inv = info_inv["ch_names"]
     n_chs_inv = len(ch_names_inv)
-    bads_inv = inverse_operator['info']['bads']
+    bads_inv = inverse_operator["info"]["bads"]
 
     # indices of bad channels
     ch_idx_bads = [ch_names_inv.index(ch) for ch in bads_inv]
@@ -494,7 +561,7 @@ def _get_matrix_from_inverse_operator(inverse_operator, forward, method='dSPM',
     id_mat = np.eye(n_chs_inv)
 
     # convert identity matrix to evoked data type (pretending it's an epoch)
-    ev_id = EvokedArray(id_mat, info=info_inv, tmin=0.)
+    ev_id = EvokedArray(id_mat, info=info_inv, tmin=0.0)
 
     # apply inverse operator to identity matrix in order to get inverse matrix
     # free orientation constraint not possible because apply_inverse would
@@ -507,11 +574,12 @@ def _get_matrix_from_inverse_operator(inverse_operator, forward, method='dSPM',
     if is_fixed_inv:
         pick_ori = None
     else:
-        pick_ori = 'vector'
+        pick_ori = "vector"
 
     # columns for bad channels will be zero
-    invmat_op = apply_inverse(ev_id, inverse_operator, lambda2=lambda2,
-                              method=method, pick_ori=pick_ori)
+    invmat_op = apply_inverse(
+        ev_id, inverse_operator, lambda2=lambda2, method=method, pick_ori=pick_ori
+    )
 
     # turn source estimate into numpy array
     invmat = invmat_op.data
@@ -538,5 +606,5 @@ def _get_matrix_from_inverse_operator(inverse_operator, forward, method='dSPM',
 
 def _check_fixed_ori(inst):
     """Check if inverse or forward was computed for fixed orientations."""
-    is_fixed = inst['source_ori'] != FIFF.FIFFV_MNE_FREE_ORI
+    is_fixed = inst["source_ori"] != FIFF.FIFFV_MNE_FREE_ORI
     return is_fixed
diff --git a/mne/minimum_norm/tests/test_resolution_matrix.py b/mne/minimum_norm/tests/test_resolution_matrix.py
index 95bb7fccc8f..7681ace402f 100755
--- a/mne/minimum_norm/tests/test_resolution_matrix.py
+++ b/mne/minimum_norm/tests/test_resolution_matrix.py
@@ -6,46 +6,58 @@
 # License: BSD-3-Clause
 
 import os.path as op
+
 import numpy as np
-from numpy.testing import (assert_equal, assert_array_almost_equal,
-                           assert_array_equal, assert_allclose)
+from numpy.testing import (
+    assert_allclose,
+    assert_array_almost_equal,
+    assert_array_equal,
+    assert_equal,
+)
 
 import mne
 from mne.datasets import testing
 from mne.minimum_norm import apply_inverse_cov, prepare_inverse_operator
-from mne.minimum_norm.resolution_matrix_tmp import (make_inverse_resolution_matrix,
-                                                    get_cross_talk,
-                                                    get_point_spread,
-                                                    _vertices_for_get_psf_ctf)
+from mne.minimum_norm.resolution_matrix_tmp import (
+    _vertices_for_get_psf_ctf,
+    get_cross_talk,
+    get_point_spread,
+    make_inverse_resolution_matrix,
+)
 
 data_path = testing.data_path(download=False)
-subjects_dir = op.join(data_path, 'subjects')
-fname_inv = op.join(data_path, 'MEG', 'sample',
-                    'sample_audvis_trunc-meg-eeg-oct-6-meg-inv.fif')
-fname_evoked = op.join(data_path, 'MEG', 'sample',
-                       'sample_audvis_trunc-ave.fif')
-fname_raw = op.join(data_path, 'MEG', 'sample', 'sample_audvis_trunc_raw.fif')
-fname_t1 = op.join(data_path, 'subjects', 'sample', 'mri', 'T1.mgz')
-fname_src = op.join(data_path, 'MEG', 'sample',
-                    'sample_audvis_trunc-meg-eeg-oct-6-fwd.fif')
-fname_src_fs = op.join(data_path, 'subjects', 'fsaverage', 'bem',
-                       'fsaverage-ico-5-src.fif')
-fname_src_3 = op.join(data_path, 'subjects', 'sample', 'bem',
-                      'sample-oct-4-src.fif')
-fname_stc = op.join(data_path, 'MEG', 'sample', 'sample_audvis_trunc-meg')
-fname_vol = op.join(data_path, 'MEG', 'sample',
-                    'sample_audvis_trunc-grad-vol-7-fwd-sensmap-vol.w')
-fname_vsrc = op.join(data_path, 'MEG', 'sample',
-                     'sample_audvis_trunc-meg-vol-7-fwd.fif')
-fname_inv_vol = op.join(data_path, 'MEG', 'sample',
-                        'sample_audvis_trunc-meg-vol-7-meg-inv.fif')
+subjects_dir = op.join(data_path, "subjects")
+fname_inv = op.join(
+    data_path, "MEG", "sample", "sample_audvis_trunc-meg-eeg-oct-6-meg-inv.fif"
+)
+fname_evoked = op.join(data_path, "MEG", "sample", "sample_audvis_trunc-ave.fif")
+fname_raw = op.join(data_path, "MEG", "sample", "sample_audvis_trunc_raw.fif")
+fname_t1 = op.join(data_path, "subjects", "sample", "mri", "T1.mgz")
+fname_src = op.join(
+    data_path, "MEG", "sample", "sample_audvis_trunc-meg-eeg-oct-6-fwd.fif"
+)
+fname_src_fs = op.join(
+    data_path, "subjects", "fsaverage", "bem", "fsaverage-ico-5-src.fif"
+)
+fname_src_3 = op.join(data_path, "subjects", "sample", "bem", "sample-oct-4-src.fif")
+fname_stc = op.join(data_path, "MEG", "sample", "sample_audvis_trunc-meg")
+fname_vol = op.join(
+    data_path, "MEG", "sample", "sample_audvis_trunc-grad-vol-7-fwd-sensmap-vol.w"
+)
+fname_vsrc = op.join(
+    data_path, "MEG", "sample", "sample_audvis_trunc-meg-vol-7-fwd.fif"
+)
+fname_inv_vol = op.join(
+    data_path, "MEG", "sample", "sample_audvis_trunc-meg-vol-7-meg-inv.fif"
+)
 rng = np.random.RandomState(0)
 
-fname_fwd = op.join(data_path, 'MEG', 'sample',
-                    'sample_audvis_trunc-meg-eeg-oct-4-fwd.fif')
-fname_cov = op.join(data_path, 'MEG', 'sample', 'sample_audvis_trunc-cov.fif')
+fname_fwd = op.join(
+    data_path, "MEG", "sample", "sample_audvis_trunc-meg-eeg-oct-4-fwd.fif"
+)
+fname_cov = op.join(data_path, "MEG", "sample", "sample_audvis_trunc-cov.fif")
 
-fname_label = op.join(data_path, 'subjects', 'sample', 'label', 'lh.V1.label')
+fname_label = op.join(data_path, "subjects", "sample", "label", "lh.V1.label")
 
 
 @testing.requires_testing_data
@@ -54,8 +66,7 @@ def test_resolution_matrix():
     # read forward solution
     forward = mne.read_forward_solution(fname_fwd)
     # forward operator with fixed source orientations
-    forward_fxd = mne.convert_forward_solution(forward, surf_ori=True,
-                                               force_fixed=True)
+    forward_fxd = mne.convert_forward_solution(forward, surf_ori=True, force_fixed=True)
 
     # noise covariance matrix
     noise_cov = mne.read_cov(fname_cov)
@@ -65,32 +76,39 @@ def test_resolution_matrix():
     # make inverse operator from forward solution
     # free source orientation
     inverse_operator = mne.minimum_norm.make_inverse_operator(
-        info=evoked.info, forward=forward, noise_cov=noise_cov, loose=1.,
-        depth=None)
+        info=evoked.info, forward=forward, noise_cov=noise_cov, loose=1.0, depth=None
+    )
     # fixed source orientation
     inverse_operator_fxd = mne.minimum_norm.make_inverse_operator(
-        info=evoked.info, forward=forward, noise_cov=noise_cov, loose=0.,
-        depth=None, fixed=True)
+        info=evoked.info,
+        forward=forward,
+        noise_cov=noise_cov,
+        loose=0.0,
+        depth=None,
+        fixed=True,
+    )
 
     # regularisation parameter based on SNR
     snr = 3.0
-    lambda2 = 1.0 / snr ** 2
+    lambda2 = 1.0 / snr**2
     # resolution matrices for free source orientation
     # compute resolution matrix for MNE with free source orientations
-    rm_mne_free = make_inverse_resolution_matrix(forward, inverse_operator,
-                                                 method='MNE', lambda2=lambda2)
+    rm_mne_free = make_inverse_resolution_matrix(
+        forward, inverse_operator, method="MNE", lambda2=lambda2
+    )
     # compute resolution matrix for MNE, fwd fixed and inv free
-    rm_mne_fxdfree = make_inverse_resolution_matrix(forward_fxd,
-                                                    inverse_operator,
-                                                    method='MNE',
-                                                    lambda2=lambda2)
+    rm_mne_fxdfree = make_inverse_resolution_matrix(
+        forward_fxd, inverse_operator, method="MNE", lambda2=lambda2
+    )
     # resolution matrices for fixed source orientation
     # compute resolution matrix for MNE
-    rm_mne = make_inverse_resolution_matrix(forward_fxd, inverse_operator_fxd,
-                                            method='MNE', lambda2=lambda2)
+    rm_mne = make_inverse_resolution_matrix(
+        forward_fxd, inverse_operator_fxd, method="MNE", lambda2=lambda2
+    )
     # compute resolution matrix for sLORETA
-    rm_lor = make_inverse_resolution_matrix(forward_fxd, inverse_operator_fxd,
-                                            method='sLORETA', lambda2=lambda2)
+    rm_lor = make_inverse_resolution_matrix(
+        forward_fxd, inverse_operator_fxd, method="sLORETA", lambda2=lambda2
+    )
     # rectify resolution matrix for sLORETA before determining maxima
     rm_lor_abs = np.abs(rm_lor)
 
@@ -109,48 +127,83 @@ def test_resolution_matrix():
     # Some arbitrary vertex numbers
     idx = [1, 100, 400]
     # check various summary and normalisation options
-    for mode in [None, 'sum', 'mean', 'maxval', 'maxnorm', 'pca']:
+    for mode in [None, "sum", "mean", "maxval", "maxnorm", "pca"]:
         n_comps = [1, 3]
-        if mode in [None, 'sum', 'mean']:
+        if mode in [None, "sum", "mean"]:
             n_comps = [1]
         for n_comp in n_comps:
-            for norm in [None, 'max', 'norm', True]:
+            for norm in [None, "max", "norm", True]:
                 stc_psf = get_point_spread(
-                    rm_mne, forward_fxd['src'], idx, mode=mode, n_comp=n_comp,
-                    norm=norm, return_pca_vars=False)
+                    rm_mne,
+                    forward_fxd["src"],
+                    idx,
+                    mode=mode,
+                    n_comp=n_comp,
+                    norm=norm,
+                    return_pca_vars=False,
+                )
                 stc_ctf = get_cross_talk(
-                    rm_mne, forward_fxd['src'], idx, mode=mode, n_comp=n_comp,
-                    norm=norm, return_pca_vars=False)
+                    rm_mne,
+                    forward_fxd["src"],
+                    idx,
+                    mode=mode,
+                    n_comp=n_comp,
+                    norm=norm,
+                    return_pca_vars=False,
+                )
                 # for MNE, PSF/CTFs for same vertices should be the same
                 assert_array_almost_equal(stc_psf.data, stc_ctf.data)
                 # with free orientations and vector source estimates
                 stc_psf_free = get_point_spread(
-                    rm_mne_free, forward_fxd['src'], idx, mode=mode,
-                    n_comp=n_comp, norm=norm, return_pca_vars=False,
-                    vector=True)
+                    rm_mne_free,
+                    forward_fxd["src"],
+                    idx,
+                    mode=mode,
+                    n_comp=n_comp,
+                    norm=norm,
+                    return_pca_vars=False,
+                    vector=True,
+                )
                 stc_ctf_free = get_cross_talk(
-                    rm_mne_free, forward_fxd['src'], idx, mode=mode,
-                    n_comp=n_comp, norm=norm, return_pca_vars=False,
-                    vector=True)
+                    rm_mne_free,
+                    forward_fxd["src"],
+                    idx,
+                    mode=mode,
+                    n_comp=n_comp,
+                    norm=norm,
+                    return_pca_vars=False,
+                    vector=True,
+                )
                 assert_array_almost_equal(stc_psf_free.data, stc_ctf_free.data)
 
     # check SVD variances
     n_comp = 3
     stc_psf, s_vars_psf = get_point_spread(
-        rm_mne, forward_fxd['src'], idx, mode=mode, n_comp=n_comp,
-        norm='norm', return_pca_vars=True)
+        rm_mne,
+        forward_fxd["src"],
+        idx,
+        mode=mode,
+        n_comp=n_comp,
+        norm="norm",
+        return_pca_vars=True,
+    )
     stc_ctf, s_vars_ctf = get_cross_talk(
-        rm_mne, forward_fxd['src'], idx, mode=mode, n_comp=n_comp,
-        norm='norm', return_pca_vars=True)
+        rm_mne,
+        forward_fxd["src"],
+        idx,
+        mode=mode,
+        n_comp=n_comp,
+        norm="norm",
+        return_pca_vars=True,
+    )
     assert_array_almost_equal(s_vars_psf, s_vars_ctf)
     # variances for SVD components should be ordered
     assert s_vars_psf[0] > s_vars_psf[1] > s_vars_psf[2]
     # all variances should sum up to 100
-    assert_allclose(s_vars_psf.sum(), 100.)
+    assert_allclose(s_vars_psf.sum(), 100.0)
 
     # Test application of free inv to fixed fwd
-    assert_equal(rm_mne_fxdfree.shape, (3 * rm_mne.shape[0],
-                 rm_mne.shape[0]))
+    assert_equal(rm_mne_fxdfree.shape, (3 * rm_mne.shape[0], rm_mne.shape[0]))
 
     # Test PSF/CTF for labels
     label = mne.read_label(fname_label)
@@ -158,82 +211,142 @@ def test_resolution_matrix():
     label = [label]
     label2 = 2 * label
     # get relevant vertices in source space
-    verts = _vertices_for_get_psf_ctf(label, forward_fxd['src'])[0]
+    verts = _vertices_for_get_psf_ctf(label, forward_fxd["src"])[0]
 
-    stc_psf_label = get_point_spread(rm_mne, forward_fxd['src'], label,
-                                     norm='max')
+    stc_psf_label = get_point_spread(rm_mne, forward_fxd["src"], label, norm="max")
     # for list of indices
-    stc_psf_idx = get_point_spread(rm_mne, forward_fxd['src'], verts,
-                                   norm='max')
-    stc_ctf_label = get_cross_talk(rm_mne, forward_fxd['src'], label,
-                                   norm='max')
+    stc_psf_idx = get_point_spread(rm_mne, forward_fxd["src"], verts, norm="max")
+    stc_ctf_label = get_cross_talk(rm_mne, forward_fxd["src"], label, norm="max")
     # For MNE, PSF and CTF for same vertices should be the same
     assert_array_almost_equal(stc_psf_label.data, stc_ctf_label.data)
     # same with free orientations and vector source estimates
     stc_psf_label_free = get_point_spread(
-        rm_mne_free, forward_fxd['src'], label, norm='norm', vector=True)
+        rm_mne_free, forward_fxd["src"], label, norm="norm", vector=True
+    )
     stc_ctf_label_free = get_cross_talk(
-        rm_mne_free, forward_fxd['src'], label, norm='norm', vector=True)
+        rm_mne_free, forward_fxd["src"], label, norm="norm", vector=True
+    )
     assert_array_almost_equal(stc_psf_label_free.data, stc_ctf_label_free.data)
 
     # test multiple labels
-    stc_psf_label2 = get_point_spread(rm_mne, forward_fxd['src'], label2,
-                                      norm='max')
+    stc_psf_label2 = get_point_spread(rm_mne, forward_fxd["src"], label2, norm="max")
     m, n = stc_psf_label.data.shape
-    assert_array_equal(
-        stc_psf_label.data, stc_psf_label2[0].data)
-    assert_array_equal(
-        stc_psf_label.data, stc_psf_label2[1].data)
-    assert_array_equal(
-        stc_psf_label.data, stc_psf_idx.data)
+    assert_array_equal(stc_psf_label.data, stc_psf_label2[0].data)
+    assert_array_equal(stc_psf_label.data, stc_psf_label2[1].data)
+    assert_array_equal(stc_psf_label.data, stc_psf_idx.data)
 
     # test noisy resolution matrix
     # resmats without noise
     rm_mne1 = make_inverse_resolution_matrix(
-        forward_fxd, inverse_operator_fxd, method='MNE', lambda2=lambda2,
-        noise_cov=None, snr=None)
+        forward_fxd,
+        inverse_operator_fxd,
+        method="MNE",
+        lambda2=lambda2,
+        noise_cov=None,
+        snr=None,
+    )
     rm_mne2 = make_inverse_resolution_matrix(
-        forward_fxd, inverse_operator_fxd, method='dSPM', lambda2=lambda2,
-        noise_cov=None, snr=None)
+        forward_fxd,
+        inverse_operator_fxd,
+        method="dSPM",
+        lambda2=lambda2,
+        noise_cov=None,
+        snr=None,
+    )
     rm_mne3 = make_inverse_resolution_matrix(
-        forward_fxd, inverse_operator_fxd, method='sLORETA', lambda2=lambda2,
-        noise_cov=None, snr=None)
+        forward_fxd,
+        inverse_operator_fxd,
+        method="sLORETA",
+        lambda2=lambda2,
+        noise_cov=None,
+        snr=None,
+    )
     rm_mne4 = make_inverse_resolution_matrix(
-        forward_fxd, inverse_operator_fxd, method='eLORETA', lambda2=lambda2,
-        noise_cov=None, snr=None)
+        forward_fxd,
+        inverse_operator_fxd,
+        method="eLORETA",
+        lambda2=lambda2,
+        noise_cov=None,
+        snr=None,
+    )
 
     # same resmats with noise but very high SNR
     rm_mne1h = make_inverse_resolution_matrix(
-        forward_fxd, inverse_operator_fxd, method='MNE', lambda2=lambda2,
-        noise_cov=noise_cov, snr=1e20)
+        forward_fxd,
+        inverse_operator_fxd,
+        method="MNE",
+        lambda2=lambda2,
+        noise_cov=noise_cov,
+        snr=1e20,
+    )
     rm_mne2h = make_inverse_resolution_matrix(
-        forward_fxd, inverse_operator_fxd, method='dSPM', lambda2=lambda2,
-        noise_cov=noise_cov, snr=1e20)
+        forward_fxd,
+        inverse_operator_fxd,
+        method="dSPM",
+        lambda2=lambda2,
+        noise_cov=noise_cov,
+        snr=1e20,
+    )
     rm_mne3h = make_inverse_resolution_matrix(
-        forward_fxd, inverse_operator_fxd, method='sLORETA', lambda2=lambda2,
-        noise_cov=noise_cov, snr=1e20)
+        forward_fxd,
+        inverse_operator_fxd,
+        method="sLORETA",
+        lambda2=lambda2,
+        noise_cov=noise_cov,
+        snr=1e20,
+    )
     rm_mne4h = make_inverse_resolution_matrix(
-        forward_fxd, inverse_operator_fxd, method='eLORETA', lambda2=lambda2,
-        noise_cov=noise_cov, snr=1e20)
+        forward_fxd,
+        inverse_operator_fxd,
+        method="eLORETA",
+        lambda2=lambda2,
+        noise_cov=noise_cov,
+        snr=1e20,
+    )
 
     # same resmats with noise and realistic/medium SNR
     rm_mne1m = make_inverse_resolution_matrix(
-        forward_fxd, inverse_operator_fxd, method='MNE', lambda2=lambda2,
-        noise_cov=noise_cov, snr=3)
+        forward_fxd,
+        inverse_operator_fxd,
+        method="MNE",
+        lambda2=lambda2,
+        noise_cov=noise_cov,
+        snr=3,
+    )
     rm_mne2m = make_inverse_resolution_matrix(
-        forward_fxd, inverse_operator_fxd, method='dSPM', lambda2=lambda2,
-        noise_cov=noise_cov, snr=3)
+        forward_fxd,
+        inverse_operator_fxd,
+        method="dSPM",
+        lambda2=lambda2,
+        noise_cov=noise_cov,
+        snr=3,
+    )
     rm_mne3m = make_inverse_resolution_matrix(
-        forward_fxd, inverse_operator_fxd, method='sLORETA', lambda2=lambda2,
-        noise_cov=noise_cov, snr=3)
+        forward_fxd,
+        inverse_operator_fxd,
+        method="sLORETA",
+        lambda2=lambda2,
+        noise_cov=noise_cov,
+        snr=3,
+    )
     rm_mne4m = make_inverse_resolution_matrix(
-        forward_fxd, inverse_operator_fxd, method='eLORETA', lambda2=lambda2,
-        noise_cov=noise_cov, snr=3)
+        forward_fxd,
+        inverse_operator_fxd,
+        method="eLORETA",
+        lambda2=lambda2,
+        noise_cov=noise_cov,
+        snr=3,
+    )
 
     # same resmats with noise and very low SNR
     rm_mne1l = make_inverse_resolution_matrix(
-        forward_fxd, inverse_operator_fxd, method='MNE', lambda2=lambda2,
-        noise_cov=noise_cov, snr=1e-20)
+        forward_fxd,
+        inverse_operator_fxd,
+        method="MNE",
+        lambda2=lambda2,
+        noise_cov=noise_cov,
+        snr=1e-20,
+    )
 
     # resmats should be very close without noise and at high SNR
     assert_array_almost_equal(np.abs(rm_mne1), rm_mne1h)
@@ -242,43 +355,73 @@ def test_resolution_matrix():
     assert_array_almost_equal(np.abs(rm_mne4), rm_mne4h)
 
     # rows of resmats (CTFs) should be similar between methods
-    row_corrs = [np.corrcoef(rm_mne1[i, :], rm_mne2[i, :])[0, 1]
-                 for i in np.arange(rm_mne1.shape[0])]
+    row_corrs = [
+        np.corrcoef(rm_mne1[i, :], rm_mne2[i, :])[0, 1]
+        for i in np.arange(rm_mne1.shape[0])
+    ]
     assert_allclose(row_corrs, 1)
-    row_corrs = [np.corrcoef(rm_mne1[i, :], rm_mne3[i, :])[0, 1]
-                 for i in np.arange(rm_mne1.shape[0])]
+    row_corrs = [
+        np.corrcoef(rm_mne1[i, :], rm_mne3[i, :])[0, 1]
+        for i in np.arange(rm_mne1.shape[0])
+    ]
     assert_allclose(row_corrs, 1)
-    row_corrs = [np.corrcoef(rm_mne1[i, :], rm_mne4[i, :])[0, 1]
-                 for i in np.arange(rm_mne1.shape[0])]
+    row_corrs = [
+        np.corrcoef(rm_mne1[i, :], rm_mne4[i, :])[0, 1]
+        for i in np.arange(rm_mne1.shape[0])
+    ]
     assert_allclose(row_corrs, 1, 0.2)  # eLORETA is different
 
     # noisy:
-    row_corrs = [np.corrcoef(rm_mne1h[i, :], rm_mne2h[i, :])[0, 1]
-                 for i in np.arange(rm_mne1.shape[0])]
+    row_corrs = [
+        np.corrcoef(rm_mne1h[i, :], rm_mne2h[i, :])[0, 1]
+        for i in np.arange(rm_mne1.shape[0])
+    ]
     assert_allclose(row_corrs, 1)
-    row_corrs = [np.corrcoef(rm_mne1h[i, :], rm_mne3h[i, :])[0, 1]
-                 for i in np.arange(rm_mne1.shape[0])]
+    row_corrs = [
+        np.corrcoef(rm_mne1h[i, :], rm_mne3h[i, :])[0, 1]
+        for i in np.arange(rm_mne1.shape[0])
+    ]
     assert_allclose(row_corrs, 1)
-    row_corrs = [np.corrcoef(rm_mne1h[i, :], rm_mne4h[i, :])[0, 1]
-                 for i in np.arange(rm_mne1.shape[0])]
+    row_corrs = [
+        np.corrcoef(rm_mne1h[i, :], rm_mne4h[i, :])[0, 1]
+        for i in np.arange(rm_mne1.shape[0])
+    ]
     assert_allclose(row_corrs, 1, 0.2)  # eLORETA is different
 
-    row_corrs = [np.corrcoef(rm_mne1m[i, :], rm_mne2m[i, :])[0, 1]
-                 for i in np.arange(rm_mne1.shape[0])]
+    row_corrs = [
+        np.corrcoef(rm_mne1m[i, :], rm_mne2m[i, :])[0, 1]
+        for i in np.arange(rm_mne1.shape[0])
+    ]
     assert_allclose(row_corrs, 1)
-    row_corrs = [np.corrcoef(rm_mne1m[i, :], rm_mne3m[i, :])[0, 1]
-                 for i in np.arange(rm_mne1.shape[0])]
+    row_corrs = [
+        np.corrcoef(rm_mne1m[i, :], rm_mne3m[i, :])[0, 1]
+        for i in np.arange(rm_mne1.shape[0])
+    ]
     assert_allclose(row_corrs, 1)
-    row_corrs = [np.corrcoef(rm_mne1m[i, :], rm_mne4m[i, :])[0, 1]
-                 for i in np.arange(rm_mne1.shape[0])]
+    row_corrs = [
+        np.corrcoef(rm_mne1m[i, :], rm_mne4m[i, :])[0, 1]
+        for i in np.arange(rm_mne1.shape[0])
+    ]
     assert_allclose(row_corrs, 1, 0.2)  # eLORETA is different
 
     # for very low SNR, resmat should be equal to noise power (except for scaling)
-    inv = prepare_inverse_operator(inverse_operator_fxd, 1, lambda2, 'MNE', copy='non-src')
+    inv = prepare_inverse_operator(
+        inverse_operator_fxd, 1, lambda2, "MNE", copy="non-src"
+    )
     stc = apply_inverse_cov(
-        noise_cov, evoked.info, inv, nave=1, lambda2=lambda2,
-        method='MNE', pick_ori=None, prepared=True, label=None,
-        method_params=None, use_cps=True, verbose=None)
+        noise_cov,
+        evoked.info,
+        inv,
+        nave=1,
+        lambda2=lambda2,
+        method="MNE",
+        pick_ori=None,
+        prepared=True,
+        label=None,
+        method_params=None,
+        use_cps=True,
+        verbose=None,
+    )
 
     x = np.sqrt(stc.data)
     y = rm_mne1l

From 09bc7e519f20b499b567e025d583b8c05fea4550 Mon Sep 17 00:00:00 2001
From: olafhauk <olaf.hauk@mrc-cbu.cam.ac.uk>
Date: Wed, 15 Nov 2023 17:01:02 +0000
Subject: [PATCH 13/27] WIP: removed prepare_inverse_operator()

---
 mne/minimum_norm/resolution_matrix.py | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/mne/minimum_norm/resolution_matrix.py b/mne/minimum_norm/resolution_matrix.py
index 754684f2f8c..0be73b9bb0c 100755
--- a/mne/minimum_norm/resolution_matrix.py
+++ b/mne/minimum_norm/resolution_matrix.py
@@ -97,7 +97,6 @@ def make_inverse_resolution_matrix(
 
     # add source noise power to columns of resolution matrix
     if noise_cov is not None:
-        inv = prepare_inverse_operator(inv, 1, lambda2, method, copy="non-src")
         info = _prepare_info(inv)
         # compute source noise power
         stc = apply_inverse_cov(
@@ -108,7 +107,7 @@ def make_inverse_resolution_matrix(
             lambda2=lambda2,
             method=method,
             pick_ori=None,
-            prepared=True,
+            prepared=False,
             label=None,
             method_params=None,
             use_cps=True,

From 04aba165946698d893dfa4a71694d1a03580ba1d Mon Sep 17 00:00:00 2001
From: "pre-commit-ci[bot]"
 <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Date: Wed, 15 Nov 2023 17:03:00 +0000
Subject: [PATCH 14/27] [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci
---
 mne/minimum_norm/resolution_matrix.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/mne/minimum_norm/resolution_matrix.py b/mne/minimum_norm/resolution_matrix.py
index 0be73b9bb0c..4ea557a2cf1 100755
--- a/mne/minimum_norm/resolution_matrix.py
+++ b/mne/minimum_norm/resolution_matrix.py
@@ -17,7 +17,7 @@
 from ..forward.forward import convert_forward_solution
 from ..io.constants import FIFF
 from ..label import Label
-from ..minimum_norm import apply_inverse, apply_inverse_cov, prepare_inverse_operator
+from ..minimum_norm import apply_inverse, apply_inverse_cov
 from ..minimum_norm.spatial_resolution import _rectify_resolution_matrix
 from ..source_estimate import _prepare_label_extraction
 from ..utils import logger, verbose

From 407c1b0d9e83e0edb2485acecf3fd30962b0aaed Mon Sep 17 00:00:00 2001
From: olafhauk <olaf.hauk@mrc-cbu.cam.ac.uk>
Date: Thu, 16 Nov 2023 11:27:49 +0000
Subject: [PATCH 15/27] WIP: minor changes to tests

---
 mne/minimum_norm/tests/test_resolution_matrix.py | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

diff --git a/mne/minimum_norm/tests/test_resolution_matrix.py b/mne/minimum_norm/tests/test_resolution_matrix.py
index 7681ace402f..1a8018f1d27 100755
--- a/mne/minimum_norm/tests/test_resolution_matrix.py
+++ b/mne/minimum_norm/tests/test_resolution_matrix.py
@@ -348,7 +348,7 @@ def test_resolution_matrix():
         snr=1e-20,
     )
 
-    # resmats should be very close without noise and at high SNR
+    # resolution matrices should be very close without noise and at high SNR
     assert_array_almost_equal(np.abs(rm_mne1), rm_mne1h)
     assert_array_almost_equal(np.abs(rm_mne2), rm_mne2h)
     assert_array_almost_equal(np.abs(rm_mne3), rm_mne3h)
@@ -371,7 +371,7 @@ def test_resolution_matrix():
     ]
     assert_allclose(row_corrs, 1, 0.2)  # eLORETA is different
 
-    # noisy:
+    # for noisy resmats:
     row_corrs = [
         np.corrcoef(rm_mne1h[i, :], rm_mne2h[i, :])[0, 1]
         for i in np.arange(rm_mne1.shape[0])
@@ -408,6 +408,7 @@ def test_resolution_matrix():
     inv = prepare_inverse_operator(
         inverse_operator_fxd, 1, lambda2, "MNE", copy="non-src"
     )
+    # compute noise power in source space (what you get at very low SNR)
     stc = apply_inverse_cov(
         noise_cov,
         evoked.info,
@@ -423,6 +424,7 @@ def test_resolution_matrix():
         verbose=None,
     )
 
+    # scale and compare noise power to all PSFs
     x = np.sqrt(stc.data)
     y = rm_mne1l
     x = x / x.max()

From ff2f3e9e8774b6bdc65d4d09e47cbfe7d3c82b03 Mon Sep 17 00:00:00 2001
From: olafhauk <olaf.hauk@mrc-cbu.cam.ac.uk>
Date: Fri, 17 Nov 2023 12:30:10 +0000
Subject: [PATCH 16/27] WIP: minor formatting

---
 mne/minimum_norm/resolution_matrix.py | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/mne/minimum_norm/resolution_matrix.py b/mne/minimum_norm/resolution_matrix.py
index 4ea557a2cf1..9adab004637 100755
--- a/mne/minimum_norm/resolution_matrix.py
+++ b/mne/minimum_norm/resolution_matrix.py
@@ -68,7 +68,7 @@ def make_inverse_resolution_matrix(
         Note: If the resolution matrix is computed with a noise covariance
         matrix then only its columns, i.e. PSFs, can meaningfully be
         interpreted. It must not be used to compute CTFs or resolution metrics
-        for CTFs!
+        for CTFs.
     """
     if noise_cov is None and snr is not None:
         msg = "snr should be None if noise_cov is None."
@@ -237,7 +237,7 @@ def _get_psf_ctf(
 
     for verts in verts_all:
         # get relevant PSFs or CTFs for specified vertices
-        if type(verts) is int:
+        if isinstance(verts, int):
             verts = [verts]  # to keep array dimensions
         funcs = resmat[:, verts]
 
@@ -293,7 +293,7 @@ def _vertices_for_get_psf_ctf(idx, src):
     """Get vertices in source space for PSFs/CTFs in _get_psf_ctf()."""
     # idx must be list
     # if label(s) specified get the indices, otherwise just carry on
-    if type(idx[0]) is Label:
+    if isinstance(idx[0], Label):
         # specify without source time courses, gets indices per label
         verts_labs, _ = _prepare_label_extraction(
             stc=None,
@@ -310,7 +310,7 @@ def _vertices_for_get_psf_ctf(idx, src):
 
         for v in verts_labs:
             # if two hemispheres present
-            if type(v) is list:
+            if isinstance(v, list):
                 # indices for both hemispheres in one list
                 this_verts = np.concatenate((v[0], v[1]))
             else:
@@ -318,7 +318,7 @@ def _vertices_for_get_psf_ctf(idx, src):
             verts.append(this_verts)
     # check if list of list or just list
     else:
-        if type(idx[0]) is list:  # if list of list of integers
+        if isinstance(idx[0], list):  # if list of list of integers
             verts = idx
         else:  # if list of integers
             verts = [idx]

From 4aa27ac08c9b741db34bb6a78b81d3b1788bf489 Mon Sep 17 00:00:00 2001
From: olafhauk <olaf.hauk@mrc-cbu.cam.ac.uk>
Date: Fri, 17 Nov 2023 17:32:02 +0000
Subject: [PATCH 17/27] WIP: added example

---
 examples/inverse/psf_empirical.py | 143 ++++++++++++++++++++++++++++++
 1 file changed, 143 insertions(+)
 create mode 100755 examples/inverse/psf_empirical.py

diff --git a/examples/inverse/psf_empirical.py b/examples/inverse/psf_empirical.py
new file mode 100755
index 00000000000..0a84f3b798b
--- /dev/null
+++ b/examples/inverse/psf_empirical.py
@@ -0,0 +1,143 @@
+# -*- coding: utf-8 -*-
+"""
+.. _ex-psd-ctf:
+
+===================================================
+Plot point-spread functions (PSFs) with added noise
+===================================================
+
+Visualise PSF with added noise at one vertex for sLORETA.
+"""
+# Authors: Olaf Hauk <olaf.hauk@mrc-cbu.cam.ac.uk>
+#          Alexandre Gramfort <alexandre.gramfort@inria.fr>
+#
+# License: BSD-3-Clause
+
+import mne
+from mne.datasets import sample
+from mne.minimum_norm import get_point_spread, make_inverse_resolution_matrix
+
+print(__doc__)
+
+data_path = sample.data_path()
+subjects_dir = data_path / "subjects"
+meg_path = data_path / "MEG" / "sample"
+fname_fwd = meg_path / "sample_audvis-meg-eeg-oct-6-fwd.fif"
+fname_cov = meg_path / "sample_audvis-cov.fif"
+fname_evo = meg_path / "sample_audvis-ave.fif"
+
+# read forward solution
+forward = mne.read_forward_solution(fname_fwd)
+# forward operator with fixed source orientations
+mne.convert_forward_solution(forward, surf_ori=True, force_fixed=True, copy=False)
+
+# noise covariance matrix
+noise_cov = mne.read_cov(fname_cov)
+
+# evoked data for info
+info = mne.io.read_info(fname_evo)
+
+# make inverse operator from forward solution
+# free source orientation
+inverse_operator = mne.minimum_norm.make_inverse_operator(
+    info=info, forward=forward, noise_cov=noise_cov, loose=0.0, depth=None
+)
+
+# regularisation parameter
+snr = 3.0
+lambda2 = 1.0 / snr**2
+method = "MNE"
+
+# %%
+# We compute noisy PSFs for two locations, one deep and one superficial.
+# For each location we compute two PSFs:
+# 1) For a high SNR, approximating the noise-less case.
+# 2) With the SNR of regularization and the added noise matched (realistic).
+
+# vertices of a deep and a superficial source
+sources = {"deep": 2146, "super": 2209}
+
+stc_psf = {"deep": {}, "super": {}}
+
+# for a deep source
+# compute resolution matrix with added noise for columns (PSFs)
+# added noise matches the SNR used for regularization
+rm_matched = make_inverse_resolution_matrix(
+    forward, inverse_operator, method=method, lambda2=lambda2, noise_cov=noise_cov
+)
+stc_psf["deep"]["matched"] = get_point_spread(
+    rm_matched, forward["src"], [sources["deep"]], norm=True
+)
+del rm_matched
+
+# compute resolution matrix with added noise for very high SNR
+rm_highsnr = make_inverse_resolution_matrix(
+    forward,
+    inverse_operator,
+    method=method,
+    lambda2=lambda2,
+    noise_cov=noise_cov,
+    snr=1000.0,
+)
+stc_psf["deep"]["highsnr"] = get_point_spread(
+    rm_highsnr, forward["src"], [sources["deep"]], norm=True
+)
+del rm_highsnr
+
+# the same for a superficial source
+# compute matching resolution matrix with added noise for columns (PSFs)
+# added noise matches the SNR used for regularization
+rm_matched = make_inverse_resolution_matrix(
+    forward, inverse_operator, method=method, lambda2=lambda2, noise_cov=noise_cov
+)
+stc_psf["super"]["matched"] = get_point_spread(
+    rm_matched, forward["src"], [sources["super"]], norm=True
+)
+del rm_matched
+
+# compute resolution matrix with added noise for very high SNR
+rm_highsnr = make_inverse_resolution_matrix(
+    forward,
+    inverse_operator,
+    method=method,
+    lambda2=lambda2,
+    noise_cov=noise_cov,
+    snr=1000.0,
+)
+stc_psf["super"]["highsnr"] = get_point_spread(
+    rm_highsnr, forward["src"], [sources["super"]], norm=True
+)
+del rm_highsnr
+
+##############################################################################
+# Visualize
+# ---------
+# PSF:
+
+vertno_lh = forward["src"][0]["vertno"]
+for source in sources:
+    # Which vertex corresponds to selected source
+    verttrue = [vertno_lh[sources[source]]]  # just one vertex
+    for psf in ["matched", "highsnr"]:
+        brain_psf = stc_psf[source][psf].plot(
+            "sample", "inflated", "lh", subjects_dir=subjects_dir
+        )
+        brain_psf.show_view("ventral")
+        title_str = "%s %s %s" % (method, source, psf)
+        brain_psf.add_text(0.1, 0.9, title_str, "title", font_size=16)
+        # mark true source location
+        brain_psf.add_foci(
+            verttrue, coords_as_verts=True, scale_factor=1.0, hemi="lh", color="green"
+        )
+
+# %%
+# The green spheres indicate the true source location.
+# For high SNR, PSFs for neither location are ideal. For the superficial
+# it is at least centred around the true source location, but still widespread
+# and with sidelobes. For the deep location the peaks occur at large distances
+# to the true location. This shows that spatial resolution is not an issue of
+# insufficient SNR but is inherently limited.
+# For matched (more realistic) SNR the PSF resembles the noise-less PSF, albeit
+# with more widely spread-out sidelobes. This indicates that moderate
+# regularization still achieves good spatial resolution while effectively
+# suppressing noise.

From ef7a152f8d71b7afc700b47eb82d1b453148a93f Mon Sep 17 00:00:00 2001
From: Marijn van Vliet <w.m.vanvliet@gmail.com>
Date: Mon, 4 Mar 2024 09:39:58 +0200
Subject: [PATCH 18/27] Fix import

---
 mne/minimum_norm/resolution_matrix.py | 16 ++++++++--------
 1 file changed, 8 insertions(+), 8 deletions(-)

diff --git a/mne/minimum_norm/resolution_matrix.py b/mne/minimum_norm/resolution_matrix.py
index 3f5ef129052..9b4f12120fc 100755
--- a/mne/minimum_norm/resolution_matrix.py
+++ b/mne/minimum_norm/resolution_matrix.py
@@ -1,4 +1,5 @@
 """Compute resolution matrix for linear estimators."""
+
 # Authors: olaf.hauk@mrc-cbu.cam.ac.uk
 #
 # License: BSD-3-Clause
@@ -7,19 +8,18 @@
 
 import numpy as np
 
-from .. import (
-    EvokedArray,
-    SourceEstimate,
-    VectorSourceEstimate,
-    pick_channels_forward,
-    pick_types,
-)
+from .._fiff.pick import pick_channels_forward, pick_types
+from ..evoked import EvokedArray
 from ..forward.forward import convert_forward_solution
 from ..io.constants import FIFF
 from ..label import Label
 from ..minimum_norm import apply_inverse, apply_inverse_cov
 from ..minimum_norm.spatial_resolution import _rectify_resolution_matrix
-from ..source_estimate import _prepare_label_extraction
+from ..source_estimate import (
+    SourceEstimate,
+    VectorSourceEstimate,
+    _prepare_label_extraction,
+)
 from ..utils import logger, verbose
 
 

From e9f90a07279864d1cc2203ea91e28962c0cb383d Mon Sep 17 00:00:00 2001
From: Marijn van Vliet <w.m.vanvliet@gmail.com>
Date: Fri, 5 Apr 2024 09:44:39 +0300
Subject: [PATCH 19/27] Put back line that was erroneously removed.

---
 mne/minimum_norm/resolution_matrix.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/mne/minimum_norm/resolution_matrix.py b/mne/minimum_norm/resolution_matrix.py
index 9b4f12120fc..7b1e17027a7 100755
--- a/mne/minimum_norm/resolution_matrix.py
+++ b/mne/minimum_norm/resolution_matrix.py
@@ -512,6 +512,7 @@ def _prepare_info(inverse_operator):
     with info._unlock():
         info["sfreq"] = 1000.0  # necessary
         info["projs"] = inverse_operator["projs"]
+        info["custom_ref_applied"] = False
     return info
 
 

From 776fe69c00657a1f7a39804cfac0eb0be0a19d5f Mon Sep 17 00:00:00 2001
From: Marijn van Vliet <w.m.vanvliet@gmail.com>
Date: Fri, 5 Apr 2024 09:50:47 +0300
Subject: [PATCH 20/27] Some cosmetic changes.

---
 mne/minimum_norm/resolution_matrix.py | 9 +--------
 1 file changed, 1 insertion(+), 8 deletions(-)

diff --git a/mne/minimum_norm/resolution_matrix.py b/mne/minimum_norm/resolution_matrix.py
index 7b1e17027a7..c1486000c4f 100755
--- a/mne/minimum_norm/resolution_matrix.py
+++ b/mne/minimum_norm/resolution_matrix.py
@@ -140,15 +140,8 @@ def make_inverse_resolution_matrix(
 
         alpha = (3 * snr) / (alphas["eeg"] + alphas["mag"] + alphas["gra"])
 
-        # # previous attempt to compute scaling via whitened gain matrix
-        # whitener = inv['whitener']
-        # gain = np.dot(whitener, leadfield)
-        # gram_gain = np.dot(gain, gain.T)
-        # alpha = np.sqrt(np.trace(gram_gain) / gram_gain.shape[0])
-
         # Add square root of source noise power to every column, scale
         # depending on SNR
-        # EDIT
         resmat = resmat + (1 / alpha) * np.sqrt(stc.data)
 
     return resmat
@@ -212,7 +205,7 @@ def _get_psf_ctf(
     ):
         msg = (
             "Number of vertices (%d) and corresponding dimension of"
-            "resolution matrix ((%d, %d) do not match" % (n_verts, n_r, n_c)
+            "resolution matrix (%d, %d) do not match" % (n_verts, n_r, n_c)
         )
         raise ValueError(msg)
 

From f5d816955dba0f5d0347c9a2f202de4dc8cbafee Mon Sep 17 00:00:00 2001
From: Marijn van Vliet <w.m.vanvliet@gmail.com>
Date: Fri, 5 Apr 2024 09:53:14 +0300
Subject: [PATCH 21/27] Fix imports

---
 mne/minimum_norm/resolution_matrix.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/mne/minimum_norm/resolution_matrix.py b/mne/minimum_norm/resolution_matrix.py
index c1486000c4f..7cce5163dc3 100755
--- a/mne/minimum_norm/resolution_matrix.py
+++ b/mne/minimum_norm/resolution_matrix.py
@@ -13,14 +13,14 @@
 from ..forward.forward import convert_forward_solution
 from ..io.constants import FIFF
 from ..label import Label
-from ..minimum_norm import apply_inverse, apply_inverse_cov
-from ..minimum_norm.spatial_resolution import _rectify_resolution_matrix
 from ..source_estimate import (
     SourceEstimate,
     VectorSourceEstimate,
     _prepare_label_extraction,
 )
 from ..utils import logger, verbose
+from .inverse import apply_inverse, apply_inverse_cov
+from .spatial_resolution import _rectify_resolution_matrix
 
 
 @verbose

From 4431eab2349d711294163f239d93af3eb530e498 Mon Sep 17 00:00:00 2001
From: Marijn van Vliet <w.m.vanvliet@gmail.com>
Date: Sat, 6 Apr 2024 23:40:42 +0300
Subject: [PATCH 22/27] style fix

---
 examples/inverse/psf_empirical.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/examples/inverse/psf_empirical.py b/examples/inverse/psf_empirical.py
index 192f23e40a1..8250220914a 100755
--- a/examples/inverse/psf_empirical.py
+++ b/examples/inverse/psf_empirical.py
@@ -124,7 +124,7 @@
             "sample", "inflated", "lh", subjects_dir=subjects_dir
         )
         brain_psf.show_view("ventral")
-        title_str = "%s %s %s" % (method, source, psf)
+        title_str = f"{method} {source} {psf}"
         brain_psf.add_text(
             0.1, 0.9, title_str, "title", font_size=16, color=(255, 255, 255)
         )

From b07ae6e5a3a3f2852de305b7e481bb56b6cd1849 Mon Sep 17 00:00:00 2001
From: Marijn van Vliet <w.m.vanvliet@gmail.com>
Date: Sat, 6 Apr 2024 23:44:28 +0300
Subject: [PATCH 23/27] more style fixes

---
 mne/minimum_norm/resolution_matrix.py | 16 +++++++++-------
 1 file changed, 9 insertions(+), 7 deletions(-)

diff --git a/mne/minimum_norm/resolution_matrix.py b/mne/minimum_norm/resolution_matrix.py
index 7cce5163dc3..040d6040e3d 100755
--- a/mne/minimum_norm/resolution_matrix.py
+++ b/mne/minimum_norm/resolution_matrix.py
@@ -76,7 +76,7 @@ def make_inverse_resolution_matrix(
 
     if noise_cov is not None and snr is None:
         snr = np.sqrt(1.0 / lambda2)
-        logger.info("Inferring snr from lamdba2 as: %f." % snr)
+        logger.info(f"Inferring snr from lamdba2 as: {snr}.")
 
     # make sure forward and inverse operator match
     inv = inverse_operator
@@ -93,7 +93,9 @@ def make_inverse_resolution_matrix(
     leadfield = fwd["sol"]["data"]
     invmat = _get_matrix_from_inverse_operator(inv, fwd, method=method, lambda2=lambda2)
     resmat = invmat.dot(leadfield)
-    logger.info("Dimensions of resolution matrix: %d by %d." % resmat.shape)
+    logger.info(
+        f"Dimensions of resolution matrix: {resmat.shape[0]} by {resmat.shape[1]}."
+    )
 
     # add source noise power to columns of resolution matrix
     if noise_cov is not None:
@@ -204,8 +206,8 @@ def _get_psf_ctf(
         (n_verts != n_c) and (n_c / 3 != n_verts)
     ):
         msg = (
-            "Number of vertices (%d) and corresponding dimension of"
-            "resolution matrix (%d, %d) do not match" % (n_verts, n_r, n_c)
+            f"Number of vertices ({n_verts}) and corresponding dimension of "
+            f"resolution matrix ({n_r}, {n_c}) do not match"
         )
         raise ValueError(msg)
 
@@ -275,10 +277,10 @@ def _get_psf_ctf(
 def _check_get_psf_ctf_params(mode, n_comp, return_pca_vars):
     """Check input parameters of _get_psf_ctf() for consistency."""
     if mode in [None, "sum", "mean"] and n_comp > 1:
-        msg = "n_comp must be 1 for mode=%s." % mode
+        msg = f"n_comp must be 1 for mode={mode}."
         raise ValueError(msg)
     if mode != "pca" and return_pca_vars:
-        msg = "SVD variances can only be returned if mode=" "pca" "."
+        msg = 'SVD variances can only be returned if mode="pca".'
         raise ValueError(msg)
 
 
@@ -592,7 +594,7 @@ def _get_matrix_from_inverse_operator(
         assert np.array_equal(v0o1, invmat[1])
         assert np.array_equal(v3o2, invmat[11])
 
-    logger.info("Dimension of Inverse Matrix: %s" % str(invmat.shape))
+    logger.info(f"Dimension of Inverse Matrix: {invmat.shape}")
 
     return invmat
 

From a1c1a92b6d398967bd87da4ee447ed5cdd8e19aa Mon Sep 17 00:00:00 2001
From: Marijn van Vliet <w.m.vanvliet@gmail.com>
Date: Fri, 12 Apr 2024 09:33:31 +0300
Subject: [PATCH 24/27] fix another import

---
 mne/minimum_norm/tests/test_resolution_matrix.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/mne/minimum_norm/tests/test_resolution_matrix.py b/mne/minimum_norm/tests/test_resolution_matrix.py
index d62da8f18cc..8f27b1a9057 100755
--- a/mne/minimum_norm/tests/test_resolution_matrix.py
+++ b/mne/minimum_norm/tests/test_resolution_matrix.py
@@ -18,7 +18,7 @@
 import mne
 from mne.datasets import testing
 from mne.minimum_norm import apply_inverse_cov, prepare_inverse_operator
-from mne.minimum_norm.resolution_matrix_tmp import (
+from mne.minimum_norm.resolution_matrix import (
     _vertices_for_get_psf_ctf,
     get_cross_talk,
     get_point_spread,

From ba99e83ed2ebd5a7cb55671307bbbb723793551c Mon Sep 17 00:00:00 2001
From: Marijn van Vliet <w.m.vanvliet@gmail.com>
Date: Fri, 12 Apr 2024 13:29:34 +0300
Subject: [PATCH 25/27] doc fixes

---
 doc/changes/devel/11092.newfeature.rst | 2 +-
 examples/inverse/psf_empirical.py      | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/doc/changes/devel/11092.newfeature.rst b/doc/changes/devel/11092.newfeature.rst
index 9fe90afaae7..c325ebcd0bd 100644
--- a/doc/changes/devel/11092.newfeature.rst
+++ b/doc/changes/devel/11092.newfeature.rst
@@ -1 +1 @@
-- Added noisy point-spread functions to mne.minimum_norm.make_inverse_resolution_matrix (:gh:`11058` by `Olaf Hauk`)
+- Added noisy point-spread functions to mne.minimum_norm.make_inverse_resolution_matrix (:gh:`11058` by `Olaf Hauk`_)
diff --git a/examples/inverse/psf_empirical.py b/examples/inverse/psf_empirical.py
index 8250220914a..fa330c8339e 100755
--- a/examples/inverse/psf_empirical.py
+++ b/examples/inverse/psf_empirical.py
@@ -1,5 +1,5 @@
 """
-.. _ex-psd-ctf:
+.. _ex-psf-empirical:
 
 ===================================================
 Plot point-spread functions (PSFs) with added noise

From b3ed287cebc22b5e35b19f6f3a6707334a2a7523 Mon Sep 17 00:00:00 2001
From: "autofix-ci[bot]" <114827586+autofix-ci[bot]@users.noreply.github.com>
Date: Thu, 29 Jan 2026 20:47:10 +0000
Subject: [PATCH 26/27] [autofix.ci] apply automated fixes

---
 examples/inverse/psf_empirical.py | 293 +++++++++++++++---------------
 1 file changed, 147 insertions(+), 146 deletions(-)

diff --git a/examples/inverse/psf_empirical.py b/examples/inverse/psf_empirical.py
index fa330c8339e..e0749ea54a2 100755
--- a/examples/inverse/psf_empirical.py
+++ b/examples/inverse/psf_empirical.py
@@ -1,146 +1,147 @@
-"""
-.. _ex-psf-empirical:
-
-===================================================
-Plot point-spread functions (PSFs) with added noise
-===================================================
-
-Visualise the point-spread of a vertex in a sLORETA source estimate when we assume the
-data to be noisy.
-"""
-# Authors: Olaf Hauk <olaf.hauk@mrc-cbu.cam.ac.uk>
-#          Alexandre Gramfort <alexandre.gramfort@inria.fr>
-#
-# License: BSD-3-Clause
-
-import mne
-from mne.datasets import sample
-from mne.minimum_norm import get_point_spread, make_inverse_resolution_matrix
-
-print(__doc__)
-
-data_path = sample.data_path()
-subjects_dir = data_path / "subjects"
-meg_path = data_path / "MEG" / "sample"
-fname_fwd = meg_path / "sample_audvis-meg-eeg-oct-6-fwd.fif"
-fname_cov = meg_path / "sample_audvis-cov.fif"
-fname_evo = meg_path / "sample_audvis-ave.fif"
-
-# read forward solution
-forward = mne.read_forward_solution(fname_fwd)
-# forward operator with fixed source orientations
-mne.convert_forward_solution(forward, surf_ori=True, force_fixed=True, copy=False)
-
-# noise covariance matrix
-noise_cov = mne.read_cov(fname_cov)
-
-# evoked data for info
-info = mne.io.read_info(fname_evo)
-
-# make inverse operator from forward solution
-# free source orientation
-inverse_operator = mne.minimum_norm.make_inverse_operator(
-    info=info, forward=forward, noise_cov=noise_cov, loose=0.0, depth=None
-)
-
-# regularisation parameter
-snr = 3.0
-lambda2 = 1.0 / snr**2
-method = "MNE"
-
-# %%
-# We compute noisy PSFs for two locations, one deep and one superficial.
-# For each location we compute two PSFs::
-# 1) For a high SNR, approximating the noise-less case.
-# 2) With the SNR of regularization and the added noise matched (realistic).
-
-# vertices of a deep and a superficial source
-sources = {"deep": 2146, "super": 2209}
-
-# point-spread functions of the deep and superficial source
-stc_psf = {"deep": {}, "super": {}}
-
-# for a deep source
-# compute resolution matrix with added noise for columns (PSFs)
-# added noise matches the SNR used for regularization
-rm_matched = make_inverse_resolution_matrix(
-    forward, inverse_operator, method=method, lambda2=lambda2, noise_cov=noise_cov
-)
-stc_psf["deep"]["matched"] = get_point_spread(
-    rm_matched, forward["src"], [sources["deep"]], norm=True
-)
-del rm_matched
-
-# compute resolution matrix with added noise for very high SNR
-rm_highsnr = make_inverse_resolution_matrix(
-    forward,
-    inverse_operator,
-    method=method,
-    lambda2=lambda2,
-    noise_cov=noise_cov,
-    snr=1000.0,
-)
-stc_psf["deep"]["highsnr"] = get_point_spread(
-    rm_highsnr, forward["src"], [sources["deep"]], norm=True
-)
-del rm_highsnr
-
-# the same for a superficial source
-# compute matching resolution matrix with added noise for columns (PSFs)
-# added noise matches the SNR used for regularization
-rm_matched = make_inverse_resolution_matrix(
-    forward, inverse_operator, method=method, lambda2=lambda2, noise_cov=noise_cov
-)
-stc_psf["super"]["matched"] = get_point_spread(
-    rm_matched, forward["src"], [sources["super"]], norm=True
-)
-del rm_matched
-
-# compute resolution matrix with added noise for very high SNR
-rm_highsnr = make_inverse_resolution_matrix(
-    forward,
-    inverse_operator,
-    method=method,
-    lambda2=lambda2,
-    noise_cov=noise_cov,
-    snr=1000.0,
-)
-stc_psf["super"]["highsnr"] = get_point_spread(
-    rm_highsnr, forward["src"], [sources["super"]], norm=True
-)
-del rm_highsnr
-
-##############################################################################
-# Visualize
-# ---------
-# PSF:
-
-vertno_lh = forward["src"][0]["vertno"]
-for source in sources:
-    # Which vertex corresponds to selected source
-    verttrue = [vertno_lh[sources[source]]]  # just one vertex
-    for psf in ["matched", "highsnr"]:
-        brain_psf = stc_psf[source][psf].plot(
-            "sample", "inflated", "lh", subjects_dir=subjects_dir
-        )
-        brain_psf.show_view("ventral")
-        title_str = f"{method} {source} {psf}"
-        brain_psf.add_text(
-            0.1, 0.9, title_str, "title", font_size=16, color=(255, 255, 255)
-        )
-        # mark true source location
-        brain_psf.add_foci(
-            verttrue, coords_as_verts=True, scale_factor=1.0, hemi="lh", color="green"
-        )
-
-# %%
-# The green spheres indicate the true source location.
-# For high SNR, PSFs for neither location are ideal. For the superficial
-# it is at least centred around the true source location, but still widespread
-# and with sidelobes. For the deep location the peaks occur at large distances
-# to the true location. This shows that spatial resolution is not an issue of
-# insufficient SNR but is inherently limited.
-# For matched (more realistic) SNR the PSF resembles the noise-less PSF, albeit
-# with more widely spread-out sidelobes. This indicates that moderate
-# regularization still achieves good spatial resolution while effectively
-# suppressing noise.
+"""
+.. _ex-psf-empirical:
+
+===================================================
+Plot point-spread functions (PSFs) with added noise
+===================================================
+
+Visualise the point-spread of a vertex in a sLORETA source estimate when we assume the
+data to be noisy.
+"""
+# Authors: Olaf Hauk <olaf.hauk@mrc-cbu.cam.ac.uk>
+#          Alexandre Gramfort <alexandre.gramfort@inria.fr>
+#
+# License: BSD-3-Clause
+# Copyright the MNE-Python contributors.
+
+import mne
+from mne.datasets import sample
+from mne.minimum_norm import get_point_spread, make_inverse_resolution_matrix
+
+print(__doc__)
+
+data_path = sample.data_path()
+subjects_dir = data_path / "subjects"
+meg_path = data_path / "MEG" / "sample"
+fname_fwd = meg_path / "sample_audvis-meg-eeg-oct-6-fwd.fif"
+fname_cov = meg_path / "sample_audvis-cov.fif"
+fname_evo = meg_path / "sample_audvis-ave.fif"
+
+# read forward solution
+forward = mne.read_forward_solution(fname_fwd)
+# forward operator with fixed source orientations
+mne.convert_forward_solution(forward, surf_ori=True, force_fixed=True, copy=False)
+
+# noise covariance matrix
+noise_cov = mne.read_cov(fname_cov)
+
+# evoked data for info
+info = mne.io.read_info(fname_evo)
+
+# make inverse operator from forward solution
+# free source orientation
+inverse_operator = mne.minimum_norm.make_inverse_operator(
+    info=info, forward=forward, noise_cov=noise_cov, loose=0.0, depth=None
+)
+
+# regularisation parameter
+snr = 3.0
+lambda2 = 1.0 / snr**2
+method = "MNE"
+
+# %%
+# We compute noisy PSFs for two locations, one deep and one superficial.
+# For each location we compute two PSFs::
+# 1) For a high SNR, approximating the noise-less case.
+# 2) With the SNR of regularization and the added noise matched (realistic).
+
+# vertices of a deep and a superficial source
+sources = {"deep": 2146, "super": 2209}
+
+# point-spread functions of the deep and superficial source
+stc_psf = {"deep": {}, "super": {}}
+
+# for a deep source
+# compute resolution matrix with added noise for columns (PSFs)
+# added noise matches the SNR used for regularization
+rm_matched = make_inverse_resolution_matrix(
+    forward, inverse_operator, method=method, lambda2=lambda2, noise_cov=noise_cov
+)
+stc_psf["deep"]["matched"] = get_point_spread(
+    rm_matched, forward["src"], [sources["deep"]], norm=True
+)
+del rm_matched
+
+# compute resolution matrix with added noise for very high SNR
+rm_highsnr = make_inverse_resolution_matrix(
+    forward,
+    inverse_operator,
+    method=method,
+    lambda2=lambda2,
+    noise_cov=noise_cov,
+    snr=1000.0,
+)
+stc_psf["deep"]["highsnr"] = get_point_spread(
+    rm_highsnr, forward["src"], [sources["deep"]], norm=True
+)
+del rm_highsnr
+
+# the same for a superficial source
+# compute matching resolution matrix with added noise for columns (PSFs)
+# added noise matches the SNR used for regularization
+rm_matched = make_inverse_resolution_matrix(
+    forward, inverse_operator, method=method, lambda2=lambda2, noise_cov=noise_cov
+)
+stc_psf["super"]["matched"] = get_point_spread(
+    rm_matched, forward["src"], [sources["super"]], norm=True
+)
+del rm_matched
+
+# compute resolution matrix with added noise for very high SNR
+rm_highsnr = make_inverse_resolution_matrix(
+    forward,
+    inverse_operator,
+    method=method,
+    lambda2=lambda2,
+    noise_cov=noise_cov,
+    snr=1000.0,
+)
+stc_psf["super"]["highsnr"] = get_point_spread(
+    rm_highsnr, forward["src"], [sources["super"]], norm=True
+)
+del rm_highsnr
+
+##############################################################################
+# Visualize
+# ---------
+# PSF:
+
+vertno_lh = forward["src"][0]["vertno"]
+for source in sources:
+    # Which vertex corresponds to selected source
+    verttrue = [vertno_lh[sources[source]]]  # just one vertex
+    for psf in ["matched", "highsnr"]:
+        brain_psf = stc_psf[source][psf].plot(
+            "sample", "inflated", "lh", subjects_dir=subjects_dir
+        )
+        brain_psf.show_view("ventral")
+        title_str = f"{method} {source} {psf}"
+        brain_psf.add_text(
+            0.1, 0.9, title_str, "title", font_size=16, color=(255, 255, 255)
+        )
+        # mark true source location
+        brain_psf.add_foci(
+            verttrue, coords_as_verts=True, scale_factor=1.0, hemi="lh", color="green"
+        )
+
+# %%
+# The green spheres indicate the true source location.
+# For high SNR, PSFs for neither location are ideal. For the superficial
+# it is at least centred around the true source location, but still widespread
+# and with sidelobes. For the deep location the peaks occur at large distances
+# to the true location. This shows that spatial resolution is not an issue of
+# insufficient SNR but is inherently limited.
+# For matched (more realistic) SNR the PSF resembles the noise-less PSF, albeit
+# with more widely spread-out sidelobes. This indicates that moderate
+# regularization still achieves good spatial resolution while effectively
+# suppressing noise.

From 18c49dcd07e1207fca62f7a463f9c6287692640c Mon Sep 17 00:00:00 2001
From: Marijn van Vliet <w.m.vanvliet@gmail.com>
Date: Thu, 29 Jan 2026 22:49:34 +0200
Subject: [PATCH 27/27] little cleanup

---
 mne/minimum_norm/resolution_matrix.py | 6 +-----
 1 file changed, 1 insertion(+), 5 deletions(-)

diff --git a/mne/minimum_norm/resolution_matrix.py b/mne/minimum_norm/resolution_matrix.py
index 44c361408e6..6c8b2869ee2 100755
--- a/mne/minimum_norm/resolution_matrix.py
+++ b/mne/minimum_norm/resolution_matrix.py
@@ -15,11 +15,7 @@
 from ..evoked import EvokedArray
 from ..forward.forward import Forward, convert_forward_solution
 from ..label import Label
-from ..source_estimate import (
-    _get_src_type,
-    _make_stc,
-    _prepare_label_extraction,
-)
+from ..source_estimate import _get_src_type, _make_stc, _prepare_label_extraction
 from ..source_space._source_space import SourceSpaces, _get_vertno
 from ..utils import _validate_type, logger, verbose
 from .inverse import apply_inverse, apply_inverse_cov