disp('----------------------------------------------------------------------')
disp('                       Iterative gap filling                          ')
disp('----------------------------------------------------------------------')

% load options script 
disp('Loading required data...')

options

% medium
load(mediumFile, 'medium')
load(LB_mediumFile, 'LB_medium')

% Universal database
load(dbFile)

% start parallel pool
delete(gcp('nocreate'))
parpool(ncpu);

for i=1:numel(experiments)

    for j = 1:numel(habitat)

        disp('----------------------------------------------------------------------')
        fprintf('%s models for OTU composition in %s dataset\n', habitat{j})
    
        % load model cell array
        model_workspace = fullfile(modelDir, [habitat{j}, '_consensus_draft_models_biomass.mat']);
    
        load(model_workspace)
    
        if ~exist('merged_models', 'var')
            merged_models = models; clear models
        end
     
        % load OTU abundances to take the subset found in the current study
        if exist('otuDir', 'var') && ~isempty(otuDir)
            otu_file = fullfile(otuDir, habitat{j}, 'otutab.txt');
            tab_merged = readAbundancesFromFile(otu_file);
            order = tab_merged;
        
            model_ids = cellfun(@(x)strtok(x.id), merged_models,...
                'UniformOutput', false);
            merged_models = merged_models(ismember(model_ids, tab_merged.Properties.RowNames));
            clear tab_merged
        elseif exist('otuDir', 'var') && ~isempty(otuDir) && ~exist('otuDir', 'dir')
            error('Given path to OTU abundance directory does not exist: %s', otuDir)
        end
    
        % create cell arrays for individual auxotrophic media
        % get the list of complementary media
        if exist('mediaDir', 'var') && ~isempty(mediaDir) && exist('mediaDir', 'dir')
            auxo_media = dir(fullfile(mediaDir, '*.tsv'));
            auxo_media = fullfile({auxo_media.folder}, {auxo_media.name})';
            auxo_media = auxo_media(contains(auxo_media, strcat(model_ids, '.tsv')));
        
            for k=1:numel(auxo_media)
                tmp_tab = importdata(auxo_media{k});
                tmp_tab = tmp_tab.textdata(2:end,1);
                auxo_media{k} = strcat(translateIDs(tmp_tab, 'met', [], 'ModelSEED',...
                    'MNXref', false), '[e]');
            end
        elseif exist('mediaDir', 'var') && ~isempty(mediaDir) && ~exist('mediaDir', 'dir')
            error('Given path to auxotrophic media directory does not exist: %s', mediaDir)
        else
            auxo_media = repelem({LB_medium}, numel(merged_models), 1);
        end
    
        % write single .mat files for every model in a subdirectory
        tmpModelDir = fullfile(modelDir, 'tmp_models', [experiments{i}, tmp_spec]);
        [s, fileList] = createModelDir(merged_models, tmpModelDir);
    
        clear merged_models
    
        if s
        
            fileList = fullfile(tmpModelDir, fileList);
        
            % run iterative gap filling
            [GF, EX, gf_order, solutions, exc, gf, bio] = ...
                iterativeGapFilling(fileList, medium, auxo_media,...
                dbModel_MNXref_balanced, weights, epsilon, include_sink,...
                order, iterations, seq_sim_workspace);

            outDir = fullfile('data/gap-filling/iterative', habitat{j},'/');

            save(fullfile(outDir, 'consensus_gap-filled.mat'),...
                'GF', 'EX', 'gf_order', 'solutions', 'exc', 'gf', 'bio')
        
            clear GF EX gf_order solutions exc gf bio
            rmdir(tmpModelDir, 's')
        end
    end
end
