    %% jobIMA.m
%   This function describes an Osprey job defined in a MATLAB script.
%
%   AUTHOR:
%       Dr. Georg Oeltzschner (Johns Hopkins University, 2019-07-15)
%       goeltzs1@jhmi.edu
%
%   HISTORY:
%       2019-07-15: First version of the code.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% 1. SPECIFY SEQUENCE INFORMATION %%%

% Specify sequence type
seqType = 'unedited';           % OPTIONS:    - 'unedited' (default)
%             - 'MEGA'
%             - 'HERMES'
%             - 'HERCULES'

% Specify editing targets
editTarget = {'none'};          % OPTIONS:    - {'none'} (default if 'unedited')
%             - {'GABA'}, {'GSH'}, {'Lac'}, {'PE322'}, {'PE398'}  (for 'MEGA')
%             - {'GABA', 'GSH'}, {'GABA', 'Lac'}, {'NAA', 'NAAG'} (for 'HERMES'and 'HERCULES')

% Specify data scenario
dataScenario = 'invivo';       % OPTIONS:    - 'invivo' (default)
%             - 'phantom'
%             - 'PRIAM'
%             - 'MRSI'
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% 2. SPECIFY DATA HANDLING AND MODELING OPTIONS %%%
% Which spectral registration method should be used? Robust spectral
% registration is default, a frequency restricted spectral registration
% method is also availaible and is linked to the fit range.
opts.SpecReg = 'RobSpecReg';                  % OPTIONS:    - 'RobSpecReg' (default) Spectral aligment with Water/Lipid removal, using simialrity meric, and weighted averaging
%             - 'ProbSpecReg' Probabilistic spectral aligment to median target and weighted averaging
%             - 'RestrSpecReg' Frequency restricted (fit range) spectral aligment, using simialrity meric, and weighted averaging
%             - 'none'

% Which algorithm do you want to align the sub spectra? L2 norm
% optimization is the default. This is only used for edited MRS!
%Perform correction on the metabolite data (raw) or metabolite
%-nulled data (mm).
opts.SubSpecAlignment.mets = 'L2Norm';          % OPTIONS:    - 'L2Norm' (default)
%             - 'L1Norm'
%             - 'none'

%Perform eddy-current correction on the metabolite data (raw) or metabolite
%-nulled data (mm). This can either be done similar for all data sets by
%supplying a single value or specified for each dataset individually by supplying
% multiple entries (number has to match the number of datasets) e.g. to perform ECC
% for the second dataset only:
% opts.ECC.raw                = [0 1];
% opts.ECC.mm                = [0 1];


opts.ECC.raw                = 0;                % OPTIONS:    - '1' (default)
opts.ECC.mm                = 0;                 %             - '0' (no)
%             - [] array

% Save LCModel-exportable files for each spectrum?
opts.saveLCM                = 0;                % OPTIONS:    - 0 (no, default)
%             - 1 (yes)

% Save jMRUI-exportable files for each spectrum?
opts.savejMRUI              = 1;                % OPTIONS:    - 0 (no, default)
%             - 1 (yes)

% Save processed spectra in vendor-specific format (SDAT/SPAR, RDA, P)?
opts.saveVendor             = 0;                % OPTIONS:    - 0 (no, default)
%             - 1 (yes)

% Save processed spectra in NIfTI-MRS format?
opts.saveNII                = 0;                % OPTIONS:    - 0 (no, default)
%             - 1 (yes)

% Save PDF output for all Osprey modules and subjects?
opts.savePDF             = 1;                % OPTIONS:    - 0 (no, default)
%             - 1 (yes)

% Choose the fitting algorithm
opts.fit.method             = 'Osprey';       % OPTIONS:    - 'Osprey' (default)
%           - 'LCModel'

% Select the metabolites to be included in the basis set as a cell array,
% with entries separates by commas.
% With default Osprey basis sets, you can select the following metabolites:
% Ala, Asc, Asp, bHB, bHG, Cit, Cr, Cystat, CrCH2, EtOH, GABA, GPC, GSH, Glc, Gln,
% Glu, Gly, H2O, mI, Lac, NAA, NAAG, PCh, PCr, PE, Phenyl, sI, Ser,
% Tau, Tyros, MM09, MM12, MM14, MM17, MM20, Lip09, Lip13, Lip20.
% If you enter 'default', the basis set will include all of the above
% except for Ala, bHB, bHG, Cit, Cystat, EtOH, Glc, Gly, Phenyl, Ser, and Tyros.
opts.fit.includeMetabs      = {'default'};      % OPTIONS:    - {'default'}
%             - {'full'}
%             - {custom}

% Choose the fitting style for difference-edited datasets (MEGA, HERMES, HERCULES)
% (only available for the Osprey fitting method)
opts.fit.style              = 'Separate';   % OPTIONS:  - 'Concatenated' (default) - will fit DIFF and SUM simultaneously)
%           - 'Separate' - will fit DIFF and OFF separately

% Determine fitting range (in ppm) for the metabolite and water spectra
opts.fit.range              = [0.5 4];        % [ppm] Default: [0.5 4]
opts.fit.rangeWater         = [2.0 7.4];        % [ppm] Default: [2.0 7.4]

% Determine the baseline knot spacing (in ppm) for the metabolite spectra
opts.fit.bLineKnotSpace     = 0.4;              % [ppm] Default: 0.4.

% Add macromolecule and lipid basis functions to the fit?
opts.fit.fitMM              = 1;                % OPTIONS:    - 0 (no)
%             - 1 (yes, default)

% How do you want to model the co-edited macromolecules at 3 ppm for GABA-edited MRS?
opts.fit.coMM3              = '3to2MM';         % OPTIONS:    - {'3to2MM'} (default)
%             - {'3to2MMsoft'}
%             - {'1to1GABA'}
%             - {'1to1GABAsoft'}
%             - {'freeGauss'}
%             - {'fixedGauss'}
%             - {'none'}

opts.fit.FWHMcoMM3              = 14;

% Optional: In case the automatic basisset picker is not working you can manually
% select the path to the basis set in the osprey/fit/basis, i.e.:
% opts.fit.basisSetFile = 'osprey/fit/basis/3T/philips/mega/press/gaba68/basis_philips_megapress_gaba68.mat';
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% 3. SPECIFY MRS DATA AND STRUCTURAL IMAGING FILES %%
% When using single-average Siemens RDA or DICOM files, specify their
% folders instead of single files!

% Clear existing files
clear files files_ref files_w files_nii files_mm

% Data folder in BIDS format
te = 30; % set to 30 or 320

bids_folder = [pwd filesep 'BIDS'];
output_folder = [pwd filesep 'OUTPUT'];
statcsv = [output_folder filesep 'statfile.csv'];

db=readtable('../SCAN_DATABASE/scan_database.xlsx');

files = {};
files_ref = {};
files_nii = {};
fid = fopen(statcsv, 'w');
fprintf(fid, 'subject, session, location, te\n');

if te==30
    locations = {'ACC', 'PCC', 'RDLPFC'};
elseif te==320
    locations = {'ACC'};
else
    error('TE not recognised.');
end

counter = 1;
for sub = 1%db.subject.'
    disp(['Assigning subject ' num2str(sub)]);

    for ses = 1%:2
        disp(['    session ' num2str(ses)]);
        % find exam dir for this visit
        exam_dir = [bids_folder filesep ...
            'sub-' num2str(sub,'%03d') filesep 'ses-' num2str(ses,'%02d')];
        if ~(exist(exam_dir)==7)
            disp('    No session 2 data found.')
            continue;
        end

        for loc = {'ACC' 'PCC'} % locations
            disp(['        location ' loc{1}])
            disp(['            TE ' num2str(te)])
            % get series numbers from database (make sure scan exists/should be processed)
            series_ws = db{db.subject==sub, ['sn_ses_' num2str(ses) '_' loc{1} '_' num2str(te)]} - 1;
            series_ref = db{db.subject==sub, ['sn_ses_' num2str(ses) '_' loc{1} '_' num2str(te)]} - 2;
            if isnan(series_ws)
                disp(['Skipping subject ' num2str(sub) ' session ' num2str(ses) ' location ' loc{1} ' te ' num2str(te) '.']);
                continue;
            end
            % check data for this acquisition exists, then add
            % locations to files and files_ref cell arrays
            if te==30
                echo_idx=1;
            elseif te==320
                echo_idx=2;
            else
                error('TE not recognised.');
            end

            scan_dir = [bids_folder filesep ...
                'sub-' num2str(sub,'%03d') filesep ...
                'ses-' num2str(ses,'%02d') filesep ...
                'mrs' filesep ...
                'sub-' num2str(sub,'%03d') ...
                '_echo-' num2str(echo_idx) '_nuc-1H' ...
                '_voi-' loc{1} '_svs'];
            assert(exist(scan_dir)==7);
            files{counter} = scan_dir;

            ref_dir = [bids_folder filesep ...
                'sub-' num2str(sub,'%03d') filesep ...
                'ses-' num2str(ses,'%02d') filesep ...
                'mrs' filesep ...
                'sub-' num2str(sub,'%03d') ...
                '_echo-' num2str(echo_idx) '_nuc-1H' ...
                '_voi-' loc{1} '_ref'];
            assert(exist(ref_dir)==7);
            files_ref{counter} = ref_dir;

            anat_file = [bids_folder filesep ...
                'sub-' num2str(sub,'%03d') filesep ...
                'ses-' num2str(ses,'%02d') filesep ...
                'anat' filesep ...
                'sub-' num2str(sub,'%03d') ...
                '_T1w.nii' ];
            assert(exist(anat_file)==2);
            files_nii{counter} = anat_file;

            fprintf(fid, '%d, %d, %s, %d\n', sub, ses, loc{1}, te);
            counter = counter + 1;
        end
    end
end
fclose(fid);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% 4. SPECIFY STAT FILE %%%
% Supply location of a csv file, which contains possible correlation
% measures and group variables. Each column must start with the name of the
% measure. For the grouping variable use 'group' and numbers between 1 and
% the number of included groups. If no group is supplied the data will be
% treated as one group. (You can always use the direct path)

file_stat = statcsv;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% 5. SPECIFY OUTPUT FOLDER %%
% The Osprey data container will be saved as a *.mat file in the output
% folder that you specify below. In addition, any exported files (for use
% with jMRUI, TARQUIN, or LCModel) will be saved in sub-folders.

% Specify output folder (you can always use the direct path)
% (MANDATORY)
outputFolder = output_folder;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%