% Make contrast maps and principal component dimensions of content dimensions 
% [Figure S3, Figure 3a]
%
% 1. extract dimension score values
% 2. Examine activity during "concept reflection" (regardless of the scores)
% 4. Contrast map with t-test (Figure 3A)
% 5. Beta map of 3 PC with t-test and pruning (Figure S3B-D)
% 6. relationship between 3 PC dimensions (Figure S3A)

%% directory setup
wh_comp = 'accumbens';

[basedir, fastdatadir, fastprojectdir, syncdir] = fast_vary_path(wh_comp);
[D, fast_sub, jth, sub_f] = fast_load_data(syncdir, 'mr-sess1');
modeldir = fullfile(fastprojectdir, 'analysis', 'imaging', 'first_level', 'model07_spm_sing_tri_remv_vw');
load(fullfile(syncdir, 'data', 'FAST_high_vif_all.mat'));

%% 1. extract dimension score values
dimension_view = get_var(D, {'valence','self-related', 'time', 'vividness', 'safe_threat'}, 'conditional', {'eventname', 'view'});
dimension_view = dimension_view(jth,:,:);

[pca_coeff, score, latent, ~, explained] = pca(reshape(dimension_view, size(dimension_view,1)*160, 5));
pca_coeff = pca_coeff(:,1:3);   % Top 3th PC

subjs = fast_sub.fastnum(jth);  % 61 people 

%% 2. Examine activity during "concept reflection" (regardless of the scores)
for subj_i = 1:numel(subjs)
    clear dim_dat*;
    t = tic;
    
    disp('======================================================')
    fprintf('WORKING ON MODEL7: %s (LAST: %s)\n', subjs{subj_i}, subjs{end});
    disp('======================================================')
    
    mask = which('gray_matter_mask.img');
    dat = fmri_data(filenames(fullfile(modeldir, subjs{subj_i}, 'view_*.nii')), mask);
    j_vifs = find(strcmp(vifs_all.subject_id, subjs{subj_i}));
    highvif_idx = vifs_all.vifs{j_vifs}'>3;
    dat.dat(:,highvif_idx) = [];
    dat = remove_empty(dat);
    dat.dat = double(dat.dat);
    
    run_idx = repmat(1:4, 40, 1);
    run_idx = run_idx(:);
    run_idx(highvif_idx) = [];
    
    dim_score = squeeze(dimension_view(subj_i,:,:));
    dim_score(highvif_idx,:) = [];
        
    % content dimensions regression
    dat.X = dim_score*pca_coeff;        % (160x3)
    temp = regress(dat, 'nodisplay');
    
    for reg_i = 1:4     % 3 beta maps and 1 intercept
        if subj_i == 1
            dat_reg{reg_i} = dat;
            dat_reg{reg_i}.dat = [];
        end
        dat_reg{reg_i}.dat(:,subj_i) = temp.b.dat(:,reg_i);     % beta value (211339x4)
    end
end

%% 4. Contrast map with t-test (Figure 3A)

load(fullfile(fastprojectdir, 'data/gray_matter_step4voxels.mat'))
mdat = fmri_data;
mdat.dat =[];
mdat.volInfo = dat_template.volInfo;
for subj_i = 1:numel(subjs)
    fprintf('WORKING ON MODEL7: %s (LAST: %s)\n', subjs{subj_i}, subjs{end});
    
    dat = fmri_data(filenames(fullfile(modeldir, subjs{subj_i}, 'view_*.nii')), 'noverbose');
    load(fullfile(modeldir, subjs{subj_i}, 'vifs.mat'));
    highvif_idx = out.allvifs(contains(out.name, 'view'))>3;
    dat.dat(:,highvif_idx) = [];
    mdat.dat = [mdat.dat mean(dat.dat,2)];
end

tdat = ttest(mdat, .05, 'bfr');
tdat.volInfo = dat.volInfo;
mas = apply_mask(tdat, which('gray_matter_mask.img'));
mas = apply_mask(tdat, which('keuken_2014_enhanced_for_underlay.img'));
orthviews(mas)

mapdir = fullfile(fastprojectdir, 'data','maps');
tdat.fullpath = fullfile(mapdir, 'view_average_model7_bfr05.nii');
write(tdat, 'thresh')

tdat.fullpath = fullfile(mapdir, 'view_average_model7_unthresholded.nii');
write(tdat);

%% 5. Beta map of 3 PC with t-test and pruning (Figure S3B-D)
dimensions = {'PC1_valence', 'PC2_self', 'PC3_time'};
for i = 1:numel(dimensions)
    tdat = ttest(dat_reg{i}, .01, 'unc');
    tdat.fullpath = fullfile(mapdir, ['view_dim' num2str(i) '_' dimensions{i} '_model7_unc01.nii']);
    write(tdat, 'thresh')
end

p = [.05 .01 .001];
k = [1 1 1];

for i = 1:numel(dimensions)
    tdat = ttest(dat_reg{i}, .01, 'unc');
    [dat, r] = pruning_img(tdat, p, k);
    dat.fullpath = fullfile(mapdir, ['view_dim' num2str(i) '_' dimensions{i} '_model7_unc001_01_05pruned.nii']);
    write(dat)
end

%% 6. relationship between 3 PC dimensions (Figure S3A)
load(which('colormap_colorbrewer_wani.mat'))
cols = col11_from_colorbrewer([8 2 1 4 10],:);
for i = 1:2
    subplot(1,2,i);
    scatter(pca_coeff(:,i), pca_coeff(:,i+1), 100, cols, 'filled');
    set(gca, 'linewidth', 1.5, 'tickdir', 'out', 'ticklength', [.015 .015], 'fontsize', 20);
    if i == 1
        set(gca, 'xlim',[0 .8], 'ylim', [-.6 .7])
    elseif i == 2
        set(gca, 'xlim',[-.6 .7], 'ylim', [-.3 1], 'xtick', -6:.3:.6)
    end
end

set(gcf, 'position', [307 1004 873 334], 'color', 'w');
