function [maticeSpojuDleCeny,VsechnySpojeDleCeny] = OptDleCeny (BW, BWcel, jednotkaBW, pocatek, konec, pocetuzlu, cenaza1m, cenaza1b, BWpl, cenaza1pa, cenaza1pb, delka, limit)

%m file urcen pro vykresleni zakladni podoby site a optimalizace v podobe
%nalezeni minimalni kostry. Zaroven k nalezene minimalni kostre dodame
%redundantni spoje pro uzly pripojene pouze jednou linkou, vybereme ty
%nejmene nakladne.
%celkovacena - vystup funkce CenaZaSpoj
%pocatek - vektor pocatecnich bodu spoju
%konec - vektor koncovych bodu spoju
%pocetuzlu - pocet uzlu site
%BWpl - kapacita spoje pristupove linky, zadavejte v Bps, Mbps, Gbps
%BW - kapacita jednotlivych linek
%BWcel - celkova kapacita spoju do uzlu, zadavejte v Bps, Mbps, Gbps
%cenaza1pa - zadejte jednotkovou cenu za prepojeny paket
%cenaza1pb - zadejte jednotkovou cenu za pametovy blok
%delka - delka spoje, zadavejte v metrech
%cenaza1m, cenaza1b - zadejte jednotkovou cenu delky a kapacity
%[maticeSpojuDleCeny,VsechnySpojeDleCeny] = OptDleCeny (BW, BWcel, jednotkaBW, pocatek, konec, pocetuzlu, cenaza1m, cenaza1b, BWpl, cenaza1pa, cenaza1pb, delka, limit)

disp('******************PARAMETRY OPTIMALIZACE DLE CENY SPOJU*************************')

if nargin ~= 13
    error('Spatny pocet vstupu')
end

delkaPocatek = length(pocatek);
delkaKonec = length(konec);
delkaDelka = length(delka);
delkaBW = length(BW');
delkaBWpl = length(BWpl');

if delkaPocatek ~= delkaKonec || delkaPocatek ~= delkaDelka || delkaPocatek ~= delkaBW
    error('Pocet prvku vektoru pocatek, konec, delka a BW musi byt stejny')
elseif delkaBWpl ~= pocetuzlu
    error('Pocet prvku vektoru BWpl musi byt stejny jako pocet uzlu')
end

[celkovacena,cenalinek,cenauzlu] = CelCenaSpoje (BW, BWcel, jednotkaBW, pocetuzlu, cenaza1m, cenaza1b, BWpl, cenaza1pa, cenaza1pb, delka);

spoje = sparse(pocatek,konec,cenalinek,pocetuzlu,pocetuzlu); %vytvoreni matice spoju
VsechnySpojeDleCeny = spoje + zeros(pocetuzlu,pocetuzlu);
grafCelk = biograph(spoje,[], 'Label', 'Zakladni navrh dle ceny','Description', 'ZakladniNavrh','EdgeType', 'straight','ShowArrows','off','LayoutScale', 2.5); %vychozi navrh
view(grafCelk)

doltrojmatice = tril(spoje + spoje'); %prevod na dolni trojuhelnikovou matici, nutnost pro dalsi vypocet
[MK,pred] = graphminspantree(doltrojmatice); %vypocet minimalni kostry vzhledem k ohodnoceni
grafOpt = biograph(MK,[], 'Label', 'Minimalni kostra navrhu dle ceny','Description', 'MinimKostra','EdgeType', 'straight','ShowArrows','off'); %optimalizovany navrh
view(grafOpt)

disp ('Spoje minimalni kostry grafu jsou:')
MKmoje = MK' %pouzite spoje v min. kostre

disp ('Nevyuzite spoje puvodniho navrhu jsou:')
nevyuzitespoje = spoje - MKmoje %nevyuzite spoje v min kostre

%Redundanci doplnim tak, ze k uzlum se stupnem 1 hledeam spolecne cesty,
%pokud takove neexistuji, hledam minima do uzlu a testuji na stupen

%Otestovat matici dodanych spoju podle Opt dle zpozdeni

redundance = MKmoje + zeros(pocetuzlu,pocetuzlu); %matice hodnot min kostry
volneSpoje = nevyuzitespoje + zeros(pocetuzlu,pocetuzlu); %nevyuzite spoje
kontrola = zeros(pocetuzlu,pocetuzlu); %pro testovani stupne uzlu
VsechnySpojeDleCenyPR = VsechnySpojeDleCeny - redundance; %rozdilova matice vsech spoju a tech pouzitych v navrhu min kostry

redundanceStupen = redundance;
redundanceStupen(redundanceStupen > 0) = ones(size(find(redundanceStupen > 0))); %nahrad prvky vetsi nez 0 jednickou

%zpusob jak doplnovat spoje, ktere patri spolecne uzlum se
%stupnem 1
for i = 1:pocetuzlu
kontrola(i,i) = sum(redundanceStupen(i,:)) + sum(redundanceStupen(:,i)); %nenulove hodnoty prevede v danem sloupci a radku na 1, vytvori jejich soucet a dle toho se urci stupen grafu
end

stupen1 = find(any(kontrola == 1)); %hledam vsechny sloupce, kde se vyskytuje stupen 1
[x,y] = size(stupen1);
nalezeneSpolSpoje = zeros(pocetuzlu,pocetuzlu);

if y > 1
for i = 1 : y-1    
    nalezeneSpolSpoje(stupen1(1,i),stupen1(1,i+1)) = sum(volneSpoje(stupen1(1,i),stupen1(1,i+1))); %matice spolecnych spoju, ktere mohou byt pridany 
end
end

redundance = redundance + nalezeneSpolSpoje; %matice spoju po pridani spolednych

% kontrola pro pripad, ze se nanasly spolecne spoje
redundanceStupen = redundance;
redundanceStupen(redundanceStupen > 0) = ones(size(find(redundanceStupen > 0))); %nahrad prvky vetsi nez 0 jednickou
for i = 1:pocetuzlu
kontrola(i,i) = sum(redundanceStupen(i,:)) + sum(redundanceStupen(:,i)); %nenulove hodnoty prevede v danem sloupci a radku na 1, vytvori jejich soucet a dle toho se urci stupen grafu
end

for i = 1:pocetuzlu    
    
    if kontrola(i,i) < 2 %pokud ma nektery uzel stupen mesi nez 2   
        
            VsechnySpojeDleCenyPR(VsechnySpojeDleCenyPR == 0) = Inf*ones(size(find(VsechnySpojeDleCenyPR == 0)));
            index1 = min(VsechnySpojeDleCenyPR(i,:)); %hledej minimum v patricnem radku matice smazanych hodnot
            index2 = min(VsechnySpojeDleCenyPR(:,i)); %hledej minimum v patricnem sloupci matice smazanych hodnot
            
            %musi zkouset ,jestli pridani minimalniho spoje nezapricini
            %pridani zbytecne mnoha dalsich spoju
            
            if index1 > index2 && index2 < Inf
            index3 = index2;                       
            elseif index1 < index2 && index1 < Inf
            index3 = index1;
   %         elseif index1 == index2 && index1 || index2 ~= Inf
    %        index3 = index1;    
            end
            
            L = VsechnySpojeDleCenyPR == index3; 
            najdi = find(L); %pozice minima v matici smazanych hodnot
            redundance(L) = index3*ones(size(najdi)); %prirazeni hodnoty do matice spoju
                        
            VsechnySpojeDleCenyPR(VsechnySpojeDleCenyPR == Inf) = 0*ones(size(find(VsechnySpojeDleCenyPR == Inf)));
                    
     end
      
end

%kontrola na zpetne spoje z horni poloviny uzlu. Takze vezmu matici
%prvku a vyberu sloupce od poloviny nahoru a prvni polovinu radku
zbyleSpoje = VsechnySpojeDleCeny - redundance;
zbyleSpojeVysec = zbyleSpoje(1:(pocetuzlu/2),(pocetuzlu/2):pocetuzlu);
redundanceStupen = redundance(1:(pocetuzlu/2),(pocetuzlu/2):pocetuzlu);
redundanceStupen(redundanceStupen > 0) = ones(size(find(redundanceStupen > 0))); %nahrad prvky vetsi nez 0 jednickou

%mozna jeste zapremyslet, jestli by nebzlo vhodne vyclenit z prohledavanych
%sloupcu ten posledni

for i = 1:pocetuzlu/2
    if sum(redundanceStupen(:,i)) > 0
        continue
    else
        zbyleSpojeVysec(zbyleSpojeVysec == 0) = Inf*ones(size(find(zbyleSpojeVysec == 0)));
        index1 = min(zbyleSpojeVysec(i,:)); %hledej minimum v patricnem radku matice smazanych hodnot
        index2 = min(zbyleSpojeVysec(:,i)); %hledej minimum v patricnem sloupci matice smazanych hodnot
            
        %musi zkouset ,jestli pridani minimalniho spoje nezapricini
        %pridani zbytecne mnoha dalsich spoju
            
        if index1 > index2 && index2 < Inf
        index3 = index2;              
        elseif index1 < index2 && index1 < Inf && index1 == index2
        index3 = index1;
       % elseif index1 == index2 && index1 || index2 ~= Inf
        %index3 = index1;    
        end
            
        L = zbyleSpoje == index3; 
        najdi = find(L); %pozice minima v matici smazanych hodnot
        redundance(L) = index3*ones(size(najdi)); %prirazeni hodnoty do matice spoju
                       
        zbyleSpojeVysec(zbyleSpojeVysec == Inf) = 0*ones(size(find(zbyleSpojeVysec == Inf)));
    end
end
        

redundance = sparse(redundance);
maticeSpojuDleCeny = redundance + zeros(pocetuzlu,pocetuzlu); %promenna pouzita v celkovem navrhu

grafOptSRedun = biograph(redundance,[], 'Label', 'Minimalni kostra s redundanci navrhu dle ceny','Description', 'MinimKostraSRedun','EdgeType', 'segmented','ShowArrows','off','LayoutScale', 3.0); %navrh s redundanci
view(grafOptSRedun)

disp ('Nevyuzite spoje dodane pro redundanci:')
dodanespoje = redundance - MKmoje; %redundantni spoje
dodanespoje(dodanespoje < 0) = 0*ones(size(find(dodanespoje < 0))) %redundantni spoje

%zpodeni-zjistit jake je a vypsat hotnoty u minimalni kostry

fixnizpozdeni = 1.4*(10^-3); % fixni zpozdeni v uzlech pro paket velikosti 1500B
zpozdenilinky = (delka./10^5)*(500*10^-6); %zpozdeni linek pri sireni signalu rychlosti 500mikro s/100km
packet = 1500*8;

switch jednotkaBW
    case 'Bps'
    serializacniZpozdeni = (BW.\packet).*2; %serializacni zpozdeni na zdrojovem a cilovem uzlu transportni site
    zpozdeniPL = BWpl.\packet; %serializacni zpozdeni vzhleem k pristupova lince
    
    case 'Mbps'
    serializacniZpozdeni = ((BW.*10^6).\packet).*2;
    zpozdeniPL = (BWpl.*10^6).\packet;
    
    case 'Gbps'
    serializacniZpozdeni = ((BW.*10^9).\packet).*2;
    zpozdeniPL = (BWpl.*10^9).\packet;
    
    otherwise
        disp('zadej spravnou jednotku BW> Bps, Mbps, Gbps')
      
end

zpTransportSite = serializacniZpozdeni + zpozdenilinky; %celkove zpozdeni linek transportni site

MKmojePZ = MKmoje;%pro potraby testu zpozdeni
sitDleZpozdeni = sparse(pocatek,konec,zpTransportSite,pocetuzlu,pocetuzlu); %vytvoreni matice spoju dle zpozdeni
NNsdz = MKmojePZ > 0;%musim najit spoje v matici spoju min kostry a dle nich vybrat ty z matice dle zpozdeni
hledejNNsdz = find(NNsdz);
MKmojePZ(hledejNNsdz) = sitDleZpozdeni(hledejNNsdz);

UG = tril(MKmojePZ + MKmojePZ');%dolni trojuhel matice

%musim urcit zpozdeni vsech cest v siti, zatim bez fixniho zpozdeni
maticeZpozdeni = graphallshortestpaths(UG,'directed',false); %matice sousednosti na zaklade doby celkoveho zpozdeni kazde cesty v transport. siti návrhu dle min kostry.
zpozdeniPL = zpozdeniPL';


for i = 1:pocetuzlu
    
    %musim pricist zpozdeni vychozi PL a vsech cilovych PL
    maticeZpozdeni(i,:) = (maticeZpozdeni(i,:) + zpozdeniPL(1,i)) + zpozdeniPL(1,:); %cyklus pro pricteni zpozdeni pristupovych linek a vytvoreni matice celkoveho zpozdeni
    maticeZpozdeni(i,i) = 0;
    
end

%musim urcit, pres kolik uzlu to jde a podle toho pricist fixni zpozdeni

for i = 1:(pocetuzlu) %opakovani pro zdrojove uzly
    
    for j = 1: pocetuzlu %opakovani pro cilove uzly
        [dist,path,pred] = graphshortestpath(UG,i,'directed',false); %hledani vsech nejkratsich tras ze zdrojoveho uzlu
        nejrychlejsiTrasa = path{j}; % nejrzchlejsi cesta k cili trasa k danemu cili
        pocetMezUzlu = length(nejrychlejsiTrasa);%pocet mezilehlych uzlu
        maticeZpozdeni(i,j) = maticeZpozdeni(i,j) + (pocetMezUzlu * fixnizpozdeni);%podle poctu mezilehlych uzlu prictu fixni zpozdeni
        
    end
end

%odstanim cesty k uzlum samotnym 1 do 1 atd
for i = 1:pocetuzlu
    maticeZpozdeni(i,i) = 0;
end

%disp('Zpozdeni jednotlivych linek mezi koncovymi ucastniky v [s] v navrhu minimalni kostry mezi koncovymi uzly:')
maticeZpozdeni = triu(maticeZpozdeni);
zpozdeni = sparse(maticeZpozdeni);

% V pripade ze dojde k prekroceni limitu zpozdeni, tak urcit co to je za spoj, o kolik a
% jak to resit
zpozdeneLinky = zeros(pocetuzlu, pocetuzlu);
Delay = zpozdeni > limit; %hledej hodnoty presahujici limit
najdiDelay = find(Delay); %hledej pozice v matici zpozdeni linek
zpozdeneLinky(najdiDelay) = zpozdeni(najdiDelay) - limit; %prirazeni hodnoty do matice zpozdenych linek a jeji zmena na hodnotu o kterou je prekrocen limit
disp('Linky, ktere prekracuji maximalni limit zpozdeni navrhu minimalni kostry a jakym rozdilem v [s]')
zpozdeneLinky = sparse(zpozdeneLinky)       

zpozdeni = zpozdeni(:);
NNzp = zpozdeni > 0;%hledam nenulove prvky a temi budu delit
hledejNNzp = find(NNzp);
[x,y] = size(hledejNNzp);
average = sum(zpozdeni)/x; %prumerna doba zpozdeni linky
maximum = max(zpozdeni);

disp(sprintf('Prumerna doba zpozdeni cest v minimalni kostre = %9.9f s \n',average(1,1)));
disp(sprintf('Maximalni doba zpozdeni cest v minimalni kostre = %9.9f s \n',maximum(1,1)));

%pro redundanci
redundancePZ = redundance; %pro potraby testu zpozdeni
NNsdz = redundancePZ > 0;%musim najit spoje v matici spoju min kostry a dle nich vybrat ty z matice dle zpozdeni
hledejNNsdz = find(NNsdz);
redundancePZ(hledejNNsdz) = sitDleZpozdeni(hledejNNsdz);

doltrojmatice2 = tril(redundancePZ + redundancePZ'); %prevod na dolni trojuhelnikovou matici, nutnost pro dalsi vypocet
maticeZpozdeniSredun = graphallshortestpaths(doltrojmatice2,'directed',false); %matice sousednosti na zaklade doby celkoveho zpozdeni kazde cesty v transport. siti návrhu dle min kostry s redundanci

for i = 1:pocetuzlu
    
    maticeZpozdeniSredun(i,:) = (maticeZpozdeniSredun(i,:) + zpozdeniPL(1,i)) + zpozdeniPL(1,:); %cyklus pro pricteni zpozdeni pristupovych linek a vytvoreni matice celkoveho zpozdeni
    maticeZpozdeniSredun(i,i) = 0;
    
end

%musim urcit, pres kolik uzlu to jde a podle toho pricist fixni zpozdeni

for i = 1:(pocetuzlu) %opakovani pro zdrojove uzly
    
    for j = 1: pocetuzlu %opakovani pro cilove uzly
        [dist,path,pred] = graphshortestpath(doltrojmatice2,i,'directed',false); %hledani vsech nejkratsich tras ze zdrojoveho uzlu
        nejrychlejsiTrasa = path{j}; % nejrychlejsi cesta k cili trasa k danemu cili
        pocetMezUzlu = length(nejrychlejsiTrasa);%pocet mezilehlych uzlu
        maticeZpozdeniSredun(i,j) = maticeZpozdeniSredun(i,j) + (pocetMezUzlu * fixnizpozdeni);%podle poctu mezilehlych uzlu prictu fixni zpozdeni
        
    end
end

%odstanim cesty k uzlum samotnym 1 do 1 atd
for i = 1:pocetuzlu
    maticeZpozdeniSredun(i,i) = 0;
end

%disp('Zpozdeni jednotlivych linek mezi koncovymi ucastniky v [s] mezi koncovymi uzly v navrhu s redundanci:')
maticeZpozdeniSredun = triu(maticeZpozdeniSredun);
zpozdeniSredun = sparse(maticeZpozdeniSredun)

% V pripade ze nastane hodnota presahujici limit zpozdeni, tak urcit co to je za spoj, o kolik a
% jak to resit
zpozdeneLinkySRedun = zeros(pocetuzlu, pocetuzlu);
Delay = zpozdeniSredun > limit; %hledej hodnoty presahujici limit
najdiDelay = find(Delay); %hledej pozice v matici zpozdeni linek
zpozdeneLinkySRedun(najdiDelay) = zpozdeniSredun(najdiDelay) - limit; %prirazeni hodnoty do matice zpozdenych linek a jeji zmena na hodnotu o kterou je prekrocen limit
disp('Linky, ktere prekracuji maximalni limit zpozdeni v navrhu s redundanci a jakym rozdilem v [s]')
zpozdeneLinkySRedun = sparse(zpozdeneLinkySRedun)       

zpozdeniSredun = zpozdeniSredun(:);
NNzp = zpozdeniSredun > 0;%hledam nenulove prvky a temi budu delit
hledejNNzp = find(NNzp);
[x,y] = size(hledejNNzp);
average = sum(zpozdeniSredun)/x; %prumerna doba zpozdeni linky
maximum = max(zpozdeniSredun);

disp(sprintf('Prumerna doba zpozdeni cest v navrhu s redundanci = %9.9f s \n',average(1,1)));
disp(sprintf('Maximalni doba zpozdeni cest v navrhu s redundanci = %9.9f s \n',maximum(1,1)));

%%%toky v navrhu s redundanci

sitDleToku = sparse(pocatek,konec,BW,pocetuzlu,pocetuzlu); %vytvoreni matice spoju dle zpozdeni

redundancePT = redundance; %pro potraby testu zpozdeni
NNsdz = redundancePT > 0;%musim najit spoje v matici spoju min kostry a dle nich vybrat ty z matice dle zpozdeni
hledejNNsdz = find(NNsdz);
redundancePT(hledejNNsdz) = sitDleToku(hledejNNsdz);

zdroje = pocetuzlu - 1;
tok = zeros(pocetuzlu,pocetuzlu); %budouci matice toku

for i = 1:zdroje
[MaxTok,TokVeSpoji,Rez] = graphmaxflow(redundancePT,i,pocetuzlu); %hleda max. toky
      tok = tok + TokVeSpoji; % a vytvari novou matici spoju na zaklade souctu nalezenych max. toku
     
end

disp(sprintf('Toky ve spojich navrhu s redundanci po vypoctu max. toku: %s',jednotkaBW))
toky = sparse(tok) %zobrazeni matice souctu max. toku

%%%%%%%Naklady

%Urceni celkove vstupni kapacity z transportni site do jednotlivych uzlu v
%navrhu bez redundance
BWCelOptNavrhu = sparse(pocatek,konec,BW,pocetuzlu,pocetuzlu);
ind = find(MKmoje <= 0);
BWCelOptNavrhu(ind) = 0*(ind);
BWCelOptNavrhu = BWCelOptNavrhu + zeros(pocetuzlu,pocetuzlu);

for i = 1:pocetuzlu
    
    BWcel(i,1) = sum(BWCelOptNavrhu(i,:)) + sum(BWCelOptNavrhu(:,i)); %nove hodnoty toku do uzlu
    
end

[cenapaketu, cenapambloku, cenacelkemuzlu] = CenaZaUzel (BWpl, BWcel,cenaza1pa, cenaza1pb, jednotkaBW); %vypocet cen uzlu v danem navrhu s minimalni kostrou

sumavychozi = sum(celkovacena);
disp (sprintf('Celkove ohodnoceni puvodniho navrhu plneho propojeni: %9.9f\n', sumavychozi));

sumaspoju = sum(MKmoje);%cena spoju
sumaspoju = sum(sumaspoju);
sumaspoju = sumaspoju(1,1);
sumauzlu = sum(cenacelkemuzlu);%cena uzlu
cenaoptnavrhu = sumaspoju + sumauzlu;
disp (sprintf('Celkove ohodnoceni optimalizovaneho navrhu dle minimalni kostry: %9.9f\n', cenaoptnavrhu));

%Urceni celkove vstupni kapacity z transportni site do jednotlivych uzlu v
%navrhu s redundanci
BWCelOptNavrhuSR = sparse(pocatek,konec,BW,pocetuzlu,pocetuzlu);
ind = find(redundance <= 0);
BWCelOptNavrhuSR(ind) = 0*(ind);
BWCelOptNavrhuSR = BWCelOptNavrhuSR + zeros(pocetuzlu,pocetuzlu);

for i = 1:pocetuzlu
    
    BWcelSR(i,1) = sum(BWCelOptNavrhuSR(i,:)) + sum(BWCelOptNavrhuSR(:,i)); %nove hodnoty toku do uzlu
    
end

[cenapaketu, cenapambloku, cenacelkemuzlu] = CenaZaUzel (BWpl, BWcelSR,cenaza1pa, cenaza1pb, jednotkaBW); %vypocet cen uzlu v danem navrhu s redundanci

sumalinekSR = sum(redundance);%cena spoju
sumalinekSR = sum(sumalinekSR);
sumalinekSR = sumalinekSR(1,1);
sumauzluSR = sum(cenacelkemuzlu);%cena uzlu
cenaoptnavrhusredun = sumalinekSR + sumauzluSR;
disp (sprintf('Celkove ohodnoceni optimalizovaneho navrhu s redundanci: %9.9f\n', cenaoptnavrhusredun));

rozdilnakladusMKsRed = cenaoptnavrhusredun - cenaoptnavrhu;
disp(sprintf('Rozdil navrhu s redundanci a minim kostry je: %9.9f\n',rozdilnakladusMKsRed));

rozdilnakladusMK = cenaoptnavrhu - sum(celkovacena);
disp(sprintf('Rozdil stavu plneho propojeni oproti minimalni kostre je: %9.9f\n',rozdilnakladusMK));

rozdilnakladusSR = cenaoptnavrhusredun - sum(celkovacena);
disp(sprintf('Rozdil stavu plneho propojeni oproti navrhu s redundanci je: %9.9f\n',rozdilnakladusSR));


end