%% This code generates Figure 4.a-f
% Citation: Saleh F, Fakharian M & Ghazizadeh A "Stimulus presentation can enhance spiking irregularity across subcortical and cortical regions." 2021

addpath(genpath('../lib'))
%%
clc;
clear;

% Two SNR Samples from Good and Bad sessions
%load Datasets/SNR_Sample_Good.mat
load Datasets/SNR_Sample_Bad.mat

raster_ = raster;

% start/end times
params.t_start = -200;
params.t_end = 800;
% sliding window and binwidth in ms
params.slide = 50;
params.bin_width = 200;

% CV2_Local
cv = var_cv(raster_,params);
params.kass_params.cv2_all_neu = cv.cv2_all_neu;
params.kass_params.K = 14;

% Var Decompostion
result = var_decom(raster_,params);
result.cv = cv;

% VarCE
result.varce = var_varce(raster_{1},params);

t_start = params.t_start;
t_end = params.t_end;
bin_width = params.bin_width;

raster_local = raster_{1};
n = numel(raster_local);

bin_size = 20;

t_rate = t_start:t_end-bin_size;

fr = zeros(n,length(t_rate));

spike_mat = cell2mat(cellfun(@(x) histcounts(x,t_start:t_end),...
    raster_local,'UniformOutput',false));

for i = 1:n
    perccount(i,numel(raster))
    fr_temp = spike_mat(i,:);
    fr(i,:) = movmean(fr_temp,bin_size,'Endpoints','discard')*1e3;
end

rate = mean(fr);
rate_sd = std(fr);
t_rate = t_start:t_end-bin_size;

result.raster = raster_;
result.rate.rate = rate;
result.rate.t_rate = t_rate;
result.rate.rate_trial = fr;
result.rate.rate_sd = rate_sd;

%% plot
fig = figure('Color', 'White','unit','normalized','position',[0 0 .9 .4]);

T = [-150,300];

red = [255, 35, 71]/255;
blue = [0, 145, 207]/255;
orange = [255, 111, 58]/255;
skyblue = [0, 188, 236]/255;
green = [0, 148, 50]/255;
lightgreen = [163, 203, 56]/255;
viol = [142, 68, 173]/255;
pink = 	[217, 128, 250]/255;

k = 20;

raster = result.raster{1};
n = numel(raster);

if n<k
    k = n;
end

rng('default');
tr_inds = randsample(n,k);

r = 4;
c = 4;

subplot(r,c,[1,5]);
plotErrorbar(result.rate.t_rate,result.rate.rate,...
    result.rate.rate_sd,'black');
xlim([t_start+bin_width/2,t_end-bin_width/2]);
xticklabels({})

subplot(r,c,[c*(2:3)+1]);

raster_ = raster(tr_inds);

plot_raster(raster_,'black',.5);
axis tight;
box off;
xlim([t_start+bin_width/2,t_end-bin_width/2]);

subplot(r/2,c,3)
yyaxis left
plotErrorbar(result.times,...
    result.nRV.kass,result.nRV.kass_SE,green);
ylabel('nRV')
xlim([t_start+bin_width/2,t_end-bin_width/2]);
yyaxis right
plotErrorbar(result.times,...
    result.FF,result.FF_SE,red);

ax = gca;
ax.YAxis(1).Color = green;
ax.YAxis(2).Color = red;

subplot(r/2,c,7)
yyaxis left
plotErrorbar(result.times,...
    result.nSI.kass,result.nSI.kass_SE,viol);
ylabel('nSI')

ax = gca;
ax.YAxis(1).Color = viol;
ax.YAxis(2).Color = 'k';

yyaxis right
plotErrorbar(result.times,...
    result.cv.cv2,result.cv.cv2_SE,'b');
ylabel('CV2')
ylim([0,.8])

xlim([t_start+bin_width/2,t_end-bin_width/2]);

subplot(r/2,c,4)
h1=plotErrorbar(result.varce.times,...
    double(result.varce.VEC.VarCE),double(result.varce.VEC.VarCE_SE),red);
hold on
h2=plotErrorbar(result.times,result.VEC.kass,result.VEC.kass_SE,blue);
ylabel('VEC')
xlim([t_start+bin_width/2,t_end-bin_width/2]);legend ([h1, h2],'VarCE','Kass','Location','best');legend('boxoff');

subplot(r/2,c,8)
h1=plotErrorbar(result.varce.times,...
    result.varce.nSI.VarCE .*result.FR*bin_width*1e-3,...
    result.varce.nSI.VarCE .*result.FR_SE*bin_width*1e-3,red);

hold on
h2=plotErrorbar(result.times,result.EVC.kass,result.EVC.kass_SE,blue);
ylabel('EVC')
xlim([t_start+bin_width/2,t_end-bin_width/2]); legend ([h1, h2],'VarCE','Kass','Location','best');legend('boxoff');