% Appies HSMM-MVPA as in Anderson 2016 to find the boundaries of cognitive
% processing stages


bumpM   = 4; % number of bumps, 4 bumps == 5 stages

nBuMax  = max(bumpM);
nModel  = length(bumpM);
nPca    = 32; % 32 PCA -> more than 95% variance, 25 -> more than 85%; 20 -> 80%  variance on PCA 
nSj     = 18;
rep     = 150;
shape   = 2; % gamma distributoin, shape parameter
roiDataR = 'nScoreAveR'; % 'nScoreAve' or 'nScorePca'

dataPath = './data';
dataRead = '/vars4hsmmBL_R.mat';

roiTypeR  = ['roi' roiDataR(end-3:end)];
roiTypeR(end) = [];
saveName = ['/hmmBu' num2str(bumpM) 'R' roiTypeR 'Map3StagBL.mat']; 


addpath('./hsmm/');


info = {['test ', num2str(bumpM),' bumps n repetitions one models: R; map 4 stage, PCA cov by subject, Baseline']};

load([dataPath, dataRead], roiDataR, 'condR','xR','yR','subjtR')

% right hand
data_   = eval(roiDataR);
data    = single(data_(:,1:nPca));
x       = xR;
y       = yR;

condR_ = condR;
condR_(logical( condR==2 | condR==3 )) = 1; % fan1-target
condR_(logical( condR==4 | condR==5 )) = 2; % fan2-target
condR_(logical( condR==6 | condR==7 )) = 3; % fan1-foil
condR_(logical( condR==8 | condR==9 )) = 4;  % fan2-foil

cond  = condR_;
subj = subjtR;
% left hand

clear data_ roiDataR xR yR condR condR_ subjtR 

% define map settings, 3th stage (retrieval) different gamma
map = [1, 1, 1, 1, 1;     
       1, 1, 2, 1, 1;
       1, 1, 3, 1, 1;
       1, 1, 4, 1, 1];  

gamIni  = repmat([2 200], nBuMax+1,1);
magIni  = -1 + 2*rand(nPca, nBuMax, rep); 

% bump kernel correlation, removed form inside HSMM
data = calcBumps(data{1,1});
lkhr = zeros(repnModel); % [repetitions,  hand response, models(number bumps),] 
parpool(30);
parfor r = 1:rep
    % do HSMM
        for m = 1:nModel % iterate models
            [lkhr(r,m),~,~,~]=hsmmEEG50CondMapAllDNew(data,cond,magIni(:,1:bumpM(m),r),gamIni(1:bumpM(m)+1,:),1,x,y,map,shape);
        end
    end
end
disp('End of random initializations')

% recompute best models
lkh     = cell(1,nModel); % hand, model(nBump)
bumpMag = cell(1,nModel);
flatGam = cell(1,nModel);
proBump = cell(1,nModel);
for m = 1:nModel
        [~,iBest]  = max(lkhr(:,m));
        [lkh{1,m},bumpMag{1,m},flatGam{1,m},proBump{1,m}] = hsmmEEG50CondMapAllDNew(data,cond,magIni(:,1:bumpM(m),iBest),gamIni(1:bumpM(m)+1,:),1,x,y,map,shape);
end

% LOOCV, leave-one-out cross validation

lkhCV = zeros(nSj,nModel);
parfor sj = 1:nSj
    for m = 1:nModel % iterate models  
        sTrain  = logical(subj ~= sj);
        sTest   = logical(subj == sj);
        [~, magTrain, gammTrain, ~] = hsmmEEG50CondMapAllDNew(data,cond(sTrain),bumpMag{1,m},flatGam{1,m}, 1 ,x{h}(sTrain),y(sTrain),map,shape);
        [ lkhCV(sj, m), ~, ~, ~]  = hsmmEEG50CondMapAllDNew(data,cond(sTest), magTrain,    gammTrain,    0 ,x(sTest),y(sTest),map,shape);
    end
end


%save(['./dataSave' saveName], 'lkhr','lkh','bumpMag','flatGam','proBump','info','nPca', 'rep', 'x','y','cond','lkhCV','subj')
%format longG
disp(char(datetime(now,'ConvertFrom','datenum')))


delete(gcp('nocreate'))

clear all
close all