-
Notifications
You must be signed in to change notification settings - Fork 2
Expand file tree
/
Copy pathsolver.py
More file actions
118 lines (95 loc) · 3.81 KB
/
Copy pathsolver.py
File metadata and controls
118 lines (95 loc) · 3.81 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
import numpy as np
import mne
from .get_options import GetHeadModel, GetMEMOptions
def _solver(eng, head_model, mem_options):
"""Run L2 penalized regression and keep 10 strongest locations.
Parameters
----------
head_model : dict
The head model, as returned by GetHeadModel.
mem_options : dict
The options for the MEM solver, as returned by GetMEMOptions.
Returns
-------
X : array, (n_active_dipoles, n_times)
The time series of the dipoles in the active set.
active_set : array (n_dipoles)
Array of bool. Entry j is True if dipole j is in the active set.
We have ``X_full[active_set] == X`` where X_full is the full X matrix
such that ``M = G X_full``.
"""
# # Transfer M and G to MATLAB workspace
eng.workspace['HeadModel'] = head_model
eng.workspace['MEMOptions'] = mem_options
# Prepare and call the MATLAB function using the variables in the workspace
(Results, OPTIONS) = eng.eval("be_main_call(HeadModel, MEMOptions)", nargout=2)
ImageGridAmp_np = np.array(Results["ImageGridAmp"])
# TODO : Extract the active set, do not create a np.ones array
return ImageGridAmp_np, np.ones(len(head_model["vertex_connectivity"][1]), dtype=bool)
def _mem_solver(eng, evoked, forward, noise_cov, PyMemOptions=None):
"""Call a custom solver on evoked data.
This function does all the necessary computation:
- to select the channels in the forward given the available ones in
the data
- to take into account the noise covariance and do the spatial whitening
- to apply loose orientation constraint as MNE solvers
- to apply a weigthing of the columns of the forward operator as in the
weighted Minimum Norm formulation in order to limit the problem
of depth bias.
Parameters
----------
eng: matlab.engine
The MATLAB engine.
evoked : instance of mne.Evoked
The evoked data
forward : instance of Forward
The forward solution.
noise_cov : instance of Covariance
The noise covariance.
head_model : dict
The head model, as returned by GetHeadModel.
memOptions : dict
The options for the MEM solver, as returned by GetMEMOptions.
Returns
-------
stc : instance of SourceEstimate
The source estimates.
"""
if PyMemOptions is None:
raise ValueError("MEM Options can't be None")
undefined = PyMemOptions.get_defined_mandatory()
if len(undefined) != 0:
raise ValueError("Some errors occured in the the MEM Options : " + ", ".join(undefined))
# Import the necessary private functions
from mne.inverse_sparse.mxne_inverse import (
_make_sparse_stc,
_prepare_gain,
_reapply_source_weighting,
is_fixed_orient,
)
all_ch_names = evoked.ch_names
# Handle depth weighting and whitening (here is no weights)
forward, gain, gain_info, whitener, source_weighting, mask = _prepare_gain(
forward,
evoked.info,
noise_cov,
pca=False,
depth=None,
loose=0,
weights=None,
weights_min=None,
rank=None,
)
# Select channels of interest
sel = [all_ch_names.index(name) for name in gain_info["ch_names"]]
M = evoked.data[sel]
####### SOLVER
vertex_connectivity_matrix = mne.spatial_src_adjacency(forward['src'])
default_cmem_pipeline_options = eng.be_cmem_pipelineoptions()
head_model = GetHeadModel(gain, vertex_connectivity_matrix, PyMemOptions)
mem_options = GetMEMOptions(default_cmem_pipeline_options, M, evoked.times, noise_cov, PyMemOptions)
X, active_set = _solver(eng, head_model, mem_options)
stc = _make_sparse_stc(
X, active_set, forward, tmin=evoked.times[0], tstep=1.0 / evoked.info["sfreq"]
)
return stc