view 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'};