function [EFM,supports,t] = shortestEFMs(network,k)

% ########################################################
% for a given metabolic network this function computes the k shortest EFMs.

% INPUT:
% struct called network with:
% S: Stoichiometric Matrix
% rxns: Reactions
% mets: Metabolites
% rev: vector of the reversible reactions, as an (0,1) vector, where 1 stands for: reaction is reversible
% description: Name of the network
% k is the number of EFMs which should be computed

% OUTPUT:
% the k-shortest EFMs of the network, saved as a mat file called "k_EFMs_description_fluxes"
% the indices of the corresponding supportsets, saved as a mat file called "k_EFMs_description_supports"
% t is the total time needed
% ########################################################

clear tic
clear toc
tic

tol = 1e-9;

%For performance reasons, it would be useful, to remove the blocked
%reactions
% first of all we compute all blocked reactions (we do this with the F2C2 Tool
% by L. David)
%[fctable, blocked,time_f2c2] = F2C2('glpk', network);

S = network.S; 
rev = network.rev;  % reversibility vector
M = 1000; % big M

% delete the blocked reactions from the network:
%S(:,logical(blocked)) = [];
%rev(logical(blocked)) = [];
%network.rxns(logical(blocked)) = [];

[m n] = size(S);

% split now the reversible reactions and define some numbers:
S_split = [S -S(:,logical(rev))];
which_rev = find(rev);
n_split = size(S_split,2);
n_vartotal = 2*n_split;
n_rev = nnz(rev);

% we have two kind of variables: v the fluxvector and a the corresponding
% binaries
v = [1:n_split];
a = [n_split+1:n_vartotal];

% first equation: Sv = 0
I = [S_split  zeros(m,n_split)];

% second equation: v - a*M <= 0
II = [eye(n_split) -M*eye(n_split)];

%third equation: a - v <= 0
III = [-eye(n_split) eye(n_split)];

%fourth equation: a_rev + a_rev' <= 1 (split reactions cannot appear in the same flux)
IVa = zeros(n_rev,n);
for i = 1:n_rev
    IVa(i,which_rev(i)) = 1;
end
IVb = eye(n_rev);
IV = [zeros(n_rev,n_split) IVa IVb];

%fifth equation: sum(a_i) >= 1 (at least one reaction should carry flux)
V = [zeros(1,n_split) ones(1,n_split)];


%now we construct the right hand side of the model:
rhsI = zeros(size(I,1),1);
rhsII = zeros(size(II,1),1);
rhsIII = zeros(size(III,1),1);
rhsIV = ones(size(IV,1),1);
rhsV = 1;

% what kind of constraints do we have:
senseI = ones(size(I,1),1);
senseI(:) = '=';
senseII = ones(size(II,1),1);
senseII(:) = '<';
senseIII = ones(size(III,1),1);
senseIII(:) = '<';
senseIV = ones(size(IV,1),1);
senseIV(:) = '<';
senseV = 1;
senseV(:) = '>';

%type of the variables:
vtype = ones(1,n_vartotal);
vtype(v) = 'C';
vtype(a) = 'B';

tosolve.obj = [zeros(1,n_split) ones(1,n_split)];
tosolve.A = sparse([I;II;III;IV;V]);
tosolve.rhs = [rhsI ; rhsII ; rhsIII ; rhsIV ; rhsV];
tosolve.sense = [char(senseI) ; char(senseII) ; char(senseIII) ; char(senseIV) ; char(senseV)];
tosolve.vtype = char(vtype);

EFM_split = compute_ith_EFM(tosolve,k,n_split,n_vartotal);

% resplit of the reversible reactions
EFM_split = EFM_split';
EFM_nonblocked = EFM_split(1:n,:);
EFM_rev = EFM_split(n+1:n_split,:);

EFM_minus = zeros(size(EFM_nonblocked));
for i = 1:n_rev
    EFM_minus(which_rev(i),:) = EFM_rev(i,:);
end

EFM_nonblocked = EFM_nonblocked - EFM_minus;

% re-insert the blocked reactions:
%EFM = zeros(size(S_orig,2),size(EFM_nonblocked, 2));
%nonblocked = setdiff(1:size(S_orig, 2), find(blocked));
%EFM(nonblocked,:) = EFM_nonblocked(1:size(S, 2),:);

%since we don't have excluded blocked reactions:
EFM=EFM_nonblocked; 

EFM(abs(EFM) <= tol) = 0;
% save the computed EFMs
filename = strcat(num2str(size(EFM_nonblocked, 2)), '_EFMs_of_', network.description,'_fluxes.mat');
save(filename, 'EFM')

disp('The EFM can be found in the file ');disp(filename);


%compute the supports:
supports=cell(1);
for i=1:size(EFM,2);
    %first the positive, then the negatie supports:
    %supports{i}=[find(EFM(:,i)>0); -find(EFM(:,i)<0)];
    
    %first, get the indices of the support, then add the sign of the
    %corresponding value
    supports_i=find(EFM(:,i));
    supports{i}= sign(EFM(supports_i,i)).*supports_i;   
end    
filename = strcat(num2str(size(EFM_nonblocked, 2)), '_EFMs_of_', network.description,'_supports.mat');
save(filename, 'supports')

disp('The supports can be found in the file ');disp(filename);


%sol = HelpFunctions.EFM_Test(EFM,S_orig,'yes');
%fprintf('%d are really EFMs.\n', nnz(sol));

t = toc;
end


function EFM_split = compute_ith_EFM(tosolve,k,n_split,n_vartotal)

% we will search now in every iteration for a new EFM. For that we need to
% tell the solver not to enumerate the EFMs which were computed in the
% previous steps (Z and rhsz)
% For that we need to put in each step one row to the matrix and one entry
% to the right hand side

opts.OutputFlag=0;  % Solver should not show its iterations
opts.ResultFile = 'shortestEFMs.lp';
opts.IntFeasTol = 1e-9;
opts.OptimalityTol = 1e-9;

Z = zeros(k,n_vartotal);
rhsZ = zeros(k,1);
senseZ = ones(size(Z,1),1);
senseZ(:) = '<';
EFM_split = zeros(k,n_vartotal);

[m_A n_A] = size(tosolve.A);
[m_rhs n_rhs] = size(tosolve.rhs);

tosolve.A = sparse([tosolve.A ; Z]);
tosolve.rhs = [tosolve.rhs ; rhsZ];
tosolve.sense = [char(tosolve.sense); char(senseZ)];

for i = 1:k
    result = gurobi(tosolve,opts);
    if strcmp(result.status,'INFEASIBLE')
        disp(' ');
        disp('#######################################');
        disp(' ');
        disp('There exist not more than ');disp(i-1);disp('EFMs. Or something else is wrong ...');
        disp(' ');
        disp('#######################################');
        disp(' ');
        break;
    end
    efm = result.x;
    EFM_split(i,:) = efm;   
    tosolve.A(m_A+i,:) = [zeros(1,n_split) efm((n_split+1):(n_vartotal))'];
    tosolve.rhs(m_rhs+i,:) = nnz(efm((n_split+1):(n_vartotal))) -1 ;    
end

% If there exists less then k EFMs the zeros need to be removed from the
% matrix:
if not(i==k)
    EFM_split(i:k,:) = [];
end
end