function [description,x,eta,it,B,N] = simplex(A_in,b_in,c_in)
%SIMPLEX takes a m x n matrix A_in, an m- vector b_in and an n-vector c_in, and returns an optimal n-vector x and the corresponding
%o.f. value eta obtained through the simplex algorithm applied to the LP
%problem associated with A_in,b_in and c_in. 

MaxIt=10000;
it_start=0;

[m,n]=size(A_in);
if (size(b_in)~=m | size(c_in)~=n)
    msg='The dimensions of the data do not match.';
    error(msg)
end
B=[n+1:n+m]';
N=[1:n]';
if (sum(b_in<0)~=0)  %The initial dictionary is not feasible. We need to use phase I
    A_0=-A_in';
    b_0=ones(n,1);
    c_0=-b_in;
    [description_0,~,~,it_0,B_0,N_0] = simplex (A_0,b_0,c_0);
    it_start = it_0;
    B=N_0;
    N=B_0;
    % x(j) <---> z(m+j)    j=1...n
    %x(n+i)<--->z(i)         i=1...m
    for i=1:m
        if B(i)>m
            B(i)=B(i)-m;
        else
            B(i)=B(i)+n;
        end
    end
    for j=1:n
        if N(j)>m
            N(j)=N(j)-m;
        else
            N(j)=N(j)+n;
        end
    end
    if (strcmp(description_0,'unbounded'))
        description = 'unfeasible';
        x=inf;
        eta=inf;
        it=it_start;
        return;
    end
end
%The initial dictionary is feasible
A=[A_in,eye(m)];
b=b_in;
c=[c_in;zeros(m,1)];
x=zeros(n+m,1);
x(B)=inv(A(:,B))*b;   %x(B)=inv(A(:,B))\b
z=zeros(n+m,1);
z(N)=(inv(A(:,B))*A(:,N))'*c(B)-c(N);
for it = it_start : MaxIt
    A_B = A( : , B);
    A_N = A( : , N);
    c_B = c(B);
    c_N = c(N);
    x_B = x(B);
    x_N = x(N);
    z_B = z(B);
    z_N = z(N);
    %STEP 1
    if (sum(z_N<0)==0)
        eta=c'*x;
        description = 'feasible and bounded';
        return
    end
    %STEP 2
    j = find(z_N<0,1,'first');  % N(j) is the variable becoming basic
    %STEP 3
    e_j=zeros(n,1);
    e_j(j)=1;
    Dx_B = inv(A_B)*A_N*e_j;
    %STEP 4 & STEP 5
    if (sum(Dx_B>0)==0)
        description= 'unbounded';
        x=inf;
        eta=inf;
        return;
    end
    u = inf*ones(m,1);
    u(Dx_B>0) = x_B(Dx_B>0)./Dx_B(Dx_B>0);
    [ t , i ] = min(u); % B(i) is the variable which leaves the basis    
   %STEP 6
   e_i=zeros(m,1);
   e_i(i)=1;
   Dz_N = - (inv(A_B)*A_N)'*e_i;
   %STEP 7
   s=z_N(j)/Dz_N(j);
   %STEP 8
   x(N(j))=t;
   x(B)=x_B-t*Dx_B;
   z(B(i)) = s;
   z(N) = z_N-s*Dz_N;
   % STEP 9 : We interchange B(i) and N(j)
   temp = B(i);
   B(i)=N(j);
   N(j)=temp;
end

msg='Maximum number of iterations exceeded.'
disp(msg);

end
   
   
   
   
   
        
        
    

