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

% model parameters
[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, 2));
VX = log(zeros(T1+1, T2+1, 2));
VY = log(zeros(T1+1, T2+1, 2));

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

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

        d1 = log(eps) + VY(i, j-1);
        d2 = log(del) + VM(i, j-1);
        [m ind] = max([d1 d2]);
        VY(i, j, 1) = log(qy(Obs2(j-1))) + m;
        VY(i, j, 2) = ind - 2;
    end
end

[~, ind] = max([VM(end, end, 1), VX(end, end, 1), VY(end, end, 1)]);

alignment = zeros(2, T1+T2);
num = 0;
t1 = T1;
t2 = T2;
num = num+1;
if ind == 1
    alignment(1, num) = Obs1(T1);
    alignment(2, num) = Obs2(T2);
    s = VM(t1, t2, 2);
    t1 = t1 - 1;
    t2 = t2 - 1;
elseif ind == 2
    alignment(1, num) = Obs1(T1);
    alignment(2, num) = -1;
    t1 = t1 - 1;
    s = VX(t1, t2, 2);
else
    alignment(1, num) = -1;
    alignment(2, num) = Obs2(T2);
    t2 = t2 - 1;
    s = VY(t1, t2, 2);
end
    
while t1 >= 1 || t2 >= 1
    num = num + 1;
    if s == 0
        alignment(1, num) = Obs1(t1);
        alignment(2, num) = Obs2(t2);
        s = VM(t1, t2, 2);
        t1 = t1 - 1;
        t2 = t2 - 1;
    elseif s == 1
        alignment(1, num) = Obs1(t1);
        alignment(2, num) = -1;
        s = VX(t1, t2, 2);
        t1 = t1 - 1;
    elseif s == -1
        alignment(1, num) = -1;
        alignment(2, num) = Obs2(t2);
        s = VY(t1, t2, 2);
        t2 = t2 - 1;
    end
end
alignment = alignment(:, 1:num);

end

