function [viterbi_path, viterbi_score] = viterbi(pi, a, b, num_states, obs)

%*****************************************************************
%Preprocessing
pi_tilda = log(pi);
b_tilda = log(b);
a_tilda = log(a);

T=size(obs,1);
delta_tilda = zeros(T,num_states);
psi = zeros(T,num_states);
%*****************************************************************
%Initialization
delta_tilda(1,:) = pi_tilda' + b_tilda(:,obs(1)+1)';
psi(1,:) = zeros(1,num_states);

%*****************************************************************
%Recursion
for t=2:T
    for j=1:num_states
        [max_val, max_index] = max(delta_tilda(t-1,:) + a_tilda(:,j)');
        delta_tilda(t,j) = max_val + b_tilda(j, obs(t)+1);
        psi(t,j) = max_index;
    end
end

%*****************************************************************
%Termination and Backtracking
[max_val, max_index] = max(delta_tilda(T,:));
viterbi_score = max_val;

viterbi_path=[max_index];

for t=1:T-1
    viterbi_path = [psi(T-t+1, viterbi_path(1)) viterbi_path];
end
