Daniel@0: function [Gs, op, nodes] = mk_nbrs_of_digraph(G0)
Daniel@0: % MK_NBRS_OF_DIGRAPH Make all digraphs that differ from G0 by a single edge deletion, addition or reversal
Daniel@0: % [Gs, op, nodes] = mk_nbrs_of_digraph(G0)
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: 
Daniel@0: debug = 0; % the vectorized version is about 3 to 10 times faster
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: % SINGLE EDGE DELETIONS
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: % 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: % if debug
Daniel@0: % % Non-vectorized version
Daniel@0: % ctr = 1;
Daniel@0: % for e=1:E
Daniel@0: %   i = I(e); j = J(e);
Daniel@0: %   Gdel2(:,:,ctr) = G0;
Daniel@0: %   Gdel2(i,j,ctr) = 0;
Daniel@0: %   ctr = ctr + 1;
Daniel@0: % end
Daniel@0: % assert(isequal(Gdel, Gdel2));
Daniel@0: % end
Daniel@0: 
Daniel@0: 
Daniel@0: % SINGLE EDGE REVERSALS
Daniel@0: 
Daniel@0: % rev_edge_ndx(k) is the scalar location of the k'th reversed edge
Daniel@0: %rev_edge_ndx = find(G0'); % different order to edge_ndx, which is bad
Daniel@0: rev_edge_ndx = subv2ind([n n], [J I]);
Daniel@0: % We set (rev_edge_ndx(k), k) to 1 for k=1:E in G1
Daniel@0: % We have already deleted i->j in the previous step
Daniel@0: ndx = subv2ind(size(Grep), [rev_edge_ndx(:) (1:E)']);
Daniel@0: G1(ndx) = 1;
Daniel@0: Grev = reshape(G1, [n n E]);
Daniel@0: 
Daniel@0: % if debug
Daniel@0: % % Non-vectorized version
Daniel@0: % ctr = 1;
Daniel@0: % for e=1:E
Daniel@0: %   i = I(e); j = J(e);
Daniel@0: %   Grev2(:,:,ctr) = G0;
Daniel@0: %   Grev2(i,j,ctr) = 0;
Daniel@0: %   Grev2(j,i,ctr) = 1;
Daniel@0: %   ctr = ctr + 1;
Daniel@0: % end
Daniel@0: % assert(isequal(Grev, Grev2));
Daniel@0: % end
Daniel@0: 
Daniel@0: 
Daniel@0: % SINGLE EDGE ADDITIONS
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: [Ibar,Jbar] = find(Gbar); 
Daniel@0: 
Daniel@0: bar_edge_ndx = find(Gbar);
Daniel@0: Ebar = length(Ibar); % num edges present in Gbar
Daniel@0: Grep = repmat(G0(:), 1, Ebar); % each column is a copy of G0
Daniel@0: ndx = subv2ind(size(Grep), [bar_edge_ndx(:) (1:Ebar)']);
Daniel@0: Grep(ndx) = 1;
Daniel@0: Gadd = reshape(Grep, [n n Ebar]);
Daniel@0: 
Daniel@0: % if debug
Daniel@0: % % Non-vectorized version
Daniel@0: % ctr = 1;
Daniel@0: % for e=1:length(Ibar)
Daniel@0: %   i = Ibar(e); j = Jbar(e);
Daniel@0: %   Gadd2(:,:,ctr) = G0;
Daniel@0: %   Gadd2(i,j,ctr) = 1;
Daniel@0: %   ctr = ctr + 1;
Daniel@0: % end
Daniel@0: % assert(isequal(Gadd, Gadd2));
Daniel@0: % end
Daniel@0: 
Daniel@0: 
Daniel@0: Gs = cat(3, Gdel, Grev, Gadd);
Daniel@0: 
Daniel@0: nodes = [I J;
Daniel@0: 	 I J;
Daniel@0: 	 Ibar Jbar];
Daniel@0: 
Daniel@0: op = cell(1, E+E+Ebar);
Daniel@0: op(1:E) = {'del'};
Daniel@0: op(E+1:2*E) = {'rev'};
Daniel@0: op(2*E+1:end) = {'add'};
Daniel@0: 
Daniel@0: 
Daniel@0: % numeric output:
Daniel@0: % op(i) = 1, 2, or 3, if the i'th neighbor was created by adding, deleting or reversing an arc.
Daniel@0: 
Daniel@0: ADD = 1;
Daniel@0: DEL = 2;
Daniel@0: REV = 3;
Daniel@0: 
Daniel@0: %op = [repmat(DEL, 1, E) repmat(REV, 1, E) repmat(ADD, 1, Ebar)];