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