Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

add the possibility to compute the MEG (forward FEM) with a block of sensors #367

Closed
wants to merge 7 commits into from
Closed

add the possibility to compute the MEG (forward FEM) with a block of sensors #367

Show file tree
Hide file tree
Changes from 1 commit
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
0e2bcbe
add the possibility to compute the MEG (with FEM) with a block of sen…
tmedani Dec 1, 2020
c5db49c
add the meg sensor divider as suggested by ftadel
tmedani Dec 7, 2020
2740fd1
add the meg sensor divider to the panel as suggested
tmedani Dec 7, 2020
bf70ee8
adapt the bst_duneuro to the new options
tmedani Dec 7, 2020
dff7613
bug fix: avoid multiple computation of the EEG with the combined opti…
tmedani Dec 7, 2020
80d24e5
add clock for fem computation per block
tmedani Dec 7, 2020
60992c8
some weired thing on the java tool
tmedani Feb 23, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
48 changes: 44 additions & 4 deletions toolbox/forward/bst_duneuro.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -523,7 +523,47 @@
disp(['DUNEURO> System call: ' callStr]);
tic;
% Call DUNEuro
[status,cmdout] = system(callStr);
cfg.SizeOfBlockOfSensor = 250; % we may ask the user to change this value from the panel? Q for Francois
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

250?
For a 256-electrode EEG channel, it means that it would run once with 250 and once with 6 sensors?
For an Elekta MEG: one execution with 250 and one with 56?
This sounds like a weird logic.

Have you tried to quantify the increase in computation time vs. the decrease in memory usage?
For example: when you divide the number of sensors by 2, by 4, by 8.
If you need to split this list, I guess the most optimal way to balance the computation time and the memory usage would be to split in equal blocks, isn't it?

I am not in favor of adding useless parameters that no one will ever modify. Maybe one check box that allows to divide in 2 or 4 could be useful for re-running when it crashes because of memory errors?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

250?
For a 256-electrode EEG channel, it means that it would run once with 250 and once with 6 sensors?
For an Elekta MEG: one execution with 250 and one with 56?
This sounds like a weird logic.

as discussed, this is not for EEG but for MEG sensors where four integrations points are used for magnetometers and 8 for gradiometers, so the total number of points for computation is around 8x250 ~= 2000

Have you tried to quantify the increase in computation time vs. the decrease in memory usage?
For example: when you divide the number of sensors by 2, by 4, by 8.

Yes, in some basic case, less block faster is the solution, with more block the process of reading/writing is multiplied.

Maybe one check box that allows to divide in 2 or 4 could be useful for re-running when it crashes because of memory errors?

ok, this is another useful option, then I can modify the code to this option

if isMeg && cfg.MegPerBlockOfSensor
% Define the group of channles
megNbOfBlock = 1: cfg.SizeOfBlockOfSensor : length(MegChannels);
groupOfSensors = cell( length(megNbOfBlock) ,1);
for iBlock = 1 : length(megNbOfBlock)
if ~(iBlock == length(megNbOfBlock))
groupOfSensors{iBlock} = megNbOfBlock(iBlock) : megNbOfBlock(iBlock+1) - 1;
else
groupOfSensors{iBlock} = megNbOfBlock(iBlock) : length(MegChannels);
end
end
% Update the sensor file
GainMeg = [];
for iBlock =1 : length(megNbOfBlock)
disp(['block ' num2str(iBlock) '/' num2str(length(megNbOfBlock)) ': sensors ' ])
% Update the channels file
% Write coil file
CoilsLocTemp = MegChannels(groupOfSensors{iBlock},2:4);
fid = fopen(CoilFile, 'wt+');
fprintf(fid, '%d %d %d \n', CoilsLocTemp');
fclose(fid);
% Write projection file
CoilsOrientTemp = MegChannels(groupOfSensors{iBlock},5:7);
fid = fopen(ProjFile, 'wt+');
fprintf(fid, '%d %d %d \n', CoilsOrientTemp');
fclose(fid);
% call the DUNEuro
[status,cmdout] = system(callStr);
GainTmp = in_duneuro_bin(fullfile(TmpDir, cfg.BstMegLfFile))';
GainMeg = [GainMeg; GainTmp];
% TODO cfg.BstSaveTransfer % not possible with this version
% solution : concatenate and saveback the transfer matrix?, or
% disable this option in this case, ask Francois how to do it from the GUI?
% when the option "use meg block sensor" is set to 1, disable the save of the
% transfer matrix.? Q for Francois
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In panel_duneuro/UpdatePanel, add something like:

if ~isempty(jCheckMegPerBlockOfSensor)
    jCheckMegPerBlockOfSensor.setEnabled(~jCheckSaveTransfer.isSelected());
    if jCheckSaveTransfer.isSelected()
        jCheckMegPerBlockOfSensor.setSelected(0);
    end
end

You could also rename jCheckMegPerBlockOfSensor into jCheckMegBlock.
These super long names are painful to read.

And remove this comment :)

end
else
[status,cmdout] = system(callStr);
end

if (status ~= 0)
disp('DUNEURO> Error log:');
disp(cmdout);
Expand All @@ -532,7 +572,6 @@
end
disp(['DUNEURO> FEM computation completed in: ' num2str(toc) 's']);


%% ===== READ LEADFIELD ======
bst_progress('text', 'DUNEuro: Reading leadfield...');
% EEG
Expand All @@ -542,8 +581,9 @@

%MEG
if isMeg
GainMeg = in_duneuro_bin(fullfile(TmpDir, cfg.BstMegLfFile))';

if ~cfg.MegPerBlockOfSensor
GainMeg = in_duneuro_bin(fullfile(TmpDir, cfg.BstMegLfFile))';
end
% === POST-PROCESS MEG LEADFIELD ===
% Compute the total magnetic field
dipoles_pos_orie = [kron(cfg.GridLoc,ones(3,1)), kron(ones(length(cfg.GridLoc),1), eye(3))];
Expand Down
21 changes: 14 additions & 7 deletions toolbox/forward/panel_duneuro.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -234,7 +234,7 @@
% high memory, or want to use the integration points
jCheckEnableCacheMemory = gui_component('checkbox', jPanelMegComputationOption, 'br', 'Enable cache memory', [], '', [], []);
% Enable the MEG Computation per block of sensors
... jCheckMegPerBlockOfSensor = gui_component('checkbox', jPanelMegComputationOption, 'br', 'Compute per block of sensors [Todo]', [], '', [], []);
jCheckMegPerBlockOfSensor = gui_component('checkbox', jPanelMegComputationOption, 'br', 'Compute per block of sensor ', [], '', [], []);
% Set jCheckUseIntegrationPoint to 1 as default option
if (OPTIONS.UseIntegrationPoint)
jCheckUseIntegrationPoint.setSelected(1);
Expand All @@ -254,6 +254,13 @@
if (OPTIONS.BstSaveTransfer)
jCheckSaveTransfer.setSelected(1);
end
% disable the save transfer option if the option per block is
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Delete this section, move it to UpdatePanel()

% activated, ask Francois for help to be interactive from the gui
% disable this option only if the MegPerBlockOfSensor is checked
if ~isempty(jCheckMegPerBlockOfSensor)
isMegPerBlockOfSensor = jCheckMegPerBlockOfSensor.isSelected();
jCheckSaveTransfer.setEnabled(isMegPerBlockOfSensor);
end
c.gridy = 5;
jPanelRight.add(jPanelOutput, c);

Expand Down Expand Up @@ -313,7 +320,7 @@
'jCheckSaveTransfer', jCheckSaveTransfer, ...
'jCheckUseIntegrationPoint', jCheckUseIntegrationPoint,...
'jCheckEnableCacheMemory', jCheckEnableCacheMemory,...
...'jCheckMegPerBlockOfSensor', jCheckMegPerBlockOfSensor,...
'jCheckMegPerBlockOfSensor', jCheckMegPerBlockOfSensor,...
'UseTensor', OPTIONS.UseTensor);
ctrl.FemNames = OPTIONS.FemNames;
% Create the BstPanel object that is returned by the function
Expand Down Expand Up @@ -469,11 +476,11 @@ function UpdatePanel(isForced)
s.EnableCacheMemory = 0;
end

% if ~isempty(ctrl.jCheckMegPerBlockOfSensor)
% s.MegPerBlockOfSensor = ctrl.jCheckMegPerBlockOfSensor.isSelected();
% else
% s.MegPerBlockOfSensor = 0;
% end
if ~isempty(ctrl.jCheckMegPerBlockOfSensor)
s.MegPerBlockOfSensor = ctrl.jCheckMegPerBlockOfSensor.isSelected();
else
s.MegPerBlockOfSensor = 0;
end
end


Expand Down