function [x,eta,dict,B,N] = simplex(A,b,c)
    % This function solves the LP problem
    % maximize eta = c' * x, subject to Ax <= b, x >= 0,
    % using the simplex algorithm.
    % Outputs:
    % x: solution vector, x = [x_1,...,x_n]'.
    % eta: optimal value of objective function
    % dict: Optimal dictionary
    % B: Indices of basic variables in dict
    % N: Indices of nonbasic variables in dict
    % Test from section 5.6:  [x,eta,dict,B,N]=simplex([-2 -1;-2 4;-1 3],[4; -8; -7],[-1;-1]);
    
    [m,n]=size(A);
    assert(length(b) == m, 'mismatch in the sizes of b and A')
    assert(length(c) == n, 'mismatch in the sizes of c and A')
  
    if sum(b<0) == 0
        % We can make an initial feasible dictionary.
        dict = [0 c';b -A];
        B=[n+1:n+m]';
        N=[1:n]';
    else
        % Replace the objective with -x_1-...-x_n, to obtain a dual feasible dictionary
        [x_0,eta_0,dict_0,B_0,N_0]=simplex(-A',ones(n,1),-b);
        assert(eta_0 ~= inf, 'The problem is infeasible') % Unbounded dual implies infeasible primal
      
        % Convert dual dictionary to primal
        dict = -dict_0'; % The negative transpose
        I=[(n+1:n+m) (1:n)]'; % Used to map dual variables to their complementary primal variables, see the notes 
        B=I(N_0); N=I(B_0);
        % Replace the objective function by the original one, see the notes
        c_ext = [c;zeros(m,1)];
        CB = c_ext(B);
        CN = c_ext(N);
        dict(1,1) = CB'*dict(2:(m+1),1);
        dict(1,2:n+1) = CB'*dict(2:(m+1),2:(n+1)) + CN';
    end

    MaxIt=10000; it = 1;
    tol = 1.0e-12; % tolerance used for sign tests.
    while it < MaxIt
        inds = find(dict(1,2:end)>tol); % Choice set for entering variable
        if isempty(inds) % Dictionary is optimal.
            break;
        else
            [t,q] = min(N(inds)); % Blands rule for choosing entering variable
            r = inds(q); % r is entering
            
            ratios = -dict(2:end,r+1)./dict(2:end,1);
            % Apply the covention 0/0=0. Matlab computes this to NaN otherwise
            inds = find(( abs(dict(2:end,r+1)) < tol) & ( abs(dict(2:end,1)) < tol));
            ratios(inds)=0;
            
            val= max(ratios);
            if val<=tol % The problem is unbounded
                eta = inf; break;
            else
                inds = find(ratios==val); % Choice set for leaving variable
                [t,q] = min(B(inds)); % Blands rule for choosing leaving variable
                s = inds(q); % s is leaving
            
                dict=pivot0(dict,r,s);
                temp=B(s);
                B(s)=N(r);
                N(r)=temp;
                it = it+1;
            end
        end
    end

    assert(it < MaxIt, 'Maximum number of iterations exceeded.')

    x_ext = zeros(n+m,1);
    x_ext(B) = dict(2:end,1);
    x = x_ext(1:n);
    eta = dict(1,1);
end

function d1 = pivot0(d0,in,out)
    % See the notes
    r = in+1;  % r enters
    s = out+1; % s leaves
    d1      = d0 - d0(:,r) * d0(s,:) / d0(s,r);
    d1(:,r) =   d0(:,r) / d0(s,r);
    d1(s,:) =  -d0(s,:) / d0(s,r);
    d1(s,r) = 1 / d0(s,r);
end