function Result = var_decom(raster, params)

% Inputs:
% raster: nx1 cell: the timing of spikes as array in n trials with 0
%             as the stim onset
% params: struct with following att.s:
%         t_start:   time of start of recording in ms
%         t_end:     time of end of recording in ms
%         bin_width: bin width in ms
%         slide:     sliding window length in ms
%         keep_indx: mean matching index
%         kass_params: parameters of Kass:
%               params
%               K [ind]: number of correlated samples
%               bin_count_kass
%               smooth_bin
%               cv_count [int] (optional, default: 100)
%         FFS_params:  parameters of FFS:
%               bin ratio(ratio) [vector in range [0,1]] (optional)
%
% Citation: Saleh F, Fakharian M & Ghazizadeh A "Stimulus presentation can enhance spiking irregularity across subcortical and cortical regions." 2021

raster = reshape(raster(:,1),size(raster,1),1);

n = length(raster);

% params
t_start = params.t_start;
t_end = params.t_end;
bin_width = params.bin_width;
slide = params.slide;

if(~isfield(params,'min_spike'))
    params.min_spike = 0;
end

% find average spike count of each neuron
valid_index = ones(n,1);

for i = 1:n
    raster_local = raster{i};
    spike_mat = cell2mat(cellfun(@(x) ...
        histcounts(x,t_start:bin_width:t_end),...
        raster_local,'UniformOutput',false));
    if mean(spike_mat(:)) < params.min_spike
        valid_index(i) = 0;
    end
end

Result.valid_index = valid_index;
Result.keep_perc = sum(valid_index)/length(valid_index)*100;
raster = raster(valid_index == 1);
n = length(raster);

% finding times
Result.times = t_start+(bin_width/2):slide:t_end-(bin_width/2);

if(~isfield(params,'keep_indx'))
    params.keep_indx = repmat(1:n,length(Result.times),1)';
end

% method_params
params.kass_params.bin_width = bin_width;
params.FFS_params.bin_width = bin_width;

i_max = (bin_width/slide) - 1;

Result.VEC.FFS = zeros(size(Result.times));
Result.VEC.kass = zeros(size(Result.times));
Result.EVC.FFS = zeros(size(Result.times));
Result.EVC.kass = zeros(size(Result.times));
Result.nSI.FFS = zeros(size(Result.times));
Result.nSI.kass = zeros(size(Result.times));
Result.nRV.FFS = zeros(size(Result.times));
Result.nRV.kass = zeros(size(Result.times));

Result.VEC.FFS_SE = zeros(size(Result.times));
Result.VEC.kass_SE = zeros(size(Result.times));
Result.EVC.FFS_SE = zeros(size(Result.times));
Result.EVC.kass_SE = zeros(size(Result.times));
Result.nSI.FFS_SE = zeros(size(Result.times));
Result.nSI.kass_SE = zeros(size(Result.times));
Result.nRV.FFS_SE = zeros(size(Result.times));
Result.nRV.kass_SE = zeros(size(Result.times));

Result.FR = zeros(size(Result.times));
Result.FR_SE = zeros(size(Result.times));
Result.FF = zeros(size(Result.times));
Result.FF_SE = zeros(size(Result.times));
Result.Var = zeros(size(Result.times));
Result.Var_SE = zeros(size(Result.times));

for i = 0:i_max
    
    perccount(i,i_max)
    
    times = t_start+(bin_width/2)+i*slide:bin_width:t_end-(bin_width/2);
    
    params.kass_params.keep_indx = params.keep_indx...
        (:,ismember(Result.times, times));
    params.FFS_params.keep_indx  = params.keep_indx...
        (:,ismember(Result.times, times));
    
    FFS = FFS_method(raster,times,params.FFS_params);
    kass = kass_method(raster,times,params.kass_params);
    
    Result.FR(ismember(Result.times, times)) = kass.fr;
    Result.FR_SE(ismember(Result.times, times)) = kass.fr_SE;
    Result.FF(ismember(Result.times, times)) = kass.FF;
    Result.FF_SE(ismember(Result.times, times)) = kass.FF_SE;
    Result.Var(ismember(Result.times, times)) = kass.Var;
    Result.Var_SE(ismember(Result.times, times)) = kass.Var_SE;
    
    Result.VEC.FFS(ismember(Result.times, times)) = FFS.VEC;
    Result.VEC.kass(ismember(Result.times, times)) = kass.VEC;
    Result.EVC.FFS(ismember(Result.times, times)) = FFS.EVC;
    Result.EVC.kass(ismember(Result.times, times)) = kass.EVC;
    Result.nSI.FFS(ismember(Result.times, times)) = FFS.nSI;
    Result.nSI.kass(ismember(Result.times, times)) = kass.nSI;
    Result.nRV.FFS(ismember(Result.times, times)) = FFS.nRV;
    Result.nRV.kass(ismember(Result.times, times)) = kass.nRV;
    
    Result.VEC.FFS_SE(ismember(Result.times, times)) = FFS.VEC_SE;
    Result.VEC.kass_SE(ismember(Result.times, times)) = kass.VEC_SE;
    Result.EVC.FFS_SE(ismember(Result.times, times)) = FFS.EVC_SE;
    Result.EVC.kass_SE(ismember(Result.times, times)) = kass.EVC_SE;
    Result.nSI.FFS_SE(ismember(Result.times, times)) = FFS.nSI_SE;
    Result.nSI.kass_SE(ismember(Result.times, times)) = kass.nSI_SE;
    Result.nRV.FFS_SE(ismember(Result.times, times)) = FFS.nRV_SE;
    Result.nRV.kass_SE(ismember(Result.times, times)) = kass.nRV_SE;

end
