%% TASK GED

%% Time-delay covariances

clearvars -except EEG

% Signal
EEG.data = double(EEG.data);
nbchans  = size(EEG.data,1);
trials   = size(EEG.data,3);
pnts     = size(EEG.data,2);

%channel sorting
chanlocs                = struct2table(EEG.chanlocs);
[chanlocs_sorted,chidx] = sortrows(chanlocs, 'X','descend');
electrodes              = chanlocs_sorted.labels;
[EEG.chanlocs,chanlocs] = deal(table2struct(chanlocs_sorted));
EEG.data = EEG.data (chidx,:,:);

Rcov_time = [-300 -100 ]; % R-covariance time 
[~,Ridx(1)] = min(abs(EEG.times-Rcov_time(1)));
[~,Ridx(2)] = min(abs(EEG.times-Rcov_time(2)));

Scov_time = [ 0 300]; % S-covariance time 
[~,Sidx(1)] = min(abs(EEG.times-Scov_time(1)));
[~,Sidx(2)] = min(abs(EEG.times-Scov_time(2)));


%% ==== GED ====


cov_R = zeros(nbchans);
for triali=1:trials
    % COVARIANCE
end 
   cov_R = cov_R ./ EEG.trials;
    
   % apply shrinkage to R 
   shr = .01;
   cov_R = (1-shr)*cov_R + shr*mean(eig(cov_R))*eye(size(cov_R));  
 
cov_S = zeros(nbchans);
for triali=1:trials
    % COVARIANCE
end
cov_S = cov_S ./ EEG.trials;

% GED
[evecs,evals] = 
[evals,sidx]  = sort( diag(evals),'descend' );
evecs = evecs(:,sidx);

% normalize eigenvectors to unit length
for v = 1:size(evecs,2)
    evecs(:,v) = evecs(:,v)/norm(evecs(:,v)); 
end

nb_comps = 10;
comp_map = zeros(nb_comps,64); 
comp_ts  = zeros(nb_comps,EEG.pnts,EEG.trials);

for compi = 1:nb_comps       
    % compute map and adjust sign
    comp_map(compi,:) = ;
    [~,idx]           = max( abs(comp_map(compi,:)) ); % if the largest is negative, then multiply by -1;
    comp_map(compi,:) = comp_map(compi,:) * sign(comp_map(compi,idx));
           
    % compute time series
    tmp_tmp = ;
    comp_ts(compi,:,:) = reshape( tmp_tmp,[pnts trials] );
end 


%%

% Topographical maps:
figure(4),clf
set(gcf, 'color','w', 'WindowState', 'maximize');
for compi = 1:nb_comps

    subplot(2,5,compi)
    topoplot(,chanlocs,'plotrad',0.537,'maplimits', 'absmax','headrad', 'rim','electrodes', 'on', 'style', 'both','whitebk','on');
    title(['Component | number:', num2str(compi)], FontSize=20)
    colormap jet; axis tight


end

% Time-series of each source:
figure(5),clf
set(gcf, 'color','w', 'WindowState', 'maximize');
for compi = 1:nb_comps

    subplot(2,5,compi)
    plot(EEG.times,,'linewidth', 3,'color','k','LineStyle','-')
    box off
    title(['Component | number:', num2str(compi)], FontSize=20)

end