% clean up for report
clc; close all;

% save figure?
savefig     =   0;
figpath     =   fullfile( paths.local, 'images');

if savefig && ~exist( figpath, 'dir')
    mkdir( figpath);
end

% prepare names
conditions   =  {'size',            {'near',  'far'};   ...
                 'shape',           {'mel',   'box'}};
names.train  =  {'Object size',     {'large', 'small'}; ...
                 'Object type',     {'melon', 'box'}};   
names.test   =  {'Search distance', {'near',  'far'};   ...
                 'Object cue',      {'melon', 'box'}}; 

coltyp.shape =  {'m', 'd-'};
coltyp.size  =  {'b', 'o-'};
metric_label =  'Parameter estimate (a.u.)';                    

% retrieve some handy variables
allrois     =   fields(RS);
exroi       =   allrois{1};

sided       =   RS.(exroi).info.stats.sided;
cutoff      =   RS.(exroi).info.stats.cutoff;
mytest      =   RS.(exroi).info.mytest;
all_cons    =   fields( RS.(exroi).data.vox);
nsubs       =   size( RS.(exroi).data.all.(all_cons{1}).(conditions{1,1}), 1);
   
% default plot settings
transp      =   .2;     % opacity of error region
transp_i    =   .5;     % opacity of individual data

nsigs       =   3;      % number of significance levels (factor 10)
sig_spac    =   8;      % space used by significance reports (fraction of graph height)
sig_dist    =   40;     % distance between  (~size of) significance reports (fraction of graph height)

jitsize     =   0.5;    % jitter width of individual data points

a_size      =   8;      % text size of asterisks (across ROIs)
a_size_a    =   20;     % text size of asterisks (across ROIs)

marksize    =   6;      % marker size (condition averages)
capsize     =   10;     % cap size (condition averages)
ind_frac    =   0.3;    % size of individual data wrt condition averages

% additional settings for hemifield graphs
bandw       =   0.05;               % kernel size of density plot (#SD in classification metric)
scalef      =   10;                 % horizontal (down)scaling of density plot
pdist       =   0.1;                % horizontal separation of density plot
transp3     =   transp;             % transparency of density plot
lwidth      =   1;                  % linewidth of 95%-CI per hemifield
jits        =   rand( nsubs, 1);    % ensure same jitter for all participants


% start plotting (per ROI category)
for nrt = 1:length( show.rois)
    
    % get ROI name
    gen_roi     =   show.rois{nrt};
    
    % retrieve handy variables (plotting)
    vcount          =   RS.(gen_roi).info.voxels;
    X               =   1 : length( vcount);
    
    for nmod = 1 : length( all_cons)
    
        % retrieve model name
        gen_mod = all_cons{nmod};
        
        % Open figure for current ROI
        figure('Name', sprintf( '%s: %s data', gen_roi, gen_mod), 'NumberTitle', 'off')
        hold on
        
        % retrieve stats and data
        voxstat =   RS.(gen_roi).stats.vox.(gen_mod);
        voxdata =   RS.(gen_roi).data.vox.(gen_mod);
       
        % retrieve available field names
        fnames      =   conditions( contains( conditions(:,1), fields( voxdata)), 1);  
        legends     =   [names.(gen_mod)(contains( conditions(:,1), fields( voxdata)),1)', '95% CIs'];
        
        % plot confidence intervals
        for nplot = 1 : length(fnames)
            Y = [ voxstat.(fnames{nplot}).ci(2,:), fliplr( voxstat.(fnames{nplot}).ci(1,:))];
            h =   fill( [X fliplr(X)], Y, coltyp.(fnames{nplot}){1});
            set( h, 'facealpha', transp, 'EdgeColor', 'none')
        end
        
        % plot chance line
        plot( [0 length(X)+1], [0 0], 'k.:');
        
        % plot means
        for nplot = 1 : length(fnames)
            
            % plot means
            hl(nplot) = plot( X, voxstat.(fnames{nplot}).mean, [coltyp.(fnames{nplot}){:}], 'MarkerFaceColor', coltyp.(fnames{nplot}){1}, 'MarkerEdgeColor', 'k', 'MarkerSize', marksize); %#ok<SAGROW>
        end
    
        % retrieve some variables reflecting graph size
        if isempty( show.fixedy)
            theselims   =   ylim;
        else
            theselims      =   show.fixedy;
            ylim( theselims);
        end
                
        % plot significance (using requested statistical method and null)
        for nplot =  1 : length(fnames)
            
            % retrieve y locations for signifance plots
            yloc = theselims(1) + diff(theselims)/sig_spac - diff(theselims)/sig_dist * nplot;
            
            switch show.multi
                
                % choose method
                case {'tfce','none'}
                    
                    % get requested p values
                    if strcmp( show.multi, 'tfce')
                        pvals   =   voxstat.(fnames{nplot}).tfce.(show.null) / (3-sided);
                        pvals   =   min( pvals, 0.5-pvals);
                    elseif strcmp( show.stats, 'ttype')
                        pvals   =   voxstat.(fnames{nplot}).pvals;
                    end
                    
                    % plot asterisks
                    for nsig = 1:nsigs
                        if nsig < nsigs
                            wsigs   =   pvals < (cutoff / 10^(nsig-1)) & pvals >= (cutoff / 10^(nsig));
                        else
                            wsigs   =   pvals < (cutoff / 10^(nsig-1));
                        end
                        sigsX   =   X( wsigs);
                        sigsY   =   ones(1, length( sigsX)) * yloc;
                        text( sigsX, sigsY, repmat('\ast',1,nsig), 'FontSize', a_size, 'Color', coltyp.(fnames{nplot}){1}, 'HorizontalAlignment', 'left');
                    end
                    
                    % cluster-based corrected significance
                case 'clus'
                    
                    % retrieve significant clusters
                    ks = voxstat.(fnames{nplot}).clus.(show.null).sig_k;
                    maxclus = max( ks);
                    if maxclus > 0
                        
                        % plot significant clusters
                        for nclus = 1:maxclus
                            xloc = find( ks==nclus);
                            plot( xloc, ones(1, length( xloc))*yloc, [coltyp.(fnames{nplot}){1} '.-'], 'LineWidth', 4)
                        end
                    end
            end
        end % end of plot significance loop
        
        % add legend
        [~, h_l] = legend( [hl h], legends, 'Location', 'east');
        lpatch   = findobj(h_l, 'type', 'patch');
        set( lpatch, 'FaceAlpha', transp);
        set( lpatch, 'FaceColor', [0 0 0]);
        
        % add ticks and labels
        xticks( X)
        xticklabels( vcount)
        ylabel( metric_label);
        xlabel( 'Voxel count');
        
        % final touch
        title( sprintf( 'Univariate contrast per %s ROI size (%s data)', gen_roi, mytest));
        pbaspect([3 1 1]);
        set(gcf, 'renderer', 'painters')
        hold off
        
        % save figure
        if savefig
            figname = sprintf( 'Exp2_%s_%s_%s', mytest, gen_roi, show.multi);
            print( gcf, fullfile( figpath, figname), '-depsc');
        end
        
        % Report cluster significance of different conditions
        if strcmp( show.multi, 'clus')
            
            % printf info
            fprintf( '\n%s - cluster-level signifcance per condition (2-sided):\n\n', gen_roi);
            
            % iterate over conditions
            for nc = 1:size( conditions, 1)
                
                % retrieve cluster indices for current condition/metric
                con_clus    =   voxstat.(conditions{nc,1}).clus.(show.null).sig_k;
                clus_sig    =   voxstat.(conditions{nc,1}).clus.(show.null).k_prob;
                nsigni      =   max( con_clus);
                
                if nsigni
                    for ncl = 1:nsigni
                        c_size = sum( con_clus==ncl);
                    end
                    fprintf( '%s [%s] cluster %g/%g: \tp = %1.3f\n', gen_roi, conditions{nc,1}, ncl, nsigni, clus_sig(c_size));
                end
            end
            fprintf('\n==================================================\n')
        end
    end % end of train/test loop
end % end of ROI loop


% hemi-field plots
if isfield( show, 'hemi') && ~isempty( show.hemi)
    
    % iterate over requested ROIs
    for nroi = 1 : length( show.hemi)
        
        % define average and compound ROIs
        m_roi = show.hemi{nroi};
        l_roi = [m_roi(1:2) 'L'];
        r_roi = [m_roi(1:2) 'R'];
        
        
        for nmod = 1 : length( all_cons)
            
            % retrieve model name
            gen_mod = all_cons{nmod};
            
            % retrieve available field names
            fnames      =   conditions( contains( conditions(:,1), fields( RS.(m_roi).data.vox.(gen_mod))), 1);
            
            % print info
            fprintf( '\n%s data: %s - Significance test for full ROI (%g-sided):\n\n', upper(gen_mod), m_roi, sided);
            
            for ncon = 1:length( fnames)
            
                % Create Figure for "all" ROI
                cname = names.(gen_mod){ncon,1};
                figure('Name', sprintf( '%s %s (%s data)', m_roi, cname, gen_mod), 'NumberTitle', 'off')
                hold on
                
                % find current conditions
                thesecons   =   [conditions{ncon,2}, conditions{ncon,1}];
                thesenames  =   [conditions{ncon,2}, 'difference'];
                plotcons    =   thesenames{show.wplots};
                ncons       =   length( show.wplots);
                
                % prepare to save min max values (for determing y length)
                mm  =   NaN( ncons, 2);
                
                % plot individual data
                for wplot = show.wplots
                    
                    % get data
                    yM      =   RS.(m_roi).data.all.(gen_mod).(thesecons{wplot})(:,:,1);
                    yR      =   RS.(r_roi).data.all.(gen_mod).(thesecons{wplot})(:,:,1);
                    yL      =   RS.(l_roi).data.all.(gen_mod).(thesecons{wplot})(:,:,1);
                    
                    mM      =   RS.(m_roi).stats.all.(gen_mod).(thesecons{wplot}).mean;
                    cM      =   RS.(m_roi).stats.all.(gen_mod).(thesecons{wplot}).ci;
                    mR      =   RS.(r_roi).stats.all.(gen_mod).(thesecons{wplot}).mean;
                    cR      =   RS.(r_roi).stats.all.(gen_mod).(thesecons{wplot}).ci;
                    mL      =   RS.(l_roi).stats.all.(gen_mod).(thesecons{wplot}).mean;
                    cL      =   RS.(l_roi).stats.all.(gen_mod).(thesecons{wplot}).ci;
                    
                    % get min max
                    mm(wplot, :)    =   [ min( [yR;yL]), max( [yR;yL])];
                    
                    % retrieve requested p-values
                    if strcmp( show.stats, 'ttype')
                        pvals   =   RS.(m_roi).stats.all.(gen_mod).(thesecons{wplot}).pvals;
                        pvals_l =   RS.(l_roi).stats.all.(gen_mod).(thesecons{wplot}).pvals;
                        pvals_r =   RS.(r_roi).stats.all.(gen_mod).(thesecons{wplot}).pvals;
                        
                    elseif strcmp( show.stats, 'prob')
                        pvals   =   RS.(m_roi).null.all.(gen_mod).(thesecons{wplot}).pvals;
                        pvals_l =   RS.(l_roi).null.all.(gen_mod).(thesecons{wplot}).pvals;
                        pvals_r =   RS.(r_roi).null.all.(gen_mod).(thesecons{wplot}).pvals;
                    end
                    
                    % compute density
                    [FR,XIR] = ksdensity( yR, 'Bandwidth', bandw);
                    [FL,XIL] = ksdensity( yL, 'Bandwidth', bandw);
                    
                    % cutoff low value
                    FR(FR<bandw/10) = 0;
                    FL(FL<bandw/10) = 0;
                    
                    % scale and position
                    FR = FR/scalef;
                    FL = FL/scalef;
                    bR = wplot+pdist;
                    bL = wplot-pdist;
                    
                    % determine jitter position of dots (within density plots)
                    xRs =   interp1( XIR, FR, yR, 'pchip') *0.9;
                    xLs =   interp1( XIL, FL, yL, 'pchip') *0.9;
                    xR  =    xRs.*jits + bR;
                    xL  =   -xLs.*jits + bL;
                    
                    % draw half violin plot
                    patch( [ FR+bR, fliplr(  FR+bR)], [zeros( size( XIR)), fliplr( XIR)      ], coltyp.(fnames{ncon}){1}, 'EdgeAlpha', 0', 'FaceAlpha', transp3);
                    patch( [-FL+bL, fliplr( -FL+bL)], [XIL,                zeros( size( XIL))], coltyp.(fnames{ncon}){1}, 'EdgeAlpha', 0', 'FaceAlpha', transp3);
                    
                    % draw conf intervals
                    ycR =   round( interp1( XIR, 1:length(XIR), cR, 'pchip'));
                    ycL =   round( interp1( XIL, 1:length(XIL), cL, 'pchip'));
                    cFR =   FR(ycR(1):ycR(2));
                    cFL =   FL(ycL(1):ycL(2));
                    cXR =   XIR(ycR(1):ycR(2));
                    cXL =   XIL(ycL(1):ycL(2));
                    
                    patch( [ cFR+bR, bR*ones(1, length( cXR))], [cXR, fliplr( cXR)], coltyp.(fnames{ncon}){1}, 'EdgeAlpha', 0', 'FaceAlpha', transp3);
                    patch( [-cFL+bL, bL*ones(1, length( cXL))], [cXL, fliplr( cXL)], coltyp.(fnames{ncon}){1}, 'EdgeAlpha', 0', 'FaceAlpha', transp3);
                    
                    % draw scattered dots (within density plots)
                    scatter( xR, yR, coltyp.(fnames{ncon}){1}, 'MarkerFaceColor', coltyp.(fnames{ncon}){1}, 'Marker', coltyp.(fnames{ncon}){2}(1), 'SizeData', round( marksize*ind_frac), 'MarkerFaceAlpha', transp_i, 'MarkerEdgeAlpha', transp_i);
                    scatter( xL, yL, coltyp.(fnames{ncon}){1}, 'MarkerFaceColor', coltyp.(fnames{ncon}){1}, 'Marker', coltyp.(fnames{ncon}){2}(1), 'SizeData', round( marksize*ind_frac), 'MarkerFaceAlpha', transp_i, 'MarkerEdgeAlpha', transp_i);
                    
                    % add 95% CI for L & R (and mean)
                    %plot( [bR bR; bL bL]', [cR;cL]', sprintf( '%s-', 'k'), 'LineWidth', lwidth);
                    %plot( [bR bR+pdist; bL bL-pdist]', [mR mR; mL mL]', sprintf( '%s-', 'k'), 'LineWidth', lwidth);
                    
                    % draw error bar of full ROI
                    errorbar( wplot, mM, mM-cM(1), cM(2)-mM, 'ko', 'MarkerFaceColor', coltyp.(fnames{ncon}){1}, 'CapSize', capsize, 'MarkerSize', marksize, 'Marker', coltyp.(fnames{ncon}){2}(1));
                    
                    % set y length
                    if ~isempty( show.fixedy)
                        ylim( fixedy);
                    end
                    theselims = ylim;
                    
                    % plot asterisks
                    for nsig = 1:nsigs
                        if nsig < nsigs
                            wsigs   =   pvals < (cutoff / 10^(nsig-1)) & pvals >= (cutoff / 10^(nsig));
                        else
                            wsigs   =   pvals < (cutoff / 10^(nsig-1));
                        end
                        if wsigs
                            sigsX   =   wplot;
                            sigsY   =   theselims(1) + diff(theselims)/sig_spac - diff(theselims)/sig_dist;
                            text( sigsX, sigsY, repmat('\ast',1,nsig), 'FontSize', a_size_a, 'Color', coltyp.(fnames{ncon}){1}, 'HorizontalAlignment', 'center');
                        end
                    end
                    
                    % report significance of hemi ROI
                    fprintf( '%s \t[%s] %s %s: \tp = %1.3f\n',    m_roi,                gen_mod, conditions{ncon,1}, thesenames{wplot}, pvals);
                    fprintf( '%s-%s \t[%s] %s %s: \tp = %1.3f\n', l_roi(1:2), l_roi(3), gen_mod, conditions{ncon,1}, thesenames{wplot}, pvals_l);
                    fprintf( '%s-%s \t[%s] %s %s: \tp = %1.3f\n', r_roi(1:2), r_roi(3), gen_mod, conditions{ncon,1}, thesenames{wplot}, pvals_r);
                end
                
                fprintf('\n==================================================\n')
                
                % plot chance line
                xlims = [show.wplots(1)-0.5 show.wplots(end)+0.5];
                plot( xlims, [0 0], 'k.:');
                
                % set y length
                if isempty( show.fixedy)
                    minmax      =   [min( mm(:,1)), max( mm(:,2))];
                    theselims   =   [minmax(1)-diff( minmax)*0.1, minmax(2)+diff( minmax)*0.1];
                else
                    theselims      =   show.fixedy;
                end
                ylim( theselims);
                
                % add ticks and labels
                xticks( show.wplots)
                xticklabels( plotcons)
                ylabel( metric_label);
                xlabel( 'Condition');
                xlim( xlims);
                
                % final touch
                title( sprintf( '%s-ROI %s (%s data)', m_roi, fnames{ncon}, gen_mod));
                pbaspect([3 4 1]);
                set( gcf, 'Renderer', 'painters');
                hold off
                
                % save figure
                if savefig
                    figname = sprintf( 'Exp2_%s_%s_%s_hemi_%s', gen_mod, mytest, m_roi, fnames{ncon});
                    print( gcf, fullfile( figpath, figname), '-depsc');
                    pause(0.5);
                end
                
            end % end of loop over shape/size
        end % end of loop over train /test
    end % end of hemi-ROI loop
end % end of hemi if-statement