% flux variablity analysis 
habitat = {'coral','seawater'};
methods = {'carveme','gapseq','kbase','consensus'};

modelDir = '~/COMMIT/data/gap-filling/iterative';
tableDir = '~/COMMIT/figures/exchanged_reactions_flux';

% medium that has been used for gap filling
mediumFile = 'media/M9-medium.mat';
LB_mediumFile = 'media/LB-medium.mat';

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

all_medium = [medium; LB_medium];

for i=1:numel(habitat)
    habitat{i}
    
    sink_max_100flux_t = table;
    sink_min_100flux_t = table;
    sink_min_100flux_v = [];
    sink_max_100flux_v = [];

    sink_max_80flux_t = table;
    sink_min_80flux_t = table;
    sink_min_80flux_v = [];
    sink_max_80flux_v = [];

    model_id = [];
       
    % load model workspace
    carveme = load(fullfile(modelDir, habitat{i},[methods{1}, '.mat']), 'GF');
    gapseq = load(fullfile(modelDir, habitat{i},[methods{2}, '.mat']), 'GF');
    kbase = load(fullfile(modelDir, habitat{i},[methods{3}, '.mat']), 'GF');
    consensus = load(fullfile(modelDir, habitat{i},[methods{4}, '.mat']), 'GF');

    for k = 1:numel(carveme.GF)

        % find export reactions        
        c_sink = carveme.GF{k}.rxns(contains(carveme.GF{k}.rxns,'sink'));
        g_sink = gapseq.GF{k}.rxns(contains(gapseq.GF{k}.rxns,'sink'));
        k_sink = kbase.GF{k}.rxns(contains(kbase.GF{k}.rxns,'sink'));
        con_sink = consensus.GF{k}.rxns(contains(consensus.GF{k}.rxns,'sink'));

        % find import reactions in common
        rxn_list = intersect(c_sink,g_sink);
        rxn_list = intersect(rxn_list,k_sink);
        rxn_list = intersect(rxn_list,con_sink);

        % remove metabolites from medium
        %rxn_list = rxn_list(~contains(rxn_list, all_medium));

        % remove sink_BIOMASS
        rxn_list = rxn_list(~contains(rxn_list, 'BIOMASS'));
   
        % FVA
        if ~isempty(rxn_list)
            [c_sink_min_flux, c_sink_max_flux] = fluxVariability(carveme.GF{k},100, 'max', rxn_list);
            [g_sink_min_flux, g_sink_max_flux] = fluxVariability(gapseq.GF{k},100, 'max', rxn_list);
            [k_sink_min_flux, k_sink_max_flux] = fluxVariability(kbase.GF{k},100, 'max', rxn_list);
            [con_sink_min_flux, con_sink_max_flux] = fluxVariability(consensus.GF{k},100, 'max', rxn_list);

            sink_min_100flux_v = [sink_min_100flux_v; c_sink_min_flux, g_sink_min_flux, k_sink_min_flux, con_sink_min_flux];
            sink_max_100flux_v = [sink_max_100flux_v; c_sink_max_flux, g_sink_max_flux, k_sink_max_flux, con_sink_max_flux];

            model_id = [model_id; cellstr(repmat(carveme.GF{k}.id,numel(rxn_list),1))];
        end
    end

    sink_max_100flux_t.CarveMe = sink_max_100flux_v(:,1);
    sink_max_100flux_t.gapseq = sink_max_100flux_v(:,2);
    sink_max_100flux_t.KBase = sink_max_100flux_v(:,3);
    sink_max_100flux_t.Consensus = sink_max_100flux_v(:,4);

    sink_min_100flux_t.CarveMe = sink_min_100flux_v(:,1);
    sink_min_100flux_t.gapseq = sink_min_100flux_v(:,2);
    sink_min_100flux_t.KBase = sink_min_100flux_v(:,3);
    sink_min_100flux_t.Consensus = sink_min_100flux_v(:,4);
    
    sink_max_100flux_t.modelID = model_id;
    sink_min_100flux_t.modelID = model_id;

    writetable(sink_max_100flux_t, fullfile(tableDir, [habitat{i}, '_', 'sink_max_100flux.csv']))
    writetable(sink_min_100flux_t, fullfile(tableDir, [habitat{i}, '_', 'sink_min_100flux.csv']))
end