%  simplx.m
%  solve the LP: 
%        max   ctilde'xtilde
%        s.t.  Atilde xtilde <= b, xtilde >= 0, 
%  by a highly inefficient but easy to understand version of the 
%  simplex method, where Atilde is m by n.   
%  Assume that b >= 0, so that xtilde = 0 is an initial feasible point.
%  Written by M. Overton for class Feb 2, 1998.
%  Warning: may not be fully debugged.
%
%  Start by introducing slack variables w, converting the LP to
%      max  ctilde'xtilde + 0'w 
%      s.t. Atilde xtilde + w = b, x >= 0, w >= 0.
%
%  Now rename x = [xtilde w], and A = [Atilde I], where I is the
%  identity matrix of order m, and so A is m by (n+m), so the LP is
%
%      max c'x 
%      s.t. A x = b, x >= 0.
%
%  Start with B, the basis matrix, set to I, and N, the
%  nonbasic part of A, set to the original Atilde.
%
%  Requires ctilde, Atilde, b and maxit (max number of simplex iterations)

[m,ntilde] = size(Atilde);
B = eye(m); N = Atilde;    % semicolon suppresses output
if size(b) ~= [m 1] | size(ctilde) ~= [ntilde,1],
   error('size mismatch')
end
cB = zeros(m,1); cN = ctilde;

for iter = 1:maxit,
% primal feasible initially by assumption and afterwards
% by construction: check this
   xB = B\b        % compute B^{-1}*b, i.e.
                   % solve system of equations B x = b
   if min(xB) < 0,
      error('x has neg component')
   end
   obj = cB'*xB
   BinvN = B\N     % the matrix B^{-1} N defines the dictionary, showing
                   % how the basic variables depend on the nonbasic variables
   y = BinvN'*cB - cN   % and the vector y shows how the objective function
                        % depends on the nonbasic variables
   
% check if optimal yet   
   [ymin, j] = min(y)
   if ymin >= 0,
      disp('stop: optimal')
      break;
   end
% j will enter the basis and increase the objective function, however
% we must see how the other basic variables  will change then.
% Ratio test: t is steplength, i is index of entering variable
   dxB = BinvN(:,j)
   [t,i] = stepsize(xB, dxB)     % function call: returns t >= 0
   if t == inf,                  % infinity is a valid number
      disp('stop: unbounded')
      break;
   end;
   if t == 0,
      disp('warning: degenerate')
   end;
% swap the roles of entering and leaving variables
   fprintf('\n nonbasic variable %d enters and basic variable %d leaves\n', j, i);
   temp_col = B(:,i);    temp_value = cB(i);
   B(:,i) = N(:,j);      cB(i) = cN(j);
   N(:,j) = temp_col;    cN(j) = temp_value;
   disp('press any key to continue')
   pause
end;
