function[A, subs, x, z] = psimplex(c, A, b, subs, varargin);
%PSIMPLEX - simplex algorithm for maximization problems in standard form
%
% [A, subs, x, z] = psimplex(c, A, b, subs, 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(4,5,nargin));

c = -c;
[m, n] = size(A);
b = b(:);
c = c(:)';

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 length(subs) ~= m
    error('psimplex: wrong size of basis');
end

if nargin == 5
   disp(sprintf('     ________________________________________________'))
   disp(sprintf('\n              Tableaux of the Simplex Algorithm\n'))
   disp(sprintf('     ________________________________________________'))
   disp(sprintf('\n      Initial problem, bringing to tableau form\n'))
   A
   disp(sprintf(' Press any key to continue ...\n'))
   pause
end

% do the initial pivots to bring to tableau form
rows = 1:m; subss = [];
for i = 1:m
    j = subs(i);
    y_j = A(1:m, j);

    ind = 1:m;
    rs = intersect(ind(y_j ~= 0), rows);
    if isempty(rs)
        error('psimplex: basis is not linearly indpendent');
    end
    r = rs(1);
    subss(r) = j;
    
    A(r,:)= A(r,:) ./ A(r,j);
    for k = 1:m+1
        if k ~= r
            A(k,:)= A(k,:)-A(k,j)*A(r,:);
        end
    end

%     disp(sprintf('     ________________________________________________'))
%     disp(sprintf('\n      Basis column %i: ready\n', j))
%     A
%     pause

    rows = setdiff(rows, r);
end

subs = subss;

if min( A(1:m, n+1) ) < 0
    error('psimplex: initial basis infeasible');
end

if nargin == 5
    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

[z, j] = Dpr(A(m+1,1:n));
while z < -eps
    y_j = A(1:m, j);
    if all(y_j <= 0)
        if nargin == 5
            disp(sprintf('     ________________________________________________'))
        end
        disp(sprintf('\n       Problem has unbounded objective function'));
        A
        x = zeros(n,1);
        z = inf;
        return;
    end
    r = MRT(A(1:m, n+1),A(1:m, j));
    
    A(r,:) = A(r,:) / A(r ,j);
    for i = 1:m+1
        if i ~= r
            A(i,:)= A(i,:)-A(i,j)*A(r,:);
        end
    end
    
    if nargin == 5
        disp(sprintf('     ________________________________________________\n'))
        disp(sprintf('         %d enters, %d leaves the basis', j, subs(r)))
        A
    end
    
    subs(r) = j;

    if nargin == 5
        disp(sprintf('         current basis:'))
        subs
        disp(sprintf(' Press any key to continue ...\n\n'))
        pause
    end
    
    [z, j] = Dpr(A(m+1, 1:n));
end

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

if nargin == 5
    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
