function varargout = exp_barumerli2021(varargin)
%EXP_BARUMERLI2021 - Experiments of Barumerli et al. (2021)
%
%   Usage: [] = exp_barumerli2021(flag) 
%
%
%   To display Fig.2 use :
%
%     exp_barumerli2021('fig2');
%
%   To display Fig.3 use :
%
%     exp_barumerli2021('fig3');
%
%   To display Fig.4 use :
%
%     exp_barumerli2021('fig4');
%
%   To display exp_middlebrooks1999 use :
%
%     exp_barumerli2021('exp_middlebrooks1999');
%
%   To display exp_macpherson2003 use :
%
%     exp_barumerli2021('exp_macpherson2003');
%
%   To display exp_best2005 use :
%
%     exp_barumerli2021('exp_best2005');
%
%   Url: http://amtoolbox.sourceforge.net/amt-0.10.0/doc/experiments/exp_barumerli2021.php

% Copyright (C) 2009-2020 Piotr Majdak and the AMT team.
% This file is part of Auditory Modeling Toolbox (AMT) version 1.0.0
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.


%   See also: barumerli2021 barumerli2021_featureextraction

%   AUTHOR: Roberto Barumerli
%   Information Engineering Dept., University of Padova, Italy, 2021

%% ------ Check input options ---------------------------------------------

definput.import = {'amt_cache'};
definput.keyvals.MarkerSize = 6;
definput.keyvals.FontSize = 9;
definput.keyvals.repetitions = 300;
definput.keyvals.sbj = [];

definput.flags.type = {'missingflag', 'fig2', 'fig3', 'fig4', 'exp_middlebrooks1999', 'exp_macpherson2003', 'exp_best2005', 'verify_calibration'};
definput.flags.plot = {'plot', 'no_plot'};

definput.flags.plot_type = {'interp','scatter'};
definput.flags.redo = {'no_redo_fast','redo_fast'};
definput.flags.test = {'no_test','test'};

[flags,kv]  = ltfatarghelper({},definput,varargin);

if flags.do_missingflag
  flagnames=[sprintf('%s, ',definput.flags.type{2:end-2}),...
             sprintf('%s or %s',definput.flags.type{end-1},...
             definput.flags.type{end})];
  error('%s: You must specify one of the following flags: %s.', ...
      upper(mfilename),flagnames);
end

%% ------ fig.2 from paper
if flags.do_fig2
    
    data_majdak = data_majdak2010('Learn_M');
    data_majdak([1:5]) = [];

    calibrations = amt_load('barumerli2021', 'barumerli2021_calibration.mat');
    calibrations = calibrations.cache.value;
    
    if size(calibrations.sigma, 1) ~= length(data_majdak)
        warning('sigma values not enough for the provided subejcts')
    end

    fig2 = amt_cache('get', 'barumerli2021_fig2',flags.cachemode);

    % Preallocation
    if isempty(fig2)
        s = 1;

        calibrations = calibrations.sigma(s,2); % select pge
        
%         sofa = SOFAload(['dtf_' lower '.sofa']);
        sofa = SOFAload(fullfile(SOFAdbPath,'barumerli2021',['ARI_' upper(data_majdak(s).id) '_hrtf_M_dtf 256.sofa']));

        mtx = data_majdak(s).mtx;

        [template, target] = barumerli2021_featureextraction(sofa, ...
                                'pge', 'monoaural',...
                                'targ_az', mtx(:, 1), ...
                                'targ_el', mtx(:, 2));


        m_noprior = barumerli2021('template', template, ...
                                 'target', target, ...
                                 'num_exp', 1, ...
                                 'sigma_lateral',calibrations.values(1), ...
                                 'sigma_lateral_bis', calibrations.values(2), ...
                                 'sigma_polar', calibrations.values(3), ...
                                 'sigma_prior', [],...
                                 'sigma_sensorimotor', []);

         m_prior = barumerli2021('template', template, ...
                                 'target', target, ...
                                 'num_exp', 1, ...
                                 'sigma_lateral',calibrations.values(1), ...
                                 'sigma_lateral_bis', calibrations.values(2), ...
                                 'sigma_polar', calibrations.values(3), ...
                                 'sigma_prior', calibrations.values(5),...
                                 'sigma_sensorimotor', []);

         m_prior_motor = barumerli2021('template', template, ...
                                 'target', target, ...
                                 'num_exp', 1, ...
                                 'sigma_lateral',calibrations.values(1), ...
                                 'sigma_lateral_bis', calibrations.values(2), ...
                                 'sigma_polar', calibrations.values(3), ...
                                 'sigma_prior', calibrations.values(5),...
                                 'sigma_sensorimotor', calibrations.values(4));
                             
        results  = {m_noprior, m_prior, m_prior_motor, mtx};
        amt_cache('set', 'barumerli2021_fig2', results);
    end

    m_noprior = fig2{1,1};
    m_prior = fig2{1,2};
    m_prior_motor = fig2{1,3};
    mtx = fig2{1,4};
    
    FontSize = 8;
    Size = 10;
    fig = figure('Units', 'points', 'Position', [1000 1000 430 150]);
    
    subplot(1,3,1)
    scatter(mtx(:, 6), mtx(:, 8), Size, 0.6*[1 1 1], 'filled');
    hold on
    m_noprior(~((m_noprior(:, 6) > -30) & (m_noprior(:, 6) < 210)),:)=[];
    scatter(m_noprior(:, 6), m_noprior(:, 8),  Size, 0.2*[1 1 1]);  
    grid on
    ax = gca;
    axis equal
    set(ax, 'YTick',[-90, 0, 90, 180, 270], 'XTick',[-90, 0, 90, 180, 270], 'XLim', [-90, 270], 'YLim', [-90, 270]);
    ax.XLabel.String  = '      Actual Angle [deg]';
    ax.TitleFontWeight = 'normal';
    ax.Title.String  = 'Maximum likelihood';
    ax.YLabel.String  = 'Estimated Angle [deg]';
    
    subplot(1,3,2)
    scatter(mtx(:, 6), mtx(:, 8), Size, 0.6*[1 1 1], 'filled');
    hold on
    m_prior(~((m_prior(:, 6) > -30) & (m_prior(:, 6) < 210)),:)=[];
    scatter(m_prior(:, 6), m_prior(:, 8), Size, 0.2*[1 1 1]); 
    grid on
    ax = gca;
    axis equal
    set(ax, 'YTick',[-90, 0, 90, 180, 270], 'XTick',[-90, 0, 90, 180, 270], 'XLim', [-90, 270], 'YLim', [-90, 270]);
    ax.XLabel.String  = '      Actual Angle [deg]';
    ax.TitleFontWeight = 'normal';
    ax.Title.String  = 'w/ prior';
    set(ax,'YTickLabel','')
    
    subplot(1,3,3)
    scatter(mtx(:, 6), mtx(:, 8), Size, 0.6*[1 1 1], 'filled');
    hold on
    m_prior_motor(~((m_prior_motor(:, 6) > -30) & (m_prior_motor(:, 6) < 210)),:)=[];
    scatter(m_prior_motor(:, 6), m_prior_motor(:, 8), Size, 0.2*[1 1 1]); 
    grid on
    ax = gca;
    axis equal
    set(ax, 'YTick',[-90, 0, 90, 180, 270], 'XTick',[-90, 0, 90, 180, 270], 'XLim', [-90, 270], 'YLim', [-90, 270]);
    ax.XLabel.String  = '      Actual Angle [deg]';
    ax.TitleFontWeight = 'normal';
    ax.Title.String  = 'w/ sensorimotor';   
    set(ax,'YTickLabel','')
end

%% -------------- fig 3 paper --------------------------
if flags.do_fig3
    data_majdak = data_majdak2010('Learn_M');
    data_majdak([1:5]) = [];

    %calibrations = amt_cache('get', 'barumerli2021_calibration',flags.cachemode);
    calibrations = amt_load('barumerli2021', 'barumerli2021_calibration.mat');
    calibrations = calibrations.cache.value;
   
    
    if size(calibrations.sigma, 1) ~= length(data_majdak)
        warning('sigma values not enough for the provided subejcts')
    end

    fig3 = amt_cache('get', 'barumerli2021_fig3',flags.cachemode);

    % Preallocation
    if isempty(fig3)
        s = 1;
        calibrations = calibrations.sigma(s,2); % select pge
        
%         sofa = SOFAload(['dtf_' lower(data_baseline(s).id) '.sofa']);
        sofa = SOFAload(fullfile(SOFAdbPath,'barumerli2021',['ARI_' upper(data_majdak(s).id) '_hrtf_M_dtf 256.sofa']));
        
        mtx = data_majdak(s).mtx;

        [template, target] = barumerli2021_featureextraction(sofa, ...
                                'pge', 'monoaural',...
                                'targ_az', mtx(:, 1), ...
                                'targ_el', mtx(:, 2));

        m_motor = barumerli2021('template', template, ...
                                 'target', target, ...
                                 'num_exp', 1, ...
                                 'sigma_lateral',calibrations.values(1), ...
                                 'sigma_lateral_bis', calibrations.values(2), ...
                                 'sigma_polar', calibrations.values(3), ...
                                 'sigma_prior', calibrations.values(5),...
                                 'sigma_sensorimotor', calibrations.values(4));

        results  = {m_motor, mtx};
        amt_cache('set', 'barumerli2021_fig3', results);
    end

    m_motor = fig3{1,1};
    mtx= fig3{1,2};

    %% FIGURE
    FontSize = 9;
    Size = 10;
    fig = figure('Units', 'Points', 'Position', [1e3 1e3 245 150]);

    %% SCATTER
    subplot(1,2,1)
    scatter(mtx(:, 5), mtx(:, 7), Size, 0.6*[1 1 1],'filled');  % shadow
    hold on
    scatter(m_motor(:, 5), m_motor(:, 7), Size, 0.2*[1 1 1]);  % shadow
    grid on
    ax = gca;
    axis equal
    set(ax, 'YTick',[-90 -45 0 45 90], 'XTick',[-90 -45 0 45 90], 'XLim', [-100, 100], 'YLim', [-100, 100], 'FontSize',FontSize)

    ax.XLabel.String  = 'Actual Angle [deg]';
    ax.YLabel.String  = 'Estimated Angle [deg]';
    ax.YLabel.FontSize = FontSize;
    ax.XLabel.FontSize = FontSize;
    ax.Title.String  = 'Lateral';
    ax.TitleFontWeight = 'normal';
    title('Lateral')
    

    subplot(1,2,2)
    scatter(mtx(:, 6), mtx(:, 8), Size, 0.6*[1 1 1], 'filled');  % shadow
    hold on
    scatter(m_motor(:, 6), m_motor(:, 8), Size, 0.2*[1 1 1]);  % shadow
    grid on
    ax = gca;
    axis equal
    set(ax, 'YTick',[-90, 0, 90, 180, 270], 'XTick',[-90, 0, 90, 180, 270], 'XLim', [-90, 270], 'YLim', [-90, 270], 'FontSize',FontSize)
    ax.XLabel.String  = 'Actual Angle [deg]';
    ax.YLabel.FontSize = FontSize;
    ax.XLabel.FontSize = FontSize;
    ax.Title.String  = 'Polar';
    ax.TitleFontWeight = 'normal';
end

%% -------------- fig 4 paper --------------------------
if flags.do_fig4
    data_majdak = data_majdak2010('Learn_M');
    data_majdak([1:5]) = [];

    %calibrations = amt_cache('get', 'barumerli2021_calibration',flags.cachemode);
    calibrations = amt_load('barumerli2021', 'barumerli2021_calibration.mat');
    calibrations = calibrations.cache.value;

    
    if size(calibrations.sigma, 1) ~= length(data_majdak)
        warning('sigma values not enough for the provided subejcts')
    end

    fig4 = amt_cache('get', 'barumerli2021_fig4',flags.cachemode);

    % compute results
    if isempty(fig4)
        for s=1:length(data_majdak)
            calib = calibrations.sigma(s,2); % select pge

%             sofa = SOFAload(['dtf_' lower(data_baseline(s).id) '.sofa']);
            sofa = SOFAload(fullfile(SOFAdbPath,'barumerli2021',['ARI_' upper(data_majdak(s).id) '_hrtf_M_dtf 256.sofa']));

            mtx = data_majdak(s).mtx;

            [template, target] = barumerli2021_featureextraction(sofa, ...
                                    'pge', 'monoaural',...
                                    'targ_az', mtx(:, 1), ...
                                    'targ_el', mtx(:, 2));

            m = barumerli2021('template', template, ...
                                     'target', target, ...
                                     'num_exp', 1, ...
                                     'sigma_lateral',calib.values(1), ...
                                     'sigma_lateral_bis', calib.values(2), ...
                                     'sigma_polar', calib.values(3), ...
                                     'sigma_prior', calib.values(5),...
                                     'sigma_sensorimotor', calib.values(4));

            results{s}  = m;
        end
        amt_cache('set', 'barumerli2021_fig4', results);
    end

    %% COMPUTE RMS 
    lat = [-90, -40, 0, 40, 90];
    lat_label = [-65, -20, 20, 65];

    pol =     [-30 30 150 210];
    pol_label = [0 90 180];

    % results
    rmsL = zeros(length(data_majdak), length(lat)-1);
    rmsL_real = zeros(length(data_majdak), length(lat)-1);

    rmsP = zeros(1, length(pol)-1);
    rmsP_real = zeros(1, length(pol)-1);

    for s = 1:length(data_majdak)
        m = fig4{s};
        mtx = data_majdak(s).mtx;

        for i=1:length(lat)-1
            % real
            mtx_temp = mtx((mtx(:, 5) > lat(i) & (mtx(:, 5) < lat(i+1))),:);
            rmsL_real(s,i) = localizationerror(mtx_temp, 'rmsL');

            % simulation
            m_temp = m((m(:, 5) > lat(i) & (m(:, 5) < lat(i+1))),:);
            rmsL(s,i) = localizationerror(m_temp, 'rmsL');
        end

        mtx(abs(mtx(7,:)) <= 30,:)=[];
        m(abs(m(7,:)) <= 30,:)=[];
        for i=1:length(pol)-1
            m_temp = mtx((mtx(:, 6) > pol(i) & (mtx(:, 6) < pol(i+1))),:);
            rmsP_real(s,i) = localizationerror(m_temp, 'rmsPmedianlocal');

            m_temp = m((m(:, 6) > pol(i) & (m(:, 6) < pol(i+1))),:);
            rmsP(s,i) = localizationerror(m_temp, 'rmsPmedianlocal');    
        end
    end

    %% FIGURE
    FontSize = 9;
    Size = 36;
    fig = figure('Units', 'Points', 'Position', [1e3 1e3 510 180]);
    
    sp = 1;
    for s = 1:length(data_majdak)
        ax(sp) = subplot(2, length(data_majdak), sp);
        ax(sp).Position(4) = ax(s).Position(4)/1.612;
    
        % lateral
        scatter(lat_label, rmsL_real(s,:), Size, 0.8*[1 1 1], 'filled')
        hold on
        scatter(lat_label, rmsL(s,:), Size, 0.2*[1 1 1])
        set(gca, 'YLim', [5 20], 'Ytick', [5:5:20], 'Xtick', lat_label, 'Xlim', [-90, 90],'FontSize',FontSize)
        title({upper(data_majdak(s).id)})
        grid on

        % polar
        ax(length(data_majdak)+sp) = subplot(2, length(data_majdak), length(data_majdak) + sp);
        scatter(pol_label, rmsP_real(s,:), Size, 0.8*[1 1 1], 'filled')
        hold on
        scatter(pol_label, rmsP(s,:), Size, 0.2*[1 1 1])
        set(gca, 'Xtick', pol_label, 'Xlim', [-90, 270], 'Ylim', [0 60], 'Ytick', [0:20:60],'FontSize',FontSize)
        grid on
        
        sp = sp + 1;
    end
    
    set(ax([2:5, 7:10]),'YTickLabel','')

    for i=6:10
        ax(i).Title.String  = [];
        ax(i).YLabel.FontSize = FontSize;
    end

    ax(8).XLabel.String  = '      Actual Angle [deg]';
    ax(8).XLabel.FontSize = FontSize;

    for i=1:5
        ax(i).TitleFontWeight = 'normal';
        ax(i).YLabel.FontSize = FontSize;
    end
    
    ax(1).YLabel.String  = 'LE [deg]';
    ax(6).YLabel.String  = 'PE [deg]';
    
end

%% ------ verify_calibration -----------------------------------------------------------
if flags.do_verify_calibration
    data_majdak = data_majdak2010('Learn_M');
    data_majdak([1:5]) = [];

    %calibrations = amt_cache('get', 'barumerli2021_calibration',flags.cachemode);
    calibrations = amt_load('barumerli2021', 'barumerli2021_calibration.mat');
    calibrations = calibrations.cache.value;
    
    if size(calibrations.sigma, 1) ~= length(data_majdak)
        warning('sigma values not enough for the provided subejcts')
    end

    verify_calibration = amt_cache('get', 'barumerli2021_verify_calibration',flags.cachemode);

    % Preallocation
    if isempty(verify_calibration)
        verify_calibration = repmat(struct('err', struct([])),length(data_majdak), size(calibrations.combination, 1));
        num_calib = size(calibrations.combination,1);
        for s = 1:size(calibrations.sigma, 1) % subjects
            amt_disp('\n %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \n SUBJECT %s\n', data_majdak(s).id)

%             sofa = SOFAload(['dtf_' lower(data_baseline(s).id) '.sofa']);
%             sofa = SOFAload(fullfile(SOFAdbPath,'barumerli2021',['dtf_' lower(calibrations.name{s}) '.sofa']));
            sofa = SOFAload(fullfile(SOFAdbPath,'barumerli2021',['ARI_' upper(calibrations.name{s}) '_hrtf_M_dtf 256.sofa']));

            
            for c = 1:num_calib% feature space
                amt_disp(sprintf('\nCOMBINATION %i', c))
                amt_disp(sprintf('%s ', calibrations.combination{c,:}))

                [template_par, target_par] = barumerli2021_featureextraction(sofa, ...
                                        calibrations.combination{c,1}, ...
                                        calibrations.combination{c,2}, ...
                                        calibrations.combination{c,3});


                sigma_l = calibrations.sigma(s,c).values(1);
                sigma_l2 = calibrations.sigma(s,c).values(2);
                sigma_p = calibrations.sigma(s,c).values(3);
                sigma_m = calibrations.sigma(s,c).values(4);
                sigma_prior = calibrations.sigma(s,c).values(5);

                if calibrations.sigma(s,c).values(3) == 0
                    sigma_p = [];
                end

                m = barumerli2021('template', template_par, ...
                                    'target', target_par, ...
                                    'num_exp', 300, ...
                                    'sigma_lateral', sigma_l, ...
                                    'sigma_lateral_bis', sigma_l2, ...
                                    'sigma_polar', sigma_p, ...
                                    'sigma_sensorimotor', sigma_m,...
                                    'MAP',...
                                    'sigma_prior', sigma_prior);

                verify_calibration(s,c).err = barumerli2021_metrics(m,'middle_metrics');
            end
        end

        amt_cache('set', 'barumerli2021_verify_calibration',verify_calibration);
    end
    % metrics
    metric_calib = zeros(size(calibrations.combination, 1),3);
    for s = 1:size(calibrations.sigma, 1)
        amt_disp(sprintf('\nSUBJECT %i', s))
        mtx = data_majdak(s).mtx;
        real = barumerli2021_metrics(mtx,'middle_metrics');
        
        for c = 1:size(calibrations.combination, 2)
            % plot
            amt_disp(sprintf('\t\t\tFEATURE SPACE'))
            amt_disp(sprintf('%s ', calibrations.combination{c,:}))
            le_tau = abs(verify_calibration(s,c).err.rmsL - real.rmsL)/real.rmsL;
            pe_tau = abs(verify_calibration(s,c).err.rmsP - real.rmsP)/real.rmsP;
            qe_tau_rau = abs(rau(verify_calibration(s,c).err.querr, 1, 'PC') - rau(real.querr, 1, 'PC'))/rau(real.querr, 1, 'PC');
            amt_disp(['le_tau ', num2str(le_tau),' pe_tau ', num2str(pe_tau),' querr_tau ', num2str(qe_tau_rau)]);%.2f\n', , pe_tau, qe_tau_rau)

            metric_calib(c,1) = metric_calib(c,1) + verify_calibration(s,c).err.rmsL;
            metric_calib(c,2) = metric_calib(c,2) + verify_calibration(s,c).err.rmsP;
            metric_calib(c,3) = metric_calib(c,3) + verify_calibration(s,c).err.querr;
        end
    end
    
    metric_calib = metric_calib./size(calibrations.sigma, 1);
    amt_disp(sprintf('\n############\nAVERAGE OVER SUBJECTS'))
    for c = 1:size(calibrations.combination, 1)
        amt_disp(sprintf('\t\t\tFEATURE SPACE'))
        amt_disp(sprintf('%s ', calibrations.combination{c,:}))
        amt_disp(['le ', num2str(metric_calib(c,1)), ' pe ', num2str(metric_calib(c,2)), ' querr ', num2str(metric_calib(c,3))]);%.2f, pe %.2f, querr %.2f \n', metric_calib(c,1), metric_calib(c,2), metric_calib(c,3))
    end
end

%% ------ middlebrooks -------------------------------------
if flags.do_exp_middlebrooks1999
    data_majdak = data_majdak2010('Learn_M');
    data_majdak([1:5]) = [];

    exp_middlebrooks = [];

    if ~flags.do_redo
        exp_middlebrooks = amt_cache('get', ...
            'exp_middlebrooks1999',flags.cachemode);
    end
    
    calibrations = amt_load('barumerli2021', 'barumerli2021_calibration.mat');
    calibrations = calibrations.cache.value;
    if size(calibrations.sigma, 1) ~= length(calibrations.name)
        error('sigma values not enough for the provided subejcts')
    end
    
    sbj_num = length(calibrations.name);
    cal_num = size(calibrations.combination, 1);
    num_exp = 50;
    
    if flags.do_redo_fast
        num_exp = 5;
        exp_middlebrooks = [];
    end
    
    if flags.do_test
        num_exp = 1;
        cal_num = 1;
        sbj_num = 1;
    end
    
    if isempty(exp_middlebrooks) || flags.do_redo_fast
        % preprocess templates for each user
        amt_disp('Processing subjects'' templates');

        for s = 1:sbj_num
            amt_disp(['Pre-processing subject #' num2str(s)],'progress');
            %sofa = SOFAload(['dtf_' lower(calibrations.name{s}) '.sofa']);

            sofa = SOFAload(fullfile(SOFAdbPath,'barumerli2021',['ARI_' upper(calibrations.name{s}) '_hrtf_M_dtf 256.sofa']));
            horpolar_coords = barumerli2021_coordinates(sofa).return_positions('horizontal-polar');

            idx = find((horpolar_coords(:, 2) >= -60 ...
                            & horpolar_coords(:, 2) <= 79) ...
                            | (horpolar_coords(:, 2) >= 101 & ...
                                    horpolar_coords(:, 2) <= 240));

            for c = 1:cal_num
                [template(c,s), target(c,s)] = ...
                    barumerli2021_featureextraction(sofa, ...
                                            calibrations.combination{c,1}, ...
                                            calibrations.combination{c,2}, ...
                                            calibrations.combination{c,3});
            end
        end
        
        % preallocation for results
        amt_disp('Allocating memory for results');
        estimations = struct('m', []);
        estimations = repmat(estimations, cal_num, ...
            sbj_num, sbj_num); % all vs all
        
        for c = 1:cal_num
            amt_disp(sprintf('Combination #%i', c));
            for s = 1:sbj_num
                amt_disp(sprintf('\tSubject #%i', s));
                
                assert(length(calibrations.sigma(s,c).values) == 5, 'something is wrong with the calibration file')
                sigma_l = calibrations.sigma(s,c).values(1);
                sigma_l2 = calibrations.sigma(s,c).values(2);
                sigma_p = calibrations.sigma(s,c).values(3);
                sigma_m = calibrations.sigma(s,c).values(4);
                sigma_prior = calibrations.sigma(s,c).values(5);
                
                for j = 1:sbj_num
                    amt_disp(num2str(j));
                    estimations(c, s, j).m = ... 
                        barumerli2021('template', template(c, s), 'target', target(c, j), 'num_exp', num_exp, ...
                                            'sigma_lateral', sigma_l, ...
                                            'sigma_lateral_bis', sigma_l2, ...
                                            'sigma_polar', sigma_p,...
                                            'sigma_sensorimotor', sigma_m, ...
                                            'sigma_prior', sigma_prior);
                end
            end
        end

        % compute metrics
        for c = 1:size(estimations, 1)
            for i = 1:size(estimations, 2)
                for j = 1:size(estimations, 3)
                    metrics(c, i, j) = barumerli2021_metrics(estimations(c, i, j).m, 'middle_metrics'); 
                end
            end
        end
        
        exp_middlebrooks = metrics;
        
        if ~flags.do_redo_fast
            amt_cache('set','exp_middlebrooks1999',exp_middlebrooks);
        end
    end
    
    metrics_all = exp_middlebrooks;
    num_calib = size(metrics_all, 1);
        quants = [0,0.05,0.25,0.5,0.75,0.95,1];
    
    % iterate over calibrations
    for c=1:num_calib
        metrics = squeeze(metrics_all(c,:,:));
        
        % aggregate metrics
        ns = size(metrics,1);
        own = logical(eye(ns));
        other = not(own);

        % code similar to baumgartner2014 - fig9
        le_own(c,1).quantiles = quantile([metrics(own).rmsL], quants);
        lb_own(c,1).quantiles = quantile([metrics(own).accL], quants);
        qe_own(c,1).quantiles = quantile([metrics(own).querr], quants);
        pe_own(c,1).quantiles = quantile([metrics(own).rmsP], quants);
        pb_own(c,1).quantiles = quantile([metrics(own).accP], quants);
        le_own(c,1).mean = mean([metrics(own).rmsL]);
        lb_own(c,1).mean = mean([metrics(own).accL]);
        qe_own(c,1).mean = mean([metrics(own).querr]);
        pe_own(c,1).mean = mean([metrics(own).rmsP]);
        pb_own(c,1).mean = mean([metrics(own).accP]);

        le_other(c,1).quantiles = quantile([metrics(other).rmsL], quants);
        lb_other(c,1).quantiles = quantile([metrics(other).accL], quants);
        qe_other(c,1).quantiles = quantile([metrics(other).querr], quants);
        pe_other(c,1).quantiles = quantile([metrics(other).rmsP], quants);
        pb_other(c,1).quantiles = quantile([metrics(other).accP], quants);
        le_other(c,1).mean = mean([metrics(other).rmsL]);
        lb_other(c,1).mean = mean([metrics(other).accL]);
        qe_other(c,1).mean = mean([metrics(other).querr]);
        pe_other(c,1).mean = mean([metrics(other).rmsP]);
        pb_other(c,1).mean = mean([metrics(other).accP]);
    end
    
    % load reference data
    data_middle = data_middlebrooks1999;

    % load actual data and use the middlebrooks grid
    for i=1:size(data_majdak, 2)
        metric_majdak(i) = barumerli2021_metrics(data_majdak(i).mtx, 'middle_metrics');
    end
    
    clear data_majdak
    data_majdak.le.quantiles = quantile([metric_majdak.rmsL], quants);
    data_majdak.lb.quantiles = quantile([metric_majdak.accL], quants);
    data_majdak.qe.quantiles = quantile([metric_majdak.querr], quants);
    data_majdak.pe.quantiles = quantile([metric_majdak.rmsP], quants);
    data_majdak.pb.quantiles = quantile([metric_majdak.accP], quants);
    data_majdak.le.mean = mean([metric_majdak.rmsL]);
    data_majdak.lb.mean = mean([metric_majdak.accL]);
    data_majdak.qe.mean = mean([metric_majdak.querr]);
    data_majdak.pe.mean = mean([metric_majdak.rmsP]);
    data_majdak.pb.mean = mean([metric_majdak.accP]);


    % baumgartner data
    data_baum_temp = exp_baumgartner2014('fig9', 'no_plot');
    data_baum.qe_pool = data_baum_temp(1).qe;
    data_baum.pe_pool = data_baum_temp(1).pe;
    data_baum.pb_pool = data_baum_temp(1).pb;

    ns = size(data_baum.pe_pool,1);
    own = eye(ns) == 1;
    other = not(own);
    data_baum.pb_pool = abs(data_baum.pb_pool);
    data_baum.qe_own.quantiles = quantile(data_baum.qe_pool(own),quants);
    data_baum.pe_own.quantiles = quantile(data_baum.pe_pool(own),quants);
    data_baum.pb_own.quantiles = quantile(data_baum.pb_pool(own),quants);
    data_baum.qe_own.mean = mean(data_baum.qe_pool(own));
    data_baum.pe_own.mean = mean(data_baum.pe_pool(own));
    data_baum.pb_own.mean = mean(data_baum.pb_pool(own));

    data_baum.qe_other.quantiles = quantile(data_baum.qe_pool(other),quants);
    data_baum.pe_other.quantiles = quantile(data_baum.pe_pool(other),quants);
    data_baum.pb_other.quantiles = quantile(data_baum.pb_pool(other),quants);
    data_baum.qe_other.mean = mean(data_baum.qe_pool(other));
    data_baum.pe_other.mean = mean(data_baum.pe_pool(other));
    data_baum.pb_other.mean = mean(data_baum.pb_pool(other));
    
    % reijniers
    data_reij_temp = exp_reijniers2014('fig2_barumerli2020forum', 'no_plot');

    ns = size(data_reij_temp,1);
    own = logical(eye(ns));
    other = not(own);
    data_reij.le_own.quantiles = quantile([data_reij_temp(own).rmsL], quants);
    data_reij.lb_own.quantiles = quantile([data_reij_temp(own).accL], quants);
    data_reij.qe_own.quantiles = quantile([data_reij_temp(own).querr], quants);
    data_reij.pe_own.quantiles = quantile([data_reij_temp(own).rmsP], quants);
    data_reij.pb_own.quantiles = quantile([data_reij_temp(own).accP], quants);
    data_reij.le_own.mean = mean([data_reij_temp(own).rmsL]);
    data_reij.lb_own.mean = mean([data_reij_temp(own).accL]);
    data_reij.qe_own.mean = mean([data_reij_temp(own).querr]);
    data_reij.pe_own.mean = mean([data_reij_temp(own).rmsP]);
    data_reij.pb_own.mean = mean([data_reij_temp(own).accP]);

    data_reij.le_other.quantiles = quantile([data_reij_temp(other).rmsL], quants);
    data_reij.lb_other.quantiles = quantile([data_reij_temp(other).accL], quants);
    data_reij.qe_other.quantiles = quantile([data_reij_temp(other).querr], quants);
    data_reij.pe_other.quantiles = quantile([data_reij_temp(other).rmsP], quants);
    data_reij.pb_other.quantiles = quantile([data_reij_temp(other).accP], quants);
    data_reij.le_other.mean = mean([data_reij_temp(other).rmsL]);
    data_reij.lb_other.mean = mean([data_reij_temp(other).accL]);
    data_reij.qe_other.mean = mean([data_reij_temp(other).querr]);
    data_reij.pe_other.mean = mean([data_reij_temp(other).rmsP]);
    data_reij.pb_other.mean = mean([data_reij_temp(other).accP]);
    
    % plot
    if flags.do_plot
        calib_plot_order = [3,1,2]; %[lat, dtf, pge]
        % spacing
        dx = 0.11;
        % multiplier for horizontal shift
        middle_off = 3;
        majdak_off = 2;
        reij_off = 1;
        baum_off = 0;
       
        Marker = 's-';
        LineColor = [[0.9290, 0.6940, 0.1250]; ...
            [0.4940, 0.1840, 0.5560]; ...
            [0.4660, 0.6740, 0.1880]; ...
            [0.3010, 0.7450, 0.9330]; ...
            [0.6350, 0.0780, 0.1840]];
        data_middle.Marker = 'ko-';
        data_middle.LineColor = 'k';%[1 1 1]*0.3;
        
        data_majdak.Marker = 'b^-';
        data_majdak.LineColor = 'b';%[0 0 1]*0.3;
        
        data_baum.Marker = 'd-';
        data_baum.LineColor = [0.8500 0.3250 0.0980];
        
        data_reij.Marker = 'v-';
        data_reij.LineColor = [0 0.4470 0.7410];
        
        mFig = figure;
        mFig.Units = 'points';
        mFig.Position = [0, 0, 510, 300];
        tile_left = 0.02;%[left bottom width height]
        tile_width = 0.32;

        %% SUBPLOT 1
        sp = 1; 
        subplot(1, 3, sp);
        ax=gca;
        ax.OuterPosition(1) = tile_left;
        ax.OuterPosition(3) = tile_width;
        ax.Position(3) = 0.26;
        
        % reference
        local_middlebroxplot(ax, 1-middle_off*dx,data_middle.le_own, data_middle.Marker, kv.MarkerSize, data_middle.LineColor, data_middle.LineColor);
        local_middlebroxplot(ax, 2-middle_off*dx,data_middle.le_other, data_middle.Marker, kv.MarkerSize, data_middle.LineColor, data_middle.LineColor);
        
        % baseline
        local_middlebroxplot(ax, 1-majdak_off*dx,data_majdak.le, data_majdak.Marker, kv.MarkerSize, data_majdak.LineColor, data_majdak.LineColor);
        
        % simulation
        cdist = 1;
        for c=calib_plot_order
            local_middlebroxplot(ax, 1+cdist*dx,le_own(c,1), Marker, kv.MarkerSize, LineColor(c,:), 'w');
            local_middlebroxplot(ax, 2+cdist*dx,le_other(c,1), Marker, kv.MarkerSize, LineColor(c,:), 'w');
            cdist = cdist+1;
        end
        
        % reijniers2014
        local_middlebroxplot(ax, 1-reij_off*dx,data_reij.le_own, data_reij.Marker, kv.MarkerSize, data_reij.LineColor, 'w');
        local_middlebroxplot(ax, 2-reij_off*dx,data_reij.le_other, data_reij.Marker, kv.MarkerSize, data_reij.LineColor, 'w');
%         
        ylabel('Lateral RMS Error [deg]','FontSize',kv.FontSize)
        set(gca,'YLim',[0 45],'YTick', 0:10:40,'XLim',[0.5 2.5],...
          'XTick',1:2,'XTickLabel',{'Own' 'Other'},'FontSize',kv.FontSize,...
            'TickLength',2*get(gca,'TickLength'))
        
        %% SUBPLOT 2
        sp = sp +1;
        subplot(1, 3, sp);
        ax=gca;
        ax.OuterPosition(1) = tile_left*2 + tile_width;
        ax.OuterPosition(3) = tile_width;
        ax.Position(3) = 0.26;

        % reference
        local_middlebroxplot(ax, 1-middle_off*dx, data_middle.pe_own, data_middle.Marker, kv.MarkerSize, data_middle.LineColor, data_middle.LineColor);
        local_middlebroxplot(ax, 2-middle_off*dx, data_middle.pe_other, data_middle.Marker, kv.MarkerSize, data_middle.LineColor, data_middle.LineColor);
        
        % baseline
        local_middlebroxplot(ax, 1-majdak_off*dx, data_majdak.pe, data_majdak.Marker, kv.MarkerSize, data_majdak.LineColor, data_majdak.LineColor);

        % simulations
        cdist = 1;
        for c=calib_plot_order
            local_middlebroxplot(ax, 1+cdist*dx, pe_own(c,1), Marker, kv.MarkerSize, LineColor(c,:),'w');
            local_middlebroxplot(ax, 2+cdist*dx, pe_other(c,1), Marker, kv.MarkerSize, LineColor(c,:),'w');
            cdist = cdist + 1;
        end
        
        % reijniers2014
        local_middlebroxplot(ax, 1-reij_off*dx, data_reij.pe_own,data_reij.Marker, kv.MarkerSize, data_reij.LineColor,'w');
        local_middlebroxplot(ax, 2-reij_off*dx, data_reij.pe_other,data_reij.Marker, kv.MarkerSize, data_reij.LineColor,'w');

        % baumgartner2014
        local_middlebroxplot(ax, 1-baum_off*dx, data_baum.pe_own,data_baum.Marker, kv.MarkerSize, data_baum.LineColor,'w');
        local_middlebroxplot(ax, 2-baum_off*dx, data_baum.pe_other,data_baum.Marker, kv.MarkerSize, data_baum.LineColor,'w');
        
        ylabel('Local Polar RMS Error [deg]','FontSize',kv.FontSize)
        set(gca,'YLim',[0 65],'YTick', 0:10:60,'XLim',[0.5 2.5],...
          'XTick',1:2,'XTickLabel',{'Own' 'Other'},'FontSize',kv.FontSize,...
            'TickLength',2*get(gca,'TickLength'))

         %% SUBPLOT 3
        sp = sp +1;
        sp_ref = subplot(1, 3, sp);
        ax=gca;
        
        ax.OuterPosition(1) = tile_left*3 + tile_width*2;
        ax.OuterPosition(3) = tile_width;
                ax.Position(3) = 0.26;

        % reference
        middle = local_middlebroxplot(ax, 1-middle_off*dx,data_middle.qe_own, data_middle.Marker, kv.MarkerSize, data_middle.LineColor, data_middle.LineColor);
        local_middlebroxplot(ax, 2-middle_off*dx,data_middle.qe_other, data_middle.Marker, kv.MarkerSize, data_middle.LineColor, data_middle.LineColor);

        % baseline
        baseline = local_middlebroxplot(ax, 1-majdak_off*dx,data_majdak.qe, data_majdak.Marker, kv.MarkerSize, data_majdak.LineColor, data_majdak.LineColor);
        
        % simulations
        cdist = 1;
        for c=calib_plot_order
            baru(1,c) = local_middlebroxplot(ax, 1+cdist*dx,qe_own(c,1), Marker, kv.MarkerSize, LineColor(c,:),'w');
            local_middlebroxplot(ax, 2+cdist*dx,qe_other(c,1), Marker, kv.MarkerSize, LineColor(c,:),'w');
            cdist = cdist + 1;
        end
        
        % reijniers2014
        reij = local_middlebroxplot(ax, 1-reij_off*dx,data_reij.qe_own, data_reij.Marker,kv.MarkerSize, data_reij.LineColor,'w');
        local_middlebroxplot(ax, 2-reij_off*dx,data_reij.qe_other, data_reij.Marker,kv.MarkerSize, data_reij.LineColor,'w');

        % baumgartner2014
        baum = local_middlebroxplot(ax, 1-baum_off*dx,data_baum.qe_own, data_baum.Marker,kv.MarkerSize, data_baum.LineColor,'w');
        local_middlebroxplot(ax, 2-baum_off*dx,data_baum.qe_other, data_baum.Marker,kv.MarkerSize, data_baum.LineColor,'w');

        ylabel('Quadrant Errors [%]','FontSize',kv.FontSize)
        set(gca,'YLim',[-5 55],'YTick', 0:10:50,'XLim',[0.5 2.5],...
          'XTick',1:2,'XTickLabel',{'Own' 'Other'},'FontSize',kv.FontSize,...
            'TickLength',2*get(gca,'TickLength'))
        
        for c=calib_plot_order
            switch lower(calibrations.combination{c,2})
                case 'dtf'
                    labels{1,c} = '$T_{DTF}^\varphi$';
                case 'pge'
                    labels{1,c} = '$T_{PGE}^\varphi$';
            end
            switch lower(calibrations.combination{c,3})
%                 case 'reijniers'
%                 labels{1,c} = 'Reijniers2014';
                case 'interaural_difference'
                    labels{1,c} = '$T_{LAT}^\varphi$';
            end
        end

        leg = legend(sp_ref, [middle, baseline, reij, baum, baru(calib_plot_order)], horzcat({'middlebrooks1999'},{'majdak2010'},{'reijniers2014'}, {'baumgartner2014'}, labels(calib_plot_order)));
%         leg = legend(sp_ref, [middle, baseline, baum, baru(calib_plot_order)], horzcat({'Middlebrooks1999'},{'Majdak2010'}, {'Baumgartner2014'}, labels(calib_plot_order)));
        leg.FontSize = kv.FontSize - 2;
        leg.Units = 'centimeters';
        leg.Interpreter = 'latex';
        leg.Units = 'normalized';
        leg.Position = [0.0996373699690041 0.679833338340123 0.183774799022303 0.21349999499321];
    end
end

%% ------ exp_best -------------------------------------
if flags.do_exp_best2005
    exp_best = [];

    if ~flags.do_redo
        exp_best = amt_cache('get', ...
            'exp_best2005',flags.cachemode);
    end
    
    calibrations = amt_load('barumerli2021', 'barumerli2021_calibration.mat');
    calibrations = calibrations.cache.value;
    if size(calibrations.sigma, 1) ~= length(calibrations.name)
        error('sigma values not enough for the provided subejcts')
    end
    
    % Settings
    num_calib = size(calibrations.combination, 1);
    num_sbj = length(calibrations.name);
    num_exp = 50;
    num_sample = 260;
    num_exp_LP = num_exp;
    
    if flags.do_redo_fast
        amt_disp('Redo fast')
        exp_best = [];
        num_exp = 5;
    end
    
    if flags.do_test
        amt_disp('Test')
        exp_best = [];
        num_exp = 1;
        num_sample = 1;
        num_exp_LP = num_exp;
        num_calib = 1;
        num_sbj = 1;
    end
    
    if isempty(exp_best)
        % extract directions
        % retrieved from figure 3 and 4 of the original paper

        dir_hor = [-160,-150,-140,-130,-120,-90,-60,-50,-40,-30,-20,0,20,30,40,50,60,90,120,130,140,149,160,180];
        dir_hor(dir_hor < -90) = dir_hor(dir_hor < -90) + 360; % adjust to interval [-90, 270]
        dir_lat = [-90,-70,-60,-50,-45,-30,-20,-10,0,10,20,30,45,50,60,70,90];
        [x,y] = meshgrid(dir_lat,dir_hor);

        %sofa_obj = SOFAload(['dtf_' lower(calibrations.name{1}) '.sofa']);
        sofa = SOFAload(fullfile(SOFAdbPath,'barumerli2021',['ARI_' upper(calibrations.name{1}) '_hrtf_M_dtf 256.sofa']));

        coords = barumerli2021_coordinates(sofa_obj);
        coords.normalize_distance();
        
        coords_cart = coords.return_positions('cartesian');

        targ_coords = barumerli2021_coordinates([dirs, ones(size(dirs,1),1)], 'horizontal-polar');
        targ_cart = targ_coords.return_positions('cartesian');
        targ_cart(targ_cart(:, 3) < min(coords_cart(:,3)), :) = []; % remove lowest positions
        targ_coords = Coordinates(targ_cart, 'cartesian');

        targ_idx = coords.find_positions(targ_coords);
        targ_idx = unique(targ_idx);
        
        coords_sph = coords.return_positions('spherical');
        
        targ_az = coords_sph(targ_idx,1);
        targ_el = coords_sph(targ_idx,2);
        targ_cart = targ_coords.return_positions('cartesian');
%         figure;
%         scatter3(coords.pos(:,1), coords.pos(:,2), coords.pos(:,3), 10, 'k', 'filled')
%         hold on
%         scatter3(targ_cart(:,1), targ_cart(:,2), targ_cart(:,3), 20, 'b')
%         scatter3(coords.pos(targ_idx,1), coords.pos(targ_idx,2), coords.pos(targ_idx,3), 50, 'r')
%         legend({'sofa', 'targ', 'chosen'})
%         
        
        % Load Data
        % Speech Samples from Harvard Word list
        speechsample = amt_cache('get','../experiments%2Fexp_baumgartner2014.m/best2005speechSamples');
        speech_fs = 48e3;
        speechsample(num_sample+1:end) = [];
        samples_num = length(speechsample);
        for s = 1:samples_num
            speechsample{s} = speechsample{s} ./ max(speechsample{s});
        end
        
        % prepare directions, stimuli combinations
        if flags.do_test 
            comb_stim = 1;
        else
            repetitions = 5; % see paper sect II.D
            num_stim_rep = 50;
            comb_stim = zeros(length(targ_idx), repetitions);
            for r = 1:repetitions
                idx = repmat(transpose(((r-1)*num_stim_rep+1):(r*num_stim_rep)), ceil(length(targ_idx)/num_stim_rep), 1);
                comb_stim(:,r) = idx(1:size(comb_stim,1),:);
            end
            % randomize everything
            comb_stim = comb_stim(randperm(size(comb_stim, 1)), :);
        end
        
        % FIR Low-pass filters at 8kHz
        % Brick-wall (aka sinc-filter): fir1(200,1/3) -> -60 dB
        x=amt_load('baumgartner2014','highfreqatten_filters.mat');
        LP{1} = [1 zeros(1,100)];
        LP{2} = x.fir20db;
        LP{3} = x.fir40db;
        LP{4} = x.fir60db;

        % Model Data
        estimations = struct('m', []);
        est_expflat = repmat(estimations, num_calib, num_sbj, 1); 

        est_expLP = repmat(estimations, num_calib, ...
            num_sbj,length(LP)); 
        
        % Simulations
        amt_disp('Processing subjects HRIR');
        
        data_majdak = data_majdak2010('Learn_M');
        data_majdak([1:5]) = [];
        data_actual = [];
        
        for s = 1:num_sbj
            amt_disp(['Pre-processing subject #' num2str(s)],'progress');
            % load data for prior
            for s_i = 1:length(data_majdak)
                if strcmpi(data_majdak(s_i).id, calibrations.name{s})
                    data_actual(s).m = data_majdak(s_i).mtx;
                end
            end
            
            %sofa = SOFAload(['dtf_' lower(calibrations.name{s}) '.sofa']);
            sofa = SOFAload(fullfile(SOFAdbPath,'barumerli2021',['ARI_' upper(calibrations.name{s}) '_hrtf_M_dtf 256.sofa']));
            for c = 1:num_calib
                amt_disp(sprintf('\tCalibration #%i', c),'progress');

                % noise stimulus
                [template, target_flat] = ...
                    barumerli2021_featureextraction(sofa, 'targ_az', targ_az, ...
                                            'targ_el', targ_el, ...
                                            calibrations.combination{c,1}, ...
                                            calibrations.combination{c,2}, ...
                                            calibrations.combination{c,3});
                
                assert(length(calibrations.sigma(s,c).values) == 5, 'something is wrong with the calibration file')
                sigma_l = calibrations.sigma(s,c).values(1);
		sigma_l2 = calibrations.sigma(s,c).values(2);
		sigma_p = calibrations.sigma(s,c).values(3);
		sigma_m = calibrations.sigma(s,c).values(4);
		sigma_prior = calibrations.sigma(s,c).values(5);
                
                % flat spectrum estimations
                est_expflat(c, s).m = ... 
                    barumerli2021('template', template, ...
                                    'target', target_flat, ...
                                    'num_exp', num_exp, ...
                                    'sigma_lateral', sigma_l, ...
                                    'sigma_lateral_bis', sigma_l2, ...
                                    'sigma_polar', sigma_p, ...
                                    'sigma_sensorimotor', sigma_m, ...
                                    'sigma_prior', sigma_prior);
    
                for f = 1:length(LP)
                    amt_disp(sprintf('\tFilter %i\n', f))
                    estimations = [];% zeros(size(comb_stim,1)*size(comb_stim,2)*num_exp_LP, 8);

                    for sa = 1:size(comb_stim, 2)
                        for d = 1:size(comb_stim,1)
                            %fprintf('%i-%i-%i\n', f, sa, d)
                            amt_disp([num2str(f), '-' , num2str(sa), '-', num2str(d)]);%'%i-%i-%i\n', f, sa, d)
                            
                            stim = filter(LP{f},1,speechsample{comb_stim(d, sa)});

                            trgt = barumerli2021_featureextraction(sofa, 'target', ...
                                                'targ_az', targ_az(d), ...
                                                'targ_el', targ_el(d), ...
                                                'source', 'source_ir', stim, ...
                                                'source_fs', speech_fs, ...
                                                calibrations.combination{c,1}, ...
                                                calibrations.combination{c,2}, ...
                                                calibrations.combination{c,3});

                            m_LP = barumerli2021('template', template, 'target', trgt, 'num_exp', num_exp_LP, ...
                                                            'sigma_lateral', sigma_l, ...
                                                            'sigma_lateral_bis', sigma_l2, ...
                                                            'sigma_polar', sigma_p, ...
                                                            'sigma_sensorimotor', sigma_m, ...
                                                            'sigma_prior', sigma_prior);
                            
                            estimations = [estimations; m_LP];                          
                        end
                    end
                    est_expLP(c, s, f).m = estimations;
                end
            end
        end

        % metrics
        % allocate memory for results
        % aggregate over different lateral angles
        ale_expflat = zeros(num_calib, num_sbj,1);
        ale_expLP = zeros(num_calib, num_sbj, length(LP));
        ape_expflat = zeros(num_calib, num_sbj,1);
        ape_expLP = zeros(num_calib, num_sbj, length(LP));
        qe_expflat = zeros(num_calib, num_sbj,1);
        qe_expLP = zeros(num_calib, num_sbj, length(LP));

        for c = 1:num_calib
            for i = 1:num_sbj
                amt_disp(['Computing metrics #' num2str(i)],'progress');

                % flat spectrum estimations
                m_flat = est_expflat(c, i, 1).m;
                ale_expflat(c, i, 1) = localizationerror(m_flat, 'accabsL');
                m_flat = m_flat((abs(m_flat(:,7))<=30),:);
                ape_expflat(c, i, 1) = localizationerror(m_flat, 'accabsP');
                qe_expflat(c, i, 1)  = localizationerror(m_flat, 'querr');

                for f = 1:length(LP)
                    m_LP = est_expLP(c, i, f).m;
                    ale_expLP(c, i, f) = localizationerror(m_LP, 'accabsL');
                    m_LP = m_LP((abs(m_LP(:,7))<=30),:);
                    ape_expLP(c, i, f) = localizationerror(m_LP, 'accabsP');
                    qe_expLP(c, i, f)  = localizationerror(m_LP, 'querr');
                end
            end
        end 
        % save cache
        exp_best.ale_expflat = ale_expflat;
        exp_best.ale_expLP   = ale_expLP;
        exp_best.ape_expflat = ape_expflat;
        exp_best.ape_expLP   = ape_expLP;        
        exp_best.qe_expflat  = qe_expflat;
        exp_best.qe_expLP    = qe_expLP;
        
        if ~flags.do_redo_fast
            amt_cache('set','exp_best2005', exp_best);
        end
    end
    
    varargout{1} = [exp_best];
    
    for c = 1:num_calib
        ale_expflat = exp_best.ale_expflat(c, :);
        ale_expLP = squeeze(exp_best.ale_expLP(c, :, :, :));
        ape_expflat = exp_best.ape_expflat(c, :);
        ape_expLP = squeeze(exp_best.ape_expLP(c, :, :, :));        
        qe_expflat = exp_best.qe_expflat(c, :);
        qe_expLP = squeeze(exp_best.qe_expLP(c, :, :, :));

        % Pool Samples
        ale_pooled = mean(ale_expLP,3);
        ape_pooled = mean(ape_expLP,3);
        qe_pooled = mean(qe_expLP,3);

        % Confidence Intervals or standard errors
        tquant_speech = 1;
        tquant_noise = 1;
        % barumerli model 
        df_speech = size(ale_expLP,1)-1;
        seale_speech = std(ale_pooled,0,1)*tquant_speech/(df_speech+1);
        df_noise = size(ale_expflat,1)-1;
        seale_noise = std(ale_expflat)*tquant_noise/(df_noise+1);
        seale = [seale_noise, seale_speech];

        df_speech = size(ape_expLP,1)-1;
        seape_speech = std(ape_pooled,0,1)*tquant_speech/(df_speech+1);
        df_noise = size(ape_expflat,1)-1;
        seape_noise = std(ape_expflat)*tquant_noise/(df_noise+1);
        seape = [seape_noise, seape_speech];

        % averages
        % baru2021
        ale(:,c) = mean([ale_expflat', ale_pooled],1,'omitnan');
        ape(:,c) = mean([ape_expflat', ape_pooled],1,'omitnan');
        qe(:,c) = mean([qe_expflat', qe_pooled],1, 'omitnan');
    end
    
    % load real and baum2014 data
    data = data_best2005;
    
    % reijniers2014
    data_reij = exp_reijniers2014('fig3_barumerli2020forum', 'no_plot');
    data_reij.ale_pooled = mean(data_reij.ale_expLP,3);
    data_reij.ape_pooled = mean(data_reij.ape_expLP,3);
    data_reij.qe_pooled = mean(data_reij.qe_expLP,3);

    df_speech = size(data_reij.ale_expLP,1)-1;
    seale_speech = std(data_reij.ale_pooled,0,1)*tquant_speech/(df_speech+1);
    df_noise = size(data_reij.ale_expflat,1)-1;
    seale_noise = std(data_reij.ale_expflat)*tquant_noise/(df_noise+1);
    data_reij.seale = [seale_noise, seale_speech];

    df_speech = size(data_reij.ape_expLP,1)-1;
    seape_speech = std(data_reij.ape_pooled,0,1)*tquant_speech/(df_speech+1);
    df_noise = size(data_reij.ape_expflat,1)-1;
    seape_noise = std(data_reij.ape_expflat)*tquant_noise/(df_noise+1);
    data_reij.seape = [seape_noise, seape_speech];

    data_reij.ale = mean([data_reij.ale_expflat, data_reij.ale_pooled],1);
    data_reij.ape = mean([data_reij.ape_expflat, data_reij.ape_pooled],1);
    data_reij.qe = mean([data_reij.qe_expflat, data_reij.qe_pooled],1);

    % baumgartner2014 model
    [baum_temp, ~] = exp_baumgartner2014('fig11', 'no_plot'); 
    % num_sbj, num_filters, num_samples
    data_baum.ape = zeros([size(baum_temp{1}, 2), size(ape_expLP, 2), size(ape_expLP, 3)]);
    data_baum.qe = zeros([size(baum_temp{1}, 2), size(qe_expLP, 2), size(qe_expLP, 3)]);
    for i = 1:size(ape_expLP,3)
        data_baum.ape(:,:,i) = transpose(baum_temp{1}(:,:,i));
        data_baum.qe(:,:,i) = transpose(baum_temp{2}(:,:,i));
    end
    data_baum.ape_expflat = baum_temp{3};
    data_baum.qe_expflat = baum_temp{4};
    data_baum.ape_pooled = mean(data_baum.ape,3);
    data_baum.qe_pooled = mean(data_baum.qe,3);
    df_speech = size(data_baum.ape,1)-1;
    seape_speech = std(data_baum.ape_pooled,0,1)*tquant_speech/(df_speech+1);
    df_noise = size(data_baum.ape_expflat,1)-1;
    seape_noise = std(data_baum.ape_expflat)*tquant_noise/(df_noise+1);
    data_baum.seape = [seape_noise, seape_speech];
    

    % baum2014
    data_baum.ape = mean([data_baum.ape_expflat', data_baum.ape_pooled],1);
    data_baum.qe = mean([data_baum.qe_expflat', data_baum.qe_pooled],1);
    
    if flags.do_plot
        num_calib = size(calibrations.combination, 1);
        calib_plot_order = [3,1,2];
        
        MarkerSize = kv.MarkerSize -1;
        FontSize = kv.FontSize;
        LineColor = [[0.9290, 0.6940, 0.1250]; ...
                    [0.4940, 0.1840, 0.5560]; ...
                    [0.4660, 0.6740, 0.1880]; ...
                    [0.3010, 0.7450, 0.9330]; ...
                    [0.6350, 0.0780, 0.1840]];
        Marker = 's-';
        data.Marker = 'ko-';
        data.LineColor = [0 0 0]; 
        
        data_reij.LineColor = [0 0.4470 0.7410];
        data_reij.Marker = 'v--';
        
        data_baum.Marker = 'd--';
        data_baum.LineColor = [0.8500, 0.3250, 0.0980];

        mFig = figure;
        mFig.Units = 'points';

        mFig.Position = [1e3,1e3,250,370];
        tile_position = [-0.025, 0, 1.1, 0.27];
        % sgtitle(sprintf('%s ', calibrations.combination{c,:}), 'Interpreter', 'none')

%         % compute delta
%         data.ale = data.ale - data.ale(1,1);
%         ale = ale - ale(1,:);
%         data.ape = data.ape - data.ape(1,1);
%         ape = ape - ape(1,:);
%         data_baum.ape = data_baum.ape - data_baum.ape(1,1);
%         data.qe([1 2 5]) = data.qe([1 2 5]) - data.qe(1,1);
%         qe = qe - qe(1,:);
%         data_baum.qe = data_baum.qe-data_baum.qe(1,1);
        
        dx = 0;
        xticks = 0:size(ale_pooled,2);
        sp1 = subplot(3,1,1);

        plot(xticks-dx,data.ale, data.Marker,'Color', data.LineColor, 'MarkerFaceColor',data.LineColor,'MarkerSize',MarkerSize)
        hold on
        for c = calib_plot_order
            plot(xticks-dx,ale(:,c), Marker,'Color', LineColor(c,:), 'MarkerFaceColor', 'w','MarkerSize',MarkerSize)
        end
        plot(xticks,data_reij.ale, data_reij.Marker,'Color', data_reij.LineColor, 'MarkerFaceColor','w','MarkerSize',MarkerSize)
%         ylabel('\Delta LE [deg]','FontSize',FontSize)
%         set(gca,'XLim',[-0.5 4.5],'YLim',[0 20],'YMinorTick','on')
        ylabel('Lat. acc. [deg]','FontSize',FontSize)
        set(gca,'XLim',[-0.5 4.5],'YLim',[0 30],'YMinorTick','on')
        set(gca,'XTick',xticks,'XTickLabel',[],'FontSize',FontSize)


        sp2 = subplot(3,1,2);
        
        plot(xticks,data.ape,data.Marker,'Color', data.LineColor, 'MarkerFaceColor',data.LineColor,'MarkerSize',MarkerSize)
        hold on
        for c=calib_plot_order
            plot(xticks-dx,ape(:,c), Marker,'Color', LineColor(c,:), 'MarkerFaceColor','w','MarkerSize',MarkerSize)
        end
        plot(xticks,data_reij.ape, data_reij.Marker,'Color', data_reij.LineColor, 'MarkerFaceColor', 'w','MarkerSize',MarkerSize)
        plot(xticks,data_baum.ape, data_baum.Marker,'Color', data_baum.LineColor, 'MarkerFaceColor','w','MarkerSize',MarkerSize)
%         ylabel('\Delta PE [deg]','FontSize',FontSize)
%         set(gca,'XLim',[-0.5 4.5],'YLim',[0 50],'YMinorTick','on')
        ylabel('Pol. Acc. [deg]','FontSize',FontSize)
        set(gca,'XLim',[-0.5 4.5],'YLim',[0 100],'YTick',10:20:100,'YMinorTick','on')
        set(gca,'XTick',xticks,'XTickLabel',[],'FontSize',FontSize)

        sp3 = subplot(3,1,3);
%         best = plot(xticks([1 2 5]),data.qe([1 2 5]), data.Marker,'Color', data.LineColor, 'MarkerFaceColor',data.LineColor,'MarkerSize',MarkerSize);
        best = plot(xticks([1 2 5]),data.qe([1 2 5]), data.Marker,'Color', data.LineColor, 'MarkerFaceColor',data.LineColor,'MarkerSize',MarkerSize);
        hold on
        for c=calib_plot_order
            baru(1,c) = plot(xticks-dx,qe(:,c), Marker,'Color', LineColor(c,:), 'MarkerFaceColor','w','MarkerSize',MarkerSize);
        end
        reij = plot(xticks-dx,data_reij.qe, data_reij.Marker,'Color', data_reij.LineColor, 'MarkerFaceColor', 'w','MarkerSize',MarkerSize);
        baum = plot(xticks,data_baum.qe, data_baum.Marker,'Color', data_baum.LineColor, 'MarkerFaceColor','w','MarkerSize',MarkerSize);
%         ylabel('\Delta QE [%]','FontSize',FontSize)
%         set(gca,'XLim',[-0.5 4.5],'YLim',[0 50],'YMinorTick','on')
        ylabel('QE [%]','FontSize',FontSize)
        set(gca,'XLim',[-0.5 4.5],'YLim',[0 70],'YTick',10:20:70,'YMinorTick','on')
        set(gca,'XTick',xticks,'XTickLabel',data.meta,'FontSize',FontSize)
                
        sp3.OuterPosition = tile_position;
        sp3.OuterPosition(2) = 0.01;
        sp3.Position(4) = 0.24;
    
        sp2.OuterPosition = tile_position;
        sp2.OuterPosition(2) = sum(sp3.OuterPosition([2,4]));
        sp2.Position(4) = 0.24;
        
        sp1.OuterPosition = tile_position;
        sp1.OuterPosition(2) = sum(sp2.OuterPosition([2,4]));
        sp1.Position(4) = 0.24;
        
        for c=calib_plot_order
            switch lower(calibrations.combination{c,2})
                case 'dtf'
                    labels{1,c} = '$T_{DTF}^\varphi$';
                case 'pge'
                    labels{1,c} = '$T_{PGE}^\varphi$';
            end
            if strcmpi(calibrations.combination{c,3}, 'interaural_difference')
                labels{1,c} = '$T_{LAT}^\varphi$';
            end
        end
        
        leg = legend([best, reij, baum, baru(calib_plot_order)], horzcat({'best2005'}, 'reijniers2014', {'baumgartner2014'}, labels(calib_plot_order)));
        leg.FontSize = FontSize - 2;
        leg.Units = 'normalized';
        leg.Interpreter = 'latex';
        leg.Position= [0.367949860921234 0.828427722881692 0.359055268858493 0.167246848265153];
    end
end

%% ------ exp_macpherson -------------------------------------
if flags.do_exp_macpherson2003
    exp_macpherson = [];

    if ~flags.do_redo
        exp_macpherson = amt_cache('get', ...
            'exp_macpherson2003',flags.cachemode);
    end
    
    calibrations = amt_load('barumerli2021', 'barumerli2021_calibration.mat');
    calibrations = calibrations.cache.value;
    if size(calibrations.sigma, 1) ~= length(calibrations.name)
        error('sigma values not enough for the provided subejcts')
    end
    
    % Settings
    num_exp = 50;
    num_sbj = length(calibrations.name);
    num_calib = size(calibrations.combination, 1);
    
    if flags.do_redo_fast
        exp_macpherson = [];
        num_exp = 10;
    end
    
    if flags.do_test
        exp_macpherson = [];
        num_exp = 1;
        num_calib = 1;
        num_sbj = 1; 
    end
    
    if isempty(exp_macpherson)
        %sofa = SOFAload(['dtf_' lower(calibrations.name{1}) '.sofa']);
        sofa = SOFAload(fullfile(SOFAdbPath,'barumerli2021',['ARI_' upper(calibrations.name{1}) '_hrtf_M_dtf 256.sofa']));
        % generate stimulus
        % copyed from exp_baumgartner2014/do_fig10
        density = [0.25, 0.5, 0.75, 1, 1.5, 2, 3, 4, 6, 8]; % ripples/oct
        depth =   10:10:40;        % ripple depth (peak-to-trough) in dB
        
        % 250-ms bursts, 20-ms raised-cosine fade in/out, flat from 0.6-16kHz
        fs = sofa.Data.SamplingRate;
        flow = 1e3;   % lower corner frequency of ripple modification in Hz
        fhigh = 16e3; % upper corner frequency of ripple modification in Hz
        Nf = 2^10;    % # Frequency bins

        f = 0:fs/2/Nf:fs/2;	% frequency bins
        id600 = find(f<=600,1,'last'); % index of 600 Hz (lower corner frequency of stimulus energy)
        idlow = find(f<=flow,1,'last'); % index of flow (ripples)
        idhigh = find(f>=fhigh,1,'first');  % index of fhigh (ripples)
        N600low = idlow - id600 +1;   % # bins without ripple modification
        Nlowhigh = idhigh - idlow +1; % # bins with ripple modification     % 
        O = log2(f(idlow:idhigh)/1e3);   % freq. trafo. to achieve equal ripple density in log. freq. scale

        % Raised-cosine '(i.e., cos^2)' ramp 1/8 octave wide
        fup = f(idlow)*2^(1/8);       % upper corner frequency of ramp upwards 
        idup = find(f<=fup,1,'last');
        Nup = idup-idlow+1;
        rampup = cos(-pi/2:pi/2/(Nup-1):0).^2;
        fdown = f(idhigh)*2^(-1/8);  % lower corner frequency of ramp downwards
        iddown = find(f>=fdown,1,'first');
        Ndown = idhigh-iddown+1;
        rampdown = cos(0:pi/2/(Ndown-1):pi/2).^2;
        ramp = [rampup ones(1,Nlowhigh-Nup-Ndown) rampdown];
        ramp = [-inf*ones(1,id600-1) zeros(1,N600low) ramp -inf*ones(1,Nf - idhigh)];

        % Ripples of Experiment I
        Sexp1 = zeros(Nf+1,length(density),2);  % 3rd dim: 1:0-phase 2:pi-phase
        Sexp1(idlow:idhigh,:,1) = (40/2* sin(2*pi*density'*O+ 0))';  % depth: 40dB, 0-phase
        Sexp1(idlow:idhigh,:,2) = (40/2* sin(2*pi*density'*O+pi))';  % depth: 40dB, pi-phase
        Sexp1 = repmat(ramp',[1,length(density),2]) .* Sexp1;
        Sexp1 = [Sexp1;Sexp1(Nf:-1:2,:,:)];
        Sexp1(isnan(Sexp1)) = -100;
        sexp1 = ifftreal(10.^(Sexp1/20),2*Nf);
        sexp1 = circshift(sexp1,Nf);  % IR corresponding to ripple modification
        sexp1 = squeeze(sexp1(:,:,1));
        % Ripples of Experiment II
        Sexp2 = zeros(Nf+1,length(depth),2);  % 3rd dim: 1:0-phase 2:pi-phase
        Sexp2(idlow:idhigh,:,1) = (depth(:)/2*sin(2*pi*1*O+ 0))';  % density: 1 ripple/oct, 0-phase
        Sexp2(idlow:idhigh,:,2) = (depth(:)/2*sin(2*pi*1*O+pi))';  % density: 1 ripple/oct, pi-phase
        Sexp2 = repmat(ramp',[1,length(depth),2]) .* Sexp2;
        Sexp2 = [Sexp2;Sexp2(Nf-1:-1:2,:,:)];
        Sexp2(isnan(Sexp2)) = -100;
        sexp2 = ifftreal(10.^(Sexp2/20),2*Nf);
        sexp2 = circshift(sexp2,Nf);  % IR corresponding to ripple modification
        sexp2 = squeeze(sexp2(:,:,1));

        if flags.do_test
             density(2:end) = []; % ripples/oct
            depth(2:end) = []; 
        end
        % preprocess templates for each user
        for i = 1:num_sbj
            amt_disp(['Processing subject ', num2str(i)]);

            for c = 1:num_calib
                amt_disp(['Pre-processing calibration #' num2str(c)],'progress');
    
                %sofa = SOFAload(['dtf_' lower(calibrations.name{i}) '.sofa']);
                sofa = SOFAload(fullfile(SOFAdbPath,'barumerli2021',['ARI_' upper(calibrations.name{i}) '_hrtf_M_dtf 256.sofa']));
                % extract directions
                % filter targets' coordinates
                % convert from spherical to horizontal-polar coordinates
                horpolar_coords = barumerli2021_coordinates(sofa).return_positions('horizontal-polar');
  
                % polar in [60, 120]
                % lateral = 0
                idx = find(((horpolar_coords(:, 2) >= -60 ...
                                & horpolar_coords(:, 2) <= 60) ...
                                | (horpolar_coords(:, 2) >= 120 & ...
                                        horpolar_coords(:, 2) <= 240)) ...
                                & (horpolar_coords(:, 1) <= 30 & horpolar_coords(:, 1) >= -30));

                amt_disp(['Pre-processing subject #' num2str(i)],'progress');
                [template(c, i), target_flat(c, i)] = ...
                    barumerli2021_featureextraction(sofa, ...
                        'targ_az', sofa.SourcePosition(idx, 1), ...
                        'targ_el', sofa.SourcePosition(idx, 2), ...
                        calibrations.combination{c,1}, ...
                        calibrations.combination{c,2}, ...
                        calibrations.combination{c,3});

                amt_disp('Densities conditions');
                for j = 1:length(density)
                    target_exp1(c, i, j) = ...
                        barumerli2021_featureextraction(sofa, 'target', 'source', ...
                        'source_ir', squeeze(sexp1(:, j)), 'source_fs', fs, ...
                        'targ_az', sofa.SourcePosition(idx, 1), ...
                        'targ_el', sofa.SourcePosition(idx, 2), ...
                                            calibrations.combination{c,1}, ...
                                            calibrations.combination{c,2}, ...
                                            calibrations.combination{c,3});
                end

                amt_disp('Depth conditions');
                for j = 1:length(depth)
                    target_exp2(c, i, j) = ...
                        barumerli2021_featureextraction(sofa, 'target', 'source', ...
                        'source_ir', squeeze(sexp2(:, j)), 'source_fs', fs, ...
                        'targ_az', sofa.SourcePosition(idx, 1), ...
                        'targ_el', sofa.SourcePosition(idx, 2), ...
                        calibrations.combination{c,1}, ...
                        calibrations.combination{c,2}, ...
                        calibrations.combination{c,3});
                end
            end
        end
        
        % preallocation for results
        amt_disp('Allocating memory for results');
        estimations = struct('m',[]);
        est_expflat = repmat(estimations, num_calib, num_sbj); 
        est_exp1 = repmat(estimations, num_calib, ...
            num_sbj,length(density)); 
        est_exp2 = repmat(estimations, num_calib, ...
            num_sbj,length(depth)); 

        % data for prior computation
        data_majdak = data_majdak2010('Learn_M');
        data_majdak([1:5]) = [];
        
        % simulations
        for i = 1:num_sbj
            for c = 1:num_calib
                amt_disp(sprintf('\tCalibration #%i', c),'progress');

	        assert(length(calibrations.sigma(i,c).values) == 5, 'something is wrong with the calibration file')
                sigma_l = calibrations.sigma(i,c).values(1);
                sigma_l2 = calibrations.sigma(i,c).values(2);
                sigma_p = calibrations.sigma(i,c).values(3);
                sigma_m = calibrations.sigma(i,c).values(4);
                sigma_prior = calibrations.sigma(i,c).values(5);
                
                amt_disp(sprintf('\tSubject #%i', i),'progress');
                % flat spectrum estimations
                est_expflat(c, i, 1).m = barumerli2021('template', template(c, i), 'target', target_flat(c, i), ...
                                    'num_exp', num_exp, ...
                                    'sigma_lateral', sigma_l, ...
                                    'sigma_lateral_bis', sigma_l2, ...
                                    'sigma_polar', sigma_p,...
                                    'sigma_sensorimotor', sigma_m, ...
                                    'sigma_prior', sigma_prior);

                % rippled estimations
                for j = 1:length(density)
                    est_exp1(c, i, j).m = barumerli2021('template', template(c, i), ...
                                    'target', target_exp1(c, i, j), ...
                                    'num_exp', num_exp, ...
                                    'sigma_lateral', sigma_l, ...
                                    'sigma_lateral_bis', sigma_l2, ...
                                    'sigma_polar', sigma_p,...
                                    'sigma_sensorimotor', sigma_m, ...
                                    'sigma_prior', sigma_prior);
                end

                for j =1:length(depth)
                    est_exp2(c, i, j).m = barumerli2021('template', template(c, i), ...
                                    'target', target_exp2(c, i, j), ...
                                    'num_exp', num_exp, ...
                                    'sigma_lateral', sigma_l, ...
                                    'sigma_lateral_bis', sigma_l2, ...
                                    'sigma_polar', sigma_p,...
                                    'sigma_sensorimotor', sigma_m, ...
                                    'sigma_prior', sigma_prior);
                end
            end        
        end

        % metrics
        % allocate memory for results
        % aggregate over different lateral angles
        pe_exp1 = zeros(num_calib, num_sbj, length(density));
        pe_exp2 = zeros(num_calib, num_sbj, length(depth));
        pe_flat = zeros(num_calib, num_sbj, 1);

        for c = 1:num_calib
            for i = 1:num_sbj
                % compute iterative regression (see Macpherson paper and localizationerror.m)   
                [f,r] = localizationerror(est_expflat(c,i).m, 'sirpMacpherson2000');
                
                pe_flat(c,i) = localizationerror(est_expflat(c,i).m, f, r, 'perMacpherson2003');

                for j = 1:length(density)
                    pe_exp1(c, i, j) = localizationerror(est_exp1(c, i, j).m, f, r, 'perMacpherson2003');
                end

                for j = 1:length(depth)
                    pe_exp2(c, i, j) = localizationerror(est_exp2(c ,i, j).m, f, r, 'perMacpherson2003');
                end
            end
        end
        
        % save cache
        exp_macpherson.pe_flat = pe_flat;
        exp_macpherson.pe_exp1 = pe_exp1;
        exp_macpherson.pe_exp2 = pe_exp2;
        
        if ~flags.do_redo_fast
            amt_cache('set','exp_macpherson2003', exp_macpherson);
        end
    end
    
    % Original data:
    data = data_macpherson2003;

    % Reijniers2014's data
    data_reij = exp_reijniers2014('fig4_barumerli2020forum', 'no_plot');

    % Baumgartner2014's data
    % varargout{1} = {pe_exp1,pe_exp2,pe_flat,noDCN};
    data_baum_temp = exp_baumgartner2014('fig10', 'no_plot');
    data_baum.pe_exp1 = data_baum_temp{1,1};
    data_baum.pe_exp2 = data_baum_temp{1,2};
    data_baum.pe_flat = data_baum_temp{1,3};

    % Phase condition handling
    % average across the phase condition
    % real data
    data.pe_exp1 = mean(data.pe_exp1,3);
    data.pe_exp2 = mean(data.pe_exp2,3);
    % baumgartner data 
    data_baum.pe_exp1 = mean(data_baum.pe_exp1,3);
    data_baum.pe_exp2 = mean(data_baum.pe_exp2,3);
    idphase = 1;

    % Increase
    % reijniers2014
    data_reij.pe_exp1 = data_reij.pe_exp1 - repmat(data_reij.pe_flat(:), 1, size(data_reij.pe_exp1, 2));
    data_reij.pe_exp2 = data_reij.pe_exp2 - repmat(data_reij.pe_flat(:), 1, size(data_reij.pe_exp2, 2));
    % baumgartner data
    data_baum.pe_exp1 = data_baum.pe_exp1 - repmat(data_baum.pe_flat(:),1,size(data_baum.pe_exp1,2));
    data_baum.pe_exp2 = data_baum.pe_exp2 - repmat(data_baum.pe_flat(:),1,size(data_baum.pe_exp2,2));

    % Statistics
    % real data
    data.quart_pe_flat = quantile(data.pe_flat,[.25 .50 .75]);
    data.quart_pe_exp1 = quantile(data.pe_exp1,[.25 .50 .75]);
    data.quart_pe_exp2 = quantile(data.pe_exp2,[.25 .50 .75]);

    % reijniers data
    data_reij.quart_pe_flat = quantile(data_reij.pe_flat,[.25 .50 .75]);
    data_reij.quart_pe_exp1 = quantile(data_reij.pe_exp1,[.25 .50 .75]);
    data_reij.quart_pe_exp2 = quantile(data_reij.pe_exp2,[.25 .50 .75]);
    
    % baumgartner data
    data_baum.quart_pe_flat = quantile(data_baum.pe_flat,[.25 .50 .75]);
    data_baum.quart_pe_exp1 = quantile(data_baum.pe_exp1,[.25 .50 .75]);
    data_baum.quart_pe_exp2 = quantile(data_baum.pe_exp2,[.25 .50 .75]);
    
    for c = 1:num_calib
        % simulations data
        pe_flat = exp_macpherson.pe_flat(c,:);
        pe_exp1 = squeeze(exp_macpherson.pe_exp1(c,:,:));
        pe_exp2 = squeeze(exp_macpherson.pe_exp2(c,:,:));

        % simulations
        pe_exp1 = pe_exp1 - repmat(pe_flat(:), 1, size(pe_exp1, 2) );
        pe_exp2 = pe_exp2 - repmat(pe_flat(:), 1, size(pe_exp2, 2) );
        
        % simulations
        quart_pe_flat(c,:) = quantile(pe_flat,[.25 .50 .75]);
        quart_pe_exp1(c,:,:) = quantile(pe_exp1,[.25 .50 .75]);
        quart_pe_exp2(c,:,:) = quantile(pe_exp2,[.25 .50 .75]);
    end


    % plot
    if flags.do_plot
        calib_plot_order = [3,1,2]; % lat, dtf, pge
        
        dx = 1.05;
        FontSize = kv.FontSize;
        MarkerSize = kv.MarkerSize;

        LineColor = [[0.9290, 0.6940, 0.1250]; ...
                    [0.4940, 0.1840, 0.5560]; ...
                    [0.4660, 0.6740, 0.1880]; ...
                    [0.3010, 0.7450, 0.9330]; ...
                    [0.6350, 0.0780, 0.1840]];
        data.Marker = 'ko-';
        data.LineColor = [1 1 1]*0;

        data_reij.Marker = 'v-';
        data_reij.LineColor = [0 0.4470 0.7410];
        
        data_baum.Marker = 'd-';
        data_baum.LineColor = [0.8500 0.3250 0.0980];

        % Exp1
        mFig = figure;
        mFig.Units = 'points';
        mFig.Position = [0 0 500 265];
%         sgtitle(sprintf('%s ', calibrations.combination{c,:}), 'Interpreter', 'none')

        subplot(2,8,1:8)
        mach = errorbar(data.density,data.quart_pe_exp1(2,:,idphase),...
        data.quart_pe_exp1(2,:,idphase) - data.quart_pe_exp1(1,:,idphase),...
        data.quart_pe_exp1(3,:,idphase) - data.quart_pe_exp1(2,:,idphase),...
        'o-','MarkerSize',MarkerSize, 'Color', data.LineColor, ...
        'MarkerFaceColor', data.LineColor);
        hold on
        for c = calib_plot_order
            baru(1,c) = errorbar(data.density/dx,squeeze(quart_pe_exp1(c, 2,:)),...
            squeeze(quart_pe_exp1(c,2,:) - quart_pe_exp1(c,1,:)),...
            squeeze(quart_pe_exp1(c,3,:) - quart_pe_exp1(c,2,:)),...
            's-','MarkerSize',MarkerSize, 'Color', LineColor(c,:),'MarkerFaceColor','w');
        end
        hold on
        reij = errorbar(data.density*dx,data_reij.quart_pe_exp1(2,:),...
        data_reij.quart_pe_exp1(2,:) - data_reij.quart_pe_exp1(1,:),...
        data_reij.quart_pe_exp1(3,:) - data_reij.quart_pe_exp1(2,:),...
        'v--','MarkerSize',MarkerSize, 'Color', data_reij.LineColor,'MarkerFaceColor','w');

        baum = errorbar(data.density*dx,data_baum.quart_pe_exp1(2,:,idphase),...
        data_baum.quart_pe_exp1(2,:,idphase) - data_baum.quart_pe_exp1(1,:,idphase),...
        data_baum.quart_pe_exp1(3,:,idphase) - data_baum.quart_pe_exp1(2,:,idphase),...
        'd--','MarkerSize',MarkerSize, 'Color', data_baum.LineColor,'MarkerFaceColor','w');

        set(gca,'XScale','log','YMinorTick','on')
        set(gca,'XLim',[0.25/1.2 8*1.2],'XTick',data.density,'YLim',[-16 59],'FontSize',FontSize)
        xlabel('Ripple Density [ripples/octave]','FontSize',FontSize)
        ylabel({'Increase in';'Polar Error Rate [%]'},'FontSize',FontSize)

        % Exp2
        sp_ref = subplot(2,8,9:13);
        errorbar(data.depth,data.quart_pe_exp2(2,:,idphase),...
        data.quart_pe_exp2(2,:,idphase) - data.quart_pe_exp2(1,:,idphase),...
        data.quart_pe_exp2(3,:,idphase) - data.quart_pe_exp2(2,:,idphase),...
        'o-','MarkerSize',MarkerSize, 'Color', data.LineColor, ...
        'MarkerFaceColor', data.LineColor);
        hold on
        for c = calib_plot_order
            errorbar(data.depth-0.5,squeeze(quart_pe_exp2(c, 2,:)),...
            squeeze(quart_pe_exp2(c,2,:) - quart_pe_exp2(c,1,:)),...
            squeeze(quart_pe_exp2(c,3,:) - quart_pe_exp2(c,2,:)),...
            's-','MarkerSize',MarkerSize, 'Color', LineColor(c,:),'MarkerFaceColor','w');
        end
        hold on
        errorbar(data.depth+1,data_reij.quart_pe_exp2(2,:),...
        data_reij.quart_pe_exp2(2,:) - data_reij.quart_pe_exp2(1,:),...
        data_reij.quart_pe_exp2(3,:) - data_reij.quart_pe_exp2(2,:),...
        'v--','MarkerSize',MarkerSize, 'Color', data_reij.LineColor,'MarkerFaceColor','w');

        errorbar(data.depth+1,data_baum.quart_pe_exp2(2,:,idphase),...
        data_baum.quart_pe_exp2(2,:,idphase) - data_baum.quart_pe_exp2(1,:,idphase),...
        data_baum.quart_pe_exp2(3,:,idphase) - data_baum.quart_pe_exp2(2,:,idphase),...
        'd--','MarkerSize',MarkerSize, 'Color', data_baum.LineColor,'MarkerFaceColor','w');

        set(gca,'XLim',[data.depth(1)-5 data.depth(end)+5],'XTick',data.depth,...
        'YLim',[-16 59],'YMinorTick','on','FontSize',FontSize)
        xlabel('Ripple Depth [dB]','FontSize',FontSize)
        ylabel({'Increase in';'Polar Error Rate [%]'},'FontSize',FontSize)
        ytick = get(gca,'YTick');
        ticklength = get(gca,'TickLength');

        % Baseline
        subplot(2,8,14:16)
        errorbar(-1,data.quart_pe_flat(2),...
        data.quart_pe_flat(2) - data.quart_pe_flat(1),...
        data.quart_pe_flat(3) - data.quart_pe_flat(2),...
        'o-','MarkerSize',MarkerSize, 'Color', data.LineColor, ...
        'MarkerFaceColor', data.LineColor);
        hold on
        for c = calib_plot_order
            errorbar(-0.6 + (c-1)*0.5,quart_pe_flat(c,2),...
                quart_pe_flat(c,2) - quart_pe_flat(c,1),...
                quart_pe_flat(c,3) - quart_pe_flat(c,2),...
                's-','MarkerSize',MarkerSize, 'Color', LineColor(c,:),'MarkerFaceColor','w');
        end
        hold on
        errorbar(1,data_reij.quart_pe_flat(2),...
        data_reij.quart_pe_flat(2) - data_reij.quart_pe_flat(1),...
        data_reij.quart_pe_flat(3) - data_reij.quart_pe_flat(2),...
        'd--','MarkerSize',MarkerSize, 'Color', data_reij.LineColor,'MarkerFaceColor','w');

        errorbar(1,data_baum.quart_pe_flat(2),...
        data_baum.quart_pe_flat(2) - data_baum.quart_pe_flat(1),...
        data_baum.quart_pe_flat(3) - data_baum.quart_pe_flat(2),...
        'd--','MarkerSize',MarkerSize, 'Color', data_baum.LineColor,'MarkerFaceColor','w');

        set(gca,'XLim',[-3 3],'XTick',0,'XTickLabel',{'Baseline'},...
        'YLim',[-15 59],'YTick',ytick,'TickLength',3*ticklength,...
        'FontSize',FontSize,'YAxisLocation','right')
        xlabel(' ','FontSize',FontSize)
        ylabel({'Polar Error Rate [%]'},'FontSize',FontSize)

        %legend
        for c = calib_plot_order
            switch lower(calibrations.combination{c,2})
                case 'dtf'
                    labels{1,c} = '$T_{DTF}^\varphi$';
                case 'pge'
                    labels{1,c} = '$T_{PGE}^\varphi$';
            end

            if strcmpi(calibrations.combination{c,3}, 'interaural_difference')
                labels{1,c} = '$T_{LAT}^\varphi$';
            end
        end
        leg = legend(sp_ref, [mach, reij, baum, baru(calib_plot_order)], horzcat({'macpherson2003'}, {'reijniers2014'}, {'baumgartner2014'}, labels(calib_plot_order)));
        leg.FontSize = FontSize - 2;
        leg.Units = 'normalized';
        leg.Interpreter = 'latex';
        leg.Orientation = 'horizontal';
        leg.Position = [0.147112074200069 0.944044062817491 0.739274312186318 0.0430594892744974];
        % Overall correlation between actual and predicted median values
        for c=calib_plot_order
            m_pe_pred = [squeeze(quart_pe_exp1(c,2,:))' squeeze(quart_pe_exp2(c,2,:))'];
            m_pe_actual = [data.quart_pe_exp1(2,:) data.quart_pe_exp2(2,:)];
            r = corrcoef(m_pe_pred,m_pe_actual);
            r_sqr = r(2);

            amt_disp('Correlation between actual and predicted median values (15 conditions):')
            amt_disp(sprintf('%s: r = %0.2f', labels{1,c}, r_sqr))
        end
    end
end

function hg = local_middlebroxplot(ax, x, data, Marker, MarkerSize, LineColor, FaceColor)
    lilen = 0.05; % length of horizontal lines
    
    hb=[];
    % Symbols
    i=1; hb(i) = plot(ax, x, data.quantiles(1),'x','MarkerSize',MarkerSize, 'MarkerEdgeColor', LineColor, 'MarkerFaceColor', LineColor); % min
    hold on
    i=i+1; hb(i) = plot(ax, x,data.quantiles(7),'x','MarkerSize',MarkerSize, 'MarkerEdgeColor', LineColor, 'MarkerFaceColor', LineColor); % max

    % Horizontal lines
    i=i+1; hb(i:(i+1)) = line(ax, x+0.5*[-lilen,lilen],repmat(data.quantiles(2),2),'Color',LineColor); % lower whisker
    i=i+2; hb(i:(i+1)) = line(ax, x+[-lilen,lilen],repmat(data.quantiles(3),2),'Color',LineColor); % 25% Quartile
    i=i+2; hb(i:(i+1)) = line(ax, x+[-lilen,lilen],repmat(data.quantiles(4),2),'Color',LineColor); % Median
    i=i+2; hb(i:(i+1)) = line(ax, x+[-lilen,lilen],repmat(data.quantiles(5),2),'Color',LineColor); % 75% Quartile
    i=i+2; hb(i:(i+1)) = line(ax, x+0.5*[-lilen,lilen],repmat(data.quantiles(6),2),'Color',LineColor); % upper whisker

    % Vertical lines
    i=i+2; hb(i:(i+1)) = line(ax, [x,x],data.quantiles(2:3),'Color',LineColor); % connector lower whisker
    i=i+2; hb(i:(i+1)) = line(ax, [x,x],data.quantiles(5:6),'Color',LineColor); % connector upper whisker
    i=i+2; hb(i:(i+1)) = line(ax, [x,x]-lilen,data.quantiles([3,5]),'Color',LineColor); % left box edge
    i=i+2; hb(i:(i+1)) = line(ax, [x,x]+lilen,data.quantiles([3,5]),'Color',LineColor); % left box edge
    
    % middle value
    i=i+1; hb(i) = plot(x,data.mean, Marker,'MarkerSize', MarkerSize, 'MarkerFaceColor', FaceColor, 'MarkerEdgeColor',LineColor);
    
    % create a group to avoid issues with the legend
    % https://stackoverflow.com/questions/12894652/matlab-how-to-make-a-custom-legend
    hg = hggroup;
    set(hb,'Parent',hg);
    set(get(get(hg,'Annotation'),'LegendInformation'),...
      'IconDisplayStyle','off');
  
function rau=rau(X,N,opt)

% RAU   rationalized arcsine transform
% RAU(X,N) transforms the number of correct responses X to the
% rationalized arcsine (rau). N gives the number of repetitions.
% 
% This function allows to use ANOVA statistics with percent correct scores
% because: 1) RAUs are normally distributed; 2) mean and variance of RAUs
% are not correlated with eachother; and 3) likelihood that a score will
% increase/decrease will remain constant over the range.
% 
% RAU=RAU(X,N,opt) defines one of the following options:
%  'Pc'  ... X is given in percent correct scores (0..100%)
%  'X'   ... X is given in the number of correct responses (default)
%
% The formula are based on Sherbecoe and Studebaker,
% Int. J. of Audiology 2004; 43; 442-448
%
% See also IRAU.

% 30.8.2007, Piotr Majdak
%

if exist('opt','var')
  if strcmp(upper(opt),'PC')
    X=X/100*N;
  elseif strcmp(upper(opt),'X')
  else
    error('OPT must be Pc (=percent correct) or X (=number of correct responses)');
  end
end
th=asin(sqrt(X/(N+1)))+asin(sqrt((X+1)/(N+1)));
rau=146/pi*(th)-23;

