function [xmin, fval, flag] = LinProg(f, AInEq, bInEq, AEq, bEq, debug)

% x = LinProg(f, AInEq, bInEq)
% x = LinProg(f, AInEq, bInEq,AEq, bEq)
% x = LinProg(f, AInEq, bInEq,AEq, bEq,debug)
% [x,fval] = LinProg(___)
% [x,fval,exitflag] = LinProg(___)

%This function implments the simplex matrix algorithm.
%It accepts row vector f defining the objective function as f*x
%It can accept only inequality constraints(as in x = LinProg(f, AInEq, bInEq)) ,
% or only equality constraints (as in x = LinProg(f, [], [], AEq, bEq).
% The debug defaults false if the use did not specify to see the stages
% during the solution process
%
%It runs both phase one and phase two automatically.
%
%The input is
%
%AInEq and bInEq: defined the inequality constrain AInEq*x <= bInEq
%
%AEq and bEq: defined the equality constrain AEq*x = bEq
%
%f: Vector. This is from minimize  F(x) = fx. As defined in
%   standard Matlab documentations.
%
%debug: flag. Set to true to see lots of internal steps.

%
% Version 7/16/2017( modification of original version by  Nasser M. Abbasi)
% by Lateef A. Kareem
% Free for use.

if nargin < 6
    debug = false;
end
if nargin < 5
    bEq = [];
end
if nargin < 4
    AEq = [];
end

InEq = ~isempty(AInEq);
Eq   = ~isempty(AEq);

if (InEq)
    validate_input(AInEq,bInEq,f);
    A = [AInEq, eye(size(AInEq, 1))];
    c = [f, zeros(1, size(AInEq, 1))];
    b = bInEq;
end

if (Eq)
    if (InEq)
        validate_input(AEq,bEq,f);
        p = size(A, 2) - size(AEq, 2);
        A = [A; [AEq, zeros(size(AEq, 1), p)]];
        b = [b; bEq];
    else
        A = AEq;
        c = f;
        b = bEq;
    end
end


for j=1:length(b)
    if b(j)<0
        A(j,:)=-A(j,:);
        b(j)=-b(j);
    end
end

[A,b,c] = make_phase_one(A,b,c,debug);
try
    [throw_away, x]  = SimplexPhase2(A,b,c,debug,'phase two');
    xmin    = x(1:numel(f));
    fval    = f*xmin;
    flag    = 'Solution was found';
catch
    xmin    = x(1:numel(f));
    fval    = f*xmin;
    flag    = 'Solution was not found: primal is unbounded';
end
end
%=================================================================
function [A,b,c]     = make_phase_one(A,b,c,debug)
[m,n]              = size(A);
tab                = zeros(m+2,n+m+1);
tab(1:m,1:n)       = A;
tab(end,n+1:end-1) = 1;
tab(1:m,end)       = b(:);
tab(m+1,1:n)       = c;
tab(1:m,n+1:n+m)   = eye(m);

if debug
    fprintf('Current tableau [phase one]\n');
    disp(tab);
end

for i = 1:m %now make all entries in bottom row zero
    tab(end,:) = tab(end,:)-tab(i,:);
end

A = tab(1:m+1,1:n+m);
b = tab(1:m+1,end);
c = tab(end,:);
tab = SimplexPhase1(A,b,c,debug,'phase one');

A = tab(1:m,1:n);
b = tab(1:m+1,end);
c = tab(m+1,1:n);
end
%=================================================================
function tab = SimplexPhase1(A,b,c,debug,phase_name)
[m,n]        = size(A);
tab          = zeros(m+1,n+1);
tab(1:m,1:n) = A;
tab(m+1,:)   = c;
tab(1:m,end) = b(:);

keep_running = true;
while keep_running
    if debug
         fprintf('***********************\n');
         fprintf('Current tableau [%s] \n',phase_name);
         disp(tab);
    end

    if any(tab(end,1:n)<0)%check if there is negative cost coeff.
         [throw_away,J] = min(tab(end,1:n)); %yes, find the most negative
         % now check if corresponding column is unbounded
         if all(tab(1:m-1,J)<=0)
           error('problem unbounded. All entries <= 0 in column %d',J);
         %do row operations to make all entries in the column 0
         %except pivot
         else
             pivot_row = 0;
             min_found = inf;
             for i = 1:m-1
                 if tab(i,J)>0
                      tmp = tab(i,end)/tab(i,J);
                      if tmp < min_found
                          min_found = tmp;
                          pivot_row = i;
                      end
                 end
             end
             if debug
                 fprintf('pivot row is %d\n',pivot_row);
                 fprintf('pivot col is %d\n',J);
             end
             %normalize
             tab(pivot_row,:) = tab(pivot_row,:)/tab(pivot_row,J);
             %now make all entries in J column zero.
             for i=1:m+1
                 if i ~= pivot_row
                      tab(i,:)=tab(i,:)-tab(i,J)*tab(pivot_row,:);
                 end
             end
         end
         if debug  %print current basic feasible solution
             fprintf('current basic feasible solution is\n');
             disp(get_current_x());
         end
    else
         keep_running=false;
    end
end

    %internal function, finds current basis vector
    function current_x = get_current_x()
         current_x = zeros(n,1);
         for j=1:n
             if length(find(tab(:,j)==0))==m
                 idx= tab(:,j)==1;
                 current_x(j)=tab(idx,end);
             end
         end
    end
end
%=================================================================
function [tab, x] = SimplexPhase2(A,b,c,debug,phase_name)
[m,n]        = size(A);
tab          = zeros(m+1,n+1);
tab(1:m,1:n) = A;
tab(m+1,1:n)   = c;
tab(1:m+1,end) = b(:);
inneriter = 0;

keep_running = true;
while keep_running
    if debug
         fprintf('***********************\n');
         fprintf('Current tableau [%s] \n',phase_name);
         disp(tab);
    end

    if any(tab(end,1:n)<0)%check if there is negative cost coeff.
         [throw_away,J] = min(tab(end,1:n)); %yes, find the most negative
         % now check if corresponding column is unbounded
         if all(tab(1:m,J)<=0)
           error('problem unbounded. All entries <= 0 in column %d',J);
         %do row operations to make all entries in the column 0
         %except pivot
         else
             pivot_row = 0;
             min_found = inf;
             for i = 1:m
                 if tab(i,J)>0
                      tmp = tab(i,end)/tab(i,J);
                      if tmp < min_found
                          min_found = tmp;
                          pivot_row = i;
                      end
                 end
             end
             if debug
                 fprintf('pivot row is %d\n',pivot_row);
                 fprintf('pivot col is %d\n',J);
             end
             %normalize
             tab(pivot_row,:) = tab(pivot_row,:)/tab(pivot_row,J);
             %now make all entries in J column zero.
             for i=1:m+1
                 if i ~= pivot_row
                      tab(i,:)=tab(i,:)-tab(i,J)*tab(pivot_row,:);
                 end
             end
         end
         x = get_current_x();
         if debug  %print current basic feasible solution
             fprintf('current basic feasible solution is\n');
             disp(x);
         end
    elseif(any(sum(abs(tab(1:m,1:n-m)))) > 1 && inneriter < m*5)
        % dealing with a stubbon artificial variable
        colmax = zeros(1,n-m);
        for i = 1:n-m
            colmax(i) = max(tab(:,i));
        end
        [throw_away,J] = max(colmax);
        pivot_row = 0;
         min_found = inf;
         for i = 1:m
             if tab(i,J)>0
                  tmp = tab(i,end)/tab(i,J);
                  if tmp < min_found
                      min_found = tmp;
                      pivot_row = i;
                  end
             end
         end
         if debug
             fprintf('pivot row is %d\n',pivot_row);
             fprintf('pivot col is %d\n',J);
         end
         %normalize
         tab(pivot_row,:) = tab(pivot_row,:)/tab(pivot_row,J);
         %now make all entries in J column zero.
         for i=1:m+1
             if i ~= pivot_row
                  tab(i,:)=tab(i,:)-tab(i,J)*tab(pivot_row,:);
             end
         end
    else
         x = get_current_x();
         keep_running=false;
    end
end

    %internal function, finds current basis vector
    function current_x = get_current_x()
         current_x = zeros(n,1);
         for j=1:n
             if length(find(tab(:,j)==0))==m
                 idx= tab(:,j)==1;
                 current_x(j)=tab(idx,end);
             end
         end
    end
end
%=================================================================
function validate_input(A,b,c)


%if ~ismatrix(A)
%    error('A must be matrix');
%end

if ~isvector(b)
    error('b must be vector');
end
if ~isvector(c)
    error('c must be vector');
end

[m,n]=size(A);
if length(b) ~= m
    error('b must have same size as number of rows of A');
end
if length(c) ~= n
    error('c must have same size as number of columns of A');
end
end