%
% Preparation
%

%%% 1) Create a structure array specifying details of data
ACEfit_Par.P_nm    = 'HCP_imgs.txt';    % A list of images or other data
                                        % specification; see FileFormats.txt
                                        
ACEfit_Par.InfMx   = 'KinInf.csv';      % Kinship information matrix of 4
                                        % columns with headers
                                        % ('SubjectID','MotherID',
                                        % 'FatherID','Zygosity'). MotherID
                                        % and FatherID must be numeric,
                                        % and Zygosity is one of 'MZ',
                                        % 'NotMZ' and 'NotTwin'.  The order
                                        % of subjects must match the order
                                        % of image paths or subjects in
                                        % ACEfit_Par.P_nm.
ACEfit_Par.ResDir  = '/my/path/ResDir';

%%% The rest are optional; omit to use default values

ACEfit_Par.Pmask   = 'HCP_mask.nii';    % Brain mask image (default: whole volume)

ACEfit_Par.Dsnmtx  = '';                % Design matrix (default: all-ones vector)

ACEfit_Par.Nlz     = 1;                 % Inverse Gaussian normalisation options 
                                        % (applied to each phenotype voxel/element)
                                        %  0 - None
                                        %  1 - Gaussian normalisation *before* 
                                        %      forming residuals
                                        %  2 - Gaussian normalisation *after* 
                                        %      forming residuals; note, though, 
                                        %      that after this Gaussian 
                                        %      normalisation the residuals
                                        %      will be formed again to ensure
                                        %      orthognality of residuals to
                                        %      the design matrix.  

ACEfit_Par.AggNlz  = 0;                 % Aggregate heritability normalisation options
                                        % (applied to each phenotype voxel/element)
                                        %  0 - Default, meaning that only
                                        %      determined by Nlz option.
                                        %      This always entails
                                        %      de-meaning and removal of any effects
                                        %      in Dsnmtx, but variance at each
                                        %      voxel/element unchanged.  
                                        %  1 - Same as 0, but with variance
                                        %      normalisation; this is full
                                        %      'studentization'.
                                        %  2 - *Undo* mean centering of
                                        %      residual formation. (If Dsnmtx
                                        %      is empty or default value of
                                        %      an all-ones vector, this is
                                        %      equivalent to using the raw
                                        %      input data.  If Dsnmtx
                                        %      contains nuisance regressors
                                        %      the residuals will be formed
                                        %      and the mean added back in. 
                                        %  3 - Same as 2, but with
                                        %      variance normalisation;
                                        %      note that for each
                                        %      voxel/element, mean is
                                        %      first added back and then
                                        %      data are divided by stdev.

ACEfit_Par.ContSel = [];                % Select a single contrast (a volume
                                        % in a 4D Nifti file supplied for
                                        % each subject, or at the last
                                        % dimension in a cifti image; NOT
                                        % compatibile with a single file
                                        % containing all subjects' data.)
                                        
ACEfit_Par.NoImg   = 0;                 % If 1, suppress image-wise inference,
                                        % and only compute summaries.
                                        
                                        
%%% 2) Updata 'ACEfit_Par' with the input data information
ACEfit_Par = PrepData(ACEfit_Par);

%%% 3) Run the original data once
ACEfit_Par.alpha_CFT = [];              % Cluster-forming threshold (default: 0.05)

ACEfit_Par = ACEfit(ACEfit_Par);

%%% 4) Add permutation and bootstrapping information, and save "ACEfit_Par.mat"
ACEfit_Par.nPerm = 1000;                % Number of permutations
ACEfit_Par.nBoot = 1000;                % Number of bootstrap replicates
nParallel        = [];                  % Number of parallel runs (default: 1, without parallelization)

PrepParallel(ACEfit_Par,nParallel);


%
% Permutations
%

if ACEfit_Par.nPerm>0
    
    %%% 1)
    %%%%% The following code can be scripted or parallelized as you wish.  
    %%%%% Please refer to "README_APACE_intro.pdf" for the example snippets
    %%%%% for parallelization.
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    load(fullfile(ACEfit_Par.ResDir,'ACEfit_Par.mat'));
    RunID = 1;
    ACEfit_Perm_Parallel(ACEfit_Par,RunID);
    % Inside ResDir, it will create result sets, ACEfit_Parallel_XXXX.mat,
    % one for each RunID.
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    %%% 2) This will merge together all the results from the specified RunID's
    load(fullfile(ACEfit_Par.ResDir,'ACEfit_Par.mat'));
    ACEfit_Perm_Parallel_Results(ACEfit_Par);
    
    %%% 3)
    FWEalpha = 0.05;
    FDRalpha = 0.05;
    load(fullfile(ACEfit_Par.ResDir,'ACEfit_Par.mat'));
    ACEfit_Results(ACEfit_Par,FWEalpha,FDRalpha);
    
end


%
% Bootstrapping
%

if ACEfit_Par.nBoot>0
    
    %%% 1)
    %%%%% The following code can be scripted or parallelized as you wish.
    %%%%% Please refer to "README_APACE_intro.pdf" for the example snippets
    %%%%% for parallelization.
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    load(fullfile(ACEfit_Par.ResDir,'ACEfit_Par.mat'));
    RunID = 1;
    ACEfit_Boot_Parallel(ACEfit_Par,RunID);
    % Inside ResDir, it will create result sets, BootCI_Parallel_XXXX.mat,
    % one for each RunID.
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    %%% 2) This will merge together all the results from the specified RunID's
    load(fullfile(ACEfit_Par.ResDir,'ACEfit_Par.mat'));
    ACEfit_Boot_Parallel_Results(ACEfit_Par);
    
    %%% 3)
    Balpha = 0.05;
    load(fullfile(ACEfit_Par.ResDir,'ACEfit_Par.mat'));
    Boot_CIs(ACEfit_Par,Balpha)
    
end


%
% Aggregate heritability (aka "Steve's method"), with P-values via
% permutation and CI's via boostrapping.
%
% Note that permuation (steps 1-3) and bootstrapping (steps 1-3) can be
% skipped by setting ACEfit_Par.nPerm=0 and ACEfit_Par.nBoot=0 separately; 
% the following code can be run immediately after "PrepParallel". 
%
load(fullfile(ACEfit_Par.ResDir,'ACEfit_Par.mat'));
Palpha = 0.05; % Significance threhold for permutations (plotting only)
Balpha = 0.05; % Confidence level, where CI's have level 100*(1-alpha)

AgHe_Method(ACEfit_Par,Palpha,Balpha)
% % Once with no variance normalisation
% ACEfit_Par.AggNlz  = 0;   % de-meaning only
% AgHe_Method(ACEfit_Par,Palpha,Balpha,'_NoNorm')
% 
% % Now, again, but with Variance normalisation
% ACEfit_Par.AggNlz  = 1;   % de-meaning and scaling to have stdev of 1.0
% AgHe_Method(ACEfit_Par,Palpha,Balpha,'_Norm')


%
% Generate summary file 
%

APACEsummary(ACEfit_Par,'ResultSummary')


