function [disc,pred,close] = graphtraverse(G,S,NVPArgs)
%GRAPHTRAVERSE performs traversal of the graph by following adjacent nodes.
% 
% [DISC,PRED,CLOSE] = GRAPHTRAVERSE(G,S) traverses graph G starting from
% the node S. G is a directed graph represented by a square matrix
% in which all nonzero entries indicate the presence of an edge. DISC is a
% list of nodes in the order in which they are discovered. PRED is a vector
% of predecessor node indices of the resulting tree, and, CLOSE is
% a list of node indices in the order in which they are closed. 
% 
% GRAPHTRAVERSE(...,'METHOD',METHOD) selects the algorithm to use, options
% are:
%    'BFS'         - Breadth First Search.
%   ['DFS']        - Depth First Search. 
% 
% The time complexity of either of these algorithms is O(n+e), where n and
% e are number of nodes and edges respectively.
% 
% GRAPHTRAVERSE(...,'DIRECTED',false) indicates that the graph G is
% undirected, upper triangle of the sparse matrix is ignored. Default is
% true.
% 
% GRAPHTRAVERSE(...,'DEPTH',D) limits the depth of the search. Default is
% Inf.
%  
% Examples:
%
%   % Create a directed graph with 10 nodes and 12 edges
%   DG = sparse([1 2 3 4 5 5 5 6 7 8 8 9],[2 4 1 5 3 6 7 9 8 1 10 2],true,10,10)
%
%   % Find the DFS discovery order starting at node 4
%   order = graphtraverse(DG,4)
%   % Label the nodes with the DFS discovery order
%   h = view(biograph(DG))
%   for i = 1:10
%       h.Nodes(order(i)).Label = sprintf('%s:%d',h.Nodes(order(i)).ID,i);
%   end
%   h.ShowTextInNodes = 'label'
%   dolayout(h)
%
%   % Find the BFS discovery order starting at node 4
%   order = graphtraverse(DG,4,'Method','BFS')
%   % Label the nodes with the BFS discovery order
%   h = view(biograph(DG))
%   for i = 1:10
%       h.Nodes(order(i)).Label = sprintf('%s:%d',h.Nodes(order(i)).ID,i);
%   end
%   h.ShowTextInNodes = 'label'
%   dolayout(h)
%
%   % Find and mark nodes that are close (2 edges) from node 4
%   h = view(biograph(DG))
%   node_idxs = graphtraverse(DG, 4, 'Depth',2, 'Method', 'BFS')
%   set(h.nodes(node_idxs),'Color',[1 0 0])
%
%
% See also: GRAPHALLSHORTESTPATHS, GRAPHCONNCOMP, GRAPHISDAG,
% GRAPHISOMORPHISM, GRAPHISSPANTREE, GRAPHMAXFLOW, GRAPHMINSPANTREE,
% GRAPHPRED2PATH, GRAPHSHORTESTPATH, GRAPHTOPOORDER.
%
% Reference: 
%  [1]	R. Sedgewick "Algorithms in C++, Part 5 Graph Algorithms"
%       Addison-Wesley, 2002.  

%   Copyright 2006-2021 The MathWorks, Inc.
arguments
    G {bioinfo.internal.ValidationHelper.mustBeAdjacencyGraph}
    S(1,1) double {mustBePositive, mustBeInteger}
    NVPArgs.Directed(1,1) logical = true;
    NVPArgs.Depth(1,1) double {mustBeNonnegative}
    NVPArgs.Method(1,1) string {bioinfo.internal.ValidationHelper.mustBeMemberCaseInsensitive(NVPArgs.Method, ["BFS", "DFS"])} = "DFS";
end

G_graph = bioinfo.internal.biograph2matlab(G, 'Directed', NVPArgs.Directed);

if strcmpi(NVPArgs.Method, "DFS") && isfield(NVPArgs, 'Depth')
    error(message('bioinfo:graphtraverse:DepthDFS'));
end

dist = distances(G_graph, S, 'Method', 'unweighted');

if strcmpi(NVPArgs.Method, "BFS")
    searchEvt = bfsearch(G_graph, S, {'discovernode', 'finishnode', 'edgetonew'});
else
    searchEvt = dfsearch(G_graph, S, {'discovernode', 'finishnode', 'edgetonew'});
end

disc = searchEvt.Node(searchEvt.Event == 'discovernode')';
close = searchEvt.Node(searchEvt.Event == 'finishnode')';
dist = dist(disc); % Reorder to match discovery

% Predecessor nodes are NOT listed in order of discovery.
predEvents = searchEvt.Event == 'edgetonew';
pred = NaN(1, numnodes(G_graph));
pred(S) = 0; % Start node has no predecessors
pred(searchEvt.Edge(predEvents, 2)) = searchEvt.Edge(predEvents, 1);


% Now filter by depth
if isfield(NVPArgs, 'Depth')
    toFilter = disc(dist > NVPArgs.Depth);
    pred(toFilter) = NaN;
    close = setdiff(close, toFilter, 'stable');
    disc = disc(dist <= NVPArgs.Depth);
end
