[ENH] improve MEM performance

best/cmem/clusters/be_scores2alpha.m

@@ -48,15 +48,6 @@
 ALPHA_METHOD = OPTIONS.model.alpha_method;
 alpha_threshold = OPTIONS.model.alpha_threshold;
 
-
-% Progress bar
-hmem = [];
-if numel(varargin)
-    hmem    = varargin{1}(1);
-    st      = varargin{1}(2);
-    dr      = varargin{1}(3);
-end
-
 for jj=1:size(SCR,2)
     clusters = CLS(:,jj);
     scores = SCR(:,jj);
@@ -94,11 +85,6 @@
             curr_cls = curr_cls + 1;
         end
 
-        % update progress bar
-         if hmem
-             prg = round( (st + dr * (jj - 1 + ii/nb_clusters) / (size(SCR,2))) * 100 );
-             waitbar(prg/100, hmem, {[OPTIONS.automatic.InverseMethod, 'Step 1/2 : Running cortex parcellization ... ' num2str(prg) ' % done']});
-         end
     end
 
 end

best/cmem/clusters/be_spatial_priorw.m

@@ -1,6 +1,6 @@
 function [OPTIONS, W] = be_spatial_priorw(OPTIONS, neighbors)
-%   This function returns the W Green matrix from which the local 
-%   covariance matrices will be obtained. 
+%   This function returns the W spatial filter matrix.
+%   The local covariance matrices can be obtained as C = W' * W
 %
 %   INUPTS:
 % 		- OPTIONS    	: structure of parameters
@@ -41,10 +41,10 @@
     OPTIONS = Def_OPTIONS;
     return
 end
+
 % Check field names of passed OPTIONS and fill missing ones with default values
 OPTIONS = be_struct_copy_fields(OPTIONS, Def_OPTIONS, {'solver'}, 0);
-clear Def_OPTIONS
-%%
+
 % Parameters:
 nb_vertices = size(OPTIONS.automatic.Modality(1).gain,2);
 rho = OPTIONS.solver.spatial_smoothing; % scalar that weight the adjacency matrix 
@@ -66,6 +66,5 @@
     A0 = rho/2*A0*A / (i+1);
 end
 W = W.*(W > exp(-8));
-W = W'*W; 
 
-return
+end

best/cmem/preprocess/be_unormalize_and_units.m

@@ -0,0 +1,34 @@
+function [obj, OPTIONS] = be_unormalize_and_units(obj, OPTIONS)
+% be_unormalize_and_units: Un-normalize the result source map.
+% Inverse function of be_normalize_and_units
+% -------------------------------------------------------------------------
+%   Author: LATIS 2012
+%
+% ==============================================
+% License
+%
+% BEst is free software: you can redistribute it and/or modify
+%    it under the terms of the GNU General Public License as published by
+%    the Free Software Foundation, either version 3 of the License, or
+%    (at your option) any later version.
+%
+%    BEst is distributed in the hope that it will be useful,
+%    but WITHOUT ANY WARRANTY; without even the implied warranty of
+%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+%    GNU General Public License for more details.
+%
+%    You should have received a copy of the GNU General Public License
+%    along with BEst. If not, see <http://www.gnu.org/licenses/>.
+% -------------------------------------------------------------------------
+
+
+if  strcmp(OPTIONS.optional.normalization,'adaptive')
+    obj.ImageGridAmp = obj.ImageGridAmp/OPTIONS.automatic.Modality(1,1).ratioAmp;
+else
+    obj.ImageGridAmp = obj.ImageGridAmp*OPTIONS.automatic.Modality(1).units_dipoles; %Modified by JSB August 17th 2015
+end
+
+OPTIONS.automatic.Modality(1).Jmne = OPTIONS.automatic.Modality(1).Jmne * OPTIONS.automatic.Modality(1).MNEAmp;
+
+
+end

best/cmem/solver/be_cmem_pipelineoptions.m

@@ -27,8 +27,9 @@
 
         % solver
         DEF.solver.NoiseCov_method           = 2;
-
+        DEF.solver.mne_use_noiseCov          = 0;
+
         % optional
         DEF.optional.normalization           = 'adaptive'; 
 
-return
+end

best/cmem/solver/be_cmem_solver.m

@@ -98,7 +98,9 @@
 if OPTIONS.optional.verbose
     fprintf('\n\n===== pipeline cMEM\n');
 end 
- obj = struct('hfig', [] , 'hfigtab', [] );
+
+obj = struct();
+[obj.hfig, obj.hfigtab] = be_create_figure(OPTIONS);
 
 %% Retrieve vertex connectivity - needed for clustering
 [OPTIONS, obj.VertConn] = be_vertex_connectivity(HeadModel, OPTIONS);
@@ -154,8 +156,6 @@
 % and the leadfields
 [OPTIONS] = be_normalize_and_units(OPTIONS);
 
-
-
 %% ===== Double precision to single  ===== %%
 % relax the double precision for the msp (leadfield and data)
 [OPTIONS] = be_switch_precision( OPTIONS, 'single' );
@@ -184,6 +184,12 @@
 %% ===== Solve the MEM ===== %%
 
 [obj.ImageGridAmp, OPTIONS] = be_launch_mem(obj, OPTIONS);
+if OPTIONS.optional.display
+    be_display_entropy_drops(obj,OPTIONS);
+end
+
+%% ===== Un-Normalization  ===== %%
+[obj, OPTIONS] = be_unormalize_and_units(obj, OPTIONS);
 
 
 %% ===== Inverse temporal data window  ===== %%
@@ -199,12 +205,13 @@
                                 'minimum_norm',OPTIONS.automatic.Modality(1).Jmne, ...
                                 'MSP',obj.SCR);     
 
-% Results (full temporal sequence)                  
+% Results                   
 Results = struct(...
-    'ImageGridAmp',     obj.ImageGridAmp, ...
-    'ImagingKernel',    [], ...
-    'MEMoptions',       OPTIONS); % ...
-     %'MEMdata',          obj);
+    'ImageGridAmp',    [], ...
+    'ImagingKernel',   [], ...
+    'nComponents',     round( length(obj.iModS) / obj.nb_sources ), ...
+    'MEMoptions',      OPTIONS); 
+Results.ImageGridAmp =  obj.ImageGridAmp;
 
 disp('Bye.')
 end

best/main/be_main_call.m

@@ -258,6 +258,7 @@
 % MEM solver
 Def_OPTIONS.solver.NoiseCov                     = [];
 Def_OPTIONS.solver.NoiseCov_method              = [];
+Def_OPTIONS.solver.mne_use_noiseCov             = 0;
 Def_OPTIONS.solver.NoiseCov_recompute           = 1;
 Def_OPTIONS.solver.spatial_smoothing            = 0.6;
 Def_OPTIONS.solver.active_var_mult              = 0.05;

best/main/be_main_clustering.m

@@ -53,15 +53,15 @@
     % rmem
     switch OPTIONS.mandatory.pipeline
         case 'cMEM'
-            [CLS, SCR, OPTIONS] = be_cmem_clusterize_multim(obj, OPTIONS); %TO DO: CHECK IF DATA SHOULD BE MINUS BSL OR NOT
+            [CLS, SCR, OPTIONS] = be_cmem_clusterize_multim(obj, OPTIONS); 
             [ALPHA, CLS, OPTIONS]   = be_scores2alpha(SCR, CLS, OPTIONS);
         case 'wMEM'
             [CLS, SCR, OPTIONS] = be_wmem_clusterize_multim(obj, OPTIONS);
 
-            if OPTIONS.model.alpha_method ~= 6
+            if OPTIONS.model.alpha_method < 6
                 [ALPHA, CLS, OPTIONS] = be_scores2alpha(SCR, CLS, OPTIONS);
             else % We compute the score using MNE
-                [ALPHA, CLS, OPTIONS] = be_gain2alpha(obj , CLS, OPTIONS);
+                [ALPHA, CLS, OPTIONS] = be_mne2alpha(obj , CLS, OPTIONS);
             end
 
         case 'rMEM'

best/main/be_main_data_preprocessing.m

@@ -43,7 +43,7 @@
 % ===== call of specific pre-processing depending on the pipeline:
 switch OPTIONS.mandatory.pipeline
 
-    case 'cMEM'  % === this is the standard MEM (Chowdhurry and Grova) =====
+    case {'cMEM','cMNE'}  % === this is the standard MEM (Chowdhurry and Grova) =====
         [noise_var] = covariance_processing(OPTIONS);
         % we nomalize the cov of each modalities (we regularize the cov matrices)
         for ii = 1 : numel( OPTIONS.automatic.Modality )

best/mem/be_estimate_active_state.m

@@ -70,12 +70,15 @@
     xi_trans = xi';                              
 
     % F0 is set to a dirac by default (omega=0).
-    if isempty(omega)
-        F0=0;
-    else
+    F0 = 0;
+    if ~isempty(omega)
         F0 = 1/2 * xi_trans * omega * xi;             
     end
-    F1 = 1/2 * xi_trans * sigma * xi + xi_trans * mu;
+
+    F1 = 1/2 * xi_trans * sigma * xi;
+    if ~isempty(mu)
+        F1 = F1 + xi_trans * mu;
+    end
 
     F = F0 - F1;
 
@@ -92,10 +95,19 @@
 
 
     % Estimating amplitudes*
-    estimated_amp{ii} = estimated_alpha * mu + ...
-        ((1 - estimated_alpha) * omega * xi) + ...
-        estimated_alpha * sigma * xi;
+    estimated_amp{ii} = estimated_alpha * sigma * xi;
+
+    if ~isempty(mu)
+        estimated_amp{ii} = estimated_amp{ii}  + ...
+                             estimated_alpha * mu;
+    end
+
+    if ~isempty(omega)
+        estimated_amp{ii}  = estimated_amp{ii}  + ...
+                             ((1 - estimated_alpha) * omega * xi);
+    end
+
 end
 
-return
+end