function A = rwe( A , pn )
%RWE  randomly rewires egdes from graph
%
%   RWE(A,pn)  rewires edges from the network (given by its
%      adjacency matrix A).
%      The no. of rewired edges depends on  pn:
%      - if 0 < pn < 1  its the *probability* rewiring (per edge)
%      - if pn >= 1  its the *number* of edges rewired
%
%  Last change:  10/2/2005 - Florian Knorn

if nargin < 2
	error('Please provide adjacency matrix A and probability p !');
else
	n = size(A,1);
% 	if n ~= size(A,2)
% 		error('Please input SQUARE matrix !');
% 	end
	if pn < 0 || (pn >= 1 && round(pn) ~= pn)
		error('pn must be a probab. (0 <= pn < 1) or an int >= 1 !');
	end
end

A = logical(sparse(triu(A,1))); % only keep upper triang part

if pn < 1 % we're dealing with probabilities
	pn = binornd(nnz(A),pn);
end

ind = find(A);      % get 1D-indexes of edges in triu(A)
lind = length(ind); % and their count
if pn > lind
	error('Cannot rewire more egdes then actually exist !');
end

indperm = ind(randperm(lind)); % randomly reorder elements of ind
curr = 1;

while curr <= pn % rewire pn edges

	% get 2D-subscr of curr edge and create rand perm for column ind
	[subi,subj] = ind2sub([n n],indperm(curr)); 
	j_attempt = randperm(n);

	% now try to find a spot for the new egde that is ...
	l = 1;%     ... nonoccupied                          AND  not on the main diag
	while A(subi,j_attempt(l)) || ( j_attempt(l) == subi ) || A(j_attempt(l),subi)
		l = l+1;
	end

	% found spot -> create new edge ...
	if subi > j_attempt(l)% edge would end up in lower triangle part
		A(j_attempt(l),subi) = true;
	else
		A(subi,j_attempt(l)) = true;
	end

	% ... and remove the old one and increase counter
	A(subi,subj) = false; curr = curr + 1;

end

A = A | A';