MRG: Fix the computation of the orientation in (..)MxNE and gamma-map #4633

yousrabk · 2017-10-02T16:33:36Z

No description provided.

yousrabk · 2017-10-02T16:35:21Z

mne/inverse_sparse/_gamma_map.py

            for ii in range(len(active_src)):
                for jj in range(3):
                    if in_pos >= len(active_set):
                        break
                    if (active_set[in_pos] + jj) % 3 == 0:
                        X_xyz[3 * ii + jj] = X[in_pos]
+                        X_ori_[3 * ii + jj] = X_ori[in_pos]


XXX: There is still a problem here, I have to check what is happening when we give xyz_same_gamma=False

codecov-io · 2017-10-02T17:28:25Z

Codecov Report

Merging #4633 into master will decrease coverage by <.01%.
The diff coverage is 96.55%.

@@            Coverage Diff             @@
##           master    #4633      +/-   ##
==========================================
- Coverage   87.69%   87.68%   -0.01%     
==========================================
  Files         351      351              
  Lines       66554    66635      +81     
  Branches    10343    11554    +1211     
==========================================
+ Hits        58362    58429      +67     
- Misses       5233     5256      +23     
+ Partials     2959     2950       -9

agramfort · 2017-10-02T19:29:57Z

mne/inverse_sparse/_gamma_map.py

@@ -277,8 +277,10 @@ def gamma_map(evoked, forward, noise_cov, alpha, loose="auto", depth=0.8,

    # Reapply weights to have correct unit
    n_dip_per_pos = 1 if is_fixed_orient(forward) else 3
+    X_ori = X.copy()
    X = _reapply_source_weighting(X, source_weighting,
                                  active_set, n_dip_per_pos)


n_dip_per_pos is not used in _reapply_source_weighting :(

agramfort · 2017-10-02T19:32:40Z

mne/inverse_sparse/_gamma_map.py

    X = _reapply_source_weighting(X, source_weighting,
                                  active_set, n_dip_per_pos)
+    X_ori /= source_weighting[active_set][:, None]


this looks bogus. _reapply_source_weighting does `X *= source_weighting[active_set][:, None]``
so I doubt it's doing the right thing but it's just a guess

agramfort · 2017-10-02T19:34:23Z

mne/inverse_sparse/_gamma_map.py


    tmin = evoked.times[0]
    tstep = 1.0 / evoked.info['sfreq']

    if return_as_dipoles:
        out = _make_dipoles_sparse(X, active_set, forward, tmin, tstep, M,
-                                   M_estimated, active_is_idx=True)
+                                   M_estimated, active_is_idx=True,
+                                   X_ori=X_ori)


it's really a mystery to me why you need to pass both X and X_ori and why internally you use X to get amplitude and X_ori to get orientation.

agramfort · 2017-10-02T19:35:09Z

mne/inverse_sparse/tests/test_gamma_map.py

@@ -124,5 +124,20 @@ def test_gamma_map_vol_sphere():

    assert_array_almost_equal(stc.times, evoked.times, 5)

+    # Compare orientation obtained using fit_dipole and mixed_norm


agramfort · 2017-10-02T19:35:44Z

mne/inverse_sparse/tests/test_gamma_map.py

+    pos2 = src[0]['rr'][stc.vertices[1]]
+    evoked_dip = mne.simulation.simulate_evoked(fwd, stc, info, cov, nave=1e9)
+
+    dip_mxne = gamma_map(evoked_dip, fwd, cov, alpha, return_as_dipoles=True)


agramfort · 2017-10-08T09:50:18Z

mne/inverse_sparse/mxne_inverse.py

@@ -24,8 +24,7 @@
 @verbose
 def _prepare_weights(forward, gain, source_weighting, weights, weights_min):
    mask = None
-    if isinstance(weights, SourceEstimate):
-        # weights = np.sqrt(np.sum(weights.data ** 2, axis=1))
+    if isinstance(weights, (SourceEstimate, VolSourceEstimate)):


if isinstance(weights, _BaseSourceEstimate):

agramfort · 2017-10-08T14:25:30Z

mne/inverse_sparse/tests/test_mxne_inverse.py

+
+    amp_max = [np.max(d.amplitude) for d in dip_mxne]
+    dip_mxne = dip_mxne[np.argmax(amp_max)]
+    assert_true(dip_mxne[0].pos[0] in pos)


this does not pass on master but the assert below on ori does.

so this test is not enough to demonstrate the bug fix on orientation
estimation.

agramfort · 2017-10-08T14:33:56Z

tutorials/plot_brainstorm_phantom_elekta.py

+
+
+def plot_pos_ori(pos, ori, color=(0., 0., 0.)):
+    mlab.points3d(pos[0], pos[1], pos[2], scale_factor=0.005,


this works if you pass pos of shape (3, n_times)

as is the example is broken

larsoner · 2017-10-09T18:42:07Z

tutorials/plot_brainstorm_phantom_elekta.py

@@ -137,3 +138,19 @@

 fig.tight_layout()
 plt.show()
+
+


Why not plot_dipole_locations?

https://mne-tools.github.io/dev/generated/mne.viz.plot_dipole_locations.html#mne.viz.plot_dipole_locations

Ahh okay.

In any case, please make this a separate section with ########### so it renders nicer in SG

agramfort · 2017-10-09T18:43:42Z

I don't think it works without an MRI...

yousrabk · 2017-10-09T19:15:32Z

Are we good here after Eric's comment?

larsoner · 2017-10-09T19:31:48Z

Good enough for me

https://6100-1301584-gh.circle-artifacts.com/0/tmp/circle-artifacts.2wP5EEF/html/auto_tutorials/plot_brainstorm_phantom_elekta.html

agramfort · 2017-10-09T20:57:52Z

let me just check tomorrow that the added tests don't work in master so it's actually obvious as a bug fix cc @joewalter

agramfort · 2017-10-10T08:10:41Z

ok I did I pass

@yousrabk have a look

FYI @joewalter we realized that the depth weighting strategy for sparse solver was not optimal for orientation estimation and we now use the Frobenius normalization you had written in the past.

larsoner

Otherwise LGTM

larsoner · 2017-10-10T13:52:32Z

examples/inverse/plot_mixed_norm_inverse.py

@@ -53,7 +53,7 @@

 ###############################################################################
 # Run solver
-alpha = 50  # regularization parameter between 0 and 100 (100 is high)
+alpha = 55  # regularization parameter between 0 and 100 (100 is high)


If you made this change because the example now otherwise looks bad, we'll have to make it clear in the whats_new (please update) that the computations have changed and that changing alpha might be necessary to partially compensate.

larsoner · 2017-10-10T13:56:40Z

mne/inverse_sparse/tests/test_gamma_map.py

+
+    dip_fit = mne.fit_dipole(evoked_dip, cov, sphere)[0]
+    assert_true(np.max(np.abs(np.dot(dip_fit.ori[0],
+                                     dip_gmap.ori.T))) > 0.99)


Why is there a max here? Aren't these single dipoles? If so you can just do np.abs(np.dot(dip_fit.ori[0], dip_gmap.ori[0])).

Something else just occurred to me -- if we guarantee our dipoles have positive amplitude values during fitting, then the abs isn't even needed, because they should both point in the "outward" direction. But we can fix this bit later if you don't want to investigate (but please add an XXX comment so we don't forget)

agramfort · 2017-10-11T13:33:06Z

mne/inverse_sparse/tests/test_gamma_map.py

+    assert_true(dip_gmap[0].pos[0] in src[0]['rr'][stc.vertices])
+
+    dip_fit = mne.fit_dipole(evoked_dip, cov, sphere)[0]
+    assert_true(np.abs(np.dot(dip_fit.ori[0], dip_gmap.ori[0])) > 0.99)  # XXX


@yousrabk why the XXX?

It was in response to my request to add a # XXX comment to signify that we should look into whether the abs is even necessary at some point (if currents are always positive it shouldn't be), but without at least a few words saying this, it is tough to know that

agramfort · 2017-10-11T13:38:12Z

yes it is necessary since there is no reason why the sign of the amplitude of the dipoles should be fixed. please remove

larsoner · 2017-10-11T13:43:02Z

yes it is necessary since there is no reason why the sign of the amplitude of the dipoles should be fixed.

My point is that we could (for fit_dipole) choose to always select the orientation that gives a positive current. Then we wouldn't need it.

larsoner · 2017-10-11T13:44:36Z

(and for gamma_map IIUC it uses the source space definition with the dipoles pointing out of cortex, which is the direction of the current, which should be positive)

agramfort · 2017-10-11T14:05:05Z

ok but gamma map can be on a volume / discrete source space so... anyhow if you keep the XXX please add a comment. I would just remove the XXX

larsoner · 2017-10-11T14:17:16Z

Ahh that's true at least some of these are volume. Sure let's just remove (maybe do it in master after merge)

agramfort · 2017-10-12T07:30:27Z

ok XXX removed and what's new updated.

good to go on my end if CIs are happy.

larsoner · 2017-10-12T13:25:02Z

Thanks @yousrabk @agramfort

yousrabk added 4 commits October 2, 2017 14:03

fix ori computation

d12add4

fix ori for gamma map

2128c94

add 3D plot (dipole pos/ori) to the tutorial

856c90e

add tests

dbee697

yousrabk commented Oct 2, 2017

View reviewed changes

yousrabk changed the title ~~Fix the computation of the orientation in (..)MxNE and gamma-map~~ WIP: Fix the computation of the orientation in (..)MxNE and gamma-map Oct 2, 2017

agramfort reviewed Oct 2, 2017

View reviewed changes

larsoner added this to the 0.15 milestone Oct 3, 2017

Fix ori computation: change method to compute depth weighting

98e640e

agramfort reviewed Oct 8, 2017

View reviewed changes

fix comment

3996ac2

agramfort reviewed Oct 8, 2017

View reviewed changes

yousrabk added 2 commits October 8, 2017 19:41

fix tutorial example

2a8a81b

tests

ec0803a

larsoner reviewed Oct 9, 2017

View reviewed changes

add section to tutorial when plotting

e62824c

agramfort added 4 commits October 10, 2017 09:32

remove code dupes

41f20ea

typo

944f80d

remove warings in gamma map tests

7fdb0d9

remove warings in mxne tests

df739d6

agramfort changed the title ~~WIP: Fix the computation of the orientation in (..)MxNE and gamma-map~~ MRG: Fix the computation of the orientation in (..)MxNE and gamma-map Oct 10, 2017

fix reg parameters to have the same images as before

cadca0c

larsoner mentioned this pull request Oct 10, 2017

ENH: plot_dipole_locations should be able to use a sphere #4655

Open

larsoner reviewed Oct 10, 2017

View reviewed changes

check only 1st ori as they are the same

98f6dbc

agramfort reviewed Oct 11, 2017

View reviewed changes

agramfort added 2 commits October 12, 2017 09:25

remove XXX

5d9df30

update what's new

93df08c

larsoner merged commit 35eb5be into mne-tools:master Oct 12, 2017

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

MRG: Fix the computation of the orientation in (..)MxNE and gamma-map #4633

MRG: Fix the computation of the orientation in (..)MxNE and gamma-map #4633

yousrabk commented Oct 2, 2017

yousrabk Oct 2, 2017

codecov-io commented Oct 2, 2017 •

edited

Loading

agramfort Oct 2, 2017

agramfort Oct 2, 2017

agramfort Oct 2, 2017

agramfort Oct 2, 2017

agramfort Oct 2, 2017

agramfort Oct 8, 2017

agramfort Oct 8, 2017

agramfort Oct 8, 2017

larsoner Oct 9, 2017

larsoner Oct 9, 2017

yousrabk Oct 9, 2017

agramfort commented Oct 9, 2017 via email

yousrabk commented Oct 9, 2017

larsoner commented Oct 9, 2017

agramfort commented Oct 9, 2017 via email

agramfort commented Oct 10, 2017

larsoner left a comment

larsoner Oct 10, 2017

larsoner Oct 10, 2017

agramfort Oct 11, 2017

larsoner Oct 11, 2017

agramfort commented Oct 11, 2017 via email

larsoner commented Oct 11, 2017

larsoner commented Oct 11, 2017

agramfort commented Oct 11, 2017 via email

larsoner commented Oct 11, 2017

agramfort commented Oct 12, 2017

larsoner commented Oct 12, 2017

		@@ -124,5 +124,20 @@ def test_gamma_map_vol_sphere():

		assert_array_almost_equal(stc.times, evoked.times, 5)

		# Compare orientation obtained using fit_dipole and mixed_norm



		def plot_pos_ori(pos, ori, color=(0., 0., 0.)):
		mlab.points3d(pos[0], pos[1], pos[2], scale_factor=0.005,

MRG: Fix the computation of the orientation in (..)MxNE and gamma-map #4633

MRG: Fix the computation of the orientation in (..)MxNE and gamma-map #4633

Conversation

yousrabk commented Oct 2, 2017

Choose a reason for hiding this comment

codecov-io commented Oct 2, 2017 • edited Loading

Codecov Report

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

agramfort commented Oct 9, 2017 via email

yousrabk commented Oct 9, 2017

larsoner commented Oct 9, 2017

agramfort commented Oct 9, 2017 via email

agramfort commented Oct 10, 2017

larsoner left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

agramfort commented Oct 11, 2017 via email

larsoner commented Oct 11, 2017

larsoner commented Oct 11, 2017

agramfort commented Oct 11, 2017 via email

larsoner commented Oct 11, 2017

agramfort commented Oct 12, 2017

larsoner commented Oct 12, 2017

codecov-io commented Oct 2, 2017 •

edited

Loading