function [mlseqs, cost] = viterbi_decode(trels_edges, n, p, inseq, m, find_metric)
% Return the maximum-likelihood sequence of complex numbers represented by
% .. the input sequence of encoded symbol values, given a trellis graph
% .. and a metric function.
% n is the number of encoder states.
% p is the number of distinct code symbols.
% The trellis is a directed bipartite graph represented by an "edge array".
% (For the trellis, the vertices are implied.)
% An "edge array" is a 3-dim [(n*p) x q x 3] array, where q is #columns.
% Each edge in the array contains 3 numbers. These are the 3 values of 3rd dim:
%    destination row, an integer
%    input, a complex number
%    output, a complex number
% Each edge has a source vertex (row, col) and a destination vertex (row, col).
% For each edge in an "edge array", the (src_col, dst_col) 
% .. are given by the col# of the "edge array".
% .. (src_col = edge_col) and (dst_col = edge_col + 1).
% For each edge in an "edge array", (src_row, dst_row)
% .. are determined by the row# of the "edge array".
% .. (See functions edge_to_vert() and vert_to_edge().)
% Each column of an "edge array" represents a state transition.
% trels_edges is an "edge array" for which q==1 (a single column).
% trels_edges represents a directed bipartite graph model of a symbol encoder.
% .. Each of the two graph partitions contain n vertices.
% .. One partition represents the "before" state of a given transition.
% .. The other partition represents the "after" state of the transition.
% .. Each vertex of the trellis graph has exactly p contingent edges.
% inseq is a sequence of complex numbers, each representing an encoded symbol value.
% m is the number of symbols contained in inseq
% find_metric names a Matlab mfile function that returns a single complex number.
%
% D. Bozarth 03-15-08 SSU Engineering Science
%
mlseqs  = zeros(1, m);  % Set of k minimum-cost symbol sequences, [k x m]
cost    = 0;

% build_verts gives the vertices of a new graph: n rows by (m + 1) columns.
% build_verts is a 3-dim array: [n x (m+1) x 2].
% .. The first value of the 3rd dim is a complex number that
% .. .. represents the metric value assigned to the vertex.
% .. The second value of the 3rd dim is, if defined, an integer that
% .. .. represents the row of the edge from the previous column
% .. .. whose destination is this vertex. [Note: If defined, this value
% .. .. points back to an edge whose first value of its 3rd dimension indicates
% .. .. this vertex - i.e., reciprocal pointers relate the edge and vertex.
% .. .. Note also that if this val defined, the pointed edge will have defined
% .. .. input (2nd value of 3rd dim) and defined output (3rd value of 3rd dim)]
% build_edges is an "edge array" for which q == m, the size of inseq.
% build_edges together with build_verts, define the new graph "Build".
[build_edges, build_verts] = initialize(trels_edges, n, p, m);

% Find all surviving paths in the Build graph.
% Invariant: Each vertex has at most one edge whose destination is that vertex.
for ndx = 1:m           % For each of m symbols, ndx is the graph column.
    s = inseq(ndx);     % the symbol value
    
    % Visit each vertex that serves as a destination for the previous column.
    if (ndx == 1)
        lim = 1;
    else
        lim = n;
    end
    
    for i = 1:lim         % each row of build_verts
        % If primary or leaf vertex, then process.
        if (ndx == 1) | (0 < build_verts(i, ndx, 2))
            [build_edges, build_verts] = ...
                do_vertex(trels_edges, build_edges, build_verts, ...
                          n, p, i, ndx, s, find_metric);
        end
    end
end

%build_edges
%build_verts

% Identify a minimum-metric path, and specify the return sequence.
[mlseqs, cost] = minimize(build_edges, build_verts, n, p, m);
%
% end function main
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [be, bv] = initialize(trels, n, p, m)
% Initialize the Build graph
% trels is the edge array representing the Trellis graph.
% be, bv are the Edge and Vertex array components of the Build graph.
%
be = zeros(n*p, m,  3);
bv = zeros(n,  m+1, 2);

% initialize the metric of each vertex in build_verts to "infinity",
% .. and its back-pointer to 0
for i = 1:n             % each row of build_verts
    for j = 1 :(m + 1)  % each col of build_verts
        bv(i, j, :) = [Inf 0]';

        % initialize each value of build_edges to "edge doesn't exist yet".
        if j <= m
            for k = 1:p
                erow            = (i - 1) * p + k;
                rd              = trels(erow, 1, 1);
                be(erow, j, :) 	= [rd NaN NaN]';
            end
        end        
    end
end

bv(1, 1, 1) = 0;
%
% end function
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [be, bv] = do_vertex(trels, be, bv, n, p, r, c, s, fun)
% Process one vertex of the Build graph.
% trels is the edge array representing the Trellis graph.
% be, bv are the Edge and Vertex array components of the Build graph.
% r, c identify the vertex by row and column.
% s is the symbol value corresponding to the transition c.
% fun names a Matlab mfile function that returns a single complex number.
%
for j = 1:p     % each trellis edge whose source vertex has row r
    % Test a new edge in the Build graph.
    erow        = (r - 1) * p + j;      % the edge row
    rd          = trels(erow, 1, 1);    % the destination vertex row
    output      = trels(erow, 1, 3);
    src_metric  = bv(r, c, 1);
    metric      = fun(s, output, src_metric);

    if metric < bv(rd, c + 1, 1)
        % The new edge wins.
        input   = trels(erow, 1, 2);

        % Be sure any edge whose dest is this vertex, gets removed.
        er_check = bv(rd, c + 1, 2);
        
        if 0 < er_check
            be(er_check, c, :) = [rd NaN NaN]';                       
        end

        % Establish the new edge, its destination & vertex metric.
        be(erow, c,   :)    = [rd input' output']';
        bv(rd, c + 1, :)    = [metric' erow]';
    end
end
% end function
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [mlseqs, cost] = minimize(be, bv, n, p, m)
% Return an array [k x m] of complex numbers representing
% .. the set of maximum-likelihood input sequences - k is number of solutions. 
% Also return the minimized cost metric of that sequence.
% be, bv are the populated edge array and vertex array of the Build graph.
% n is the number of encoder states.
% p is the number of distinct symbol values.
% m is the number symbols of inseq.
%
% D.Bozarth 03-15-08 SSU Engineering Science
%
cost        = Inf;
min_row     = 0;
min_list    = 0;
min_ndx     = 0;
mlseqs      = zeros(1, m);
vc_max      = m + 1;        % highest-numbered col of bv

for i = 1:m
    mlseqs(1, i) = NaN;
end

% Find the terminal vertex of a minimum-cost path.
for j = 1:n
    this_metric = bv(j, vc_max, 1);
    
    if this_metric < cost
        cost = this_metric;
    end
end

% Enumerate the minimum-cost paths.
for j = 1:n
    if bv(j, vc_max, 1) == cost
        min_ndx                 = min_ndx + 1;
        min_list(1, min_ndx)    = j;
    end
end

% Build the return values.
for j = 1:min_ndx
    min_row = min_list(1, j);
    
    for i = m:(-1):2
        erow            = bv(min_row, i + 1, 2);
        mlseqs(j, i)    = be(erow, i, 2);
        min_row         = floor((erow - 1) / p) + 1;
    end

    erow            = bv(min_row, 2, 2);
    mlseqs(j, 1)    = be(erow, 1, 2);
end
% end function