function[A, subs, x, z] = psimplex(c, A, b, varargin);
%PSIMPLEX - simplex algorithm for maximization problems in standard form
%
% [A, subs, x, z] = psimplex(c, A, b, varargin)
%
% ---------------------------------------------------------------------------
% DESCRIPTION
% ---------------------------------------------------------------------------
% 
% The simplex algorithm for the LP problem
%
%                         max z = c*x
%                 Subject to: Ax  = b
%                              x >= 0
%
% The fifth input parameter is optional. If it is set to 'y', then the
% initial and final tableaux are displayed to the screen.
% ---------------------------------------------------------------------------
% INPUT
% ---------------------------------------------------------------------------
%
% A, b, c  - parameters of the LP
%
% subs     - initial basis indices
%
% varargin - if defined, display procession of the simplex algorithm
%
% ---------------------------------------------------------------------------
% OUTPUT
% ---------------------------------------------------------------------------
% A       - final tableau of the simplex method
% subs    - indices of the basic variables in the final tableau
% x       - optimal solution
% z       - value of the objective function at x.
%
% Copyright is with the following author(s):
%
% The primal simplex version is taken almost literally from the MATLAB 
% tutorial at 
%
% http://www.hkn.umn.edu/resources/files/matlab/MatlabTutorial.htm
%
% Copyright �2005 Eta Kappa Nu. All Right Reserved. 
%
% Improvements to the original code and the dual simplex is
%
% (C) 2009, Gabor Retvari, BMT-TMIT
%
% ---------------------------------------------------------------------------

error(nargchk(3,4,nargin));

c = -c;
[m, n] = size(A);
b = b(:);
c = c(:)';
% parameter to handle numerical uncertainity
tolerance = 0.000000001;
% variable to direct the procession display of the simplex algorithm
process = false;
if nargin == 4 
    process = true;
end

if size(c, 2) ~= n
    error('psimplex: wrong size of objective');
end

if size(b, 1) ~= m
    error('psimplex: wrong size of RHS');
end

A = [A b; c 0];


if process
   disp(sprintf('     ________________________________________________'))
   disp(sprintf('\n              Tableaux of the Simplex Algorithm\n'))
   A
   disp(sprintf(' Press any key to continue ...\n'))
   pause
end


% Here starts the two-phase simplex
    if process
        disp(sprintf('     ________________________________________________'))
        disp(sprintf('\n           Finding initial basis\n'))
    end
% First step is to eliminate the negative values from the RHS(if any)
    for i = 1:m
        if A(i, n+1) < 0 
            A(i, :) =  A(i, :) * -1;
        end
    end
% Next step is to introduce artificial variables and the new objective function
% AP1 will be initial simplex tableau of the first phase
    AP1 = [A(1:m, 1:n) eye(m) A(1:m, n+1); zeros(1, n) ones(1, m) 0];
    for i = 1:m
        AP1(m + 1, :) =  AP1(m + 1, :) - AP1(i, :);
    end
% bv contains the columns of the basis variables of the first phase
% at the begining they are the columns of the artificial variables
    bv = size(AP1, 2)-m : size(AP1, 2)-1;
% SIMPLEX function is used to calculate the optimal simplex tableau of the first phase
    [AP1, bv] = SIMPLEX(AP1, m, size(AP1, 2)-1, bv, process);
% if the optimal objective value is not zero then the original problem has no solution 
    if abs(AP1(size(AP1, 1), size(AP1, 2))) > tolerance
        error('psimplex: problem has no feasible solution');
    end
    A = [AP1(1:m, 1:n) AP1(1:m, size(AP1, 2)); c 0];
% do pivot operations until A become feasible simplex tableau
    for i = 1:size(bv,2)
        if c(bv(i)) ~= 0 
            col = A(1:m, bv(i))';
            prow = col * A(1:m, :);
            A(m+1, :) = A(m+1, :) - A(m+1, bv(i)) * prow;
        end
    end
    subs = bv;
    if process
        disp(sprintf('     ________________________________________________'))
        disp(sprintf('\n                First phase finished'))
    end


if process
    disp(sprintf('     ________________________________________________'))
    disp(sprintf('\n                  Initial tableau\n'))
    A
    disp(sprintf('\n                  Initial basis\n'))
    subs
    disp(sprintf(' Press any key to continue ...\n\n'))
    pause
end

[A, subs] = SIMPLEX(A, m, n, subs, process);

z = A(m+1, n+1);
x = zeros(1, n);
x(subs) = A(1:m, n+1);
x = x(1:n)';

if process
    disp(sprintf('     ________________________________________________'))
    disp(sprintf('                 Optimal tableau, end'))
end

end

function [m, j] = Dpr(d)
%DPR - Dantzig's steepest descent pivot rule
[m, j] = min(d);

end


function [row, mi] = MRT(a, b)

% The Minimum Ratio Test (MRT) performed on vectors a and b.
% Output parameters:
% row - index of the pivot row
% mi - value of the smallest ratio. 

m = length(a);
c = 1:m;
a = a(:);
b = b(:);
l = c(b > 0);
[mi, row] = min(a(l)./b(l));
row = l(row);

end

function [IST, bvs] = SIMPLEX(IST, row, col, bvs, show_process)
% The core simplex calculation function applied in both phases of the simplex algoritm
% Input parameters:
% IST - the initial simplex tableau
% row - number of equations
% col - number of variables
% bvs - the columns of the basis variables in the initial basis
% show_process - if true, then display procession of the simplex algorithm
% Output parameters:
% FST - the final simplex tableau
% bvs - the columns of the basis variables in the final basis
[z, j] = Dpr(IST(row+1,1:col));
while z < -eps
    y_j = IST(1:row, j);
    if all(y_j <= 0)
        if show_process
            disp(sprintf('     ________________________________________________'))
        end
        disp(sprintf('\n       Problem has unbounded objective function'));
        IST
        x = zeros(col,1);
        z = inf;
        return;
    end
    r = MRT(IST(1:row, col+1),IST(1:row, j));
    IST(r,:) = IST(r,:) / IST(r ,j);
    for i = 1:row+1
        if i ~= r
            IST(i,:)= IST(i,:)-IST(i,j)*IST(r,:);
        end
    end
    
    if show_process
        disp(sprintf('     ________________________________________________\n'))
        disp(sprintf('         %d enters, %d leaves the basis', j, bvs(r)))
        IST
    end
    
    bvs(r) = j;

    if show_process
        disp(sprintf('         current basis:'))
        bvs
        disp(sprintf(' Press any key to continue ...\n\n'))
        pause
    end
    
    [z, j] = Dpr(IST(row+1, 1:col));
end


end




