Mercurial > hg > camir-aes2014

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

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
first hg version after svn
author wolffd
date Tue, 10 Feb 2015 15:05:51 +0000
parents
children
rev   line source
wolffd@0 1 function [Gs, op, nodes, A] = my_mk_nbs_of_digraph(G0,A)
wolffd@0 2 % MY_MK_NBRS_OF_DIGRAPH Make all digraphs that differ from G0 by a single edge deletion, addition or reversal, subject to acyclicity
wolffd@0 3 % [Gs, op, nodes, A] = my_mk_nbrs_of_digraph(G0,<A>)
wolffd@0 4 %
wolffd@0 5 % G0 is an adj matrix s.t. G0(i,j)=1 iff i->j in graph
wolffd@0 6 % A is the ancestor matrix for G0 (opt, creates if necessary)
wolffd@0 7 %
wolffd@0 8 % Gs(:,:,i) is the i'th neighbor
wolffd@0 9 % op{i} = 'add', 'del', or 'rev' is the operation used to create the i'th neighbor.
wolffd@0 10 % nodes(i,1:2) are the head and tail of the operated-on arc.
wolffd@0 11 % Modified from mk_nbrs_of_digraph by Sonia Leach
wolffd@0 12 %
wolffd@0 13 % Modified by Sonia Leach Feb 02
wolffd@0 14
wolffd@0 15 if nargin ==1, A = reachability_graph(G0');, end
wolffd@0 16
wolffd@0 17 n = length(G0);
wolffd@0 18 [I,J] = find(G0); % I(k), J(k) is the k'th edge
wolffd@0 19 E = length(I); % num edges present in G0
wolffd@0 20
wolffd@0 21
wolffd@0 22 % SINGLE EDGE DELETIONS
wolffd@0 23 % all deletions are valid wrt acyclity
wolffd@0 24
wolffd@0 25 Grep = repmat(G0(:), 1, E); % each column is a copy of G0
wolffd@0 26 % edge_ndx(k) is the scalar location of the k'th edge
wolffd@0 27 edge_ndx = find(G0);
wolffd@0 28
wolffd@0 29 % edge_ndx = subv2ind([n n], [I J]); % equivalent
wolffd@0 30 % We set (ndx(k), k) to 0 for k=1:E in Grep
wolffd@0 31 ndx = subv2ind(size(Grep), [edge_ndx(:) (1:E)']);
wolffd@0 32 G1 = Grep;
wolffd@0 33 G1(ndx) = 0;
wolffd@0 34 Gdel = reshape(G1, [n n E]);
wolffd@0 35
wolffd@0 36
wolffd@0 37 % SINGLE EDGE REVERSALS
wolffd@0 38
wolffd@0 39 % SML: previously Kevin had that legal structure was if
wolffd@0 40 % A(P,i)=1 for any P = { p | p in parents(j), p~=i}
wolffd@0 41 % specifically he said
wolffd@0 42 % "if any(A(ps,i)) then there is a path i -> parent of j -> j
wolffd@0 43 % so reversing i->j would create a cycle"
wolffd@0 44 % Thus put in another way:
wolffd@0 45 % for each i,j if sum(G0(:,j)' * A(:,i)) > 0, reversing i->j
wolffd@0 46 % is not legal.
wolffd@0 47 %
wolffd@0 48 % Ex. Suppose we want to check if 2->4 can be reversed in the
wolffd@0 49 % following graph:
wolffd@0 50 % G0 = A =
wolffd@0 51 % 0 0 1 0 0 0 0 0
wolffd@0 52 % 0 0 1 1 0 0 0 0
wolffd@0 53 % 0 0 0 1 1 1 0 0
wolffd@0 54 % 0 0 0 0 1 1 1 0
wolffd@0 55 %
wolffd@0 56 % Then parents(4) = G0(:,4) = [0 1 1 0]'
wolffd@0 57 % and A(:,2) = [0 0 1 1]. Thus G0(:,4)'*A(:,2) = 1 b/c 3 is
wolffd@0 58 % an ancestor of 4 and a child of 2. Note that this works b/c
wolffd@0 59 % matrix multiplication has the effect of ANDing the two vectors
wolffd@0 60 % and summing up the result (equiv. to the any(A(ps,i)) in kevin's code)
wolffd@0 61 %
wolffd@0 62 % So, we vectorize and check for all i,j pairs by looking for
wolffd@0 63 % 1's in L = (G0'*A)' which has L(i,j)=1 if rev(i,j) not legal
wolffd@0 64 % Note that this will give 1's where there are none in the G0
wolffd@0 65 % so we do a L=max(0, G0-L) to cancel out only the existing edges that
wolffd@0 66 % aren't legal (subtracting where both are 1 and setting where
wolffd@0 67 % G0=0 and A=1 back to 0).
wolffd@0 68
wolffd@0 69 L = max(0, G0-(G0'*A)');
wolffd@0 70 [IL, JL] = find(L); % I(k), J(k) is the k'th legal edge to rev.
wolffd@0 71 EL = length(IL);
wolffd@0 72
wolffd@0 73
wolffd@0 74 % SML: First we have to DELETE THE EDGES WE ARE REVERSING
wolffd@0 75 % We can't use G1 w/ reversed edges already deleted (as
wolffd@0 76 % Kevin did) b/c the space of possible deletions are different
wolffd@0 77 % now (some reverses aren't legal)
wolffd@0 78
wolffd@0 79 Grep = repmat(G0(:), 1, EL); % each column is a copy of G0
wolffd@0 80 % edge_ndx(k) is the scalar location of the k'th edge
wolffd@0 81 edge_ndx = subv2ind([n n], [IL JL]);
wolffd@0 82 % We set (ndx(k), k) to 0 for k=1:E in Grep
wolffd@0 83 ndx = subv2ind(size(Grep), [edge_ndx(:) (1:EL)']);
wolffd@0 84 G1 = Grep;
wolffd@0 85 G1(ndx) = 0;
wolffd@0 86
wolffd@0 87 % SML: Now we add in our REVERSED EDGES
wolffd@0 88 % rev_edge_ndx(k) is the scalar location of the k'th legal reversed edge
wolffd@0 89 rev_edge_ndx = subv2ind([n n], [JL IL]);
wolffd@0 90
wolffd@0 91 % We set (rev_edge_ndx(k), k) to 1 for k=1:EL in G1
wolffd@0 92 % We have already deleted i->j in the previous step
wolffd@0 93 ndx = subv2ind(size(Grep), [rev_edge_ndx(:) (1:EL)']);
wolffd@0 94 G1(ndx) = 1;
wolffd@0 95 Grev = reshape(G1, [n n EL]);
wolffd@0 96
wolffd@0 97 % SINGLE EDGE ADDITIONS
wolffd@0 98
wolffd@0 99 % SML: previously Kevin had that any addition was legal if A(i,j)=0
wolffd@0 100 % however, you can not add i->j if j is a descendent of i.
wolffd@0 101 % Thus, we create all possible additions in Gbar and then
wolffd@0 102 % subtract the descendants of each edge as possible parents
wolffd@0 103 % This means the potential parents of i (i.e. Gbar(:,i))
wolffd@0 104 % can not also be descendants if i i.e. (A(:,i)) which is accomplished
wolffd@0 105 % by subtracting (Gbar-A == 1 iff Gbar=1 & A=0)
wolffd@0 106
wolffd@0 107 Gbar = ~G0; % Gbar(i,j)=1 iff there is no i->j edge in G0
wolffd@0 108 Gbar = setdiag(Gbar, 0); % turn off self loops
wolffd@0 109
wolffd@0 110 GbarL = Gbar-A;
wolffd@0 111 [IbarL, JbarL] = find(GbarL); % I(k), J(k) is the k'th legal edge to add
wolffd@0 112 EbarL = length(IbarL);
wolffd@0 113
wolffd@0 114 bar_edge_ndx = find(GbarL);
wolffd@0 115
wolffd@0 116 Grep = repmat(G0(:), 1, EbarL); % each column is a copy of G0
wolffd@0 117 ndx = subv2ind(size(Grep), [bar_edge_ndx(:) (1:EbarL)']);
wolffd@0 118 Grep(ndx) = 1;
wolffd@0 119 Gadd = reshape(Grep, [n n EbarL]);
wolffd@0 120
wolffd@0 121
wolffd@0 122 Gs = cat(3, Gdel, Grev, Gadd);
wolffd@0 123
wolffd@0 124 nodes = [I J;
wolffd@0 125 IL JL;
wolffd@0 126 IbarL JbarL];
wolffd@0 127
wolffd@0 128 op = cell(1, E+EL+EbarL);
wolffd@0 129 op(1:E) = {'del'};
wolffd@0 130 op(E+(1:EL)) = {'rev'};
wolffd@0 131 op((E+EL+1):end) = {'add'};
wolffd@0 132