function [alignment vpath D Si] = viterbiPathPair(file1, file2)

% reading the parameters of the GPHMM
I = setInitialStateProbabilitiesGPHMM;
a = setTransitionProbabilitiesGPHMM;
b = setEmissionProbabilitiesGPHMM;
[pm qx qy] = setParametersPHMM;


Obs1 = readGenes(file1);
Obs2 = readGenes(file2);
T1 = size(Obs1, 2);
T2 = size(Obs2, 2);

% number of states
N = 4;

% the viterbi forward variable
D = log(zeros(N, T1+1, T2+1));  % the durations are shifted by 1, so as to include the zeroth timestamp

% variables for backtracking, used to store information for state 1
align = zeros(T1+1, T2+1);

% previous state - only for state 1, intergenic region
Si = zeros(T1+1, T2+1);

% the duration information (of first sequence) for state 3, coding region
Si1 = zeros(T1+1, T2+1);

% the duration information (of second sequence) for state 3, coding region
Si2 = zeros(T1+1, T2+1);



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

% statistical parameters of the coding region states 3/6
% D = 7155; % maximum duration 
Dg = 7155;
% minDg = 30; % minimum duration

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

% state 1, intergenic region
D(1, 2, 2) = log(pm(Obs1(1), Obs2(1))) + log(I(1)) + log(1 - 2*del - tau);
align(2, 2) = 0;
Si(2, 2) = 0;

% state 2, direct start codon
D(2, 2:4, 2:4) = -Inf;
D(4, 2:4, 2:4) = -Inf;
D(2, 5, 5) = log(D(1, 2, 2)) + log(a(1, 2)) + log(b(2, codonIndex(Obs1(2:4)))) + log(b(2, codonIndex(Obs2(2:4))));

D(3, 2, 2) = -Inf;
% state 4, direct stop codon
D(4, 2:4, 2:4) = -Inf;

%% recursive computation of the rest of the values..
for t1 = 2:T1+1
    disp('******');
    disp(t1);
    disp('******');
    for t2 = 2:T2+1
        disp(t2);

        % intergenic coding region
        d11 = D(1, t1-1, t2-1) + log(1-2*del-tau) + log(pm(Obs1(t1-1), Obs2(t2-1))) + log(a(1,1));
        d12 = D(1, t1-1, t2) + log(del) + log(qx(Obs1(t1-1))) + log(a(1,1));
        d13 = D(1, t1, t2-1) + log(del) + log(qy(Obs2(t2-1))) + log(a(1,1));
        
        d41 = D(4, t1-1, t2-1) + log(1-2*del-tau) + log(pm(Obs1(t1-1), Obs2(t2-1))) + log(a(4,1));
        d42 = D(4, t1-1, t2) + log(del) + log(qx(Obs1(t1-1))) + log(a(4,1));
        d43 = D(4, t1, t2-1) + log(del) + log(qy(Obs2(t2-1))) + log(a(4,1));
        ar = [d11 d12 d13 d41 d42 d43];
        [m ind] = max(ar);
        D(1, t1, t2) = m;
        if ind == 1
            Si(t1, t2) = 1;
            align(t1, t2) = 0;
        elseif ind == 2
            Si(t1, t2) = 1;
            align(t1, t2) = 1;
        elseif ind == 3
            Si(t1, t2) = 1;
            align(t1, t2) = -1;
        elseif ind == 4
            Si(t1, t2) = 4;
            align(t1, t2) = 0;
        elseif ind == 5
            Si(t1, t2) = 4;
            align(t1, t2) = 1;
        elseif ind == 6
            Si(t1, t2) = 4;
            align(t1, t2) = -1;
        end
        
        
        % coding region state (state 3)
        ar1 = [t1-1, Dg];
        ar2 = [t2-1, Dg];
        [m1 i1] = min(ar1);
        [m2 i2] = min(ar2);
        maxv = -Inf;
        maxi = i1;
        maxj = i2;
        for i = 1:m1
            for j = 1:m2
                score = alignScore(Obs1(t1-i:t1-1), Obs2(t2-j:t2-1));
                val = log(pr(i)) + log(pr(j)) + log(a(2, 3)) + score + D(2, t1-i, t2-j);
                if val > maxv
                    maxv = val;
                    maxi = i;
                    maxj = j;
                end
            end
        end
        D(3, t1, t2) = maxv;
        Si1(t1, t2) = maxi;
        Si2(t1, t2) = maxj;

        % codon states
        if t1 > 4 && t2 > 4
            D(2, t1, t2) = D(1, t1-3, t2 - 3) + log(a(1, 2)) + log(b(2, codonIndex(Obs1(t1-3:t1-1)))) + log(b(2, codonIndex(Obs2(t2-3:t2-1))));
            D(4, t1, t2) = D(3, t1-3, t2 - 3) + log(a(3, 4)) + log(b(4, codonIndex(Obs1(t1-3:t1-1)))) + log(b(4, codonIndex(Obs2(t2-3:t2-1))));
        end
    end
end

%% backtracking and computing the path
D1 = D(:, 2:T1+1, 2:T2+1);
S = Si(2:T1+1, 2:T2+1);
align1 = align(2:T1+1, 2:T2+1);
S1 = Si1(2:T1+1, 2:T2+1);
S2 = Si2(2:T1+1, 2:T2+1);
[alignment vpath] = computePath(D1, S, align1, S1, S2, Obs1, Obs2);

%{
a1 = find(alignment(1,:) ~= 0);
b1 = alignment(1,a1);
b1 = b1(end:-1:1);

a2 = find(alignment(2,:) ~= 0);
b2 = alignment(2, a2);
b2 = b2(end:-1:1);

s = size(b1, 2);
alignment = zeros(2, s);
for i = 1:s
    if b1(i) == 1
        alignment(1, i) = 'A';
    elseif b1(i) == 2
        alignment(1, i) = 'C';
    elseif b1(i) == 3
        alignment(1, i) = 'G';
    elseif b1(i) == 4
        alignment(1, i) = 'T';
    elseif b1(i) == -1
        alignment(1, i) = '-';
    end
    
    if b2(i) == 1
        alignment(2, i) = 'A';
    elseif b2(i) == 2
        alignment(2, i) = 'C';
    elseif b2(i) == 3
        alignment(2, i) = 'G';
    elseif b2(i) == 4
        alignment(2, i) = 'T';
    elseif b2(i) == -1
        alignment(2, i) = '-';
    end


end
%}
end
