%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ESERCITAZIONE 1
%
% 2 esercizio
%
% Alberton Riccardo
% Ausserer Markus
% Barazzuol Andrea
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%% Inizializzazione:

clear all
clc

% Numero delle coppie di bit trasmesse
T=20;

% Combinazioni degli stati nascosti
v1=[0 0];
v2=[0 1];
v3=[1 0];
v4=[1 1];

v=[v1' v2' v3' v4'];

% Osservazioni possibili
w=v;

% Matrice transizione di stato
A=[0.25 0.25 0.5 0; 0.5 0.5 0 0; 0.5 0 0.25 0.25; 0.25 0 0 0.75];

% Probabilità dello stato iniziale
P0=[0.25 0.25 0.25 0.25]';

% Rumore del canale
epsilon = [0.05 0.25 0.45];

% Numero di realizzazioni
num_real = 100; % da consegna si devono effettuare 100 realizzazioni

% Ciclo al variare del rumore epsilon di canale
for u=1:length(epsilon)

    e_seq_fil_viterbi_tmp =0;
    e_stima_viterbi_tmp =0;
    e_alg_alphabeta_tmp=0; 
    e_stringa_tmp=0;

    num_prob_non_0=zeros(1,4);
    coeff_z=zeros(1,4);

    % Matrice transizione stato/uscita
    % Calcolo della distanza di Hamming e della matrice C
    for l=1:4
       for k=1:4
          d=haming(v(:,l),w(:,k),1);
          C(l,k)=(1-epsilon(u))^(2-d)*epsilon(u)^d; 
       end
    end

    % Generazione di una sequenza di Markov e delle sue osservazioni
    [oss,stat] = hmmgenerate(T,A,C); %funzione matlab

    % Ciclo per effettuare num_real (=100 da consegna) realizzazioni dei 4 stimatori 
    for n = 1:num_real


        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % Calcolo del filtro alla viterbi: 
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        seq_fil_vit = [];

        for t=1:T

            seq_fil_vit = [seq_fil_vit, viterbone(t,P0,oss,A,C)'];
            [maxProb(t) q_ott(t)]=max(viterbone(t,P0,oss,A,C));

        end

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % Stima di viterbi
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        % Si calcola la stima di viterbi a partire dai dati ottenuti dal filtro
        % di viterbi. L'algoritmo parte da T e si muove all'indietro

        prob_stimaVit = maxProb;
        oss_stima = q_ott;
        
        for t=T:-1:2
            for i=1:4
                vett_stima_temp(i,t) = seq_fil_vit(i,t-1)*A(i,oss_stima(t));
            end
            [prob_stimaVit(t-1) oss_stima(t-1)] = max(vett_stima_temp(:,t));

        end

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %  Algoritmo alpha-beta
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        global betaTemp;
        betaTemp=zeros(4,T);

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % Passata all'indietro: Beta
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        for t=T:-1:1
        if t==T
            for i=1:4
                betaTemp(i,t)=1;
            end
        else
            for i=1:4
                 for j=1:4
                 betaTemp(i,t)=betaTemp(i,t)+ (A(i,j)*C(j,oss(t+1))*betaTemp(j,t+1));
                 end
            end
          end
        end

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
        % Passata in avanti: Alpha
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        global alpha
        alpha=zeros(4,T);

        for i=1:4
                   alpha(i,1)=P0(i)*C(i,oss(1)); 
        end 

        for t=1:T-1
            for j=1:4 
                for i=1:4
                  alpha(j,t+1) = alpha(j,t+1)+(A(i,j)*C(j,oss(t+1))*alpha(i,t));
                end
            end    
        end

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % Unione delle due passate: alpha_beta
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        alpha_beta=zeros(4,T);

        for t=1:T
           for i=1:4
              den=0;
              num=alpha(i,t)*betaTemp(i,t);
              for j=1:4
                 den=den + alpha(j,t)*betaTemp(j,t);
              end
                   [maxProb(t) q_ott(t)]=max(viterbone(t,P0,oss,A,C));
                   alpha_beta(i,t)=num/den;
           end
        end

        % Ricerca sequenza ottima:
        for k=1:T

            [alg_beta(k) alg_alpha(k)]=max(alpha_beta(:,k));

        end

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%% Stringa delle osservazioni
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%

        stringa=oss;


        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % Analisi delle prestazioni
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        % Riconversione degli stati in bit (necessario per la dist di hamming):
        statB = conv_bit(stat,T);
        q_ottB = conv_bit(q_ott,T);
        oss_stimaB = conv_bit(oss_stima,T);
        alg_alphaB =  conv_bit(alg_alpha,T);
        ossB = conv_bit(oss,T);


        % Valutazione Bit error rate

        % Calcolo la distanza di Hamming fra gli stati e gli stimatori

        d_ott = haming(statB,q_ottB,T);
        d_stima = haming(statB,oss_stimaB,T);
        d_alpha = haming(statB,alg_alphaB,T);
        d_oss = haming(statB,ossB,T);

        d_ott_T = 0;
        d_stima_T = 0;
        d_alpha_T = 0;
        d_oss_T = 0;

        % Sommo le distanze di Hamming fino a T

        for i =1:T
            d_ott_T = d_ott_T + d_ott(i);
            d_stima_T = d_stima_T + d_stima(i);
            d_alpha_T = d_alpha_T + d_alpha(i);
            d_oss_T = d_oss_T + d_oss(i);    
        end

        % Calcolo bit error rate
        e_seq_fil_viterbi(n)   = d_ott_T/(2*T);
        e_stima_viterbi(n)   = d_stima_T/(2*T);
        e_alg_alphabeta(n)  = d_alpha_T/(2*T);
        e_stringa(n) =  d_oss_T/(2*T);

        % Sommo la bit error rate per tutte le 100 realizzazioni
        
        e_seq_fil_viterbi_tmp   = e_seq_fil_viterbi_tmp + e_seq_fil_viterbi(n);
        e_stima_viterbi_tmp   = e_stima_viterbi_tmp + e_stima_viterbi(n);
        e_alg_alphabeta_tmp  = e_alg_alphabeta_tmp + e_alg_alphabeta(n);
        e_stringa_tmp =  e_stringa_tmp + e_stringa(n);

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %  Valuto il logaritmo della probabilità della sequenza stimata
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

         log_p=P0(1)*ones(1,4);

        for t=2:T-1
            % Probabilità delle sequenze per ognuno dei 4 stimatori
            
           log_p(1)=log_p(1)*A(q_ott(t),q_ott(t+1))*C(q_ott(t),oss(t));
           log_p(2)=log_p(2)*A(oss_stima(t),oss_stima(t+1))*C(oss_stima(t),oss(t));
           log_p(3)=log_p(3)*A(alg_alpha(t),alg_alpha(t+1))*C(alg_alpha(t),oss(t));
           log_p(4)=log_p(4)*A(oss(t),oss(t+1))*C(oss(t),oss(t));
        end

        % Calcolo dei logaritmi
        for i=1:4
          log_p(i)=log(log_p(i));
        end

        % Conto il numero di elementi di probabilità molto piccola
        for i=1:4
            if log_p(i) ~= -Inf
                num_prob_non_0(i)=num_prob_non_0(i)+1;
                coeff_z(i)=coeff_z(i) + log_p(i);
            end
        end

    end


    for i=1:4
        coeff_z(i)=coeff_z(i)/num_prob_non_0(i);
    end

    

    % media degli indici sul numero di realizzazioni:
    e_seq_fil_viterbi_media = e_seq_fil_viterbi_tmp/num_real;
    e_stima_viterbi_media = e_stima_viterbi_tmp/num_real;
    e_alg_alphabeta_media = e_alg_alphabeta_tmp/num_real;
    e_stringa_media = e_stringa_tmp/num_real;
    
    % stampo gli stati
    disp('*********************************************')
    disp(['********* Stati per epsilon = ' num2str(epsilon(u)) ' *********'])

    disp('Stati Filtro Stimatore Alfa-Beta Osservazioni')
    for iii=1:length(stat)
        str_stat = num2str(stat(iii));
        str_q_ott = num2str(q_ott(iii));
        str_oss_stima = num2str(oss_stima(iii));
        str_alg_alpha = num2str(alg_alpha(iii));
        str_oss = num2str(oss(iii));
        str_stati = ['  ' str_stat '      ' str_q_ott '        ' str_oss_stima '        ' str_alg_alpha '           ' str_oss];
        disp(str_stati)
    end
    disp('*********************************************')
    disp(' ')
    disp(' ')
    
    % salvo il vettore degli indici al variare di epsilon per la stampa    
    % Bit_e = struct([]);
    % Pr_s = struct([]);
    
    tmp1 = [e_seq_fil_viterbi_media e_stima_viterbi_media e_alg_alphabeta_media e_stringa_media];
    tmp2 =  [coeff_z(1) coeff_z(2) coeff_z(3) coeff_z(4)];
    Bit_e(u).e = tmp1; 
    Pr_s(u).s = tmp2;

end

% stampa indici:
for i=1:length(epsilon)
    disp('*********************************************')
    disp(['********* Indici per epsilon = ' num2str(epsilon(i)) ' *********'])
    disp('        Bit error rate           ')
    disp(['filtro    : ' num2str(Bit_e(i).e(1))])
    disp(['stimatore : ' num2str(Bit_e(i).e(2))])
    disp(['aplha_beta: ' num2str(Bit_e(i).e(3))])
    disp(['Stringa   : ' num2str(Bit_e(i).e(4))])
    
    disp(' ')
    disp('        Probabilità sequenza stimatore           ')
    disp(['filtro    : ' num2str(Pr_s(i).s(1))])
    disp(['stimatore : ' num2str(Pr_s(i).s(2))])
    disp(['aplha_beta: ' num2str(Pr_s(i).s(3))])
    disp(['Stringa   : ' num2str(Pr_s(i).s(4))])   
    
%     disp(' ')
%     disp('               filtro | stimatore | alpha-beta | Stringa')
%     disp(['Bit error rate: ' num2str(Bit_e(i).e)])
%     
%     disp(' ')
%     %disp()
%      disp('                       filtro |  stimatore |  alpha-beta |  Stringa')
%    disp(['Prob. sequenza stim.: ' num2str(Pr_s(i).s)])
   disp(' ')
   disp('*********************************************')
      disp(' ')

end





