function [description,x,eta,it,B,N,dict] = simplex(A,b,c)
% This routine takes a matrix A (m x n), a vector b (m) and a vector c (n) and
% returns an optimal vector x (n) and the corresponding optimal value eta,
% to the LP problem
%
% maximize eta = c' * x,
% subject to Ax <= b, x >= 0,
%
% using the simplex algorithm.
%
% Outputs:
%
% If you just want the answer:
%
% description: the problem either has a solution or is unfeasible or unbounded. 
% x: solution vector, x = [x_1,...,x_n]'.
% eta: optimal value of objective function
%
% For experts:
%
% it: number of iterations used (total of both phases if 2 phases required)
% B: Index set of basic variables in optimal dictionary.
% N: Index set of non-basic variables in optimal dictionary.
% dict: optimal dictionary, dimension (m+1)*(n+1)
% with basic variables x_{B(1)},...,x_{B(m)}, and
% non-basic variables x_{N(1)},...,x_{N(n)}.
%
% Author: Michael Floater
% e-mail: michaelf@math.uio.no
% Created: 04 March 2020, using Matlab R2019b
% Modified: 11 March 2020. Added tolerance of 1.0e-12 for the two sign tests.
% The code now works on all data sets I have tried,
% including linear L1 regression to 1000 points (where m=2000, n=1004).
% Copyright 2020

MaxIt=10000;
[m,n]=size(A);
if length(b) ~= m || length(c) ~= n
  disp('Error. The dimensions of your data do not match.');
  return;
end
it=0;
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
  % We apply the 2-phase method.
  % In Phase I we solve a dual problem to get a feasible dictionary
  % for the primal problem.
  A_0=-A';
  b_0=ones(n,1);
  c_0=-b;
  [description_0,x_0,eta_0,it_0,B_0,N_0,dict_0]=simplex(A_0,b_0,c_0);
  it=it_0
  if strcmp(description_0,'unbounded')
    description='unfeasible';
    x=inf;
    eta=inf;
    dict_0
    return
  end
  % Convert dual dictionary to primal
  dict = -dict_0'; % The negative transpose
  I=[[n+1:n+m],[1:n]]';
  B=I(N_0);
  N=I(B_0);
  % Replace the objective function by the original one.
  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

tol = 1.0e-12; % tolerance used for sign tests.
while it < MaxIt
  l = 0;
  cmax = tol;
  for j=1:n
    if dict(1,j+1) > cmax
      cmax = dict(1,j+1);
      l = j;
    end
  end
  if l == 0
  % Dictionary is optimal.
    description='feasible and bounded';
    break;
  end
  k = 0;
  rmin = 1.0e50;
  for i=1:m
    if dict(i+1,l+1) < -tol
      r = - dict(i+1,1) / dict(i+1,l+1);
      if r < rmin
        rmin = r;
        k = i;
      end
    end
  end
  if k == 0
  % The LP problem is unbounded.
    description='unbounded';
    dict
    break;
  end

  dict=pivot(dict,l,k);
  temp=B(k);
  B(k)=N(l);
  N(l)=temp;
  it = it+1;
end

if it == MaxIt
  description='Maximum number of iterations exceeded.';
  return;
end

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