function [x0,Ac,ac,Nc,nc]=initial_condition(A,b)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   [x0,Ac,ac,Nc,nc]=initial_condition(A,b)                               %
%                                                                         %
%  It returns an initial feasible point and an adequate sets of active    %
% and inactive constraints. The function solve an auxiliary problem       % 
%                                                                         %
%                       min  J(x)=sum (bi-ai'*x), for i belong to V(k)    %
%                       subject to  ai'*x>=bi                             %
%                                                                         %
%   INPUTS:                                                               %
%       A  - matrix of the constraints (left side)                        %
%       v  - column vector of the constraints (right side)                %
%                                                                         %
%   OUTPUTS:                                                              %
%       x0 - initial feasible point                                       %
%       Ac - matrix of active constraints                                 %
%       ac - column vector of active constraints                          %
%       Nc - matrix of inactive constraints                               %
%       nc - column vector of inactive constraints                        %
%                                                                         %
%                                                                         %
%                                                                         %
%                                                 Last update: 25.10.2009 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

nx=size(A,1);                   % Number of variables x

Nc=A;                           % Creating of inactive set
nc=b;

Ac=eye(nx);                     % Creating of initial active set
ac=zeros(nx,2);                 % not. - first culomn contains right side od constraints
                                %        second culomn denotes type of constraint
                                %                      0 - pseudoconstraint
                                %                      1 - inequality

% --------- Prepairing of the active and inactive set to using ------------                                
%  If some constraint from active set is saved also in inactive set, it is
% removed from the inactive set and is change the atribut in active set
i=1;                            
while (i<=length(nc))
    if (nc(i)==0) && (abs(max(Nc(:,i))==1)) && (sum(abs(Nc(:,i))==1))
        for j=1:length(ac)
            if Ac(:,j)==A(:,i)
                ac(j,2)=1;
            end
        end
        Nc=remove_column(Nc,i);
        nc=remove_element(nc,i);
    else
        i=i+1;
    end
end
% -------------------------------------------------------------------------

xk=zeros(nx,1);                 % Creating of initial x - x(1)

z=1;                            % If z=2, xk is feasible, else the
                                % has to be executed.

if (ac(:,2)==ones(size(ac,1),1))
        x0=xk;
        ac=ac(:,1);
        z=2;
end

% ------------------------- Optimization algorithm ------------------------
while (z==1)
    if (length(nc)~=0)          % It there is some constraint in inactive set
        fi=nc-Nc'*xk;           % a filter for objective function calculation
        filter=ones(size(Nc,1),1)*((fi+abs(fi))./(2*fi))';% is prepared
        lambda=-inv(Ac)*abs((filter.*Nc)*ones(size(Nc,2),1));% and index q
        [lambdaq,q]=min(lambda);% is chosen
    else
        lambdaq=0;              % else lambda is 0.
    end
    
    
    if (lambdaq>=0)             % If there isn't any feasible descent direction
        x0=xk;                  % the algorithm terminate.
        [Ac,ac]=artific_remove(Ac,ac);
        break
    else
        Acp=inv(Ac)';
        sq=Acp(:,q);     
        as=[];                  % This matrix will contain all possible 
                                % alfa.  
        j=1;                    % The index pointer for as is set to 1 on 
                                % begin.
        for i=1:size(Nc,2)      % All inactive constraints are tested
%             if (filter(1,i)==0)
                ai=Nc(:,i);
                if (ai'*sq~=0)       % and if this condition holds alfa is                
                as(j,:)=[(nc(i,1)-ai'*xk)/(ai'*sq) i];

                                    % calculated and saved to as.
                    
                    j=j+1;          % and index pointer of as is increased. 
                end                 % The variable i is index pointer of Nc and
%             end
        end                     % and is stored in second column of matrix
                                % as.
        [alfa,p]=min(as(:,1));  % suitable alfa is chosen. 
        p=as(p,2);              % index of a constraint, which will be 


        xk=xk+alfa*sq;          % New feasible solution is calculated. 
        if (p~=0)               % If some inactive constraint was used
                                % The p-th column is removed for inactive
                                % (Nc,nc) set to the active one (Ac,ac).
            [Ac,ac,Nc,nc]=set_modificationLP(Ac,ac,Nc,nc,q,p);
        end
    end
end

function [Xn,xn]=artific_remove(X,x)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   [Xn,xn]=artific_remove(X,x)                                           %
%                                                                         %
%   The function remove all artificial constraints from variables X and x.%
%                                                                         %
%   INPUTS:                                                               %
%       X  - matrix created from left site of constraints                 %
%       x  - matrix - first column is made from right site  of constraints%
%            second one containts indicators:                             %
%                                   0 - artificial constraint             %
%                                   1 - inequality constraint             %               
%                                                                         %
%   OUTPUTS:                                                              %
%       Xn - modified matrix X                                            %
%       xn - vector created from first column of x after artificial const.%
%            removing.
%                                                                         %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
sx=size(x,2);
j=1;

Xn=[];
xn=[];

for i=1:sx
    if (x(i,2)~=0)
        Xn(:,j)=X(:,i);
        xn(j,1)=x(i,1);     
        j=j+1;
%     else
%         [c,s]=max(X(:,i));
%         X=remove_row(X,s);
%         Xn=remove_row(Xn,s);
    end
end

function [Ac,ac,Nc,nc]=set_modificationLP(Ac,ac,Nc,nc,q,p)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   [Ac,ac,Nc,nc]=set_modificationLP(Ac,ac,Nc,nc,q,p)                     %
%                                                                         %
%   The function change q-th culomn in matrix Ac and p-th culomn in matrix%
%   Nc. It does same with vector ac and nc. Pseudoconstraints are         %
%   removed from the problem once they come inactive.                     %
%                                                                         %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Xq=Ac(:,q);
xq=ac(q,1);
Xp=Nc(:,p);
xp=nc(p);

if (ac(q,2)~=0)                 % If it isn't pseudoconstraint
    Nc(:,p)=Xq;
    nc(p)=xq;
else
    Nc=remove_column(Nc,p);
    nc=remove_element(nc,p);
end

Ac(:,q)=Xp;
ac(q,:)=[xp 1];

function X=remove_row(X,i)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   X=remove_column(X,i)                                                  %
%                                                                         %
%   The function remove i-th row from the matrix X.                       %
%                                                                         %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
sx=size(X,1);           % Number of column in matrix X

if (i==1)
    if (sx>1)
        X=X(2:end,:);
    else
        X=[];
    end
elseif (i==sx)
    X=X(1:i-1,:);
else
    X=[X(1:i-1,:) X(i+1:end,:)];
end

function X=remove_column(X,i)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   X=remove_column(X,i)                                                  %
%                                                                         %
%   The function remove i-th culomn from the matrix X.                    %
%                                                                         %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
sx=size(X,2);           % Number of column in matrix X

if (i==1)
    if (sx>1)
        X=X(:,2:end);
    else
        X=[];
    end
elseif (i==sx)
    X=X(:,1:i-1);
else
    X=[X(:,1:i-1) X(:,i+1:end)];
end
    
function x=remove_element(x,i)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   x=remove_element(x,i)                                                 %
%                                                                         %
%   The function remove i-th member of the vector x.                      %
%                                                                         %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
sx=length(x);             % Number of member in vector x

if (i==1)
    if (sx>1)
        x=x(2:end);
    else
        x=[];
    end
elseif (i==sx)
    x=x(1:i-1);
else
    x=[x(1:i-1);x(i+1:end)];
end