Picture-Word-N1/04 Analysis/analyse_14_plot_sample_picture_word_topo.m

%% Setup

addpath('matlab_extensions');

% import eeglab (assumes eeglab has been added to path), e.g.
addpath('C:/EEGLAB/eeglab2020_0')
[ALLEEG, EEG, CURRENTSET, ALLCOM] = eeglab;

erps = {'sample', 'sample_picture', 'sample_response'};

% settings for x ticks
time_mins = [-200, -200, -1000];
time_steps = [100, 200, 200];

% special settings for y range of fixed effect slopes
slope_y_lims = {[], [-4, 4], []};
slope_y_ticks = {[], -4:2:4, []};

% settings for the event lines and topo plots
event_times = {[500], [1000], []};
event_rects = {[], [1150, 1650], []};

erp_topo_times = {[105, 185, 235, 300, 400, 660], [110, 150, 225, 600, 800, 1150], [-800, -600, -400, -200, 0, 200]};

for e_nr = 1:numel(erps)

    e = erps{e_nr};

    % import LMM estimates
    res_file = sprintf('analyse_13_%s_res.csv', e);
    fe_res = readtable(res_file);
    fe_res.ch_name = upper(fe_res.ch_name);  % upper case for case insensitivity in joins

    % get summary variables
    channels = table2array( unique(fe_res(:, 'ch_name')) );
    times = table2array( sortrows(unique(fe_res(:, 'time'))) );
    fes = table2array( unique(fe_res(:, 'fe_lab')) );
    
    % get channel locations
    chanlocs = pop_chanedit([], 'load', {'BioSemi64.loc', 'filetype', 'loc'});
    
    % get channel order
    ch_name = upper({chanlocs.labels}');
    ch_order = cell2mat({chanlocs.urchan}');
    chan_orders = table(ch_name, ch_order);
    
    % only keep channels in the data in the table, and the chanlocs struct
    ch_in_fe = false(numel(chan_orders.ch_name), 0);
    
    for i = 1:numel(chan_orders.ch_name)
        ch = chan_orders.ch_name{i};
        ch_in_fe(i) = sum( strcmp(unique(fe_res.ch_name), ch) );
    end
    
    chan_orders = chan_orders(ch_in_fe, :);
    chanlocs = pop_chanedit(chanlocs, 'delete', find(~ch_in_fe));
    
    % join channel order to fe_res
    fe_res = join(fe_res, chan_orders, 'Keys', 'ch_name');
    
    %% get xyz colours
    
    % note: order of x and y are flipped, and y is reversed, to approximate the MNE coordinates
    rgb = [-[chanlocs.Y]; [chanlocs.X]; [chanlocs.Z]];
    
    for colour_chan = 1:3
        rgb(colour_chan, :) = normalize(rgb(colour_chan, :), 'range');
    end
    
    %% get fixed effects in matrices
    
    % coerce fixed effects into matrices of time x channel
    mats = struct();
    
    for fe = fes'
        % subset relevant rows of table
        fe_res_i = fe_res(strcmp(fe_res.fe_lab, fe), :);
        % sort rows of table
        fe_res_i = sortrows(fe_res_i, {'ch_order', 'time'});
        % extract fixed effect estimates and put in matrix
        % (matlab reshape uses fortan order)
        field_name = regexprep(fe{:}, '*', '');  % asterisk disallowed in field name
        field_name = regexprep(field_name, '\s\s', '_');  % space disallowed in field name
        mats.(field_name) = reshape( fe_res_i.Estimate, [numel(times), numel(channels)] );
    end
    
    %% Join all plots together
    
    fig = figure;
    set(0, 'DefaultLineLineWidth', 1);
    set(0, 'DefaultAxesLineWidth', 1);
    vhline_width = 1.25;
    
    fe_names = {'Intercept', 'Congruency', 'Predictability', 'Congruency_Predictability'};
    
    % get tidy labels
    fe_labels = {};
    max_nchar = 0;
    for fe = 1:numel(fe_names)
        fe_labels{fe} = split( strcat( regexprep(fe_names{fe}, '_', ' ×nl'), ' (µV)' ), 'nl' );
        if strlength(fe_labels{fe}) > max_nchar
            max_nchar = strlength(fe_labels{fe});
        end
    end
    
    % pad and justify with spaces
    for fe = 1:numel(fe_names)
        fe_labels{fe} = strjust(pad(fe_labels{fe}, max_nchar), 'center');
    end
    
    % get topoplot settings
    y_lims = [floor(min(fe_res.Estimate)), ceil(max(fe_res.Estimate))];
    map_lims = [-max(abs(y_lims)), max(abs(y_lims))];
    map_range = map_lims(2) - map_lims(1);
    map_ticks = map_lims(1):(map_range/4):map_lims(2);
    
    cmap_resolution = 101;
    topo_resolution = 101;
    bwr = [0, 0, 1; 1, 1, 1; 1, 0, 0];
    %cmap_bwr = interp1([-cmap_resolution/2; 0; cmap_resolution/2], bwr, (-cmap_resolution/2:cmap_resolution/2));
    
    cmap_bwr = twilight_shifted(cmap_resolution);
    
    % linear interpolation to focus on specific timepoints
    topo_times = erp_topo_times{e_nr};
    
    % x axis for timecourse
    time_lims = [round(min(times), -1), round(max(times), -1)];
    time_ticks = time_mins(e_nr):time_steps(e_nr):time_lims(2);
    
    % parameters for adjusting plotting space
    x_adj_l = 0.025;  % offset space to be added to the left
    x_adj_r = 0.01;  % space to be added to width
    
    y_adj_top = -0.05; % space to be added (or removed if negative) from the top
    
    % define the shape of the subplots
    sp_ncol = numel(topo_times);
    sp_nrow = 2 * numel(fe_names);
    
    % generate subplot coordinates in normalised units
    % xmin = .175;
    xmin = .255;
    xmax = .975;
    ymin = -0.05;
    ymax = 0.89;
    
    x_width = (xmax-xmin)/sp_ncol;
    x_props = xmin:x_width:xmax;
    y_height = (ymax-ymin)/sp_nrow;
    y_props = flip(ymin:y_height:ymax);
    
    y_props_topo_adj = y_props + y_height*0.075;
    topo_yheight = y_height * 0.8;
    
    % plot!
    
    % plot colours for timecourse legend (must come first to preserve colorbar)
    colour_legend_pos = [xmin/2 - x_width/1.5, y_props_topo_adj(1), x_width, topo_yheight];  % x position may need adjusting if fig size changes
    subplot('Position', colour_legend_pos);
    hold on;
    topoplot(zeros(numel(chan_orders.ch_name), 0), chanlocs, 'electrodes', 'off');
    % set line width
    set(findall(gca, 'Type', 'Line'), 'LineWidth', 1);
    for i = 1:numel(chan_orders.ch_name)
        topoplot(zeros(numel(chan_orders.ch_name), 0), chanlocs, 'colormap', [0,0,0], 'emarker', {'.', rgb(:, i)', 11, 1}, 'plotchans', i, 'headrad', 0);
    end
    hold off
    
    % now, plot everything else
    for fe = 1:numel(fe_names)
        fe_name = fe_names{fe};
    
        mat_fe = mats.(fe_name);
    
        
        if strcmp(fe_name, 'Intercept')
            fe_topo_lims = map_lims;
        else
            if numel(slope_y_lims{e_nr})>0
                fe_topo_lims = slope_y_lims{e_nr};
            else
                fe_topo_lims = map_lims;
            end
        end

        topo_nfus = zeros(numel(topo_times), 2);
    
        % plot topoplots
        for time_nr = 1:numel(topo_times)
            subplot('Position', [x_props(time_nr), y_props_topo_adj( fe * 2 - 1 ), x_width, topo_yheight]);   
            mat_fe_topo_t = interp1(times, mat_fe, topo_times(time_nr));
            topoplot(mat_fe_topo_t, chanlocs, 'colormap', cmap_bwr, 'gridscale', topo_resolution, 'electrodes', 'off', 'maplimits', fe_topo_lims);
            % get lines
            lines = findall(gca, 'Type', 'Line');
            % set line width
            set(lines, 'LineWidth', 1);
            % get surfaces to extract y locations
            surf = findall(gca, 'Type', 'Surface');
            % get topo location to draw line to
    %         [tx, ty] = ds2nfu(mean(lines(4).XData), min(lines(4).YData));
            [tx, ty] = ds2nfu(0, min(surf.YData(:)));
            topo_nfus(time_nr, :) = [tx, ty];
        end
        
        % plot timecourse
        subplot('Position', [x_props(1), y_props( fe * 2 ), x_width*sp_ncol, y_height])
        hold on

        % plot rectangle for jittered period
        rect = event_rects{e_nr};
        if ~isempty(rect)
            rectangle('Position', [rect(1) map_lims(1) rect(2)-rect(1) map_range], 'FaceColor', [0.4, 0.4, 0.4], 'EdgeColor', 'none')
        end
        
        if ~strcmp(fe_name, 'Intercept')
            plot(times, mats.('Intercept'), 'Color', [0.8, 0.8, 0.8]);
        end
        
        for ch = 1:size(mat_fe, 2)
            plot(times, mat_fe(:, ch), 'Color', rgb(:, ch), 'LineWidth', 0.55);
        end
    
        xlim(time_lims)
        xticks(time_ticks)
        set(gca,'FontSize', 8)
    
        ylim(fe_topo_lims)

        xline(0, '--', 'LineWidth', vhline_width)
        yline(0, '--', 'LineWidth', vhline_width)

        for t = event_times{e_nr}
            xline(t, '--', 'LineWidth', vhline_width)
        end
    
        tc_nfus = zeros(numel(topo_times), 2);
    
        for time_nr = 1:numel(topo_times)
            xline(topo_times(time_nr), 'LineWidth', vhline_width, 'Color', [0.4, 0.4, 0.4])
            % get the location of the top of the line
            [tcx, tcy] = ds2nfu(topo_times(time_nr), max(fe_topo_lims));
            tc_nfus(time_nr, :) = [tcx, tcy];
        end
    
        hold off
    
        % draw annotation lines between lines and topoplots
        for time_nr = 1:numel(topo_times)
            % colour is set to .6 intensity rather than .4 in vertical lines above - hacky solution for bug(?)
            annotation('line', [tc_nfus(time_nr, 1), topo_nfus(time_nr, 1)], [tc_nfus(time_nr, 2), topo_nfus(time_nr, 2)], 'Color', [0.6, 0.6, 0.6], 'LineWidth', vhline_width);
        end
    
        % plot styling
    
        % add colorbar to the y axis to show corresponding y and colour values
        ax = gca;
        ax_pos = get(ax, 'Position');
        cb = colorbar(ax);

        if strcmp(fe_name, 'Intercept')
            cb.Ticks = map_ticks;
        else
            if numel(slope_y_ticks{e_nr})>0
                cb.Ticks = slope_y_ticks{e_nr};
            else
                cb.Ticks = map_ticks;
            end
        end

%         cb.Ticks = map_ticks;
        cb.TickDirection = 'out';
        cb.Position = [ax_pos(1)-0.015, ax_pos(2), .015, ax_pos(4)];
    
        cb.Label.String = fe_labels{fe};
        cb.Label.Rotation = 0;
        cb.Label.VerticalAlignment = 'middle';
        cb.Label.HorizontalAlignment = 'right';
        cb.FontSize = 8;
        cb.Label.FontSize = 11;
        
        ax.Colormap = cmap_bwr;
        ax.CLim = fe_topo_lims;
        ax.YTickLabels = [];
        ax.YLabel = [];
        ax.YAxis.Visible = 'off';
    
        % ticks outside plot
        set(gca,'TickDir','out');
    
        if strcmp(fe_name, 'Congruency_Predictability')
            % add x label
            xlabel('Time (ms)', 'FontSize', 12)
        else
            % remove x tick labels
            set(gca,'XTickLabel', {[]});
        end
    
        % put white line over the default axes to remove them while keeping the ticks and labels (hacky)    
        ax2 = axes('Position',gca().Position,...
        'XColor',[1 1 1],...
        'YColor',[1 1 1],... 
        'Color','none',...
        'XTick',[],...
        'YTick',[]);
    
    end
    
    set(fig, 'Units', 'Inches', 'Position', [0, 0, 6.5, 7], 'PaperUnits', 'Inches', 'PaperSize', [6.5, 7])
    
    plot_file = sprintf('figs/14_topo_timecourse_%s.pdf', e);
    exportgraphics(fig, plot_file, 'BackgroundColor','none')

    if strcmp(e, 'sample')

        % plot the N400 model predictions
        N400_t = 400;
    
        pred_levels_uncoded = 0.1:0.1:1;
        pred_levels = interp1([0.07, 1], [0, 1], pred_levels_uncoded);  % because predictability was normalised between 0 and 1
        %cong_levels = [-0.5, 0.5];
        cong_levels = [-0.4938494, 0.5061506];  % [incong, cong] - copied from the unique deviation values in the script that fits the models
    
        B = struct();
        for fe_name = fe_names
    	    B.(fe_name{:}) = interp1(times, mats.(fe_name{:}), N400_t);
        end
    
        model_preds = zeros(numel(pred_levels), numel(cong_levels), height(chan_orders));

%         topoplot(B.Congruency, chanlocs, 'colormap', cmap_bwr, 'gridscale', topo_resolution, 'electrodes', 'off');
    
        for pred_lvl_nr = 1:numel(pred_levels)
            p = pred_levels(pred_lvl_nr);
    
            for cong_lvl_nr = 1:numel(cong_levels)
                c = cong_levels(cong_lvl_nr);
    
                model_preds(pred_lvl_nr, cong_lvl_nr, :) = B.Intercept +...
                    B.Congruency * c +...
                    B.Predictability * p +...
                    B.Congruency_Predictability * p * c;
    
            end
        end
    
        abs_max = max(abs([floor(min(model_preds, [], 'all')*2)/2, ceil(max(model_preds, [], 'all')*2)/2]));
        preds_topo_lims = [-abs_max, abs_max];
        N400_fig = figure;

        t = tiledlayout(numel(cong_levels), numel(pred_levels) + 2, 'TileSpacing', 'none', 'Padding', 'none');

        for cong_lvl_nr = [2, 1] %1:numel(cong_levels)
            
            nexttile;
            
            if cong_lvl_nr == 1
                cong_lab = 'Incongruent';
            else
                cong_lab = 'Congruent';
            end

            text(0.33, 0.5, cong_lab, 'FontWeight', 'normal', 'HorizontalAlignment', 'center', 'FontSize', 9)
            ax = gca;
            ax.Visible = 0;

            for pred_lvl_nr = 1:numel(pred_levels)
                t_i = nexttile;
                topoplot(model_preds(pred_lvl_nr, cong_lvl_nr, :), chanlocs, 'colormap', cmap_bwr, 'gridscale', topo_resolution, 'electrodes', 'off', 'maplimits', preds_topo_lims);
                set(findall(gca, 'Type', 'Line'), 'LineWidth', 0.5);

                if cong_lvl_nr == 2
                    title( sprintf('%g%%', pred_levels_uncoded(pred_lvl_nr)*100), 'FontWeight', 'normal', 'FontSize', 9)
                end
            end

            if cong_lvl_nr == 2
                nexttile([2 1])
                ax = gca;
                ax.Visible = 0;

                cb = colorbar(t_i);
                cb.Ticks = preds_topo_lims(1):2.5:preds_topo_lims(2);
                cb.TickDirection = 'out';
                cb.Label.String = 'N400 Amplitude (µV)';
                cb.FontSize = 8;
                cb.Label.FontSize = 9;
                cb.Position = [0.92,0.13,0.0125,0.75];
                cb.Label.Position(1) = 0.8 * cb.Label.Position(1);
            end

        end

        set(N400_fig, 'Units', 'Inches', 'Position', [0, 0, 6.5, 1.4], 'PaperUnits', 'Inches', 'PaperSize', [6.5, 1.4])

        N400_plot_file = sprintf('figs/14_N400_%s.pdf', e);
        exportgraphics(N400_fig, N400_plot_file, 'BackgroundColor','none')

    end

end