%%
clear all; clc

Rootpath='Y:\EEG_Data\PDDys\EEG\ODDBALL\';

datalocation=[Rootpath,'Processed Data\'];   % Data are here
cd(datalocation);

% Matrix for medication randomization (ON=on meds);
load([Rootpath,'ONOFF_4Oddball.mat'],'ONOFF_4Oddball','ONOFF_4Oddball_HDR')

PDsx=[804 805 806 807 808 809 810 811 813 814 815 816 817 818 819 820 821 822 823 824 825 826 827 828 829];
CTLsx=[890 891 892 893 895 896 897 898 899 900 901 902 903 904 905 906 907 909 910 911 912 913 914 8060 8070];


for subj=1:length(PDsx)
    for sessioni=1:3
        if sessioni==1 % PD ON
            temp=ONOFF_4Oddball(find(ONOFF_4Oddball(:,1)==PDsx(subj)),:);
            ONsession=temp(temp(:,3)==1,2);
            clc; disp([num2str(PDsx(subj)),'_',num2str(ONsession),'_PDDys_ODDBALL.mat']);  load([num2str(PDsx(subj)),'_',num2str(ONsession),'_PDDys_ODDBALL.mat'],'EEG','bad_chans','bad_epochs','bad_ICAs');
        elseif sessioni==2 % PD OFF
            temp=ONOFF_4Oddball(find(ONOFF_4Oddball(:,1)==PDsx(subj)),:);
            OFFsession=temp(temp(:,3)==0,2);
            clc; disp([num2str(PDsx(subj)),'_',num2str(OFFsession),'_PDDys_ODDBALL.mat']);  load([num2str(PDsx(subj)),'_',num2str(OFFsession),'_PDDys_ODDBALL.mat'],'EEG','bad_chans','bad_epochs','bad_ICAs');
        elseif sessioni==3 % CTL ON
            clc; disp([num2str(CTLsx(subj)),'_1_PDDys_ODDBALL.mat']);  load([num2str(CTLsx(subj)),'_1_PDDys_ODDBALL.mat'],'EEG','bad_chans','bad_epochs','bad_ICAs');
        end
        clear temp ONsession OFFsession;

        % -------------
        % Remove the (presumptive) bad ICAs:
        bad_ICAs_To_Remove=bad_ICAs{2};
        EEG = pop_subcomp( EEG, bad_ICAs_To_Remove, 0);
        
        % Get the good info out of the epochs
        for ai=1:size(EEG.epoch,2)
            for bi=1:size(EEG.epoch(ai).eventlatency,2)
                % Get STIMTYPE
                if EEG.epoch(ai).eventlatency{bi}==0 % If this bi is the event
                    % Get StimType
                    FullName=EEG.epoch(ai).eventtype{bi};
                    EEG.epoch(ai).STIM=str2num(FullName(2:end)) ;
                    VECTOR(ai)=str2num(FullName(2:end));
                    clear FullName
                end
            end
        end
        
        dims=size(EEG.data);
        EEG.data=eegfilt(EEG.data,500,[],20);
        EEG.data=eegfiltfft(EEG.data,500,.1,[]);
        EEG.data=reshape(EEG.data,dims(1),dims(2),dims(3));
        
        % -------------- Downsample
        EEG=pop_resample(EEG,250);
        
        % Let's See
        dims=size(EEG.data);
        tx=-2000:4:1998;
        Nb1=find(tx==-200);  Nb2=find(tx==0);
        STb1=find(tx==-200+452);  STb2=find(tx==0+452);
        Nt1=find(tx==-248);  Nt2=find(tx==748); % for Novel
        STt1=find(tx==-248+452);  STt2=find(tx==748+452); % for Std & Target
        tx2disp=-248:4:748;
        
        % Basecor your ERPs here so they are pretty. ------- diff bases for diff event offsets
        N_BASE=squeeze(  mean(EEG.data(:,Nb1:Nb2,VECTOR==202),2)  );
        ST_BASE=squeeze(  mean(EEG.data(:,STb1:STb2,VECTOR<202),2)  );
        for ai=1:dims(1)
            EEG.data(ai,:,VECTOR==202)=squeeze(EEG.data(ai,:,VECTOR==202))-repmat( N_BASE(ai,:),dims(2),1 );
            EEG.data(ai,:,VECTOR<202)=squeeze(EEG.data(ai,:,VECTOR<202))-repmat( ST_BASE(ai,:),dims(2),1 );
        end
        
        % Select 30 random standards
        tempnum=find(VECTOR==201);
        tempnum=shuffle(tempnum);
        STD_VECTOR=tempnum(1:30);
        for ai=1:30
            if ~any(STD_VECTOR(ai)==tempnum); BOOM; end
        end
        clear tempnum
        
        % Save as {subj}{session} then repeat for all condis
        STANDARD{subj}{sessioni}=single(EEG.data(:,STt1:STt2,STD_VECTOR));
        TARGET{subj}{sessioni}=single(EEG.data(:,STt1:STt2,VECTOR==200));
        NOVELTY{subj}{sessioni}=single(EEG.data(:,Nt1:Nt2,VECTOR==202));
        
        clearvars -except Rootpath datalocation PDsx CTLsx subj sessioni STANDARD TARGET NOVELTY tx2disp ONOFF_4Oddball
    end
end

save([datalocation,'SINGLETRIAL_ODDBALL.mat'],'STANDARD','TARGET','NOVELTY','PDsx','CTLsx','tx2disp')

% Trial count match all condis by session - - then agg subjs
clear S T N;
for ai=1:3
    for bi=1:3
        S{bi}=[]; T{bi}=[]; N{bi}=[];
    end
end
for subj=1:length(PDsx)
    
    SIZE(1,1)=size(STANDARD{subj}{1},3); % ON
    SIZE(1,2)=size(STANDARD{subj}{2},3); % OFF
    SIZE(1,3)=size(STANDARD{subj}{3},3); % CTL
    
    SIZE(2,1)=size(TARGET{subj}{1},3);
    SIZE(2,2)=size(TARGET{subj}{2},3);
    SIZE(2,3)=size(TARGET{subj}{3},3);
    
    SIZE(3,1)=size(NOVELTY{subj}{1},3);
    SIZE(3,2)=size(NOVELTY{subj}{2},3);
    SIZE(3,3)=size(NOVELTY{subj}{3},3);
    
    ToTake(subj)=min(min(SIZE));
    
    for ai=1:3
        for bi=1:3
            if SIZE(ai,bi)>ToTake(subj)
                if ai==1
                    Exist=1:size(STANDARD{subj}{bi},3);  Exist=shuffle(Exist);
                    STANDARD{subj}{bi}=STANDARD{subj}{bi}(:,:,Exist(1:ToTake(subj))); clear Exist;
                    S{bi}=cat(3,S{bi},STANDARD{subj}{bi});
                elseif ai==2
                    Exist=1:size(TARGET{subj}{bi},3);  Exist=shuffle(Exist);
                    TARGET{subj}{bi}=TARGET{subj}{bi}(:,:,Exist(1:ToTake(subj))); clear Exist;
                    T{bi}=cat(3,T{bi},TARGET{subj}{bi});
                elseif ai==3
                    Exist=1:size(NOVELTY{subj}{bi},3);  Exist=shuffle(Exist);
                    NOVELTY{subj}{bi}=NOVELTY{subj}{bi}(:,:,Exist(1:ToTake(subj))); clear Exist;
                    N{bi}=cat(3,N{bi},NOVELTY{subj}{bi});
                end
            elseif SIZE(ai,bi)==ToTake(subj)
                if ai==1
                    S{bi}=cat(3,S{bi},STANDARD{subj}{bi});
                elseif ai==2
                    T{bi}=cat(3,T{bi},TARGET{subj}{bi});
                elseif ai==3
                    N{bi}=cat(3,N{bi},NOVELTY{subj}{bi});
                end
            end
        end
    end
    clear SIZE;
    
end

% Display AggERPs
site=36;      
figure;
subplot(3,1,1)
hold on
plot(tx2disp,squeeze(mean(S{1}(site,:,:),3)),'r','linewidth',2)
plot(tx2disp,squeeze(mean(S{2}(site,:,:),3)),'b','linewidth',2)
plot(tx2disp,squeeze(mean(S{3}(site,:,:),3)),'g','linewidth',2)
set(gca,'xlim',[-250 750],'xtick',[-250:250:750],'ylim',[-4 4]);
legend('ON','OFF','CTL')
title('Standard')
subplot(3,1,2)
hold on
plot(tx2disp,squeeze(mean(T{1}(site,:,:),3)),'r','linewidth',2)
plot(tx2disp,squeeze(mean(T{2}(site,:,:),3)),'b','linewidth',2)
plot(tx2disp,squeeze(mean(T{3}(site,:,:),3)),'g','linewidth',2)
set(gca,'xlim',[-250 750],'xtick',[-250:250:750],'ylim',[-4 4]);
title('Target')
subplot(3,1,3)
hold on
plot(tx2disp,squeeze(mean(N{1}(site,:,:),3)),'r','linewidth',2)
plot(tx2disp,squeeze(mean(N{2}(site,:,:),3)),'b','linewidth',2)
plot(tx2disp,squeeze(mean(N{3}(site,:,:),3)),'g','linewidth',2)
set(gca,'xlim',[-250 750],'xtick',[-250:250:750],'ylim',[-4 4]);
title('Novelty')


% re-save
save([datalocation,'SINGLETRIAL_ODDBALL_AGG.mat'],'S','T','N','PDsx','CTLsx','tx2disp','ToTake')

%%
load([Rootpath,'BV_Chanlocs_60.mat']);

site=36;      
figure; hold on
plot(tx2disp,squeeze(mean(S{2}(site,:,:),3)),'r:','linewidth',2)
plot(tx2disp,squeeze(mean(S{3}(site,:,:),3)),'k:','linewidth',2)
plot(tx2disp,squeeze(mean(N{2}(site,:,:),3)),'r','linewidth',2)
plot(tx2disp,squeeze(mean(N{3}(site,:,:),3)),'k','linewidth',2)
plot([0 0],[-4 4],'k:')
set(gca,'xlim',[-250 750],'xtick',[-250:250:750],'ylim',[-4 4]);
legend({'PD Standard','CTL Standard','PD Novelty','CTL Novelty'})

figure; 
site=36;      
subplot(2,1,1);
hold on
plot(tx2disp,squeeze(mean(S{3}(site,:,:),3)),'k','linewidth',2)
plot(tx2disp,squeeze(mean(T{3}(site,:,:),3)),'r','linewidth',2)
plot(tx2disp,squeeze(mean(N{3}(site,:,:),3)),'b','linewidth',2)
plot([0 0],[-4 4],'k:')
set(gca,'xlim',[-250 750],'xtick',[-250:250:750],'ylim',[-4 4]);
 
site=11;      
subplot(2,1,2);
hold on
plot(tx2disp,squeeze(mean(S{3}(site,:,:),3)),'k','linewidth',2)
plot(tx2disp,squeeze(mean(T{3}(site,:,:),3)),'r','linewidth',2)
plot(tx2disp,squeeze(mean(N{3}(site,:,:),3)),'b','linewidth',2)
plot([0 0],[-4 4],'k:')
set(gca,'xlim',[-250 750],'xtick',[-250:250:750],'ylim',[-4 4]);


topo1=find(tx2disp==300);
topo2=find(tx2disp==400);
figure; topoplot(squeeze(mean(mean(N{3}(:,topo1:topo2,:),2),3))-squeeze(mean(mean(S{3}(:,topo1:topo2,:),2),3)),BV_Chanlocs_60,'maplimits',[-4 4])
topo1=find(tx2disp==400);
topo2=find(tx2disp==600);
figure; topoplot(squeeze(mean(mean(T{3}(:,topo1:topo2,:),2),3))-squeeze(mean(mean(S{3}(:,topo1:topo2,:),2),3)),BV_Chanlocs_60,'maplimits',[-3 3]);  


%%
clear all; clc
Rootpath='Y:\EEG_Data\PDDys\EEG\ODDBALL\';

datalocation=[Rootpath,'Processed Data\'];   % Data are here
cd(datalocation);

load([datalocation,'SINGLETRIAL_ODDBALL_AGG.mat'],'S','T','N','ON_sx','OFF_sx','CTL_sx','tx2disp','ToTake')

% S, T, N to classify separately
% {1} = ON
% {2} = OFF
% {3} = CTL

RUN_CLASSIFY(N{1},N{3},'NOVELTY_ON_vs_CTL');
RUN_CLASSIFY(N{2},N{3},'NOVELTY_OFF_vs_CTL');
RUN_CLASSIFY(N{1},N{2},'NOVELTY_ON_vs_OFF');

RUN_CLASSIFY(S{1},S{3},'STANDARD_ON_vs_CTL');
RUN_CLASSIFY(S{2},S{3},'STANDARD_OFF_vs_CTL');
RUN_CLASSIFY(S{1},S{2},'STANDARD_ON_vs_OFF');

RUN_CLASSIFY(T{1},T{3},'TARGET_ON_vs_CTL');
RUN_CLASSIFY(T{2},T{3},'TARGET_OFF_vs_CTL');
RUN_CLASSIFY(T{1},T{2},'TARGET_ON_vs_OFF');




%%

load('NOVELTY_OFF_vs_CTL_LASSO.mat');


% ----- if there are more than 1 run...
TRN=squeeze(mean(LASSO_Probability_Train,1));
TST1=squeeze(mean(LASSO_Probability_Tst1,1));
TST2=squeeze(mean(LASSO_Probability_Tst2,1));
BETAS=squeeze(mean(LASSO_Betas,1));

% Setup time, unwrap BETAS
Nelec=60;
samples=12;
setlength=samples*Nelec; % 12 samples * 60 electrodes
bins=floor(250/samples + (250/samples)-1); % + 50% overlap bins
STstart=1;  STend=samples;
for STi=1:bins
    RANGE(STi,:)=(STstart:STend).*4;
    STstart = STstart+samples/2;  STend=STend+samples/2;
end
RANGE_TIMES=RANGE-254;
tx=RANGE_TIMES(:,6);
% Unwrap BETAS
for ti=1:40
    temp=BETAS(ti,:);
    Unwrap_BETAS(ti,:,:)=reshape(temp,60,12); clear temp;
end



figure;
subplot(2,1,1); hold on
plot(tx,TST2(:,2),'b','Linewidth',2);
plot(tx,TST2(:,3),'r','Linewidth',2);
plot([tx(1) tx(40)],[.5 .5],'k:','Linewidth',2);
plot([0 0],[0 1],'k:','Linewidth',2);
set(gca,'ylim',[0  1],'xlim',[-250 750]);
legend({'specificity','sensitivity'},'location','SouthEast');
xlabel('Time Bin Center (48 ms bins)');
ylabel('LASSO  Accuracy');
% % title('Oddball Targets vs. NULL');

MAX_BIN=25;
TO_PLOT=squeeze(mean(Unwrap_BETAS(MAX_BIN,:,:),3));

subplot(2,1,2); hold on
topoplot(TO_PLOT,BV_Chanlocs_60);
title(['Beta weights at:  ',num2str(tx(MAX_BIN)),' ms']); cbar;


%%
clear all; clc

Nelec=60;
samples=12;
setlength=samples*Nelec; % 12 samples * 60 electrodes
bins=floor(250/samples + (250/samples)-1); % + 50% overlap bins
STstart=1;  STend=samples;
for STi=1:bins
    RANGE(STi,:)=(STstart:STend).*4;
    STstart = STstart+samples/2;  STend=STend+samples/2;
end
RANGE_TIMES=RANGE-254;
tx=RANGE_TIMES(:,6);


figure; hold on
for ai=1:6
    if     ai==1, load('STANDARD_ON_vs_CTL_LASSO.mat'); xloc=1;  
    elseif ai==2, load('STANDARD_OFF_vs_CTL_LASSO.mat');xloc=2;  
    elseif ai==3, load('TARGET_ON_vs_CTL_LASSO.mat');   xloc=3;  
    elseif ai==4, load('TARGET_OFF_vs_CTL_LASSO.mat');  xloc=4;  
    elseif ai==5, load('NOVELTY_ON_vs_CTL_LASSO.mat');  xloc=5; 
    elseif ai==6, load('NOVELTY_OFF_vs_CTL_LASSO.mat'); xloc=6; 
    end
    
%    MAX_BIN=25;
   MAX_BIN=40;
    
    TST2=squeeze(mean(LASSO_Probability_Tst2,1));

    bar(xloc-.2,TST2(MAX_BIN,2),.4,'b');
    bar(xloc+.2,TST2(MAX_BIN,3),.4,'r');
    
    clear LASSO*;
end
plot([0 7],[.5 .5],'k:');
set(gca,'ylim',[.5  1],'xlim',[0 7],'xtick',[1:1:6],'xticklabel',{'S ON','S OFF','T ON','T OFF','N ON','N OFF'});
legend({'specificity','sensitivity'},'location','SouthEast');
xlabel('Condition');
ylabel('LASSO  Accuracy');
title('3 Auditory Oddball Fixed Effect Classification');

%%

figure; hold on
for ai=1:3
    if ai==1, load('STANDARD_OFF_vs_CTL_LASSO.mat');xloc=1;  
    elseif ai==2, load('TARGET_OFF_vs_CTL_LASSO.mat');  xloc=2;  
    elseif ai==3, load('NOVELTY_OFF_vs_CTL_LASSO.mat'); xloc=3; 
    end
    
%    MAX_BIN=25;
   MAX_BIN=40;
    
    TST2=squeeze(mean(LASSO_Probability_Tst2,1));

    bar(xloc-.2,TST2(MAX_BIN,2),.4,'b');
    bar(xloc+.2,TST2(MAX_BIN,3),.4,'r');
    
    clear LASSO*;
end
plot([0 4],[.5 .5],'k:');
set(gca,'ylim',[.5  1],'xlim',[0 4],'xtick',[1:1:3],'xticklabel',{'Standard','Target','Novelty'});
legend({'specificity','sensitivity'},'location','SouthEast');
xlabel('Condition');
ylabel('LASSO  Accuracy');
title('3 Auditory Oddball Fixed Effect Classification');



