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