function [alignment vpath] = computePath(D, Si, align, Si1, Si2, Obs1, Obs2)
% computePath returns the viterbi path after processing the viterbi
% variables

% size of observations
% N = size(D, 1);
T1 = size(D, 2);
T2 = size(D, 3);

% allocating memory to the path variable
alignment = zeros(2, T1+T2);    % stores the alignment, '-1' stands for insertion or deletion
vpath = zeros(1, T1+T2);    % vpath(i) stores the state of the aligned pair alignment(1, i):alignment(2, i)
num = 0;
t1 = T1;
t2 = T2;

% evaluating the last state in the path
[~, ind] = max(D(:, T1, T2));
if ind == 1 % state 1, intergenic region
   s = Si(T1, T2);  % previous state
   a = align(T1, T2);
   num = num+1;
   vpath(num) = 1;
   vpath(num+1) = s;
   
   if a == 0    % match pair
       alignment(1, num) = Obs1(t1);
       alignment(2, num) = Obs2(t2);
       t1 = t1 - 1;
       t2 = t2 - 1;
   elseif a == 1    % insertion in 1st sequence
       alignment(1, num) = Obs1(t1);
       alignment(2, num) = -1;
       t1 = t1 -1;
   elseif a == -1   % deletion in 1st sequence
       alignment(1, num) = -1;
       alignment(2, num) = Obs2(t2);
       t2 = t2 -1;
   end
elseif ind == 2 % start codon
    vpath(num+1:num+3) = 2;
    alignment(1, num+1:num+3) = Obs1(t1:-1:t1-2);
    alignment(2, num+1:num+3) = Obs2(t2:-1:t2-2);
    t1 = t1-3;
    t2 = t2-3;
    num = num+3;
    vpath(num+1) = 1;
elseif ind == 3
    i = Si1(t1, t2);
    j = Si2(t1, t2);
    tAlignment = alignPair(Obs1(t1-i+1 : t1), Obs2(t2-j+1 : t2));
    d = size(tAlignment, 2);
    vpath(num+1:num+d) = 3;
    alignment(1, num+1:num+d) = tAlignment(1, 1:end);
    alignment(1, num+1:num+d) = tAlignment(2, 1:end);
    num = num+d;
    vpath(num+1) = 2;
    t1 = t1 - i;
    t2 = t2 - j;
elseif ind == 4
    vpath(num+1:num+3) = 4;
    alignment(1, num+1:num+3) = Obs1(t1 : -1 : t1-2);
    alignment(2, num+1:num+3) = Obs2(t2 : -1 : t2-2);
    t1 = t1-3;
    t2 = t2-3;
    num = num + 3;
end
ind
t1
t2
num
while t1 > 1 || t2 > 1
    if vpath(num+1) == 1
        s = Si(T1, T2);  % previous state
        a = align(T1, T2);
        num = num+1;
        vpath(num) = 1;
        vpath(num+1) = s;

        if a == 0    % match pair
            alignment(1, num) = Obs1(t1);
            alignment(2, num) = Obs2(t2);
            t1 = t1 - 1;
            t2 = t2 - 1;
        elseif a == 1    % insertion in 1st sequence
            alignment(1, num) = Obs1(t1);
            alignment(2, num) = -1;
            t1 = t1 -1;
        elseif a == -1   % deletion in 1st sequence
            alignment(1, num) = -1;
            alignment(2, num) = Obs2(t2);
            t2 = t2 -1;
        end
    elseif vpath(num+1) == 2
        vpath(num+1:num+3) = 2;
        alignment(1, num+1:num+3) = Obs1(t1:-1:t1-2);
        alignment(2, num+1:num+3) = Obs2(t2:-1:t2-2);
        t1 = t1-3;
        t2 = t2-3;
        num = num+3;
        vpath(num+1) = 1;
    elseif vpath(num+1) == 3
        i = Si1(t1, t2);
        j = Si2(t1, t2);
        tAlignment = alignPair(Obs1(t1-i+1 : t1), Obs2(t2-j+1 : t2));
        d = size(tAlignment, 2);
        vpath(num+1:num+d) = 3;
        d
        num
        size(alignment)
        size(tAlignment(2,:))
        alignment(1, num+1:num+d) = tAlignment(1, 1:end);
        alignment(2, num+1:num+d) = tAlignment(2, 1:end);
        num = num+d;
        vpath(num+1) = 2;
        t1 = t1 - i;
        t2 = t2 - j;
    elseif vpath(num+1) == 4
        vpath(num+1:num+3) = 4;
        num
        t1
        t2
        alignment(1, num+1:num+3) = Obs1(t1 : -1 : t1-2);
        alignment(2, num+1:num+3) = Obs2(t2 : -1 : t2-2);
        num = num+3;
        vpath(num+1) = 3;
        t1 = t1-3;
        t2 = t2-3;
    end

end
    
end
