function [tValues, pValues] = icatb_compare_frequency_bins(param_file, group1, group2)
%% Determine PSD and perform group differences
%
% Inputs:
%
% 1. param_file - ICA Parameter file
% 2. group1 - Group 1 subject numbers
% 3. group2 - Group 2 subject numbers
%
% Outputs:
%
% 1. tValues - T values of size components by (NPOINT_FFT_GROUP_COMPARISON/2).
% 2. pValues - p values of size components by (NPOINT_FFT_GROUP_COMPARISON/2).
%

%% Modality type
modalityType = icatb_get_modality;

if ~strcmpi(modalityType, 'fmri')
    error('Group comparison utility works only with the GIFT');
end

%% Load  defaults
icatb_defaults;
global PARAMETER_INFO_MAT_FILE;

% colors
global FONT_COLOR;
global UI_FONTUNITS;
global UI_FONTNAME;
global UI_FS;

% FFT defaults for group comparison
global NPOINT_FFT_GROUP_COMPARISON;
global NUM_BINS_GROUP_COMPARISON;

if isempty(NPOINT_FFT_GROUP_COMPARISON)
    NPOINT_FFT_GROUP_COMPARISON = 300;
end

if isempty(NUM_BINS_GROUP_COMPARISON)
    NUM_BINS_GROUP_COMPARISON = 6;
end

%% N point FFT
nPointFFT = NPOINT_FFT_GROUP_COMPARISON;

%% Number of frequency bins
numBins = NUM_BINS_GROUP_COMPARISON;

%% Check bins
checkBin = mod(nPointFFT, 2*numBins);
if (checkBin ~= 0)
    error('Error:CheckBin', 'NPOINT_FFT_GROUP_COMPARISON (%d) must be exactly divisible by 2*NUM_BINS_GROUP_COMPARISON (%d).\nPlease check icatb_defaults.m file', ...
        nPointFFT, 2*numBins);
end

%% load parameters file
filterP = ['*', PARAMETER_INFO_MAT_FILE, '*.mat'];

if ~exist('param_file', 'var')
    param_file = icatb_selectEntry('typeEntity', 'file', 'title', 'Select a valid parameter file', 'filter', filterP);
end

drawnow;

[outputDir, param_file, extn] = fileparts(param_file);

if isempty(outputDir)
    outputDir = pwd;
end

param_file = fullfile(outputDir, [param_file, extn]);

load(param_file);

if ~exist('sesInfo', 'var')
    error('Error:sesInfo', 'Selected file %s is not a valid ICA parameter file', param_file);
end

%% Get information from parameter file
numOfSub = sesInfo.numOfSub;
numOfSess = sesInfo.numOfSess;
numComp = sesInfo.numComp;
timePoints = sesInfo.diffTimePoints;
flagTimePoints = sesInfo.flagTimePoints;

if ~strcmpi(flagTimePoints, 'same_time_points')
    error('Error:TimePoints', 'Timepoints are not the same between datasets. This utility works with the datasets that have same timepoints.');
end

%% Select groups
group1Name = 'Group 1';
group2Name = 'Group 2';

subjectString = cell(numOfSub, 1);
for nSub = 1:numOfSub
    subjectString{nSub} = ['Subject ', num2str(nSub)];
end

subjectString = str2mat(subjectString);

if ~exist('group1', 'var')
    % Group Name and datasets
    [group1Name, group1] = icatb_select_groups_gui(subjectString, 'Group 1', 'select_subjects_psd');
end

if isempty(group1)
    error('Group 1 subjects are not selected');
end

if (size(subjectString, 1) < max(group1))
    error('Error:Group1', 'Maximum of group1 (%d) exceeds the maximum of subjects (%d)\n',  max(group1), size(subjectString, 1));
end

if ~exist('group2', 'var')
    % Group2 name and datasets
    [group2Name, group2] = icatb_select_groups_gui(subjectString, 'Group 2', 'select_subjects_psd');
end

if isempty(group2)
    error('Group 2 subjects are not selected');
end

if (size(subjectString, 1) < max(group2))
    error('Error:Group2', 'Maximum of group2 (%d) exceeds the maximum of subjects (%d)\n',  max(group2), size(subjectString, 1));
end

drawnow;

%% Print some information
disp(sprintf('N point FFT used for group comparison is %d', nPointFFT));
disp(sprintf('Number of frequency bins used is %d', numBins));
disp(sprintf('Selected subjects for group 1 are %s ', num2str(group1(:)')));
disp(sprintf('Selected subjects for group 2 are %s ', num2str(group2(:)')));

%% End for selecting groups

%nPointFFT = timePoints(1);
nPoints = ceil(nPointFFT / 2);

fprintf('\n');

disp('Computing fft of subject time courses ...');

%% Initialise variables
count = 0;
% fft plots for all subjects
fftPlots = zeros(numOfSub, numComp, nPoints);
% bin count
fftBins = zeros(numOfSub, numComp, numBins);

%% Loop over subjects
for nSub = 1:numOfSub
    % Initialise mean of time courses over sessions
    meanTC = zeros(timePoints(1), numComp);
    %% Loop over sessions
    for nSess = 1:numOfSess
        count = count + 1;
        % Load subject MAT file
        subjectMATFileName = fullfile(outputDir, [sesInfo.calibrate_components_mat_file, num2str(nSub), '-', num2str(nSess), '.mat']);
        load(subjectMATFileName, 'tc');
        if (size(tc, 1) ~= timePoints(count))
            tc = tc';
        end
        meanTC = meanTC + tc;
    end
    %% End loop over sessions
    meanTC = meanTC ./ numOfSess;
    
    %% Loop over components
    for nComp = 1:numComp
        tmp = abs(fft(icatb_normdet(meanTC(:, nComp), 1), nPointFFT));
        tmp = tmp(1:ceil(length(tmp) / 2));
        binSize = length(tmp) / numBins;
        fftPlots(nSub, nComp, :) = tmp;
        %% Loop over bins
        for nBin = 1:numBins
            % Take the average for the frequencies within each bin
            fftBins(nSub, nComp, nBin) = mean(tmp((1:binSize) + binSize*(nBin - 1)));
        end
        %% End loop over bins
    end
    %% End loop over components
    
    clear meanTC;
    
end
%% End loop over subjects

disp('Done computing fft of subject time courses');

%% Compute mean power of each group
mean_power_group1 = squeeze(mean(fftPlots(group1,:,:)));
mean_power_group2 = squeeze(mean(fftPlots(group2,:,:)));

clear fftPlots;

fprintf('\n');

%% Create two matrices for between-group comparisons.

disp(['Comparing frequency bins of groups (', group1Name, ' & ', group2Name, ') ...']);

group1FFTBins = fftBins(group1,:,:);
group2FFTBins = fftBins(group2,:,:);

%% Total number of subjects
numSubjects = length(group1) + length(group2);

%% Initialise tValues and pValues
tValues = zeros(numComp, size(group1FFTBins, 3));
pValues = zeros(numComp, size(group1FFTBins, 3));

%% Loop over components
for nComp = 1:numComp
    
    %% Group 1 mean and std
    group1Mean = squeeze(mean(group1FFTBins(:, nComp,:)));
    group1Std = squeeze(std(group1FFTBins(:, nComp,:)));
    
    %% Group2 mean and std
    group2Mean = squeeze(mean(group2FFTBins(:, nComp, :)));
    group2Std = squeeze(std(group2FFTBins(:, nComp, :)));
    
    %% T and p values
    temp = (group1Mean - group2Mean) ./ ((group1Std + group2Std)/(2*sqrt(numSubjects - 1)));
    temp = temp(:)';
    tValues(nComp, :) = temp;
    pValues(nComp, :) = 1 - icatb_spm_Tcdf(tValues(nComp, :), (numSubjects - 1));
end
%% End loop over components

disp('Done comparing frequency bins');

fprintf('\n');

%% XTICK Label
xtickVal = (1:numBins)'./binSize;
freqValues = round(xtickVal*100)/100;
freqValues = freqValues(:);

xtickLabel = cellstr(num2str(freqValues));
%xtickLabel = repmat(xtickLabel, size(tValues, 2), 1);

%% Plot Group 1 - Group 2
%Multi bar plot
plotTitle = {[group1Name, ' - ', group2Name]};

if (numComp > 4)
    elemPerAxes = 4;
else
    elemPerAxes = numComp;
end

legendStr = cellstr([repmat('Comp ', numComp, 1), num2str((1:numComp)')]);

figTitle = 'Power Spectrum Analysis of ICA timecourses';
icatb_multi_bar_plot(tValues', 'XTickLabel', xtickLabel, 'xlabel', 'Frequency (Hz)', 'ylabel', 'T-Values', 'title', figTitle, ...
    'axesTitle', plotTitle, 'elemPerAxes', elemPerAxes, 'legend', legendStr);

%% Save results
outFile = fullfile(outputDir, [sesInfo.userInput.prefix, '_comparison_frequency_bins.mat']);
icatb_save(outFile, 'tValues', 'pValues', 'group1Name', 'group2Name', 'mean_power_group1', 'mean_power_group2');
disp(['Results are saved in file ', outFile]);

fprintf('\n');