load('dAccuracyData.mat');
load('tunValues.mat');
load('autolatencySummary.mat');

lstLength = length(dAccuracy);

time = dAccuracy{1}.time;
dir = dAccuracy{1}.dir;

tLength = length(time);

numTr = [];

divRatio = 6/13;

tBounds = [-40:5:200];
tWin = [-15 15];

dirSTD_haccu = nan(lstLength, length(tBounds));
dirSTD_laccu = nan(lstLength, length(tBounds));

eegPred_laccu = nan(length(dir), length(time), lstLength);
eegPred_haccu = nan(length(dir), length(time), lstLength);

corrValues = nan(lstLength, length(tBounds));
mAccu = [];
aSTD = nan(lstLength, length(tBounds));

sTime = 2; % time window 50 ~ 100 ms from stimulus onset

nTrials = [];

for n = 1:lstLength
	tQual = tuningValues(n).ev(:,sTime);
	select = find(tQual > 0); % removing trials with bad tunings.
	trNum = length(select);
	nTrials = [nTrials;trNum];
	divNum = floor(divRatio*trNum);
	currAccu = tuningValues(n).ev(select, sTime);
	mAccu = [mAccu;mean(currAccu)];
	accuVec = [currAccu transpose([1:trNum])];
	accuVec = sortrows(accuVec);
	lIdx = accuVec(1:divNum, 2);
	hIdx = accuVec(end-divNum+1:end, 2);
	for m = 1:length(tBounds)
		tBehav = find(time > tBounds(m)+tWin(1)+centerLat(n) & time <= tBounds(m)+tWin(2)+centerLat(n));
		bMode = rad2deg(circ_mean(deg2rad(dAccuracy{n}.behavDir(select, tBehav)), [], 2));
		behavValue = bMode;
		sC = find(behavValue >= mean(behavValue)-2.576*std(behavValue) & behavValue <= mean(behavValue)+2.576*std(behavValue));
		behavValue = behavValue(sC); % Removing outliers from behavioral data for more robust estimation of stadard deviation
		tmpAccu = currAccu(sC);
		selectIdx = select(sC);
		tmpVec = [tmpAccu transpose([1:length(sC)])];
		tmpVec = sortrows(tmpVec);
		divNum = floor(divRatio*length(sC));
		lIdx = tmpVec(1:divNum, 2);
		hIdx = tmpVec(end-divNum+1:end, 2);
		aSTD(n, m) = std(behavValue);
		dirSTD_haccu(n, m) = std(behavValue(hIdx));
		dirSTD_laccu(n, m) = std(behavValue(lIdx));
		mb = mean(behavValue);
		corrValues(n, m) = corr(abs(behavValue-mb), tmpAccu, 'type', 'spearman');
	end
	eegPred_laccu(:, :, n) = squeeze(mean(dAccuracy{n}.eegPred(:, :, selectIdx(lIdx)), 3));
	eegPred_haccu(:, :, n) = squeeze(mean(dAccuracy{n}.eegPred(:, :, selectIdx(hIdx)), 3));
end


% Figure 4B
figure(1)
subplot(2,1,1)
pcolor(time, dir, mean(eegPred_haccu, 3));
shading flat
caxis([-2 2]);
colormap(jet);
title('high accuracy');
xlabel('Time(ms)');
ylabel('Direction');
subplot(2,1,2)
pcolor(time, dir, mean(eegPred_laccu, 3));
shading flat
caxis([-2 2]);
colormap(jet);
title('low accuracy');
xlabel('Time(ms)');
ylabel('Direction');

% Figure 4D
figure(2)
errorbar(tBounds, mean(corrValues), std(corrValues)./sqrt(15), 'k');
title('Correlations');
xlabel('Time(ms):from pursuti latency');
ylabel('Rho');
% Figure 4F
selected = find(tBounds == 50);
figure(3)
stdH = dirSTD_haccu(:, selected);
stdL = dirSTD_laccu(:, selected);
scatter(stdH, stdL, 'ko');
axis square
set(gca, 'xlim', [9 21], 'ylim', [9 21]);
hold on
plot([9 21], [9 21], 'k-.');
xlabel('STD:high accuracy');
ylabel('STD:low accuracy');

% Figure 4G
figure(4)
behavCorrs = corr(mAccu, aSTD, 'type', 'spearman');
plot(tBounds, behavCorrs, 'k');
xlabel('time(ms):from pursuit latency');
ylabel('Rho');
