%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                        %
%                              WiMAX_logMAP                              %
%                                                                        %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%logMAP DEKODOVACI ALGORITMUS
%Vstupne premenne:
%Lc_yk - prenesene data [double]
%g - vytvaracie polynomy kodera [bin]
%L_uk - a-priori pravdepodobnost z predchadzajuceho dekodera [double]
%index - urcuje poradie dekodera [int]
%naOut - nasledujuci vystup kodera [int]
%naStav - nasledujuci stav [int]
%preOut - predchadzajuci vystup [int]
%preStav - predchadzajuci stav [int]
%Vystupne premenne:
%L_uky - a-posteriori vystup dekodera [double]
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [ L_uky ] = WiMAX_logMAP( Lc_yk,g,L_uk,index, naOut, ... 
    naStav, preOut, preStav )

%%INICIALIZACIA HODNOT
dlz = length(Lc_yk)/3;
[n,K] = size(g); 
m = K - 1;
max_stav = 2^m;
Inf = 1e10;
Alpha=0;
Beta=0;
gamma=0;
stav2=0;
tempmax=0;
temp0=0;
temp1=0;
stav1=0;

    %%INICIALIZACIA ALPHA
    Alpha(1,1) = 0; 
    Alpha(1,2:max_stav) = -Inf*ones(1,max_stav-1);

    %%INICIALIZACIA BETA
    if index==1
       Beta(dlz,1) = 0;
       Beta(dlz,2:max_stav) = -Inf*ones(1,max_stav-1); 
    elseif index==2
       Beta(dlz,1:max_stav) = zeros(1,max_stav);
    else
       error('CHYBA!!!');
    end

    %%VYPOCET ALPHA PRECHOD TRELISOM V DOPREDNOM SMERE
    for k = 2:dlz+1
        for stav2 = 1:max_stav
            gamma = -Inf*ones(1,max_stav);
            a = -Lc_yk(3*k-5);
            b = Lc_yk(3*k-4:3*k-3);
            c = preOut(stav2,2:3)';
            %GAMMA PRE 0
            gamma(preStav(stav2,1)) = (a+b*c)....     
                -log(1+exp(L_uk(k-1)));
            d = preOut(stav2,5:6)';
            %GAMMA PRE 1
            gamma(preStav(stav2,2)) = (Lc_yk(3*k-5)+Lc_yk(3*k-4:3*k-3)*d)...
                +L_uk(k-1)-log(1+exp(L_uk(k-1)));

            %VYPOCET ALPHA
            if(sum(exp(gamma+Alpha(k-1,:)))<1e-300)
                Alpha(k,stav2)=-Inf;
            else
                Alpha(k,stav2) = log( sum( exp( gamma+Alpha(k-1,:) ) ) );  
            end   
        end
        tempmax(k) = max(Alpha(k,:));
        Alpha(k,:) = Alpha(k,:) - tempmax(k);
    end     

    %%VYPOCET BETA PRECHOD TRELLISOM V SPATNOM SMERE
    for k = dlz-1:-1:1
        for stav1 = 1:max_stav
            gamma = -Inf*ones(1,max_stav);
            a = -Lc_yk(3*k+1);
            b = Lc_yk(3*k+2:3*k+3);
            c = naOut(stav1,2:3)';
            %GAMMA PRE 0
            gamma(naStav(stav1,1)) = (a+b*c)....
                -log(1+exp(L_uk(k+1)));
            %GAMMA PRE 1
            gamma(naStav(stav1,2)) = (Lc_yk(3*k+1)+Lc_yk(3*k+2:3*k+3)*naOut(stav1,5:6)')....
                +L_uk(k+1)-log(1+exp(L_uk(k+1)));
            
            %VYPOCET BETA 
            if(sum(exp(gamma+Beta(k+1,:)))<1e-300)
                Beta(k,stav1)=-Inf;
            else
                Beta(k,stav1) = log(sum(exp(gamma+Beta(k+1,:))));
            end   
        end
        Beta(k,:) = Beta(k,:) - tempmax(k+1);
    end

    %%VYPOCET A-POSTERIORI PRAVDEPODOBNOSTI
    for k = 1:dlz
        for stav2 = 1:max_stav
            a = -Lc_yk(3*k-2);
            b = Lc_yk(3*k-1:3*k);
            c = preOut(stav2,2:3)';
            gamma0 = (a+b*c)....
                -log(1+exp(L_uk(k)));
            gamma1 = (Lc_yk(3*k-2)+Lc_yk(3*k-1:3*k)*preOut(stav2,5:6)')...
                +L_uk(k)-log(1+exp(L_uk(k)));
            temp0(stav2) = exp(gamma0 + Alpha(k,preStav(stav2,1))...
                + Beta(k,stav2));
            temp1(stav2) = exp(gamma1 + Alpha(k,preStav(stav2,2))...
                + Beta(k,stav2));
        end
        L_uky(k) = log(sum(temp1)) - log(sum(temp0));                       %a-posteriori pravdepodobnost          
    end
end
