function [sumStats] = eegFwdModelLCoutDecodeAccuracy(outMat, dPad, Elc, btNum, ev, outLier)
%
% A routine for forward encoding/decoding model analysis
%
% outMat: output of prepareEEGfwdModel.m function
% Elc: electrodes indices
% 
% eegFwdModelLCoutDecodeAccuracy
% 07/11/2018 Joonyeol Lee
% Leave-one-condition out and bootstrap for comparing encoding accuracy and behavior
%


rng('shuffle', 'twister');

timeTick = outMat.eegTimeTick;
tDur = length(timeTick);
NoT = length(outMat.rIdx);
eegTmp = outMat.binEEG;
cTmp = outMat.fwdModel;
stimTmp = outMat.stim;
modesTmp = outMat.behavior;
beSpdTmp = outMat.behavSpd;
beDirTmp = outMat.behavDir;
fitQualTmp = outMat.fitQual;

accSpdTmp = outMat.accelSpd;
accDirTmp = outMat.accelDir;

randIdx = transpose(randperm(NoT));

eegData = eegTmp(:, :, randIdx);
cMat = cTmp(:, randIdx);
stimInfo = stimTmp(randIdx,:);
modes = modesTmp(randIdx,:);
beSpd = beSpdTmp(randIdx,:);
beDir = beDirTmp(randIdx,:);
accSpd = accSpdTmp(randIdx,:);
accDir = accDirTmp(randIdx,:);
fitQual = fitQualTmp(randIdx);

oIdx = find(fitQual < ev);

eegData(:, :, oIdx) = [];
cMat(:, oIdx) = [];
stimInfo(oIdx,:) = [];
modes(oIdx,:) = [];
beSpd(oIdx,:) = [];
beDir(oIdx,:) = [];
accSpd(oIdx,:) = [];
accDir(oIdx,:) = [];
fitQual(oIdx) = [];

wMat = cell(tDur, 1);

dirCond = sort(unique(stimInfo(:,3)));

numCond = length(dirCond);

for c = 1:numCond
	currCondDir = dirCond(c);
	select = find(stimInfo(:,3) == currCondDir);
		
	if exist('outLier') & outLier == 1
		cr = 2.576;

		tMode = modes(select, 3);
		dMode = modes(select, 2);
		sMode = modes(select, 1);
		
		oLiers = find(tMode > mean(tMode)+cr*std(tMode) | tMode < mean(tMode)-cr*std(tMode) | dMode > mean(dMode)+cr*std(dMode) | dMode < mean(dMode)-cr*std(dMode) | sMode > mean(sMode)+cr*std(sMode) | sMode < mean(sMode)-cr*std(sMode));
		
		eegData(:, :, select(oLiers)) = [];
		cMat(:, select(oLiers)) = [];
		stimInfo(select(oLiers),:) = [];
		modes(select(oLiers),:) = [];
		beSpd(select(oLiers),:) = [];
		beDir(select(oLiers),:) = [];
		accSpd(select(oLiers),:) = [];
		accDir(select(oLiers),:) = [];
		fitQual(select(oLiers)) = [];
	end
	
end

NoT = size(modes,1);

select = [];

currCondDir = dirCond(3);

select = find(stimInfo(:,3) == currCondDir);

tsData = eegData(:,:,select);
trData = eegData;
trData(:,:,select) = [];
trC = cMat;
trC(:,select) = [];
cMatPredOri = nan(size(cMat,1), tDur, length(select));

for i = 1:tDur
	trainData = squeeze(trData(Elc, i, :));
	testData = squeeze(tsData(Elc, i, :));
	tmpCov = trC*transpose(trC);
	CCinv = pinv(tmpCov);
	B1 = trainData;
	W = B1*transpose(trC)*CCinv;
	wMat{i} = W;
	tmpCovW = transpose(W)*W;
	WWinv = pinv(tmpCovW);
	for j = 1:length(select)
		B2 = testData(:, j);
		cMatPredOri(:, i, j) = WWinv*transpose(W)*B2;
	end
end
	
btCond = nan(length(select), btNum);
for bt = 1:btNum
	currTr = randperm(NoT);
	btCond(:, bt) = currTr(1:length(select));
end

pred_boot = nan(size(cMat,1), tDur, btNum);

parfor bt = 1:btNum
	cChoice = btCond(:, bt);
	tsData = eegData(:,:,cChoice);

	cMatPred = nan(size(cMat,1), tDur, length(cChoice));

	for i = 1:tDur
		testData = squeeze(tsData(Elc, i, :));
		W = wMat{i};
		tmpCov = transpose(W)*W;
		WWinv = pinv(tmpCov);
		for j = 1:length(cChoice)
			B2 = testData(:, j);
			cMatPred(:, i, j) = WWinv*transpose(W)*B2;
		end
	end
	pred_boot(:, :, bt) = squeeze(mean(cMatPred,3));
end

mBoot = squeeze(mean(pred_boot, 3));
sBoot = squeeze(std(pred_boot, [], 3));

cMatPredCorrected = nan(size(cMat,1), tDur, length(select));

for tr = 1:length(select)
	cMatPredCorrected(:, :, tr) = (squeeze(cMatPredOri(:, :, tr))-mBoot)./sBoot;
end

selectedPred = cMatPredCorrected;
selectedSpd = beSpd(select, :);
selectedDir = beDir(select, :);
selectedModes = modes(select,:);
selectedSpd_acc = accSpd(select, :);
selectedDir_acc = accDir(select, :);

switch dPad
case 0
	dirAngles = dirCond;
case 1
	dirAngles = [dirCond(1)-30;dirCond;dirCond(end)+30];
case 2
	dirAngles = [dirCond(1)-2*30;dirCond(1)-1*30;dirCond;dirCond(end)+1*30;dirCond(end)+2*30];
case 3
	dirAngles = [dirCond(1)-3*30;dirCond(1)-2*30;dirCond(1)-1*30;dirCond;dirCond(end)+1*30;dirCond(end)+2*30;;dirCond(end)+3*30];
end

decodeAccuracy = nan(size(selectedDir));

condIdx = find(dirAngles == currCondDir);

for q = 1:size(decodeAccuracy, 2)
	parfor p = 1:size(decodeAccuracy, 1)
		currAvg = squeeze(selectedPred(:, q, p));
		pre = currAvg(condIdx-4:condIdx-1);
		post = currAvg(condIdx+4:-1:condIdx+1);
		mv = (pre+post)./2;
		mv = [mv;currAvg(condIdx)];
		dirV = dirAngles(condIdx-4:condIdx);
		pfit = polyfit(dirV, mv, 1);
		decodeAccuracy(p, q) = pfit(1);
	end
end

sumStats.mode = selectedModes;
sumStats.behavDir = selectedDir;
sumStats.behavSpd = selectedSpd;
sumStats.accDir = selectedDir_acc;
sumStats.accSpd = selectedSpd_acc;
sumStats.decodeAccuracy = decodeAccuracy;
sumStats.time = timeTick;
sumStats.dir = dirAngles;
sumStats.eegPred = selectedPred;
sumStats.eegOri = cMatPredOri;