% The goal is to minimize z=c*x-d, subject to constraints Ax=b and x>=0.
% The inputs are c (a row vector), d (a scalar), A (a matrix), and b (a
% column vector).
% The outputs are x and z_min.

% Complementary program needed is pivot.m.


% Written by helga schaffrin, oct. 2003


function [x,z]=simplex(c, d, A, b)

[m,n]=size(A);

%if (length(c) ~= n || length(b) ~= m)
 %  error('Dimensions do not agree.')
%end

% Put the system into simplex canonical form.

i=1;
while(i<=m)
    [value, position]=max(sign(A(i,:)));
	% position records the column with the first non-zero entry in row i.
        % value lets us determine whether the row is identically 0.

    if value==0
       if b(i)~=0
          error('Infeasible constraint equation in row %.0f',i)
       else
          A_new=A([1:i-1 i+1:m],:);
              % We are deleting the i'th row, since it is redundant.
          b_new=b([1:i-1 i+1:m]);
          A=zeros(m-1,n);
          A=A_new;
          b=zeros(1,m-1);
          b=b_new;
          m=length(b);
       end
    else
       % If there is a first non-zero value in row i, we will use it to pivot.
       [d,c,A,b]=pivot(d,c,A,b,m,n,i,position);

       i=i+1;
    end
end


% Now we need to ensure that b has only non-negative entries.

% If all entries in b are negative, we need to pivot around a 
% negative entry in A to get a row with positive b entry.
% If this is not possible, there is no solution to our problem.

if max(sign(b))==-1
   % We find the first negative entry in A(1,:) and pivot there.
   [sgn_A, col]=min(sign(A(1,:)));
        % sgn_A indicates the sign of the least entry in A(1,:).
        % col gives the column with the first negative entry in A(1,:)
        % (or the first entry >=0 if no neg entries).

   if(sgn_A~=-1)
     error('There is no solution with x>=0.')
   end

   [d,c,A,b]=pivot(d,c,A,b,m,n,1,col);
end

% If some entries of b are still negative, we need to consider subproblems
% to remedy the situation.

while(min(sign(b))==-1)
   [small, row]=min(b);
                % small is the smallest entry in b (must be negative).
                % row is where this entry occurs; the corresponding row in A
                % becomes the objective function for the subproblem.

   % Find the first negative entry in A(row,:) to determine the pivoting 
   % column.

   [sgn, k]=min(sign(A(row,:)));
           % sgn indicates the sign of the least entry in A(row,:).
           % k gives the column with the first negative entry in A(row,:)
           % (or the first 0 entry if no neg entries).

   while(sgn<0)     %i.e. while there is a neg. entry in A(row,:)
        % Find the minimum ratio of b(i)/A(i,k), where i is such that
        % b(i)>=0, to determine the pivoting row.

        ratios=zeros(1,m);

        % If A(i,k)<=0 or b(i)<0, we record 10^10 for the ratio, so that it 
        % will later be disregarded (being too large).

        ratios = 10^10*(b<0) + 10^10*(A(:,k)<=0).*(b>=0) ...
                    + (A(:,k)>0).*(b>=0).*(b./(A(:,k)+(A(:,k)==0)));

        [minim, h]=min(ratios);
                  % minim is the smallest ratio.
                  % h is the pivoting row.

        if minim==10^10    %unbounded form -> pivot about A(row,k) if b(row)<0
           if b(row)<0
              [d,c,A,b]=pivot(d,c,A,b,m,n,row,k);
           end
           sgn=1;  %Because we are done with this row; b(row) is now >=0.
        else
            % Pivoting around A(h,k).
            [d,c,A,b]=pivot(d,c,A,b,m,n,h,k);
            
	    % Continuing with the optimizing...
            % Find the first negative entry in A(row,:) to determine the 
            % pivoting column.

            [sgn, k]=min(sign(A(row,:)));
        end
   end
end    



% And now we can start minimizing!

% Find the first negative entry in c to determine the pivoting column.

[sgn, k]=min(sign(c));
        % sgn indicates the sign of the least entry in c.
        % k gives the column with the first negative entry in c
        % (or the first 0 entry if no neg entries).

while(sgn<0)     %i.e. while there is a neg. entry in c
     % Find the minimum ratio of b(i)/A(i,k) to determine the pivoting row.

     ratios=zeros(1,m);

     % If A(i,k)<=0, we record 10^10 for the ratio, so that it will later be
     % disregarded (being too large).

     ratios=10^10*(A(:,k)<=0) + (A(:,k)>0).*(b./(A(:,k) + (A(:,k)==0)));

     [minim, h]=min(ratios);
               % minim is the smallest ratio.
               % h is the pivoting row.

     if minim==10^10          %unbounded form
        error('There is no minimum.')
     end


     % Pivoting around A(h,k).

     [d,c,A,b]=pivot(d,c,A,b,m,n,h,k);

     % The process starts over...
     % Find the first negative entry in c to determine the pivoting column.

     [sgn, k]=min(sign(c));
end
    
% The above loop ends when all entries in c are non-negative.  Now we can
% read off the solution.

x=zeros(1,n);

x=(c==0).*(max(A)==1).*(sum(A)==1).*(b'*A);

z=-d;

if max(abs(A*x'-b))>10^(-10)
    error('simplex is messing up again')
end