% Before running the pipeline, EEGLAB must be launched (by typing 'eeglab' in the command window)

% Clear EEGLAB
STUDY=[];CURRENTSTUDY=0;ALLEEG=[];EEG=[];CURRENTSET=[];
LASTCOM='STUDY=[];CURRENTSTUDY=0;ALLEEG=[];EEG=[];CURRENTSET=[];'; ALLCOM=[ALLCOM LASTCOM];

%% To change
% PAT=30; % marker S1-02
% where='lennox';
% filename='170228a-a_0021.edf';
% remove_bool=0; % whether some electrodes should be rejected
% remove_elec={}; % which electrode(s) to reject
% freq_oi='delta'; % frequency band of the ictal activity

PAT=14;
where='lennox';
filename='150722--d_0011.edf';
remove_bool=1; % whether some electrodes should be rejected
remove_elec={'Cz','Fz'}; % which electrode(s) to reject
freq_oi='theta-alpha'; % frequency band of the ictal activity

filter_HP=2; % high-pass filter cutoff frequency
pre = 10; % how many seconds before the marker
post = 20; % how many seconds after the marker
ica_type='runica'; % algorithm

name_pipeline='pipeline_final';

% Threshold for the component selection
thres_disp=600;
thres_freq_oi=0.2;
thres_var=0.01;

scalpmaps=0;

%% STEP 1 : Data acquisition and preprocessing

filepath=sprintf('/Volumes/LACIE_SHARE/EEG-Epilog/patients_deltamed-%s/patient_%d/PAT_%d_edf/%s',where,PAT,PAT,filename);
EEG = pop_biosig(filepath);
EEG.setname=filename;

% Change to double precision
EEG.data=double(EEG.data);

% Change the name of the electrodes (from 'EEG Fp1' to 'Fp1')
for i=1:length(EEG.chanlocs)
    EEG.chanlocs(i).labels = erase(EEG.chanlocs(i).labels,'EEG '); 
end

% Select only the 19 standard electrodes (10-20 system)
EEG=pop_select(EEG,'channel',{'Fp1','Fp2','F3','F4','C3','C4','P3','P4','O1','O2','F7','F8','T3','T4','T5','T6','Fz','Cz','Pz'});

% Add locations of the electrodes based on the standard BEM model
locfile='/Users/deborman/Documents/MATLAB/MasterThesis/eeglab2020_0/plugins/dipfit/standard_BEM/elec/standard_1005.elc';
EEG = pop_chanedit( EEG, 'lookup', locfile);

% Filter the data
filter_LP=30; % low-pass filter cutoff frequency
fprintf('Filtering data ... \n')
Fs=EEG.srate;
[b,a]=butter(2,[filter_HP filter_LP]/(Fs/2),'bandpass');
[f,e]=butter(2,[49 51]/(Fs/2),'stop');
for i=1:EEG.nbchan
    EEG.data(i,:)=filtfilt(b,a,double(EEG.data(i,:)));
    EEG.data(i,:)=filtfilt(f,e,double(EEG.data(i,:)));  
end

% Find the EPILOG marker
time=[];
markers=[];
for i=1:length(EEG.urevent)
    if strncmpi(EEG.urevent(i).type,'EPILOG',6)
        time=[time; EEG.urevent(i).latency/256];
        markers=[markers; EEG.urevent(i).type];
    end
end

% Ask the user which marker to choose
n_markers=length(time);
fprintf('%d EPILOG markers found \n',n_markers)
for i=1:n_markers
    fprintf('%s \n',markers(i,:))
end
if n_markers==1
    marker_chosen=markers(1,:);
    fprintf('First one chosen\n')
else
    choice=input('@@@ Which marker to choose? ');
    marker_chosen=markers(choice,:);
end
marker_chosen2=strrep(marker_chosen,'.','-');

if n_markers==1
    fprintf('The marker is at %4.2f \n',time(1))
else
    fprintf('The marker is at %4.2f \n',time(choice))
end
   
% Ask if some electrodes should be removed
if remove_bool==1
    EEG=pop_select(EEG,'nochannel',remove_elec);
end
N=EEG.nbchan;

% Re-reference the data
EEG.nbchan = EEG.nbchan+1;
EEG.data(end+1,:) = zeros(1, EEG.pnts);
EEG.chanlocs(1,EEG.nbchan).labels = 'initialReference';
EEG = pop_reref(EEG, []);
EEG = pop_select( EEG,'nochannel',{'initialReference'});

% Select data one minute before and after the seizure marker
EEG = pop_rmdat(EEG,{marker_chosen},[-60 60],0);

% Plot the data
eegplot(EEG.data,'srate',EEG.srate,'eloc_file',EEG.chanlocs,'winlength',15,'events',EEG.event,'title',strcat('PAT',num2str(PAT),'_',marker_chosen,' - Channels'))

%% STEP 2 : Time window selection

EEG=pop_rmdat(EEG,{marker_chosen},[-pre post],0);

% Name of the ICA file
if remove_bool==1
    no_elec=strcat(remove_elec{1:end});
    icaset=sprintf('ica_-%d+%d_no%s_%dHz_%s',pre,post,no_elec,filter_HP,freq_oi);
    icaset=strrep(icaset,'.','-');
    name_ica=sprintf('PAT%d_%s_ica_-%d+%d_no%s_%dHz_%s',PAT,marker_chosen2,pre,post,no_elec,filter_HP,freq_oi);
    name_ica=strrep(name_ica,'.','-');
else
    icaset=sprintf('ica_-%d+%d_%dHz_%s',pre,post,filter_HP,freq_oi);
    icaset=strrep(icaset,'.','-');
    name_ica=sprintf('PAT%d_%s_ica_-%d+%d_%dHz_%s',PAT,marker_chosen2,pre,post,filter_HP,freq_oi);
    name_ica=strrep(name_ica,'.','-');
end

%% STEP 3 : ICA decomposition

if ica_type=='runica'
    EEG = pop_runica(EEG,'icatype', 'runica','extended',1,'rndreset','no');
elseif ica_type=='picard'
    EEG = pop_runica(EEG,'icatype', 'picard','maxiter',200);
else
    fprintf('Unknown ICA algorithm');
end
EEG.icaact = (EEG.icaweights*EEG.icasphere)*EEG.data(EEG.icachansind,:);

%% eLORETA and performance measure parameters

load('mesh_head2')
load('sa_v2'); % head model parameters from the ICBM-NY

epsilon=0.1;
gamma = 0.01; % regulation parameter for eLORETA

% Conversion of the electrodes
elec=zeros(1,EEG.nbchan);
elec108=sa_v2.EEG_clab_electrodes;
for i=1:EEG.nbchan
    found=false;
    j=1;
    current=EEG.chanlocs(i).labels;
    while found==false && j<=length(elec108)
        if strcmp(elec108(j),current)==1 || ...
                strcmp(elec108(j),'T7')==1 && strcmp(current,'T3')==1 ||...
                strcmp(elec108(j),'T8')==1 && strcmp(current,'T4')==1 ||...
                strcmp(elec108(j),'P7')==1 && strcmp(current,'T5')==1 ||...
                strcmp(elec108(j),'P8')==1 && strcmp(current,'T6')==1

            elec(i)=j;
            found=true;
        end
        j=j+1;
    end
end

nK_sources=5;
factor_sources=round(74/nK_sources);

locs_nK=downsample(sa_v2.cortex75K.vc,factor_sources); % downsample the source space
fprintf('The source space contains %d potential sources. \n', size(locs_nK,1))
cortexnK_in_from_cortex75K=1:factor_sources:74382;

LF=sa_v2.cortex75K.EEG_V_fem_normal(:,cortexnK_in_from_cortex75K);  % lead field matrix
LF=LF(elec,:); 
P = mkfilt_eloreta2(LF,gamma); % size(P) = N_elec x N_sources

% Performance measure

tol=10; % 10 mm tolerance
results.tol=tol;
coord_voi_path=sprintf('/Volumes/LACIE_SHARE/MRI_MT/PAT%d_MRI/MAT_voi/PAT%d_coord_voi.mat',PAT,PAT);
coord_voi_down_path=sprintf('/Volumes/LACIE_SHARE/MRI_MT/PAT%d_MRI/MAT_voi/PAT%d_coord_voi_down.mat',PAT,PAT);
load(coord_voi_path);
load(coord_voi_down_path);
all_ind=1:size(locs_nK);

%% STEP 4 : Components selection

% Frequency analysis
if strcmp(freq_oi,'alpha')
    freq_range=[8,12];
elseif strcmp(freq_oi,'beta')
    freq_range=[12,20];
elseif strcmp(freq_oi,'theta')
    freq_range=[4,8];
elseif strcmp(freq_oi,'delta')
    freq_range=[1,4];
elseif strcmp(freq_oi,'alpha-beta')
    freq_range=[10,15];
elseif strcmp(freq_oi,'theta-alpha')
    freq_range=[6,10];
elseif strcmp(freq_oi,'delta-theta')
    freq_range=[2,6];
else
    fprintf(2,'WARNING: Unknown frequency\n');
end

% 1 ..(pre seconds).. onset ..(pre_seconds).. end_f
% --------BEFORE------- | ----------AFTER----------

onset=Fs*pre;
end_f=2*Fs*pre;
before=1:onset;
after=onset:end_f;
full=1:end_f;

power_freq_oi=zeros(N,1);
power_freq_phys=zeros(N,1);
var_expl=zeros(N,1);
dispersion=zeros(N,1);

ind_accepted=[];
ind_rejected=[];

ind_rej_disp=[];
ind_rej_var=[];
ind_rej_freq=[];

for j=1:N
    
    icaact_j=EEG.icaact(j,:);
    icawinv_j=EEG.icawinv(:,j);
    
    % Power in the frequency band of interest
    % (during the pre seconds after the marker)
    l=length(icaact_j(after));
    power_freq_oi(j) = bandpower(icaact_j(after),Fs,freq_range)/bandpower(icaact_j(after),Fs,[0 (l-1)*Fs/(2*l)-1]);
    
    % Power in the frequency band [2 20] Hz
    % (during the pre seconds after the marker)
    power_freq_phys(j) = bandpower(icaact_j(after),Fs,[2 20])/bandpower(icaact_j(after),Fs,[0 (l-1)*Fs/(2*l)-1]);

    % Variance explained
    back_proj=icawinv_j*icaact_j(after);
    var_expl(j) = 1-mean(var((EEG.data(:,after) - back_proj)'))/mean(var(EEG.data(:,after)'));
    
    % Active cloud
    
    sources=P'*EEG.icawinv(:,j);
    sources=sources.^2;
    sources=sources/max(sources);

    all_ind=1:size(locs_nK);
    ind_act=all_ind(sources>epsilon);
    [~,ind_most]=max(sources);

    locs_act=locs_nK(ind_act,:);

    threshold=0.8;
    ind_nmost=all_ind(sources>threshold);
    n_most=length(ind_nmost);

    locs_nmost=locs_nK(ind_nmost,:);
    C=cov(locs_nmost);
    dispersion(j)=norm(C,'fro');
    
    if dispersion(j)<thres_disp && power_freq_oi(j)>thres_freq_oi && var_expl(j)>thres_var
        ind_accepted=[ind_accepted; j];
    else
        ind_rejected=[ind_rejected; j];
        
        if power_freq_oi(j)<thres_freq_oi
            ind_rej_freq=[ind_rej_freq; j];
        end
        if dispersion(j)>thres_disp
            ind_rej_disp=[ind_rej_disp; j];
        end
        
    end
    
end

[~,ind_order_acc]=sort(power_freq_oi(ind_accepted),'descend');
[~,ind_order_rej]=sort(power_freq_oi(ind_rejected),'descend');

comp_order=[ind_accepted(ind_order_acc); ind_rejected(ind_order_rej)];

% code segment from pop_runica()
EEG.icaweights = EEG.icaweights(comp_order,:);
EEG.icawinv    = pinv( EEG.icaweights *  EEG.icasphere );
if ~isempty(EEG.icaact)
    EEG.icaact = EEG.icaact(comp_order,:,:);
end

colors = cell(1,19);
colors(:)={[0.5 0.7 1]}; % light blue
%colors(:)={[0.4    0.7    0.4]}; % light green

for k=ind_rej_disp'
    index=find(comp_order==k);
    colors(index)={[0.5 0.7 1]}; % light blue
    %colors(index)={[0.6 0.6 1]}; % light purple
end

for k=ind_rej_freq'
    index=find(comp_order==k);
    colors(index)={[0.5 0.7 1]}; % light blue
end

colors(1:length(ind_accepted)) = {[0 0 0.4]};

eegplot(EEG.icaact,'srate',EEG.srate,'winlength',15,'events',EEG.event,'title',strcat(name_ica,' - Components'),'color',colors)

%% STEP 5 : Scalp maps

if scalpmaps==1
    pop_topoplot(EEG,0,[1:EEG.nbchan]);
end

%% Ask the user which component to choose

ictal_comp=input('@@@ Select ictal components : ');
while(ictal_comp(1)>0)

    % STEP 6 : Source imaging
    
    if length(ictal_comp)==1
     
        % Project the ictal component on the scalp and compute the source activity
        sources=P'*EEG.icawinv(:,ictal_comp);
        sources=sources.^2;
        
    else
        % Project the ictal component on the scalp
        EEG_ictal = pop_subcomp(EEG,ictal_comp,0,1);
        data=EEG_ictal.data';

        % Inverse problem: brain sources reconstruction
        sources_t=data*P;

        [pxx,freq]=periodogram(sources_t,[],[],Fs,'power');
        sources=mean(pxx,1); % mean over all the frequencies
    end
    
    sources=sources/max(sources);
    
    all_ind=1:size(locs_nK);
    ind_act=all_ind(sources>epsilon);
    [~,ind_most]=max(sources);

    locs_act=locs_nK(ind_act,:);
    loc_most=locs_nK(ind_most,:);

    % Plot eLORETA
    figure('Color',[1,1,1])

    oldcmap = colormap('bone');
    newcmp=colormap( flipud(oldcmap) );
    rgb_1 = vals2colormap(sources(ind_act),newcmp,[0 1]);

    scatter3(locs_act(:,1),locs_act(:,2),locs_act(:,3),30,rgb_1,'filled'); hold on %alpha 0.9
    scatter3(loc_most(:,1),loc_most(:,2),loc_most(:,3),30,'red','filled'); hold on

    ft_plot_mesh(mesh_h,'facealpha',.1,'edgealpha',.1)
    ft_plot_mesh(locs_nK,'facealpha',.1,'vertexsize',2,'vertexcolor',[1 1 1]);
    view([0 90])
    colormap(newcmp)
    colorbar

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Performance measure

    % the most active dipole
    [~,ind_most]=max(sources);
    loc_most=locs_nK(ind_most,:);
    [nearest_ind_most,dist_most]=dsearchn(coord_voi,loc_most);

    % n most active dipoles
    threshold=0.8;
    ind_nmost=all_ind(sources>threshold);
    n_most=length(ind_nmost);

    locs_nmost=locs_nK(ind_nmost,:);
    C=cov(locs_nmost);
    norm(C,'fro');

    % Average position of the n_most active dipoles
    loc_avg_nmost=[mean(locs_nmost(:,1)) mean(locs_nmost(:,2)) mean(locs_nmost(:,3))];

    % Distance from the average position of the n_most active dipoles to the lesion
    [nearest_ind_avg_nmost,dist_avg_nmost]=dsearchn(coord_voi,loc_avg_nmost);

    % Distance from the n_most active dipoles to the lesion
    outside=0;
    for i=1:n_most
        [nearest_ind_nmost,dist_nmost]=dsearchn(coord_voi,locs_nmost(i,:));
        if dist_nmost>tol
            outside=outside+1;
        end
    end

    % Percentage of the n_most active dipoles in the lesion
    percentage_nmost=100-outside/n_most*100;
    
    % Plot results
    figure('Color',[1,1,1])

    oldcmap = colormap('bone');
    newcmp=colormap( flipud(oldcmap) );
    rgb_1 = vals2colormap(sources(ind_act),newcmp,[0 1]);

    scatter3(locs_nmost(:,1),locs_nmost(:,2),locs_nmost(:,3),30,'green','filled'); hold on
    scatter3(loc_most(:,1),loc_most(:,2),loc_most(:,3),30,'red','filled'); hold on
    scatter3(loc_avg_nmost(:,1),loc_avg_nmost(:,2),loc_avg_nmost(:,3),30,'blue','filled'); hold on
    scatter3(coord_voi_down(:,1),coord_voi_down(:,2),coord_voi_down(:,3),30,[1 0.55 0],'filled','MarkerEdgeAlpha', 0.15,'MarkerFaceAlpha',0.15); hold on
    % Plot the electrodes
    %x_list=[EEG.chanlocs.X];
    %y_list=[EEG.chanlocs.Y];
    %z_list=[EEG.chanlocs.Z];
    %scatter3(y_list,x_list,z_list,30,[0.5 0.55 0.5],'filled','s','MarkerEdgeAlpha', 0.3,'MarkerFaceAlpha',0.3); hold on
    %scatter3(y_list,x_list,z_list,30,[0.5 0.55 0.5],'s','MarkerEdgeAlpha', 0.5,'MarkerFaceAlpha',0.5); hold on

    ft_plot_mesh(mesh_h,'facealpha',.1,'edgealpha',.1); hold on
    view([0 90])

    if length(ictal_comp)>1
        str1=sprintf('%g, ', ictal_comp(1:end-1));
        str2=sprintf('%g', ictal_comp(end));
        name_ica_comp=sprintf('Projection of components: %s%s',str1,str2);
    else
        name_ica_comp=sprintf('Projection of component: %d',ictal_comp);
    end

    text_max='\color{red}max\color{black}';
    text_avg='\color{blue}avg\color{black}';
    text_nmost1='\color{green}';
    text_nmost2=sprintf('%d max',n_most);
    text_nmost3='\color{black}';
    title_fig=sprintf('%s\nDistance to resection zone (%s) : %.2f mm\nDistance to resection zone (%s) : %.2f mm\nPercentage inside (%s%s%s) : %.2f %%\n\n\n\n',name_ica_comp,text_max,dist_most,text_avg,dist_avg_nmost,text_nmost1,text_nmost2,text_nmost3,percentage_nmost);
    text(-100,120,100,title_fig,'fontsize',14,'FontWeight','Normal')
    %title(title_fig,'fontsize',12,'FontWeight','Normal')


    colormap(newcmp)

    fprintf('Distance to lesion (max) : %.2f mm\nDistance to lesion (avg) : %.2f mm\nPercentage inside the lesion : %.2f %%\n',dist_most,dist_avg_nmost,percentage_nmost)

    ictal_comp=input('@@@ Select ictal components : ');
end

%% SAVE

EEG.setname=name_ica;

eeglab redraw