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 %   Copyright (c) 2013-2016, Power Systems Engineering Research Center (PSERC)
0030 %   by Daniel Molzahn, PSERC U of Wisc, Madison
0031 %
0032 %   This file is part of MATPOWER.
0033 %   Covered by the 3-clause BSD License (see LICENSE file for details).
0034 %   See http://www.pserc.cornell.edu/matpower/ for more info.
0035 
0036 %% Setup
0037 
0038 nbus = size(Aadj,1);
0039 nline = nnz(Aadj)/2;
0040 
0041 %% Create a perfect elimination ordering
0042 
0043 % Store in s{i} all unnumbered vertices adjacent to exactly i numbered
0044 % vertices.
0045 s{nline} = [];
0046 s{1} = 1:nbus; % This corresponds to s{0} in the paper's notation
0047 
0048 % Maintain the largest index j such that s{j} is nonempty
0049 jidx = 1;
0050 
0051 % Candidate perfect elimination ordering, with vertex numbering in vnum.
0052 % Index of vnum corresponds to the original vertex number, and the value of
0053 % vnum indicates the new vertex number.
0054 alpha = zeros(nbus,1);
0055 
0056 % v2s stores which set contains a given vertex?
0057 v2s = ones(nbus,1);
0058 
0059 i = nbus;
0060 
0061 while i >= 1
0062     
0063     % To carry out a step of the search, we remove a vertex v from s(jidx)
0064     % and number it.
0065     v = s{jidx}(1);
0066     if length(s{jidx}(:)) > 1
0067         s{jidx} = s{jidx}(2:end);
0068     else
0069         s{jidx} = [];
0070     end
0071     alpha(v) = i;
0072     v2s(v) = 0; % This vertex is no longer in a set
0073     
0074     % For each unnumbered vertex w adjacent to v, we move w from the set
0075     % containing it to the next set
0076     vadj = find(Aadj(:,v));
0077     for k=1:length(vadj)
0078         % If this node isn't already numbered, remove it from the original set
0079         if v2s(vadj(k)) ~= 0
0080             s{v2s(vadj(k))} = s{v2s(vadj(k))}(s{v2s(vadj(k))}(:) ~= vadj(k));
0081         
0082             % Add it to the new set
0083             s{v2s(vadj(k))+1} = [s{v2s(vadj(k))+1}(:); vadj(k)];
0084 
0085             % Update v2s
0086             v2s(vadj(k)) = v2s(vadj(k)) + 1;
0087         end
0088     end
0089     
0090     % Add one to jidx and then decrease jidx until reaching a non-empty s(jidx)
0091     jidx = jidx + 1; 
0092     if jidx > length(s)
0093         jidx = jidx - 1;
0094     end
0095     while jidx >= 1 && isempty(s{jidx})
0096         jidx = jidx - 1;
0097     end
0098     
0099     i = i-1;
0100     
0101 end
0102 
0103 
0104 %% Check for chordality
0105 
0106 f = zeros(nbus,1);
0107 index = zeros(nbus,1);
0108 ischordal = 1;
0109 for i=1:nbus
0110     w = find(alpha == i);
0111     f(w) = w;
0112     index(w) = i;
0113     
0114     valid_v = find(Aadj(:,w));
0115     valid_v = valid_v(alpha(valid_v) < i);
0116     for vidx = 1:length(valid_v)
0117         v = valid_v(vidx);
0118         index(v) = i;
0119         if f(v) == v
0120             f(v) = w;
0121         end
0122     end
0123     
0124     for vidx = 1:length(valid_v)
0125         v = valid_v(vidx);
0126         if index(f(v)) < i
0127             ischordal = 0;
0128             break;
0129         end
0130     end
0131     
0132     if ~ischordal
0133         break;
0134     end
0135 end
0136 
0137 
0138 %% Determine maximal cliques
0139 
0140 % According to https://en.wikipedia.org/wiki/Chordal_graph
0141 % "To list all maximal cliques of a chordal graph, simply find a perfect
0142 % elimination ordering, form a clique for each vertex v together with the
0143 % neighbors of v that are later than v in the perfect elimination ordering,
0144 % and test whether each of the resulting cliques is maximal."
0145 
0146 % alpha is a candidate perfect elimination ordering. First form all
0147 % cliques. Then determine which cliques are maximal. Put maximal cliques in
0148 % the columns of the variable maxclique.
0149 
0150 if ischordal
0151     % Form a matrix representation of the cliques
0152     clique = speye(nbus,nbus);
0153     for i=1:nbus
0154         % neighbors of node i that are later in the ordering
0155         neighbors = find(Aadj(i,:)); 
0156         neighbors = neighbors(alpha(neighbors) > alpha(i));
0157 
0158         clique(i,neighbors) = 1;
0159     end
0160 
0161     % Test whether each clique is maximal. A clique is not maximal if a node
0162     % neighboring any of the nodes in the clique, but not in the clique, is
0163     % connected to all nodes in the clique.
0164     i=0;
0165     nclique = size(clique,1);
0166     while i < nclique
0167         i = i+1;
0168 
0169         neighbors = [];
0170         cliquei = find(clique(i,:));
0171         cliquei_bool = clique(i,:).';
0172         nnzi = sum(cliquei_bool);
0173 
0174         for k = 1:length(cliquei) % Check all nodes adjacent to nodes in the clique, but not included in the clique
0175             neighbors = [neighbors find(Aadj(cliquei(k),:))];
0176         end
0177         neighbors = unique(setdiff(neighbors,cliquei));
0178 
0179         not_maximal = 0;
0180         for m=1:length(neighbors)
0181             % If this neighbor is connected to all other nodes in the clique,
0182             % this is not a maximal clique
0183             if Aadj(neighbors(m),:) * cliquei_bool >= nnzi
0184                 not_maximal = 1;
0185                 break;
0186             end
0187         end
0188 
0189         if not_maximal
0190             clique = clique([1:i-1 i+1:nclique],:); % Delete this clique
0191             i = i-1; % After deletion, the next clique is really numbered i
0192             nclique = nclique - 1;
0193         end
0194     end
0195     maxclique = clique.'; % put maximal cliques in the columns
0196 else
0197     maxclique = nan; % maximal clique algorithm only works for chordal graphs
0198 end
maxCardSearch

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
