
%plot((0:size(EEG.data,2)-1)/5,abs(fft(squeeze(mean(EEG.data,3))'/(size(EEG.data,2)))))

function main4(subjectName,repeats,gap,savePath)
%%
% repeats = 1000;
% gap =2;
% savePath = '../nGoggle3_private/';
% subjects = {'sbj50','sbj51','sbj52','sbj53','sbj54','sbj55','sbj58','sbj59','sbj60','sbj61','sbj62','sbj63'};
subjects = subjectName;
for sbjs = 1:size(subjects,2)
    name = subjects{sbjs};
    subjectName =  {[name '_R_A_raw_asred.set'],[name '_L_A_raw_asred.set']};%,[name '_R_cali_raw_asred.set']};
    subjectName2 = {[name '_R_B_raw_asred.set'],[name '_L_B_raw_asred.set']};%,[name '_L_cali_raw_asred.set']};
    mypath = '../nGoggle3_private/bp550h4/';
    if(isempty(savePath))
        savePath = '../nGoggle3_private/';
    end
    for s = 1:size(subjectName,2)
        sbjName = subjectName{s};
        warning('off','all')
        numHarmonics = 2;
        channels = [3:8];
        freq = 8:0.2:11.8;
        resolution = 5;
        condition =1;%block 8.8 hz
        shift1 = 0.135;
        shift2 = 5.135;
        
        EEG = pop_loadset('filename',subjectName{s},'filepath',mypath);
        EEG = pop_select( EEG,'time',[shift1 shift2] );
        data = EEG.data;
        
        EEG3 = pop_loadset('filename',subjectName2{s},'filepath',mypath);
        EEG3 = pop_select( EEG3,'time',[shift1 shift2] );
        data2 = EEG3.data;
        
        minTrials = floor(min(size(data2,3),size(data,3)));
        %%
        correlation = 0.96;
        myCorr1 = tril(corr(squeeze(data(9,:,:)),squeeze(data(9,:,:))),-1);
        [rr1,cc1,jj] = find(myCorr1>correlation);
        relatedIndex1 = union(rr1,cc1);
        myCorr2 = tril(corr(squeeze(data2(9,:,:)),squeeze(data2(9,:,:))),-1);
        [rr2,cc2,jj] = find(myCorr2>correlation);
        relatedIndex2 = union(rr2,cc2);
        minTrials = min(length(relatedIndex1),length(relatedIndex2));
        fprintf('----min(%i,%i) trials with %i correlation----\n',length(relatedIndex1),length(relatedIndex2),correlation);
        %-------
        
        fprintf('==START==\n');
        jump = floor(linspace(1,minTrials,floor(minTrials/gap)));
        
        
        
        %         jump = 10:gap:minTrials;%floor(minTrials*1/4):gap:floor(minTrials*3/4)
        for trials = jump
            for rpt = 1:repeats
                if(rpt==repeats)
                    fprintf('Trials: %i/%i, Repeats: %i\n',trials,minTrials, rpt);
                end
                tt = relatedIndex1(randperm(minTrials,trials));%round(minTrials/2));
                tt2 = relatedIndex2(randperm(minTrials,trials));
                [out1Corr(:,trials,rpt),out1Power(:,trials,rpt),afterCCA] = my_FFT_or_CCA(squeeze(mean(data(channels,:,tt),3))',resolution,freq,numHarmonics,'CCA');
                [out2Corr(:,trials,rpt),out2Power(:,trials,rpt),afterCCA2] = my_FFT_or_CCA(squeeze(mean(data2(channels,:,tt2),3))',resolution,freq,numHarmonics,'CCA');
                [outTCCACorr(:,trials,rpt),...
                    outTCCAPower(:,trials,rpt),...
                    outTCCACorr2(:,trials,rpt),...
                    outTCCAPower2(:,trials,rpt)] = myCCA4_twoDataSet_predict(data(channels,:,:),relatedIndex1,tt,...
                    data2(channels,:,:),relatedIndex2,tt2,...
                    numHarmonics,freq,resolution);
                trialsValue{trials,rpt} = tt2;
                
                %                 tt2 = randperm(minTrials,trials);
                %                 [out1Corr(:,trials,rpt),out1Power(:,trials,rpt),afterCCA] = my_FFT_or_CCA(squeeze(mean(data(channels,:,tt2),3))',resolution,freq,numHarmonics,'CCA');
                %                 [out2Corr(:,trials,rpt),out2Power(:,trials,rpt),afterCCA2] = my_FFT_or_CCA(squeeze(mean(data2(channels,:,tt2),3))',resolution,freq,numHarmonics,'CCA');
                %                 [outTCCACorr(:,trials,rpt),...
                %                     outTCCAPower(:,trials,rpt),...
                %                     outTCCACorr2(:,trials,rpt),...
                %                     outTCCAPower2(:,trials,rpt)] = myCCA4_twoDataSet_predict(data(channels,:,:),1:size(data,3),tt2,...
                %                     data2(channels,:,:),1:size(data2,3),tt2,...
                %                     numHarmonics,freq,resolution);
                %                 trialsValue{trials,rpt} = tt2;
            end
        end
        
        %% plotting ([], red(sick), black(normal))
        %         plot_by_ring_4(jump, out1Corr, out2Corr, freq,['CCA' ', ' subjectName{s} subjectName2{s}]);
        %         plot_by_ring_4(jump, out1Power, out2Power, freq,['CCA' ', ' subjectName{s} subjectName2{s}]);
        %         plot_by_ring_4(jump, outTCCACorr, outTCCACorr2, freq,['CCA-T' ', ' subjectName{s} subjectName2{s}]);
        %         plot_by_ring_4(jump, outTCCAPower, outTCCAPower2, freq,['CCA-T' ', ' subjectName{s} subjectName2{s}]);
        if(strcmpi(subjectName{s}(7:9),'L_A'))
            EEG.laCorr = out1Corr;
            EEG.lbCorr = out2Corr;
            EEG.laPower = out1Power;
            EEG.lbPower = out2Power;
            EEG.la2 = outTCCACorr;
            EEG.lb2 = outTCCACorr2;
            EEG.la2Power = outTCCAPower;
            EEG.lb2Power = outTCCAPower2;
        else
            EEG.raCorr = out1Corr;
            EEG.rbCorr = out2Corr;
            EEG.raPower = out1Power;
            EEG.rbPower = out2Power;
            EEG.ra2 = outTCCACorr;
            EEG.rb2 = outTCCACorr2;
            EEG.ra2Power = outTCCAPower;
            EEG.rb2Power = outTCCAPower2;
        end
        EEG.trivalValue = trialsValue;
        EEG.jump = jump;
        sbjName(findstr('.set',sbjName):length(sbjName)+1)='2.set';
        EEG = pop_saveset( EEG, 'filename',sbjName,'filepath',savePath);
        %% plotting (normal, sick, name, condition,algorithm)
        %         plot2D(out1Corr,out2Corr,sbjName,3,'CCA')
        %% plotting (cca result,fft result)
        %         plot_psd_5sec_together(out1Corr,out2Corr,freq,['CCA corr' ', ' EEG.filename]);
        %         plot_psd_5sec_together(out1Power,out2Power,freq,['CCA power' ', ' EEG.filename]);
        plot_psd_5sec_together(outTCCACorr,outTCCACorr2,freq,['TCCA corr' ', ' EEG.filename]);
        %         plot_psd_5sec_together(outTCCAPower,outTCCAPower2,freq,['TCCA power' ', ' EEG.filename]);
        %% plotting scaled SAP
        %     name = EEG3.filename;
        %     name = name(1:5)
        %     plot_SAP(out1Corr,['CCA-s-A, ' name],1);
        %     plot_SAP(out2Corr,['CCA-s-B, ' name],2);
        %         plot_SAP(EEG.sap(:,1),'SAP-left',4);
        %         plot_SAP(EEG.sap(:,2),'SAP-right',4);
        %         plot_SAP3(out1Corr,out2Corr,['(A+B)CCA-std, ' subjectName{s}],'std');
        %         plot_SAP3(out1Corr,out2Corr,['(A+B)CCA-mean, ' subjectName{s}],'mean');
        %         plot_SAP3(outTCCACorr,outTCCACorr2,['(A+B)T-std, ' subjectName{s}],'std');
        %         plot_SAP3(outTCCACorr,outTCCACorr2,['(A+B)T-mean, ' subjectName{s}],'mean');
        %         plot_SAP_3test(outTCCACorr,outTCCACorr2,['(A+B)T-PSD-std, ' name]);
        
        clear out1Corr out2Corr out1Power out2Power outTCCACorr outTCCACorr2 outTCCAPower outTCCAPower2;
    end
end
end