MAXCARDSEARCH Determine graph chordality and maximal cliques. [MAXCLIQUE,ISCHORDAL] = MAXCARDSEARCH(AADJ) Determine the maximal cliques for a chordal graph described by the adjacency matrix Aadj. Also test the graph for chordality. The algorithms for determining the maximal cliques and testing for chordality are described in [1]. Inputs: AADJ : The adjacency matrix of a graph. If the graph is chordal, the maximal cliques for this graph are calculated. Outputs: MAXCLIQUE : If the graph described by Aadj is chordal, maxclique returns a matrix describing the maximal cliques of the graph. Each column of maxclique represents a clique with nonzero entries indicating the nodes included in the maximal clique. If Aadj is not chordal, maxclique is set to NaN. ISCHORDAL : Returns 1 if the graph described by Aadj is chordal, otherwise returns 0. [1] Tarjan, R. E., and M. Yannakakis. "Simple Linear-Time Algorithms to Test Chordality of Graphs, Test Acyclicity of Hypergraphs, and Selectively Reduce Acyclic Hypergraphs." SIAM Journal on computing 13 (1984): 566.
0001 function [maxclique,ischordal] = maxCardSearch(Aadj) 0002 %MAXCARDSEARCH Determine graph chordality and maximal cliques. 0003 % [MAXCLIQUE,ISCHORDAL] = MAXCARDSEARCH(AADJ) 0004 % 0005 % Determine the maximal cliques for a chordal graph described by the 0006 % adjacency matrix Aadj. Also test the graph for chordality. The 0007 % algorithms for determining the maximal cliques and testing for 0008 % chordality are described in [1]. 0009 % 0010 % Inputs: 0011 % AADJ : The adjacency matrix of a graph. If the graph is chordal, 0012 % the maximal cliques for this graph are calculated. 0013 % 0014 % Outputs: 0015 % MAXCLIQUE : If the graph described by Aadj is chordal, maxclique 0016 % returns a matrix describing the maximal cliques of the graph. 0017 % Each column of maxclique represents a clique with nonzero 0018 % entries indicating the nodes included in the maximal clique. If 0019 % Aadj is not chordal, maxclique is set to NaN. 0020 % ISCHORDAL : Returns 1 if the graph described by Aadj is chordal, 0021 % otherwise returns 0. 0022 % 0023 % [1] Tarjan, R. E., and M. Yannakakis. "Simple Linear-Time Algorithms to 0024 % Test Chordality of Graphs, Test Acyclicity of Hypergraphs, and 0025 % Selectively Reduce Acyclic Hypergraphs." SIAM Journal on computing 13 0026 % (1984): 566. 0027 0028 % MATPOWER 0029 % $Id: maxCardSearch.m 2272 2014-01-17 14:15:47Z ray $ 0030 % by Daniel Molzahn, PSERC U of Wisc, Madison 0031 % Copyright (c) 2013-2014 by Power System Engineering Research Center (PSERC) 0032 % 0033 % This file is part of MATPOWER. 0034 % See http://www.pserc.cornell.edu/matpower/ for more info. 0035 % 0036 % MATPOWER is free software: you can redistribute it and/or modify 0037 % it under the terms of the GNU General Public License as published 0038 % by the Free Software Foundation, either version 3 of the License, 0039 % or (at your option) any later version. 0040 % 0041 % MATPOWER is distributed in the hope that it will be useful, 0042 % but WITHOUT ANY WARRANTY; without even the implied warranty of 0043 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 0044 % GNU General Public License for more details. 0045 % 0046 % You should have received a copy of the GNU General Public License 0047 % along with MATPOWER. If not, see <http://www.gnu.org/licenses/>. 0048 % 0049 % Additional permission under GNU GPL version 3 section 7 0050 % 0051 % If you modify MATPOWER, or any covered work, to interface with 0052 % other modules (such as MATLAB code and MEX-files) available in a 0053 % MATLAB(R) or comparable environment containing parts covered 0054 % under other licensing terms, the licensors of MATPOWER grant 0055 % you additional permission to convey the resulting work. 0056 0057 %% Setup 0058 0059 nbus = size(Aadj,1); 0060 nline = nnz(Aadj)/2; 0061 0062 %% Create a perfect elimination ordering 0063 0064 % Store in s{i} all unnumbered vertices adjacent to exactly i numbered 0065 % vertices. 0066 s{nline} = []; 0067 s{1} = 1:nbus; % This corresponds to s{0} in the paper's notation 0068 0069 % Maintain the largest index j such that s{j} is nonempty 0070 jidx = 1; 0071 0072 % Candidate perfect elimination ordering, with vertex numbering in vnum. 0073 % Index of vnum corresponds to the original vertex number, and the value of 0074 % vnum indicates the new vertex number. 0075 alpha = zeros(nbus,1); 0076 0077 % v2s stores which set contains a given vertex? 0078 v2s = ones(nbus,1); 0079 0080 i = nbus; 0081 0082 while i >= 1 0083 0084 % To carry out a step of the search, we remove a vertex v from s(jidx) 0085 % and number it. 0086 v = s{jidx}(1); 0087 if length(s{jidx}(:)) > 1 0088 s{jidx} = s{jidx}(2:end); 0089 else 0090 s{jidx} = []; 0091 end 0092 alpha(v) = i; 0093 v2s(v) = 0; % This vertex is no longer in a set 0094 0095 % For each unnumbered vertex w adjacent to v, we move w from the set 0096 % containing it to the next set 0097 vadj = find(Aadj(:,v)); 0098 for k=1:length(vadj) 0099 % If this node isn't already numbered, remove it from the original set 0100 if v2s(vadj(k)) ~= 0 0101 s{v2s(vadj(k))} = s{v2s(vadj(k))}(s{v2s(vadj(k))}(:) ~= vadj(k)); 0102 0103 % Add it to the new set 0104 s{v2s(vadj(k))+1} = [s{v2s(vadj(k))+1}(:); vadj(k)]; 0105 0106 % Update v2s 0107 v2s(vadj(k)) = v2s(vadj(k)) + 1; 0108 end 0109 end 0110 0111 % Add one to jidx and then decrease jidx until reaching a non-empty s(jidx) 0112 jidx = jidx + 1; 0113 if jidx > length(s) 0114 jidx = jidx - 1; 0115 end 0116 while jidx >= 1 && isempty(s{jidx}) 0117 jidx = jidx - 1; 0118 end 0119 0120 i = i-1; 0121 0122 end 0123 0124 0125 %% Check for chordality 0126 0127 f = zeros(nbus,1); 0128 index = zeros(nbus,1); 0129 ischordal = 1; 0130 for i=1:nbus 0131 w = find(alpha == i); 0132 f(w) = w; 0133 index(w) = i; 0134 0135 valid_v = find(Aadj(:,w)); 0136 valid_v = valid_v(alpha(valid_v) < i); 0137 for vidx = 1:length(valid_v) 0138 v = valid_v(vidx); 0139 index(v) = i; 0140 if f(v) == v 0141 f(v) = w; 0142 end 0143 end 0144 0145 for vidx = 1:length(valid_v) 0146 v = valid_v(vidx); 0147 if index(f(v)) < i 0148 ischordal = 0; 0149 break; 0150 end 0151 end 0152 0153 if ~ischordal 0154 break; 0155 end 0156 end 0157 0158 0159 %% Determine maximal cliques 0160 0161 % According to https://en.wikipedia.org/wiki/Chordal_graph 0162 % "To list all maximal cliques of a chordal graph, simply find a perfect 0163 % elimination ordering, form a clique for each vertex v together with the 0164 % neighbors of v that are later than v in the perfect elimination ordering, 0165 % and test whether each of the resulting cliques is maximal." 0166 0167 % alpha is a candidate perfect elimination ordering. First form all 0168 % cliques. Then determine which cliques are maximal. Put maximal cliques in 0169 % the columns of the variable maxclique. 0170 0171 if ischordal 0172 % Form a matrix representation of the cliques 0173 clique = speye(nbus,nbus); 0174 for i=1:nbus 0175 % neighbors of node i that are later in the ordering 0176 neighbors = find(Aadj(i,:)); 0177 neighbors = neighbors(alpha(neighbors) > alpha(i)); 0178 0179 clique(i,neighbors) = 1; 0180 end 0181 0182 % Test whether each clique is maximal. A clique is not maximal if a node 0183 % neighboring any of the nodes in the clique, but not in the clique, is 0184 % connected to all nodes in the clique. 0185 i=0; 0186 nclique = size(clique,1); 0187 while i < nclique 0188 i = i+1; 0189 0190 neighbors = []; 0191 cliquei = find(clique(i,:)); 0192 cliquei_bool = clique(i,:).'; 0193 nnzi = sum(cliquei_bool); 0194 0195 for k = 1:length(cliquei) % Check all nodes adjacent to nodes in the clique, but not included in the clique 0196 neighbors = [neighbors find(Aadj(cliquei(k),:))]; 0197 end 0198 neighbors = unique(setdiff(neighbors,cliquei)); 0199 0200 not_maximal = 0; 0201 for m=1:length(neighbors) 0202 % If this neighbor is connected to all other nodes in the clique, 0203 % this is not a maximal clique 0204 if Aadj(neighbors(m),:) * cliquei_bool >= nnzi 0205 not_maximal = 1; 0206 break; 0207 end 0208 end 0209 0210 if not_maximal 0211 clique = clique([1:i-1 i+1:nclique],:); % Delete this clique 0212 i = i-1; % After deletion, the next clique is really numbered i 0213 nclique = nclique - 1; 0214 end 0215 end 0216 maxclique = clique.'; % put maximal cliques in the columns 0217 else 0218 maxclique = nan; % maximal clique algorithm only works for chordal graphs 0219 end