1function [v, parent] = bfs (A, s, varargin)
2%GRB.BFS breadth-first search of a graph, using its adjacency matrix.
3% v = GrB.bfs (A, s) performs the breadth-first search of the directed
4% graph represented by the square adjacency matrix A.  The breadth-first
5% search starts at node s.  The output v is a sparse vector of size
6% n-by-1, with the level of each node, where v(s)=1, and v(i)=k if the
7% path with the fewest edges from from s to i has k-1 edges.  If i is not
8% reachable from s, then v(i) is implicitly zero and does not appear in
9% the pattern of v.
10%
11% [v, parent] = GrB.bfs (A, s) also computes the parent vector,
12% representing the breadth-first search tree.  parent(s)=s denotes the
13% root of the tree, and parent(c)=p if node p is the parent of c in the
14% tree.  The parent vector is sparse, and parent (i) is not present if i
15% is not found in the breadth-first search.
16%
17% Optional string arguments can be provided, after A and s:
18%
19%   'undirected' or 'symmetric':  A is assumed to be symmetric, and
20%       represents an undirected graph.  Results are undefined if A is
21%       unsymmetric, and 'check' is not specified.
22%
23%   'directed' or 'unsymmetric':  A is assumed to be unsymmetric, and
24%       presents a directed graph.  This is the default.
25%
26%   'check': with the 'undirected' or 'symmetric' option, A is checked to
27%       ensure that it is symmetric.  The default is not to check.
28%
29% For best performance, if A represents a directed graph, it should be a
30% GraphBLAS matrix stored by row on input.  That is, GrB.format (A)
31% should report 'by row'.  (If A represents a directed graph but is
32% stored 'by col' on input, it is first converted to 'by row', which is
33% costly).  If A is an undirected graph, then it can be stored in either
34% format ('by row' or 'by col').
35%
36% A must be square.  Only the pattern, spones (A), is considered; the
37% values of its entries (the edge weights of the graph) are ignored.
38%
39% Example:
40%
41%   A = bucky ;
42%   s = 1 ;
43%   [v pi] = GrB.bfs (A, s)
44%   figure (1) ;
45%   subplot (1,2,1) ; plot (graph (A)) ;
46%   pi2 = full (double (pi)) ;
47%   pi2 (s) = 0 ;
48%   subplot (1,2,2) ; treeplot (pi2) ; title ('BFS tree') ;
49%   n = size (A,1) ;
50%   for level = 1:n
51%       level
52%       inlevel = find (v == level)
53%       parents = full (double (pi (inlevel)))
54%       if (isempty (inlevel))
55%           break ;
56%       end
57%   end
58%
59% See also graph/bfsearch, graph/shortestpathtree, treeplot.
60
61% SuiteSparse:GraphBLAS, Timothy A. Davis, (c) 2017-2021, All Rights Reserved.
62% SPDX-License-Identifier: GPL-3.0-or-later
63
64% NOTE: this is a high-level algorithm that uses GrB objects.
65
66%-------------------------------------------------------------------------
67% initializations
68%-------------------------------------------------------------------------
69
70[m, n] = size (A) ;
71if (m ~= n)
72    error ('A must be square') ;
73end
74
75% get the string options
76kind = 'directed' ;
77check = false ;
78for k = 1:nargin-2
79    arg = lower (varargin {k}) ;
80    switch arg
81        case { 'undirected', 'symmetric' }
82            kind = 'undirected' ;
83        case { 'directed', 'unsymmetric' }
84            kind = 'directed' ;
85        case { 'check' }
86            check = true ;
87        otherwise
88            error ('unknown option') ;
89    end
90end
91
92% set the descriptors
93desc_rc.out  = 'replace' ;
94desc_rc.mask = 'complement' ;
95desc_s.mask = 'structural' ;
96
97% determine the method to use, and convert A if necessary
98if (isequal (kind, 'undirected'))
99    if (check && ~issymmetric (A))
100        error ('A must be symmetric') ;
101    end
102    if (GrB.isbycol (A))
103        % A is stored by column but undirected, so use q*A' instead of q*A
104        desc_rc.in1 = 'transpose' ;
105    end
106else
107    if (GrB.isbycol (A))
108        % this can be costly
109        A = GrB (A, 'by row') ;
110    end
111end
112
113% determine the integer type to use
114int_type = 'int64' ;
115if (n < intmax ('int32'))
116    int_type = 'int32' ;
117end
118
119% initialize v as a full integer vector
120v = full (GrB (1, n, int_type)) ;
121
122%-------------------------------------------------------------------------
123% do the BFS
124%-------------------------------------------------------------------------
125
126if (nargout == 1)
127
128    %---------------------------------------------------------------------
129    % just compute the level of each node
130    %---------------------------------------------------------------------
131
132    q = GrB (1, n, 'logical') ;                  % q = sparse (1,n)
133    q = GrB.subassign (q, { s }, true) ;         % q (s) = 1
134    for level = 1:n
135        % assign the current level: v<q> = level
136        v = GrB.subassign (v, q, level, desc_s) ;
137        % quit if q is empty
138        if (~any (q)), break, end
139        % move to the next level:  q<~v,replace> = q*A
140        q = GrB.mxm (q, v, 'any.pair.logical', q, A, desc_rc) ;
141    end
142
143else
144
145    %---------------------------------------------------------------------
146    % compute both the level and the parent
147    %---------------------------------------------------------------------
148
149    parent = full (GrB (1, n, int_type)) ;       % parent = zeros (1,n)
150    parent = GrB.subassign (parent, { s }, s) ;  % parent (s) = s
151    q = GrB (1, n, int_type) ;                   % q = sparse (1,n)
152    q = GrB.subassign (q, { s }, s) ;            % q (s) = s
153    id = GrB (1:n, int_type, 'by row') ;         % id = 1:n
154    semiring = ['any.1st.' int_type] ;           % any.1st.integer semiring
155    for level = 1:n
156        % assign the current level: v<q> = level
157        v = GrB.subassign (v, q, level, desc_s) ;
158        % quit if q is empty
159        if (~any (q)), break, end
160        % move to the next level:  q<~v,replace> = q*A,
161        % using the any-first-integer semiring (int32 or int64)
162        q = GrB.mxm (q, v, semiring, q, A, desc_rc) ;
163        % assign parents: parent<q> = q
164        parent = GrB.assign (parent, q, q, desc_s) ;
165        % q(i) = i for all entries in q, using q<q>=1:n
166        q = GrB.assign (q, q, id, desc_s) ;
167    end
168    % remove zeros from parent
169    parent = GrB.prune (parent) ;
170
171end
172
173% remove zeros from v
174v = GrB.prune (v) ;
175
176