%
%   Generate Figure 2D (Drift curve) for FaceSpace paper
%   _______________________________________________
%   by Marcelo G Mattar (10/21/2013)
%   mattar@sas.upenn.edu

% Parameters
datasets = {'MDS_RGI'}; % Choose between 'MDS_RGI' and 'MDS_Multi'
hemispheres_to_include = {'rh'}; % Choose {'rh'}, {'lh'}, or {'rh','lh'}
masksize = 800; % Number of vertices to form ROI (within fusiform mask)
outputDir = './Figures/'; % Directory where output images will be saved
saveBool = false;

% Data location
dataDir = '~/Data/FaceSpace/';
% Subjects to include
subjectList_RGI = {...
    % [list of subject ID strings redacted]
    };
subjectList_Multi = {...
    % [list of subject ID strings redacted]
    };

% Include auxiliar functions
addpath('./auxFunctions');

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

for d = 1:length(datasets)
    % Determine which dataset we'll work on
    dataset = datasets{d};
    switch dataset
        case 'MDS_RGI'
            subjectList = subjectList_RGI;
            dataset_lowercase = '-mds_rgi';
        case 'MDS_Multi'
            subjectList = subjectList_Multi;
            dataset_lowercase = '-mds_multi';
    end
    numSubjects = length(subjectList);
    
    % Load masks
    mask_lh = load_mgh([dataDir 'ROIs/binaryMask_' dataset '_' num2str(masksize) '_lh.mgh']);
    mask_rh = load_mgh([dataDir 'ROIs/binaryMask_' dataset '_' num2str(masksize) '_rh.mgh']);
    mask_lh = logical(mask_lh);
    mask_rh = logical(mask_rh);
    
    % Initialize variables
    DriftBetas_lh = zeros(length(mask_lh),numSubjects,numCoef);
    DriftBetas_rh = zeros(length(mask_rh),numSubjects,numCoef);
    averageBetasROI = zeros(numSubjects,length(filename_endings));
    
    % Load data
    for k=1:numCoef
        DriftBetas_lh(:,:,k) = squeeze(load_mgh([dataDir 'group_glms/' dataset '/Drift_' filename_endings{k} '/lh.Drift_' filename_endings{k} '.fwhm10.mgh']));
        DriftBetas_rh(:,:,k) = squeeze(load_mgh([dataDir 'group_glms/' dataset '/Drift_' filename_endings{k} '/rh.Drift_' filename_endings{k} '.fwhm10.mgh']));
        
        % Calculate the average effect within ROI
        for s=1:numSubjects
            for h=1:length(hemispheres_to_include)
                eval(['averageBetasROI(s,k) = averageBetasROI(s,k) + mean(DriftBetas_' hemispheres_to_include{h} '(mask_' hemispheres_to_include{h} ',s,k));']);
            end
            % Now, take the mean across hemispheres (if looking at both together)
            averageBetasROI(s,k) = averageBetasROI(s,k) / length(hemispheres_to_include);
        end
    end
    
    % Calculate means and standard errors of the mean
    meanBetas = mean(averageBetasROI,1);
    steBetas = std(averageBetasROI,1)/sqrt(numSubjects);
    
    % Calculate statistics
    [~,driftPeak_index] = max(meanBetas);
    driftPeak = coefficient_values(driftPeak_index);
    %{
    [~,p_vsAdapt,ci_vsAdapt,stats_vsAdapt] = ttest(averageBetasROI(:,driftPeak_index),averageBetasROI(:,21),'tail','right');
    [~,p_vsPrototype,ci_vsPrototype,stats_vsPrototype] = ttest(averageBetasROI(:,driftPeak_index),averageBetasROI(:,1),'tail','right');
    %}
    % Display curve peak info
    display(['The peak of the drift curve is located at EC = ' num2str(driftPeak)]);
    
    % Plot the results
    figure(1);
    hold on;
    plot(coefficient_values,meanBetas,'k');
    plot(coefficient_values,meanBetas-steBetas,'c');
    plot(coefficient_values,meanBetas+steBetas,'c');
    ylim([0 0.4*10^(-3)])
    xlim([-0.1 1.1])
    if saveBool; print('-depsc',[outputDir 'Figure3_driftcurve_' dataset '.eps']); end;
    if saveBool; print('-dpng',[outputDir 'Figure3_driftcurve_' dataset '.png']); end;
end