function [ A , reglatt ] = gen_sw( N , k_init , p , quiet)
%GEN_SW  Generates a small world network
%
%  GEN_SW(N,k_init,p) uses the method proposed by Watts and
%      Strogatz (1998) to generate a network with the small
%      world character (if the parameters are chosen 
%      appropriately):
%      N  is total number of nodes, k_init is the initial 
%      node degree of the network (k_init must be pair!, so 
%      that each node is connected to its 1st ... k_init-th
%      neighbours), and  p  is the rewiring probability
%
%  Last change:  10/2/2005 - Florian Knorn

if nargin < 3
	error('Please provide enough input arguments !');
end

if p < 0 || p > 1
	error('p must be a probability, so 0 <= p <= 1 !');
end

if mod(k_init,2) || k_init==0
	error('k must be pair and nonzero !');
end

if mod(N,2) && N-k_init < 1
	error('k too large');
end

if ~mod(N,2) && N-k_init < 2
	error('k too large');
end

if nargin ~= 4
	quiet = 0;
end

rand('state',sum(100*clock));		% reset random number generator to a different state

% step 1 : generate regular 1D-lattice (upper triangle part of A only -> reduces work)
A = logical(sparse(N,N));

for k=1:k_init/2
	A = A | spdiags(ones(N,1),k,N,N); A = A | spdiags(ones(N,1),(N-k),N,N);
end

% step 2 : rewire by first looking at each edge and then deciding whether to rewire
[i,j] = find(A); rewire_counter = 0;

for curr=1:length(i)
	
	if rand <= p       % go ahead, rewire !
		rewire_counter = rewire_counter+1;
		
		% now try to find a spot for the new egde ...
		j_attempt = randperm(N);  % get a random permuation for the column indices
		l = 1;%    ... which is nonoccupied   AND  not on the main diagonal
		while A(i(curr),j_attempt(l)) || A(j_attempt(l),i(curr)) ...
				|| ( j_attempt(l) == i(curr) )
			l = l+1;
		end

		% found candidate -> new one
		if i(curr) > j_attempt(l) % edge would end up in low triangle part -> "transpose"
			A(j_attempt(l),i(curr)) = true;
		else
			A(i(curr),j_attempt(l)) = true;
		end

		% and removing old edge
		A(i(curr),j(curr)) = false;
	end
end

% create full adjacency matrix by adding the transposed of the (up to now) upper triangle
A = A | A';

if rewire_counter
	reglatt = false;
else
	reglatt = true;
end

if ~quiet
	disp([ ...
		sprintf('Small-World-Network of %i nodes and %i edges created',N,nnz(A)/2) ...
		sprintf('\nusing k_init = %i and p = %f.',k_init,p) ...
		sprintf('  %i edges rewired.',rewire_counter)]);
end