function Results = cv2_method(raster,times,params)
%%%%%%%%%%%%%%%%%%%% input

%%%%%%%%%%%% raster(ms) [nx1 cell]:
% n session cells of trials containing spike time array,
% aligned by time of stimulus on zero

%%%%%%%%%%%% times(ms) [1xt]
% time vector of evaluation times.
% in range [bin_width/2, end-bin_width/2]
% This function does not support sliding window, please use var_decom.
% example: 100:100:1000

%%%%%%%%%%%% params
%%% bin_width(ms) [double]
% number of cross-validations
% histogram bin count
%%% keep_indx(session number) [mxt] (optional):
% specifies selected sessions/neurons for each time point

%%%%%%%%%%%%%%%%%%%% OUTPUT
%%%% shape_K:         Shape parameter      [tx1]
%%%% shape_K_std:     SE of Shape parameter [tx1]
%%%% cv2_est:         cv2_est [tx1]
%%%% cv2_est_std:     SE of cv2_est [tx1]
%%%% cv_n_dist:       histogram of cv_n [txhist_bin_count] 

%%%% note: cv_n is within [0,2] so each bin in cv_n_dist represents a 2/hist_bin_count
%%%% length

% Citation: Saleh F, Fakharian M & Ghazizadeh A "Stimulus presentation can enhance spiking irregularity across subcortical and cortical regions." 2021


n = length(raster);

% Check Inputs
try
    bin_width = params.bin_width;
catch
    error('bin_width is required!');
end
if(~isfield(params,'hist_bin_count'))
    params.hist_bin_count = 10;
end
if(~isfield(params,'keep_indx'))
    params.keep_indx = repmat(1:n,length(times),1)';
end
if(~isfield(params,'cv_count'))
    params.cv_count = 100;
end

keep_indx = params.keep_indx;
cv_count = params.cv_count;
hist_bin_count = params.hist_bin_count;
hist_edge = linspace(0,2,hist_bin_count+1);

shape_K_tot = zeros(n,length(times));
cv2_est_tot = zeros(n,length(times));
cv_n_dist_tot =  zeros(n,length(times),hist_bin_count);

shape_K_tot_SE = zeros(1,length(times));
cv2_est_tot_SE = zeros(1,length(times));

min_k_trials = 20;

parfor i = 1:n
    raster_local_ = raster{i};
    raster_local = cellfun(@(x) reshape(x,length(x),1),...
        raster_local_,'UniformOutput',false);
    
    % Cross Validation
    shape_K_cv = nan(cv_count,length(times));
    dist_cv_n = zeros(length(times),hist_bin_count);
    for tim_bin=1:length(times)
        period = [times(tim_bin)-bin_width/2,times(tim_bin)+bin_width/2];
        isi_cell_temp = cellfun(@(x) ...
            diff(x((x>period(1)) & (x<period(2)))),...
            raster_local,'UniformOutput',false);
        
        isi_cell_temp_count = cellfun(@(x) length(x),isi_cell_temp,...
            'UniformOutput',false);
        isi_mat_temp_count = cell2mat(isi_cell_temp_count);
        % remove less than one isi in each bin
        isi_cell_not_zero = isi_cell_temp(isi_mat_temp_count>1);
        cv_n = cellfun(@(x) 4 * abs((x./(x+[x(2:end);nan])) - 0.5),...
                        isi_cell_not_zero,'UniformOutput',false);
        k_trials = cell2mat(cv_n);
        k_trials(isnan(k_trials)) = [];
        dist_cv_n(tim_bin,:) = histcounts(k_trials,hist_edge,...
            'Normalization', 'probability');
        if length(k_trials)> min_k_trials
            shape_K_cv(:,tim_bin) =  bootstrp(cv_count,...
                @(x) 2/mean(x.^2) - 0.5 ,k_trials);
        end
    end
    
    shape_K_cv(shape_K_cv == inf) = NaN;
    shape_K_cv(shape_K_cv == -inf) = NaN;
    
    if n == 1
        shape_K_tot_SE(i,:) = nanstd(shape_K_cv,[],1);
        cv2_est_tot_SE(i,:) = nanstd(1./shape_K_cv,1);        
    end
    
    shape_K_tot(i,:) = nanmean(shape_K_cv,1);
    cv2_est_tot(i,:) = nanmean(1./shape_K_cv,1);
    
    cv_n_dist_tot(i,:,:) = dist_cv_n;

end

if n == 1
    Results.shape_K = shape_K_tot;
    Results.cv2_est = cv2_est_tot;

    
    Results.shape_K_SE = shape_K_tot_SE;
    Results.cv2_est_SE = cv2_est_tot_SE;
    Results.cv_n_dist  = cv_n_dist_tot;
else
    for t=1:length(times)
        nn = size(keep_indx,1);
        
        Results.shape_K(t) = nanmean(shape_K_tot(keep_indx(:,t),t),1);
        Results.cv2_est(t) = nanmean(cv2_est_tot(keep_indx(:,t),t),1);
        
        Results.cv_n_dist(t,:) = nanmean(cv_n_dist_tot(...
            keep_indx(:,t),t,:),1);

        Results.shape_K_SE(t) = nanstd(shape_K_tot(...
            keep_indx(:,t),t),0,1)/sqrt(nn);
        Results.cv2_est_SE(t) = nanstd(cv2_est_tot(...
            keep_indx(:,t),t),0,1)/sqrt(nn);

    end
    
end

Results.cv2_est_tot = cv2_est_tot;
end
