[TEST] Playing around with fixed dipole orientation in a discrete source space

[wmvanvliet]

Trying to see what it takes to perform MxNE using a source space that consists of two dipoles with fixed orientation. Currently, MNE doesn't support forward models with fixed orientation that are not of type 'surface'. In this PR I'm disabling these checks and see what happens.

Here is some example code that creates such a forward model and attempts to perform MxNE source estimation with it:

import mne
import numpy as np
import matplotlib.pyplot as plt

# Read phantom data
data_path = mne.datasets.brainstorm.bst_phantom_elekta.data_path(verbose=True)
raw = mne.io.read_raw_fif(f'{data_path}/kojak_all_200nAm_pp_no_chpi_no_ms_raw.fif')
events = mne.find_events(raw, 'STI201')

# Extract evoked data for 2 dipoles
epochs = mne.Epochs(raw, events, [1, 17], tmin=-0.1, tmax=0.1, baseline=(-0.1, -0.05))
evoked_dip1 = epochs['1'].average()
evoked_dip2 = epochs['17'].average()

# Create an evoked with both dipole active at the same time
evoked = evoked_dip1.copy()
evoked.data += evoked_dip2.data

# Fit two dipoles
cov = mne.compute_covariance(epochs, tmax=-0.05)
bem = mne.make_sphere_model(r0=[0, 0, 0], head_radius=0.08)
dip1, _ = mne.fit_dipole(evoked_dip1.copy().crop(0.035, 0.035), cov, bem)
dip2, _ = mne.fit_dipole(evoked_dip2.copy().crop(0.035, 0.035), cov, bem)

# Collect the dipoles in a single `Dipole` object
dips = mne.Dipole(
    times=np.hstack((dip1.times, dip2.times)),
    pos=np.vstack((dip1.pos, dip2.pos)),
    ori=np.vstack((dip1.ori, dip2.ori)),
    amplitude=np.hstack((dip1.amplitude, dip2.amplitude)),
    gof=np.hstack((dip1.gof, dip2.gof)),
    name=None,
    conf=None,
    khi2=np.hstack((dip1.khi2, dip2.khi2)),
    nfree=dip1.nfree,
)

# Compute dipole timecourses using MxNE
fwd, _ = mne.make_forward_dipole(dips, bem, evoked.info)
stc = mne.inverse_sparse.mixed_norm(evoked, fwd, cov)

# Plot the timecourses
plt.plot(evoked.times, stc.data.T)

[larsoner]

Currently, MNE doesn't support forward models with fixed orientation that are not of type 'surface'.

This works for me:

>>> import mne
>>> info = mne.io.read_info(mne.datasets.sample.data_path() / 'MEG' / 'sample' / 'sample_audvis_raw.fif')
>>> src = mne.setup_volume_source_space(pos=dict(rr=[[0, 0, 0.055]], nn=[[1, 0, 0]]))
>>> bem = mne.make_sphere_model(info=info)
>>> fwd = mne.make_forward_solution(info, src=src, bem=bem, trans=None)
>>> fwd_fixed = mne.convert_forward_solution(fwd, force_fixed=True)
>>> fwd_fixed['sol']['data'].shape
(366, 1)
>>> fwd['sol']['data'].shape
(366, 3)

But indeed with depth weighting, loose etc. maybe we do need to change our checks to be in ('surface', 'discrete') rather than == 'surface'. I think it's okay to shift to fixed ori for surf and discrete, but not volume.

[wmvanvliet]

Ah, I needed to clarify better: constructing a volume or discrete forward model with fixed dipole orientations is no problem at all. It's when you attempt to use it to perform source modeling. the forward._prepare_forward function has a series of checks.

[wmvanvliet]

I think it's okay to shift to fixed ori for surf and discrete, but not volume.

What makes volume special in this regard?

[larsoner]

What makes volume special in this regard?

When doing force_fixed=True:

• For a surface src, things are fine (it means use the surface normal).
• For a discrete src, I think we can assume we're in an expert use case and trust they set their normals properly (and suffer the consequences if they did not)
• For volume src, 99%+ of the time it's going to be a mistake. The volume normals in setup_volume_source_space are just +Z. If you really want this, you should make a discrete source space with these normals.

[agramfort]

note that with so few dipoles mne.inverse_sparse.mixed_norm should be used with alpha almost set to 0.

[wmvanvliet]

In that case, would plain MNE also work here? What is the added benefit of the sparse solver?

[wmvanvliet]

(@larsoner recommended MxNE in the forum thread)

[agramfort]

you will get a set of dipoles instances with MxNE code

…

Message ID: ***@***.***>

[larsoner]

With explained variance per dipole (though you could figure it out for MNE the same way I think). Otherwise, assuming you use a very small lambda2 for MNE, the total exp var should be similar. Would make a nice little unit test!

[wmvanvliet]

Would make a nice little unit test!

A test that passes:

import mne
import numpy as np
import matplotlib.pyplot as plt

# Read phantom data
data_path = mne.datasets.brainstorm.bst_phantom_elekta.data_path(verbose=True)
raw = mne.io.read_raw_fif(f'{data_path}/kojak_all_200nAm_pp_no_chpi_no_ms_raw.fif')
events = mne.find_events(raw, 'STI201')

# Extract evoked data for 2 dipoles
epochs = mne.Epochs(raw, events, [1, 17], tmin=-0.1, tmax=0.1, baseline=(-0.1, -0.05))
evoked_dip1 = epochs['1'].average()
evoked_dip2 = epochs['17'].average()

# Create an evoked with both dipole active at the same time
evoked = mne.combine_evoked([evoked_dip1, evoked_dip2], [1, 1])

# Fit two dipoles
cov = mne.compute_covariance(epochs, tmax=-0.05)
bem = mne.make_sphere_model(r0=[0, 0, 0], head_radius=0.08)
dip1, _ = mne.fit_dipole(evoked_dip1.copy().crop(0.035, 0.035), cov, bem)
dip2, _ = mne.fit_dipole(evoked_dip2.copy().crop(0.035, 0.035), cov, bem)

# Collect the dipoles in a single `Dipole` object
dips = mne.Dipole(
    times=np.hstack((dip1.times, dip2.times)),
    pos=np.vstack((dip1.pos, dip2.pos)),
    ori=np.vstack((dip1.ori, dip2.ori)),
    amplitude=np.hstack((dip1.amplitude, dip2.amplitude)),
    gof=np.hstack((dip1.gof, dip2.gof)),
    name=None,
    conf=None,
    khi2=np.hstack((dip1.khi2, dip2.khi2)),
    nfree=dip1.nfree,
)

# Compute dipole timecourses using both MxNE and plain MNE
fwd, _ = mne.make_forward_dipole(dips, bem, evoked.info)
inv = mne.minimum_norm.make_inverse_operator(evoked.info, fwd, cov, fixed=True, depth=0)
stc_mxne = mne.inverse_sparse.mixed_norm(evoked, fwd, cov)
stc_mne = mne.minimum_norm.apply_inverse(evoked, inv, method='MNE')

# Plot the timecourses
plt.plot(evoked.times, stc_mxne.data.T)
plt.plot(evoked.times, stc_mne.data.T)

[agramfort]

pretty reassuring ! thx @wmvanvliet for the code snippet

[wmvanvliet]

Closing in favor of #10464