Description of maxCardSearch

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
maxCardSearch

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
