%    
% Priklad reseni linearni optimalizace simplexovou metodou
% 
% Volne podle Numerical Recipes in Fortran a Numerical Recipes in C,
% 2. vydani
%
% Bez zaruky!    
%
% Pokud nemate v Matlabu funkci goto, pouzijte soubor goto.m
%    

%%% Hlavni rutina
function dlp = demo_linprog()
%
% Zada ulohu, zavola simplexovou metodu a vypise vysledek 
%    
    
%% Zadani ulohy
priklad = 1;
    
if (priklad==1)
    % 2 promenne, 2 podminky
    % U =  x +  y
    %     3x +  y  <= 3      
    %      x + 3y  >= 1
    n=2;
    m1=1; m2=1; m3=0;
    a = [ 0   1   1;
          3  -3  -1;
          1  -1  -3 ];
    
elseif (priklad==2)
    % 3 promenne, 3 podminky, po 1 od kazdeho typu
    % U =   x +  y +  z
    %      3x +  y +  z  <= 3
    %       x + 3y -  z  >= 1
    %      2x      + 2z  =  1
    n=3;
    m1=1; m2=1; m3=1;
    a = [ 0   1   1   1;
          3  -3  -1  -1;
          1  -1  -3   1;
          1  -2   0  -2 ];

elseif (priklad==3)
    % 4 promenne, 4 podminky, priklad z Numerical Recipes
    % U =   x +  y + 3z - (1/2)w
    %       x      + 2z           <= 740
    %           2y      - 7w      <= 0
    %            y -  z + 2w      >= 1/2
    %       x +  y +  z +  w      =  9
    n=4;
    m1=2; m2=1; m3=1;
    a = [   0   1   1   3  -0.5;
          740  -1   0  -2   0;
            0   0  -2   0   7;
          0.5   0  -1   1  -2;
            9  -1  -1  -1  -1 ];

end

%% Konec zadani


% celkovy pocet podminek
m=m1+m2+m3;
disp(['*** Vstup ***']);
disp(['  Vstupni matice A:']);
disp(a);
% pridame radek pro pomocnou ucinkovou funkci U'
a = [a ; zeros(1,n+1)];
pause 


%% Volani simplexove metody

[a icase izrov iposv] = simplex(a,m,n,m1,m2,m3);

%% Vypis vystupu rutiny simplex()

disp(['*** Vystup funkce simplex ***']);
disp(['  case= ' num2str(icase)]);
disp(['  vystupni matice A:']);
disp(a(1:m+2,1:n+1))
disp(['  vektor izrov']);
disp(izrov(1:n))
disp(['  vektor iposv']);
disp(iposv(1:m))

%% Vypis vysledku 

disp(['*** Vypis vysledku ***']);
disp(['  Maximalni hodnota ucinkove funkce = ' num2str(a(1,1))]);
disp(['  -------------------------------------------- ']);
for j=1:n
    if (izrov(j)<=n) 
        disp(['  x(' int2str(izrov(j)) ') = 0'])
    end
end
for j=1:n
    if (iposv(j)<=n) 
        disp(['  x(' int2str(iposv(j)) ') = ' num2str(a(j+1,1))])
    end
end
    

disp([' -------------------------------------------- ']);
disp(['*** Vypis rozdilu v jednotlivych podminkach ***']);
for j=1:n
    if (izrov(j)>n) 
        disp(['rozdil v podm. ' int2str(izrov(j)-n) ' je 0'])
    end
end
for j=1:n
    if (iposv(j)>n) 
        disp(['rozdil v podm. ' int2str(iposv(j)-n) ' je ' num2str(a(j+1,1))])
    end
end

%pause
end
 

%%% Rutina s implementaci simplexove metody

function [a icase izrov iposv] = simplex(a,m,n,m1,m2,m3);

% Simplexova metoda pro reseni ulohy linearniho programovani.
% Vstupni parametry
%    n je pocet promennych
%    m1 je pocet podminek <=, m2 je pocet podminek >=, m3 pocet podm. =,
%    m je pocet vsech podminek, m = m1 + m2 + m3. 
%    Matice a(1:m+2,1:n+1)
%      V 1. radku zadavame ucinek, na 1. miste absolutni clen. 
%      V dalsich radcich jsou postupne podm. v poradi typ 1(<=), 2(>=), 3(=).  
%      V radku podminky na vstupu je vzdy nejdrive prava str. na 1. miste
%        a pak nasleduje N koeficientu podminky,
%        koeficienty se udavaji prevedene na pravou stranu, tj. vynasobene -1.
%      Radek m+2 neni na vstupu zadan a je vyhrazen pro pomocnou ucinkovou fci.
% Vystup
%   Prvek a(1,1) je maximum ucinku.
%   Vektor iposv:
%      pokud iposv(j) (1..m) = i<=n, pak promenna x(i) je levostranna
%                                    a jeji hodnota v optimu je a(j+1,1).
%      pokud iposv(j) (1..m) = k>n,  tak rozdil v (k-n)-te podmince je a(j+1,1).
%   Vektor izrov:
%      pokud izrov(j) (1..n) = i<=n, pak promenna x(i)=0 (pravostranna).
%      pokud izrov(j) (1..n) = k>n,  pak (k-n)-ta podminka je splnena presne.
%      reprezentovane nyni sloupcem j+1 matice a.
%   icase: = 0 pokud je nalezeno konecne reseni,
%          = 1 pokud je ucinek neomezeny,
%          =-1 pokud nelze vyhovet podminkam.
% Pokud iterace nekonverguje, lze zvetsit pozadovanou presnost EPS. 
    eps = 1e-6;
    
    l1 = zeros(1,n+1);
    l2 = zeros(1,m);
    l3 = zeros(1,m);
    nl1 = n;
    %% na zacatku udelame ze vsech promennych pravostranne
    for k=1:n
        l1(k)=k;
        izrov(k)=k;
    end
    nl2=m;
    %% a vsechny umele promenne levostranne
    for i=1:m
        if (a(i+1,1)<0.0)
            disp('Spatna vstupni tabulka');
            %% konstanty b_i museji byt kladne
            return;
        end
        l2(i) = i;
        iposv(i) = n+i;
    end
    for i=1:m2
        l3(i) = 1;
    end
    ir=0;
    %% ir znaci zda jsme ve fazi 2, tedy zda mame mozny vektor.
    %% Pokud uz pocatek je mozny vektor, jdeme na fazi 2.
    if (m2+m3~=0) 
        ir=1;
        %% Priznak ir=1 znamena ze musime provest fazi 1.
        for k=1:n+1
            %% Spocitame pomocnou ucinkovou funkci
            q1 = 0.0;
            for i=m1+1:m
                q1 = q1 + a(i+1,k);
            end
            a(m+2,k) = -q1;
        end

        while (ir~=0)
             %% Najdeme nejvetsi koeficient pomocne ucinkove fce
             [kp, bmax] = simp1(a,m+1,l1,nl1,0);
             if ((bmax<=eps)&&(a(m+2,1)<-eps))
                 icase = -1;
                 %% Pomocna u.f. je stale zaporna a nelze ji
                 %% vylepsit, tedy neexistuje mozny vektor
                 return;
             elseif ((bmax<=eps)&&(a(m+2,1)<=eps))
                 m12 = m1+m2+1;
                 %% Pomocna u.f. je nulova a nelze ji vylepsit.
                 %% Mame tedy mozny vektor. Pomoci "goto one"
                 %% zahodime vsechny umele promenne a jdeme na fazi 2.
                 if (m12 <= m) 
                     for ip=m12:m 
                         if (iposv(ip)==ip+n) 
                             [kp, bmax] = simp1(a,ip,l1,nl1,1);
                             if (bmax > 0.0)
                                 goto('one');
                             end
                         end
                     end
                 end
                 ir=0;
                 %% Nastavime priznak ze jsme dosahli faze 2
                 m12=m12-1;
                 if (m1+1<=m12)
                     for i=m1+1:m12
                         if (l3(i-m1)==1)
                             for k=1:n+1
                                 a(i+1,k) = -a(i+1,k);
                             end
                         end 
                     end    
                 end
                 break;
                 %% Jdeme na fazi 2
             end
             [ip,q1] = simp2(a,n,l2,nl2,kp);
             %% Najdeme pivota (hlavni prvek) - faze 1
             if (ip == 0) 
                 %% Maximum pomocne u.f. je neomezene, takze
                 %% neexistuje mozne reseni
                 icase = -1;
                 return;
             end
             %% Navratovy bod pro goto
             % Label one:
             a = simp3(a,m+1,n,ip,kp);
             %% Prohodime levostrannou a pravostrannou promennou
             %% a aktualizujeme seznamy
             if (iposv(ip)>=(n+m1+m2+1)) 
                 for k=1:nl1
                     if (l1(k)==kp), break, end
                 end
                 nl1 = nl1-1;
                 for is=k:nl1
                     l1(is) = l1(is+1);
                 end
                 a(m+2,kp+1) = a(m+2,kp+1) + 1.0;
                 for i=1:m+2
                     a(i,kp+1) = -a(i,kp+1);
                 end
             else 
                 if (iposv(ip)>=(n+m1+1)) 
                     kh=iposv(ip)-m1-n;
                     if (l3(kh)~=0) 
                         l3(kh) = 0;
                         a(m+2,kp+1) = a(m+2,kp+1) + 1.0;
                         for i=1:m+2
                             a(i,kp+1) = -a(i,kp+1);
                         end
                     end
                 end
             end 
             is = izrov(kp);
             izrov(kp) = iposv(ip);
             iposv(ip) = is;
        end % end of "while (ir~=0)"
        %% Pokud jsme stale ve fazi 1, vratime se na while.
    end  % end of "if (m2+m3~=0)"

    %% Konec faze 1, tedy hledani pocatecniho mozneho vektoru.
    %% Zacina faze 2, tedy optimalizace reseni.

    disp('zacinam fazi 2 (optimalizaci reseni)') 	

    while 1
        [kp, bmax] = simp1(a,0,l1,nl1,0);
        %% otestujeme prvni radek (ucinek) zda jsme hotovi
        if (bmax <= 0.0) 
            %% ano. tradaa!
            icase = 0;
            return;
        end
        [ip,q1] = simp2(a,n,l2,nl2,kp);
        %% najdeme pivota (hlavni prvek) - faze 2
        if (ip == 0) 
            %% ucinkova funkce neomezena. Koncime.
            icase = 1;
            return;
        end
        %% Prohodime levostrannou a pravostrannou promennou
        a = simp3(a,m,n,ip,kp);
        is = izrov(kp);
        izrov(kp) = iposv(ip);
        iposv(ip) = is;
        %% Jdeme na dalsi iteraci
    end
end
   
%%% Pomocne rutiny pro rutinu simplex
        
function [kp,bmax] = simp1(a, mm, ll, nll, iabf)
%
% Najde maximum z prvku s indexem obsazenym v seznamu l1.
% Pokud iabf==1, uvazuje absolutni hodnoty.
%    
    kp  = ll(1);
    bmax= a(mm+1,kp+1);
    for k=2:nll 
        if (iabf==0)
            test= a(mm+1,ll(k)+1) - bmax;
        else
            test= abs( a(mm+1,ll(k)+1) ) - abs( bmax );
        end
        if (test>0.0) 
            bmax= a(mm+1,ll(k)+1);
            kp  = ll(k);
        end    
    end
end
        

function [ip q1] = simp2(a, n, l2, nl2, kp)
%
% Najde pivota (hlavni prvek), uvazuje i degenerovane pripady
%
    eps = 1e-6;
    ip = 0;
    for i=1:nl2
        if (a(l2(i)+1,kp+1) < -eps), break, end
        %% nejaci mozni pivoti?
    end
    if (i>nl2), return, end
    q1 = -a(l2(i)+1,1) / a(l2(i)+1,kp+1);
    ip = l2(i);
    for i=i+1:nl2 
        ii = l2(i);
        if (a(ii+1,kp+1) < -eps) 
            q = -a(ii+1,1) / a(ii+1,kp+1);
            if (q < q1) 
                ip = ii;
                q1 = q;
            elseif (q == q1) 
                %% mame degenerovany pripad
                for k=1:n
                    qp = -a(ip+1,k+1) / a(ip+1,kp+1);
                    q0 = -a(ii+1,k+1) / a(ii+1,kp+1);
                    if (q0~=qp), break, end
                end			
                if (q0<qp)
                    ip = ii;
                end
            end		
        end
    end
end

    
function  [a] = simp3(a, i1,k1,ip,kp);
%
% Maticove operace k prohozeni levostranne a pravostranne promenne    
%    
    piv = 1.0 / a(ip+1,kp+1);
    for ii =1:i1+1
        if (ii-1~=ip)
            a(ii,kp+1) = a(ii,kp+1)*piv;
            for kk=1:k1+1 
                if (kk-1~=kp)
                    a(ii,kk) = a(ii,kk) - a(ip+1,kk)*a(ii,kp+1);
                end
            end
        end
    end
    for kk=1:k1+1
        if (kk-1 ~=kp)
            a(ip+1,kk) = -a(ip+1,kk)*piv;
        end
    end
    a(ip+1,kp+1) = piv;
end
        
%%% 

