function Results = FFS_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]
%%% bin ratio(ratio) [vector in range [0,1]] (optional):
% linear zone as a ratio of bin_width (default: .7:.05:1)
%%% keep_indx(session number) [mxt] (optional):
% specifies selected sessions/neurons for each time point

%%%%%%%%%%%%%%%%%%%% OUTPUT
%%%% VEC:       FFS VEC [tx1]
%%%% VEC_std :  SE of FFS VEC [tx1]
%%%% nRV:      FFS [tx1]
%%%% nRV_std:  SE of FFS [tx1]
%%%% EVC:       FFS EVC [tx1]
%%%% EVC_std:   SE of FFS EVC [tx1]
%%%% nSI:      FF0 [tx1]
%%%% nSI_std:  SE of FF0 [tx1]

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


n = length(raster);

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

keep_indx = params.keep_indx;
bin_ratio = params.bin_ratio;

% Assign all mothod results to this struct
Results = struct();

% Linear Zone
X = (bin_width * bin_ratio)' * 1e-3;


% Estimate FF
FF = zeros(n,length(bin_ratio),length(times));
spike_count = zeros(n,length(bin_ratio),length(times));
% spike_count_var = zeros(n,length(bin_ratio),length(times));
for i = 1:n
    raster_local = raster{i};
    parfor j = 1:length(bin_ratio)
        binwidth_local = bin_width * bin_ratio(j);
        bin_start  = times - binwidth_local/2;
        bin_end = times + binwidth_local/2;
        bin_edge = sort([bin_start,bin_end]);
        spike_cell_local = cellfun(@(x) histcounts(x,bin_edge),...
            raster_local,'UniformOutput',false);
        spike_mat_local = cell2mat(spike_cell_local)';
        spike_mat_local = spike_mat_local(1:2:end,:);
        FF(i,j,:) = nanvar(spike_mat_local,0,2)./...
            nanmean(spike_mat_local,2);
        spike_count(i,j,:) = nanmean(spike_mat_local,2);
        %         spike_count_var(i,j,:) = nanvar(spike_mat_local,0,2);
    end
end

% Keep index for mean-matching if provided
mean_FF = zeros(length(bin_ratio),length(times));
se_FF = zeros(length(bin_ratio),length(times));
for t_num = 1:size(keep_indx,2)
    mean_FF(:,t_num) = squeeze(nanmean(FF(keep_indx(:,t_num),:,t_num),1));
    se_FF(:,t_num) = squeeze(nanstd...
        (FF(keep_indx(:,t_num),:,t_num),[],1))/sqrt(size(keep_indx,1));
end

n_ = size(keep_indx,1);

% Fit line to average FF in linear zone
nSI_all = zeros(size(keep_indx,1),length(times));
nSI_all_SE = nSI_all;
nRV_all = nSI_all;
nRV_all_SE = nSI_all;
parfor t = 1:length(times)
    for nue_num = 1:n_
        y_neu  = squeeze(FF(nue_num,:,t));
        mdl = fitlm(X,y_neu);
        
        nSI_all(nue_num,t)= mdl.Coefficients.Estimate(1)
        nSI_all_SE(nue_num,t)= mdl.Coefficients.SE(1);
        
        nRV_all(nue_num,t) = mdl.Coefficients.Estimate(2);
        nRV_all_SE(nue_num,t) = mdl.Coefficients.SE(2);
    end
end

Results.FF = mean_FF(end,:);
Results.FF_SE = se_FF(end,:);


if n_ == 1
    Results.nSI = nSI_all;
    Results.nRV= nRV_all;
    mean_count = squeeze(spike_count(:,end,:));
    Results.VEC = Results.nRV .* mean_count' * bin_width * 1e-3;
    Results.EVC = Results.nSI .* mean_count';
    
    Results.nSI_SE = nSI_all_SE;
    Results.nRV_SE = nRV_all_SE;
    Results.EVC_SE = nSI_all_SE .* mean_count';
    Results.VEC_SE = nRV_all_SE .* mean_count' * bin_width * 1e-3;
    
else
    Results.nSI = nanmean(nSI_all);
    Results.nRV= nanmean(nRV_all);
    mean_count = nanmean(squeeze(spike_count(:,end,:)));
    Results.VEC = nanmean(Results.nRV .* mean_count * bin_width * 1e-3);
    Results.EVC = nanmean(Results.nSI .* mean_count);
    
    Results.nSI_SE = nanstd(nSI_all_SE);
    Results.nRV_SE = nanstd(nRV_all_SE);
    Results.EVC_SE = nanstd(nSI_all_SE .* mean_count);
    Results.VEC_SE = nanstd(nRV_all_SE .* mean_count * bin_width * 1e-3);
end

end



