function adj = random_modular_graph(d, c, p, mu, sigma)
%function adj = random_modular_graph(d, c, p, mu, sigma)   
% Generate random modular network introduced in Songdechakraiwut, T. Chung, 
% M.K. 2020 Topological learning for brain networks, arXiv: 2012.00675.
%
%
% INPUT
% d : number of nodes
% c : number of clusters/modules
% p : probability of connection within module
% mu, sigma : used for random edge weights
%
% OUTPUT
% adj : weighted adjacency matrix
%
%
% 
% %If you are using any part of the code, please reference the above paper.
%
% (C) 2020 Tananun Songdechakraiwut
%          University of Wisconsin-Madison
%
%  Contact tananun@cs.wisc.edu or mkchung@wisc.edu
%  for support/permission with the codes 
%
% Update history
%     November 11, 2020 created - Tananun Songdechakraiwut
%     December 10, 2020 edited Moo Chung
%


% adjacency matrix
adj = zeros(d);

% nodes are even distributed among modules
modules = cell(c, 1);
for k = 1:c
    modules{k} = round((k-1)*d/c+1):round(k*d/c);
end

for i = 1:d
    for j = i + 1:d
        module_i = ceil(c*i/d); % the module of node i
        module_j = ceil(c*j/d); % the module of node j

        % check if nodes i and j belongs to the same module
        if module_i == module_j
            % probability of attachment within module
            if rand <= p
                w = normrnd(mu, sigma);
                adj(i, j) = max(w, 0);
                adj(j, i) = max(w, 0);
            else
                % noise
                w = normrnd(0, sigma);
                adj(i, j) = max(w, 0);
                adj(j, i) = max(w, 0);
            end

        else
            % probability of attachment between modules
            if rand <= 1 - p
                w = normrnd(mu, sigma);
                adj(i, j) = max(w, 0);
                adj(j, i) = max(w, 0);
            else
                % noise
                w = normrnd(0, sigma);
                adj(i, j) = max(w, 0);
                adj(j, i) = max(w, 0);
            end

        end
    end
end

end