comparison toolboxes/FullBNT-1.0.7/graph/mk_nbrs_of_digraph.m @ 0:e9a9cd732c1e tip

first hg version after svn

+function [Gs, op, nodes, A] = my_mk_nbs_of_digraph(G0,A)
+% MY_MK_NBRS_OF_DIGRAPH Make all digraphs that differ from G0 by a single edge deletion, addition or reversal, subject to acyclicity
+% [Gs, op, nodes, A] = my_mk_nbrs_of_digraph(G0,<A>)
+%
+% G0 is an adj matrix s.t. G0(i,j)=1 iff i->j in graph
+% A is the ancestor matrix for G0  (opt, creates if necessary)
+%
+% Gs(:,:,i) is the i'th neighbor
+% op{i} = 'add', 'del', or 'rev' is the operation used to create the i'th neighbor. 
+% nodes(i,1:2) are the head and tail of the operated-on arc.
+% Modified from mk_nbrs_of_digraph by Sonia Leach
+%
+% Modified by Sonia Leach Feb 02
+
+if nargin ==1, A = reachability_graph(G0');, end
+
+n = length(G0);
+[I,J] = find(G0); % I(k), J(k) is the k'th edge
+E = length(I);    % num edges present in G0
+
+
+% SINGLE EDGE DELETIONS
+% all deletions are valid wrt acyclity
+
+Grep = repmat(G0(:), 1, E); % each column is a copy of G0
+% edge_ndx(k) is the scalar location of the k'th edge 
+edge_ndx = find(G0);
+
+% edge_ndx = subv2ind([n n], [I J]); % equivalent
+% We set (ndx(k), k) to 0 for k=1:E in Grep
+ndx = subv2ind(size(Grep), [edge_ndx(:) (1:E)']);
+G1 = Grep;
+G1(ndx) = 0;
+Gdel = reshape(G1, [n n E]);
+
+
+% SINGLE EDGE REVERSALS
+
+% SML: previously Kevin had that legal structure was if
+% A(P,i)=1 for any P = { p | p in parents(j), p~=i}
+% specifically he said 
+%  "if any(A(ps,i)) then there is a path i -> parent of j -> j
+%   so reversing i->j would create a cycle"
+% Thus put in another way:
+%    for each i,j if sum(G0(:,j)' * A(:,i)) > 0, reversing i->j
+% is not legal.
+%
+% Ex. Suppose we want to check if 2->4 can be reversed in the 
+% following graph: 
+% G0 =                               A =
+%     0     0     1     0               0     0     0     0
+%     0     0     1     1               0     0     0     0
+%     0     0     0     1               1     1     0     0
+%     0     0     0     0               1     1     1     0
+% 
+% Then parents(4) = G0(:,4) = [0 1 1 0]'
+% and A(:,2) = [0 0 1 1]. Thus G0(:,4)'*A(:,2) = 1 b/c 3 is
+% an ancestor of 4 and a child of 2. Note that this works b/c
+% matrix multiplication has the effect of ANDing the two vectors 
+% and summing up the result (equiv. to the any(A(ps,i)) in kevin's code)
+%
+% So, we vectorize and check for all i,j pairs by looking for
+% 1's in L = (G0'*A)' which has L(i,j)=1 if rev(i,j) not legal
+% Note that this will give 1's where there are none in the G0
+% so we do a L=max(0, G0-L) to cancel out only the existing edges that 
+% aren't legal (subtracting where both are 1 and setting where
+% G0=0 and A=1 back to 0).
+
+L = max(0, G0-(G0'*A)');
+[IL, JL] = find(L);  % I(k), J(k) is the k'th legal edge to rev.
+EL = length(IL);
+
+
+% SML: First we have to DELETE THE EDGES WE ARE REVERSING
+% We can't use G1 w/ reversed edges already deleted (as
+% Kevin did) b/c the space of possible deletions are different 
+% now (some reverses aren't legal)
+
+Grep = repmat(G0(:), 1, EL); % each column is a copy of G0
+% edge_ndx(k) is the scalar location of the k'th edge 
+edge_ndx = subv2ind([n n], [IL JL]); 
+% We set (ndx(k), k) to 0 for k=1:E in Grep
+ndx = subv2ind(size(Grep), [edge_ndx(:) (1:EL)']);
+G1 = Grep;
+G1(ndx) = 0;
+
+% SML: Now we add in our REVERSED EDGES
+% rev_edge_ndx(k) is the scalar location of the k'th legal reversed edge
+rev_edge_ndx = subv2ind([n n], [JL IL]);
+
+% We set (rev_edge_ndx(k), k) to 1 for k=1:EL in G1
+% We have already deleted i->j in the previous step
+ndx = subv2ind(size(Grep), [rev_edge_ndx(:) (1:EL)']);
+G1(ndx) = 1;
+Grev = reshape(G1, [n n EL]);
+
+% SINGLE EDGE ADDITIONS
+
+% SML: previously Kevin had that any addition was legal if A(i,j)=0
+% however, you can not add i->j  if j is a descendent of i.
+% Thus, we create all possible additions in Gbar and then
+% subtract the descendants of each edge as possible parents
+% This means the potential parents of i (i.e. Gbar(:,i))
+% can not also be descendants if i i.e. (A(:,i)) which is accomplished
+% by subtracting (Gbar-A == 1 iff Gbar=1 & A=0)
+
+Gbar = ~G0;  % Gbar(i,j)=1 iff there is no i->j edge in G0
+Gbar = setdiag(Gbar, 0); % turn off self loops
+
+GbarL = Gbar-A;
+[IbarL, JbarL] = find(GbarL);  % I(k), J(k) is the k'th legal edge to add
+EbarL = length(IbarL);
+
+bar_edge_ndx = find(GbarL);
+
+Grep = repmat(G0(:), 1, EbarL); % each column is a copy of G0
+ndx = subv2ind(size(Grep), [bar_edge_ndx(:) (1:EbarL)']);
+Grep(ndx) = 1;
+Gadd = reshape(Grep, [n n EbarL]);
+
+
+Gs = cat(3, Gdel, Grev, Gadd);
+
+nodes = [I J;
+     IL JL;
+   IbarL JbarL];
+
+op = cell(1, E+EL+EbarL);
+op(1:E) = {'del'};
+op(E+(1:EL)) = {'rev'};
+op((E+EL+1):end) = {'add'};
+