%function M=Lrandnetwork(N,p,pow,rotlen)
%
%  Creates a network, where directed loops of length rotlen are promoted.
%
%  Input:
%    N - number of nodes
%    p - connection probability OR a vector of length N denoting the
%        probabilities of each number inputs
%    pow - the strength of the loop promotion [0,inf]
%    rotlen - the promoted loop length
%
%  Output:
%    M - the connectivity matrix such that M(i,j) denotes the existence of
%        an edge from i to j
%
%  Tuomo Mäki-Marttunen
%  Last modified 8.1.2013

function M=Lrandnetwork(N,p,pow,rotlen)

  M = zeros(N,N);
  if nargin < 3 || isempty(pow)
    pow = 1;
  end
  if nargin < 4 || isempty(rotlen)
    rotlen = 3;
  end
  if numel(p)==1  
    nns = binornd(N-1,p,N,1);  
  elseif numel(p)==N
    nns = zeros(N,1);
    %p(i) gives the probability that number of in-neighbours is i-1 (i=1,...,N)
    p = p/sum(p);
    pcs = cumsum(p);
    for i=1:N
      [vain,nns(i)] = max(rand() <= pcs); %[~,nns(i)] not supported in all versions 
    end
    nns = nns - 1;
  end
  
  pinds = cell(N);
  for i=1:N
      pinds{i} = [1:i-1,i+1:N];
  end
  nplaced = zeros(N,1);
  
  lastinput = 0; %this tells which node was last chosen to be an input
  for ie=1:sum(nns) %go through all the edges in a non-obvious order
    
    if lastinput %if some node has already been chosen as an input
        if nplaced(lastinput) < nns(lastinput) %if this node has to have more inputs
            i = lastinput; %choose that one
        else %otherwise pick by random
            w = nns - nplaced;
            w = w/sum(w);
            r = rand(1) < cumsum(w);
            i = find(r,1);
        end
    else %otherwise pick by random
        w = nns - nplaced;
        w = w/sum(w);
        r = rand(1) < cumsum(w);
        i = find(r,1);
    end
    
    pind = pinds{i};
    Npaths = nminpath(M,i,rotlen); %a matrix that tells the number of paths from i to other nodes

    % The points are given such that the highest points (>4) will be for
    % those nodes, to which there is a path from i of length rotlen-1, but no path
    % shorter than that. Choosing a such node as input will close the loop
    % of length rotlen. The second largest points (>3) is given to those
    % nodes, that are not reachable from i with a path of length rotlen-1,
    % but on the other hand not with any shorter path either. Choosing such
    % a node will not close a loop of length <= rotlen. The second lowest
    % points (>2) is given to nodes that are reachable from i with rotlen-1
    % edges, but which are also reachable by shorter paths. Choosing such a
    % node will close a loop of length < rotlen. The least points (1) are given
    % to the remaining nodes.
    points = (Npaths(end,:)&~any(Npaths(1:end-1,:),1)).*(4+Npaths(end,:)./((N-2)/(rotlen-2))^(rotlen-2)) +...
             (~Npaths(end,:)&~any(Npaths(1:end-1,:),1))*3 +...
             (Npaths(end,:)&any(Npaths(1:end-1,:),1)).*(2+Npaths(end,:)./((N-1)^(rotlen-2))) +...
             (~Npaths(end,:)&any(Npaths(1:end-1,:),1));
    if ~isinf(pow)
      weights = points(pind).^pow;
    else
      weights = points(pind) == max(points(pind));
    end
    if any(isinf(weights)) %if the weight of any node is inf, the choice will be made by random between such nodes
      weights = isinf(weights);
    end
    if all(weights==0)
      weights = weights+1;
    end
    weights = weights/sum(weights);
    r = rand(1) < cumsum(weights');
    rind = find(r,1);
    M(pind(rind),i) = 1;
    lastinput = pind(rind);
    nplaced(i) = nplaced(i) + 1;
    pind(rind) = [];
    pinds{i} = pind;
  end
end

%function nminpath returns a matrix that tells by how many paths you can
% reach nodes 1...N (horizontal) using 1...rotlen-1 (vertical) edges from
% node i
function Npaths = nminpath(M,i,rotlen) 
  N = size(M,1);
  Npaths = zeros(rotlen-1,N);
  Npaths(1,:) = M(i,:);
  
  for j=2:rotlen-1
    Npaths(j,:) = Npaths(j-1,:)*M;
  end
end