function [Leaks, varargout] = bktr_Run(irrevGSM,GSM,matchRev,map, scale)
% function run the FBA on the prepared model and calculates Leak fluxes
% in:
% irrevGSM:   prepared irreversible model
% GSM:        prepared reversible model
% matchRev:   matcihng of the two models (cmp CBRA: convertToIrreversible)
% map:        Mapping structure (cmp bktr_mapping)
%
% out:
% Leaks:      Calculated leak fluxes
% AddUptakes: Secretion and Uptake rates

% debug
% table(irrevGSM.rxns, irrevGSM.lb, irrevGSM.ub,  irrevGSM.lb<= irrevGSM.ub)
% pause


%% Run FBA
MinimizedFlux = optimizeCbModel(irrevGSM,'min');

% Convert result back to reversible model
vRev = convertIrrevFluxDistribution(MinimizedFlux.x,matchRev);

if isempty(vRev)
    disp('No FBA solution found')
    Leaks=nan(size(map.DMets));
    varargout={[]};
    return
end

%% Collect Leaks

% find Excretion rxns and their excr metabolites
IdxExcrRxns=find(findExcRxns(GSM,1));
[IdxExcrMets,~]=find(GSM.S(:,IdxExcrRxns));

Leaks = zeros(length(map.DMets),1);
isLeak= false(length(IdxExcrRxns),1);

% all mapped rxns
mapped_rxns=~cellfun(@ischar, map.GSMRxns);
try
    mapped_rxns=unique([map.GSMRxns{mapped_rxns}]);
catch
    trp_list=cellfun(@transpose,map.GSMRxns(mapped_rxns),'UniformOutput',false);
    mapped_rxns=unique([trp_list{:}]);
end
% add lumped rxns if applicable
if isfield(map,'GSMLumpedRxns')
    mapped_rxns = unique([mapped_rxns map.GSMLumpedRxns(map.GSMLumpedRxns~=0)]);
end

for k=1:length(map.DMets)
    if ~ischar(map.GSMMets{k})
        % find Rxn excreting the metabolites mapped to the current dyn metabolite
        mapped_mets=[map.GSMMets{k}];
        Sw=full(GSM.S).*repmat(vRev,1,size(GSM.S,1))';
        Sw(:,mapped_rxns)=[];
        Leaks(k)=sum(sum(Sw(mapped_mets,:),2))*scale;
    else
        Leaks(k)=0;
    end
    
end

%% Other required Uptakes (optional)
idx_others =  abs(vRev(IdxExcrRxns))>1e-10; % exclude numerical errors
AddUptakes =  table(GSM.rxnNames(IdxExcrRxns(idx_others)),vRev(IdxExcrRxns(idx_others))*scale);

%% find "waste" terms of lumped reactions
mapped_mets=[map.GSMMets{~cellfun(@ischar, map.GSMMets)}];
waste=cell(1,length(map.DRxns));
for k=1: length(map.DRxns)
    if ~ischar(map.GSMRxns{k})
        idx_rxns = [map.GSMRxns{k}];
        if length(idx_rxns)>1
            weights=vRev(idx_rxns(:))';
            
            % calculte summed & weighted stoichiometry
            sumstoich = sum(full(GSM.S(:,idx_rxns)).*repmat(weights,length(GSM.mets),1),2);
     
            % which of the mets are part of the DDM?
            [reactants,~]=find(full(GSM.S(:,idx_rxns)));
            reactants = unique(reactants);
            [isinDM,order] = ismember(reactants,mapped_mets);
            [~,order]=sort(order);
            
            % complie waste table
            waste{k}= table(GSM.metNames(reactants(order)),sumstoich(reactants(order))*scale,...
                isinDM(order));
        else
            waste{k}=[];
        end
        
    end
end

varargout={AddUptakes, waste};
