function [Acc,ConfMat] = Compute_accuracyNN_matrix(AmpResp,AllSamples,LOO)

%% This function calculates the similarity matrix which estimates the correlation between pairs
%% of sounds after correction for trial to trial noise
%% INPUT
% AmpResp = the responses to be correlations : ROI * stim * trials
% AllSamples = 1 (default) calculates all possible pairs of trials, 0
% instead does this for a random 50 combiations which is a lot quicker
% LOO = 1 train using leave one out cross validation
% LOO = 0 train using two halves of trials for cross validation

%% OUTPUT
% Acc : accuracy for identifying each sound using nearest neighbour
% ConfMat : confusion matrix

%% intialize
if not(exist('AllSamples'))
    AllSamples = 1;
end

if not(exist('LOO'))
    LOO = 0;
end
if LOO ==1
    AllSamples = 1;
end
% clean up NaN cells
NaN_cells = find(isnan(sum(squeeze(nanmean(AmpResp,3)'))));
AmpResp(NaN_cells,:,:) = [];
disp([num2str(length(NaN_cells)) ' ROIs out of ' num2str(size(AmpResp,1)) 'excluded because they had NaNs'])

% Not all data sets have all of the trials
maxgoodtrials = find(squeeze(sum(sum(isnan(AmpResp),1),2))==0,1,'last');
disp([num2str(maxgoodtrials) ' trials without NaNs'])
maxgoodtrials = floor(maxgoodtrials/2)*2;

for neur = 1:size(AmpResp,1)
    AmpResp(neur,:,:) = AmpResp(neur,:,:) - nanmean(AmpResp(neur,:));
end

if LOO
    % Use leave one out
    AllTrialPairs = combnk(1:maxgoodtrials,maxgoodtrials-1);
    TrialsTrain = AllTrialPairs;
    TrialsTest = [maxgoodtrials:-1:1];
    TrialsTest = reshape(TrialsTest,maxgoodtrials,1);
else
    % Get the trial pairs
    AllTrialPairs = combnk(1:maxgoodtrials,floor(maxgoodtrials/2));
    TrialsTrain = AllTrialPairs;
    TrialsTest = fliplr(TrialsTrain')';
end


if AllSamples
    
else
    % a subset of trial pairs
    RandTrials = randperm(size(TrialsTrain,1),10);
    TrialsTrain = TrialsTrain(RandTrials,:);
    TrialsTest = TrialsTest(RandTrials,:);
    
end


%% Get accuracy
for perm = 1:size(TrialsTest,1)
    
    TrainVect = squeeze(nanmean(AmpResp(:,:,TrialsTrain(perm,:)),3));
    if length(TrialsTest(perm,:))>1
        TestVect = squeeze((AmpResp(:,:,TrialsTest(perm,:))));
    else
        TestVect = ((AmpResp(:,:,TrialsTest(perm,:))));
    end
    
    
    for k = 1:size(TestVect,3)
        
        C = corr(TrainVect,squeeze(TestVect(:,:,k)));
        
        [val,ind(k,:)] = max(C);
        
    end
    
    Acc(perm,:) = nanmean((ind-repmat([1:size(TestVect,2)],[size(TestVect,3),1]))==0,1);
    
    for sd = 1:size(ind,2)
        for sd2 = 1:size(ind,2)
            ConfMat(perm,sd,sd2) = sum(ind(:,sd) ==sd2);
        end
    end
    
    for sd = 1:size(ind,2)
        for sd2 = 1:size(ind,2)
            ConfMat(perm,sd,sd2) = mean(ind(:,sd) ==sd2);
        end
    end
end