%function M=FFrandnetwork(N,p,pow)
%
% Creates a (recurrent) network, where feedforward loops 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 feedforward loop promotion [0,inf]
%
% Output:
% M  the connectivity matrix such that M(i,j) denotes the existence of
% an edge from i to j
%
% Tuomo MäkiMarttunen
% Last modified 8.1.2013
function M=FFrandnetwork(N,p,pow)
M = zeros(N,N);
if nargin < 3  isempty(pow)
pow = 1;
end
if numel(p)~=1
if numel(p) ~= N, error('The vector p should be of size N!'); end
%p(i) gives the probability that number of inneighbours is i1 (i=1,...,N)
p = p/sum(p);
pcs = cumsum(p);
end
for i=1:N %go through all nodes in order
if numel(p)==1
nn = binornd(N1,p);
else
[vain,nn] = max(rand() <= pcs); %[~,nns] not supported in all versions
nn = nn  1;
end
pind = [1:i1,i+1:N]; %possible inputs are all but the node itself
points = ones(N,1); % initially 1 point for each node
lastinput = 0; %this tells which node was last chosen to be an input
for j=1:nn %go through all possible inputs in order
% instead of recalculating the weight for all nodes, the
% recalculation can be done only for those nodes that project to the
% last chosen input and do not yet project to i
if lastinput && any(M(:,lastinput) & ~M(:,i))
points = points + (M(:,lastinput) & ~M(:,i)); % points added to the nodes that are now
% reaching i with two hops
end
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 = [0; rand(1) < cumsum(weights)];
rind = ~~(r(2:end)r(1:end1));
M(pind(rind),i) = 1;
lastinput = pind(rind);
pind(rind) = [];
end
end
