currDir = cd;
bDir = [currDir '\data\'];

dirList = {'SSK','SMK','SJH','MJS','LYM','KTJ','KSH','KKM','KDH','KCH','JJH','BSW','ABJ','LCM','LJJ'}; 

lstLength = length(dirList);

allChannel = [1:64];

dPad = 3;
bkOpt = 1;

popEncoding = [];
popEncodingSlopes = [];
stepSize = 5;
wSize = 10;
wholeDirs = [-180:30:150];
bkWin = [-99:-50];
bootNum = 500;
tWin = [-100 600];

evV = 0.5;

for n = 1:lstLength
	load([bDir dirList{n} '_organizedEEGdata.mat']);
	outMat = prepareEEGfwdModel(outData, tWin, wSize, stepSize, dPad, {'Cosine', 15}, bkOpt, bkWin);
	[wMat, eegData, cMat, tInfo] = eegFwdModelCondOutWgtOnly(outMat, allChannel, evV, 1);	
	tLength = length(outMat.eegTimeTick);
	encodedData = nan(12, tLength, tLength);
	slopeData = nan(tLength, tLength);

	for m = 1:tLength
		currW = wMat(m, :);	
		[stimSortedPred, zScoredPred, dirCon, sumStats] = eegFwdModelTemporalGenCondOut_boot(currW, eegData, cMat, tInfo, dPad, allChannel, bootNum);	
		encodedData(:, :, m) = sumStats.boot;
		zeroIdx = find(sumStats.dir == 0);
		for l = 1:tLength
			currAvg = sumStats.boot(:, l);
			pre = currAvg(zeroIdx-4:zeroIdx-1);
			post = currAvg(zeroIdx+4:-1:zeroIdx+1);
			mv = (pre+post)./2;
			mv = [mv;currAvg(zeroIdx)];
			dirV = sumStats.dir(zeroIdx-4:zeroIdx);
			p = polyfit(dirV, mv, 1);
			slopeData(m, l) = p(1);
		end
	end
	popEncoding = cat(4, popEncoding, encodedData);
	popEncodingSlopes = cat(3, popEncodingSlopes, slopeData);
	n
end

% popEncoding: direction tuning responses
% popEncodingSlopes: cross-temporal generalization map of slope of direction tuning curve (encoding accuracy)
timetick = outMat.eegTimeTick;
dir = wholeDirs;

encodingDiag = nan(lstLength, length(timetick));
tuningDiag = nan(length(dir), length(timetick), lstLength);
selectTuning = nan(length(dir), length(timetick), lstLength);
selectTime = 60;
tIdx = find(timetick == selectTime);
for i = 1:length(timetick)
	encodingDiag(:, i) = popEncodingSlopes(i, i, :);
	tuningDiag(:,i,:) = squeeze(popEncoding(:, i, i, :));
	selectTuning(:,i,:) = squeeze(popEncoding(:, i, tIdx, :));
end
% encodingDiag: Diagonal of the cross-temporal generalization map

% Testing the significance of encoding accuracy in the diagonal of the CTG.
[Sig, P_test] = CBPtest(encodingDiag, repmat(mean(encodingDiag(:, 12:20), 2), 1, length(timetick)), timetick, 0.01, 0.005, 10000);

% Testing the significant of encoding accuracy in the CTG 
[Perm, SigArea] = CBPtest_FW(permute(popEncodingSlopes, [3 1 2]), mean(encodingDiag(:, 12:20), 2), 0.01, 0.005, 10000);

% Figure 2A
figure(1);
pcolor(timetick, dir, mean(tuningDiag,3));
shading flat
caxis([-3 3]);
colormap(jet);

% Figure 2B
figure(2)
errorbar(timetick, mean(encodingDiag), std(encodingDiag)./sqrt(lstLength-1));
hold on
scatter(Sig, zeros(length(Sig), 1), 'ro', 'filled');
grid on

% Figure 3A
figure(3)
pcolor(timetick, timetick, mean(popEncodingSlopes, 3));
shading flat;
caxis([-0.03 0.03]);
colormap(jet);
figure(4) % plot for significance
pcolor(timetick, timetick, SigArea);
shading flat

% Figure 3B
figure(5)
pcolor(timetick, dir, mean(selectTuning,3));
shading flat
caxis([-2 2]);
colormap(jet);

save('baseData.mat', 'popEncoding', 'popEncodingSlopes', 'encodingDiag', 'tuningDiag', 'selectTuning');