1function iset = mis (A, check)
2%GRB.MIS variant of Luby's maximal independent set algorithm.
3%
4%   iset = GrB.mis (A) ;
5%
6% Given an n-by-n symmetric adjacency matrix A of an undirected graph,
7% GrB.mis (A) finds a maximal set of independent nodes and returns it as a
8% logical vector, iset, where iset(i) of true implies node i is a member of
9% the set.
10%
11% The matrix A must not have any diagonal entries (self edges), and it must
12% be symmetric.  These conditions are not checked by default, and results
13% are undefined if they do not hold.  In particular, diagonal entries will
14% cause the method to stall.  To check these conditions, use:
15%
16%   iset = GrB.mis (A, 'check') ;
17%
18% Reference: M Luby. 1985. A simple parallel algorithm for the maximal
19% independent set problem. In Proceedings of the seventeenth annual ACM
20% symposium on Theory of computing (STOC '85). ACM, New York, NY, USA,
21% 1-10.  DOI: https://doi.org/10.1145/22145.22146
22%
23% See also GrB.offdiag.
24
25% SuiteSparse:GraphBLAS, Timothy A. Davis, (c) 2017-2021, All Rights Reserved.
26% SPDX-License-Identifier: GPL-3.0-or-later
27
28% NOTE: this is a high-level algorithm that uses GrB objects.
29
30[m, n] = size (A) ;
31if (m ~= n)
32    error ('A must be square') ;
33end
34
35% convert A to logical
36A = GrB.apply ('1.logical', A) ;
37
38if (nargin < 2)
39    check = false ;
40else
41    if (isequal (check, 'check'))
42        check = true ;
43    else
44        error ('unknown option') ;
45    end
46end
47
48if (check)
49    if (nnz (diag (A)) > 0)
50        error ('A must not have any diagonal entries') ;
51    end
52    if (~issymmetric (A))
53        error ('A must be symmetric') ;
54    end
55end
56
57neighbor_max = GrB (n, 1) ;
58new_neighbors = GrB (n, 1, 'logical') ;
59candidates = GrB (n, 1, 'logical') ;
60
61% Initialize independent set vector
62iset = GrB (n, 1, 'logical') ;
63
64% descriptor: C_replace
65r_desc.out = 'replace' ;
66
67% descriptor: C_replace + structural complement of mask
68sr_desc.mask = 'complement' ;
69sr_desc.out  = 'replace' ;
70
71% compute the degree of each nodes
72degrees = GrB.vreduce ('+.double',  A) ;
73
74% Singletons require special treatment.  Since they have no neighbors, their
75% prob is never greater than the max of their neighbors, so they never get
76% selected and cause the method to stall.  To avoid this case they are removed
77% from the candidate set at the begining, and added to the iset.
78
79% candidates (degree != 0) = true
80candidates = GrB.assign (candidates, degrees, true) ;
81
82% add all singletons to iset
83% iset (degree == 0) = 1
84iset = GrB.assign (iset, degrees, true, sr_desc) ;
85
86% Iterate while there are candidates to check.
87ncand = GrB.entries (candidates) ;
88last_ncand = ncand ;
89
90while (ncand > 0)
91
92    % compute a random probability scaled by inverse of degree
93    % FUTURE: this is slower than it should be; rand may not be parallel,
94    % See GraphBLAS/Demo/Source/mis.c and the prand_* functions for a better
95    % approach using user-defined types and operators.
96    prob = 0.0001 + rand (n,1) ./ (1 + 2 * degrees) ;
97    prob = GrB.assign (prob, candidates, prob, r_desc) ;
98
99    % compute the max probability of all neighbors
100    neighbor_max = GrB.mxm (neighbor_max, candidates, ...
101        'max.second.double', A, prob, r_desc) ;
102
103    % select node if its probability is > than all its active neighbors
104    new_members = GrB.eadd (prob, '>', neighbor_max) ;
105
106    % add new members to independent set.
107    iset = GrB.eadd (iset, '|', new_members) ;
108
109    % remove new members from set of candidates
110    candidates = GrB.apply (candidates, new_members, 'identity', ...
111        candidates, sr_desc) ;
112
113    ncand = GrB.entries (candidates) ;
114    if (ncand == 0)
115        break ;                    % early exit condition
116    end
117
118    % Neighbors of new members can also be removed from candidates
119    new_neighbors = GrB.mxm (new_neighbors, candidates, ...
120        '|.&.logical', A, new_members) ;
121
122    candidates = GrB.apply (candidates, new_neighbors, 'identity', ...
123        candidates, sr_desc) ;
124
125    ncand = GrB.entries (candidates) ;
126
127    % this will not occur, unless the input is corrupted somehow
128    if (last_ncand == ncand)
129        error ('method stalled; rerun with ''check'' option') ;
130    end
131    last_ncand = ncand ;
132end
133
134% drop explicit false values
135iset = GrB.prune (iset) ;
136
137