%
%   Generate Figure 3A (Whole-brain surface maps with best EC) for FaceSpace paper
%   _______________________________________________
%   by Marcelo G Mattar (10/21/2013)
%   mattar@sas.upenn.edu

% Parameters
datasets = {'MDS_Multi'}; % Choose between 'MDS_RGI' and 'MDS_Multi'
saveBool = false;
sigThresh = 0.5; % minimum p-value required for the inclusion of a vertex in the figure

% Directories
baseDir = '/jet/aguirre/FaceSpace/';
dataDir = '~/Data/FaceSpace/';
outputDir = './Figures/'; % Directory where output images will be saved

filename_endings = {'000' '005' '010' '015' '020' '025' '030' '035' '040' '045' '050' '055' '060' '065' '070' '075' '080' '085' '090' '095' '100'};
numCoef = length(filename_endings);

for d = 1:length(datasets)
    
    dataset = datasets{d};
    groupBetas_lh = zeros(163842,numCoef);
    groupBetas_rh = zeros(163842,numCoef);
    groupSigs_lh = zeros(163842,numCoef);
    groupSigs_rh = zeros(163842,numCoef);
    
    for k = 1:numCoef
        disp(['Loading surface maps for Coef = ' filename_endings{k}]);
        this_coef = filename_endings{k};
        
        [groupBetas_lh(:,k), M_lh, mr_parms_lh] = load_mgh([dataDir 'group_glms/' dataset '/Drift_' this_coef '/lh.Drift_' this_coef '.glmdir/beta.mgh']);
        [groupBetas_rh(:,k), M_rh, mr_parms_rh] = load_mgh([dataDir 'group_glms/' dataset '/Drift_' this_coef '/rh.Drift_' this_coef '.glmdir/beta.mgh']);
        groupSigs_lh(:,k) = load_mgh([dataDir 'group_glms/' dataset '/Drift_' this_coef '/lh.Drift_' this_coef '.glmdir/osgm/sig.mgh']);
        groupSigs_rh(:,k) = load_mgh([dataDir 'group_glms/' dataset '/Drift_' this_coef '/rh.Drift_' this_coef '.glmdir/osgm/sig.mgh']);
    end
    
    % Load the main effect group maps for masking
    groupMainEffect_lh = load_mgh([dataDir 'group_glms/' dataset '/DriftMainEffect/lh.DriftMainEffect.glmdir/beta.mgh']);
    groupMainEffect_rh = load_mgh([dataDir 'group_glms/' dataset '/DriftMainEffect/rh.DriftMainEffect.glmdir/beta.mgh']);

    % Extract the best coefficients at each vertex from the group maps
    [max_lh,groupBestCoef_lh] = max(groupBetas_lh,[],2);
    [max_rh,groupBestCoef_rh] = max(groupBetas_rh,[],2);
    
    % Don't include any vertices where the max beta is negative
    groupBestCoef_lh(max_lh<=0) = nan;
    groupBestCoef_rh(max_rh<=0) = nan;
    
    % Also don't include any vertices where main effect is negative
    groupBestCoef_lh(groupMainEffect_lh<=0) = nan;
    groupBestCoef_rh(groupMainEffect_rh<=0) = nan;
    
    % Determine significance levels
    % The groupBestCoefSig arrays are the same as the groupSigs arrays,
    % except that it has nan on any vertex where groupBestCoef has nan
    groupBestCoefSig_lh = zeros(size(groupBestCoef_lh));
    groupBestCoefSig_rh = zeros(size(groupBestCoef_rh));
    for i=1:163842
        if isnan(groupBestCoef_lh(i))
            groupBestCoefSig_lh(i) = nan;
        else
            groupBestCoefSig_lh(i) = groupSigs_lh(i,groupBestCoef_lh(i));
        end
        if isnan(groupBestCoef_rh(i))
            groupBestCoefSig_rh(i) = nan;
        else
            groupBestCoefSig_rh(i) = groupSigs_rh(i,groupBestCoef_rh(i));
        end
    end
    
    % Create a masked image of the best coefficients
    % Start as the full (unmasked) group map
    groupBestCoefMasked_lh = groupBestCoef_lh;
    groupBestCoefMasked_rh = groupBestCoef_rh;
    % Then, remove vertices that are below the significance threshold
    groupBestCoefMasked_lh(groupBestCoefSig_lh < -log10(sigThresh)) = nan;
    groupBestCoefMasked_rh(groupBestCoefSig_rh < -log10(sigThresh)) = nan;
    
    % Finally, save surface maps
    %if saveBool; save_mgh(groupBestCoef_lh,[outputDir 'groupBestCoef_' dataset '_lh.mgh'],M_lh,mr_parms_lh); end;
    %if saveBool; save_mgh(groupBestCoef_rh,[outputDir 'groupBestCoef_' dataset '_rh.mgh'],M_rh,mr_parms_rh); end;
    %if saveBool; save_mgh(groupBestCoefSig_lh,[outputDir 'groupBestCoefSig_' dataset '_lh.mgh'],M_lh,mr_parms_lh); end;
    %if saveBool; save_mgh(groupBestCoefSig_rh,[outputDir 'groupBestCoefSig_' dataset '_rh.mgh'],M_rh,mr_parms_rh); end;
    if saveBool; save_mgh(groupBestCoefMasked_lh,[outputDir 'groupBestCoefMasked_' sprintf('%2.2f',sigThresh) '_' dataset '_lh.mgh'],M_lh,mr_parms_lh); end;
    if saveBool; save_mgh(groupBestCoefMasked_rh,[outputDir 'groupBestCoefMasked_' sprintf('%2.2f',sigThresh) '_' dataset '_rh.mgh'],M_rh,mr_parms_rh); end;

end