function [score] = alignScore(Obs1, Obs2)
% this function computes the alignment score of the two *coding* sequences
% using the viterbi algorithm for a Pair HMM

[pm qx qy] = setParametersPHMM;

% parameters for the Pair HMM
del = 0.05;
tau = 0.01;
eps = 0.05;

T1 = size(Obs1, 2);
T2 = size(Obs2, 2);

VM = log(zeros(T1+1, T2+1));
VX = log(zeros(T1+1, T2+1));
VY = log(zeros(T1+1, T2+1));

% initialization
VM(1, 1) = 0;
VM(2:T1+1, 1) = -Inf;
VM(1, 2:T2+1) = -Inf;

% recursive computation of the rest of the viterbi values
for i = 2:T1+1
    for j = 2:T2+1
        d1 = log(1-2*del-tau) + VM(i-1, j-1);
        d2 = log(1-eps-tau) + VX(i-1, j-1);
        d3 = log(1-eps-tau) + VY(i-1, j-1);
        VM(i, j) = log(pm(Obs1(i-1), Obs2(j-1))) + max([d1 d2 d3]);
        
        d1 = log(del) + VM(i-1, j);
        d2 = log(eps) + VX(i-1, j);
        VX(i, j) = log(qx(Obs1(i-1))) + max([d1 d2]);
        
        d1 = log(del) + VM(i, j-1);
        d2 = log(eps) + VY(i, j-1);
        VY(i, j) = log(qy(Obs2(j-1))) + max([d1 d2]);
    end
end

score = log(tau) + max([VM(end, end), VX(end, end), VY(end, end)]);

end
