%% CODE DESCRIPTION
% This code provides a comparison between three models: micro model, macro
% model (Markov process) and mean-field approzimation.

% List of parameters
% n: number of nodes
% N: number of possible edges
% c_1: birth constant
% c_2: death constant
% c_3: tradic closure constant
% c_3_hat: size dependent rate for micro
% c_3_tilde: size dependent rate for macro
% T: final time of the simulations
% Dt: stepsize for mean-field solution
% tot_edges: number of existing edges at initial time
% p: E--R parameter


%Initial conditions: E--R graph
%tot_edges edges.

clear all
clc

rng('shuffle')

%%Parameter
n=100;
N=n*(n-1)*0.5;
%c_1 = 0.25; % parameter for monostable regime
c_1 = 0.025; %parameter for bistable regime
c_2 = 0.25;
c_3 = 0.91;


%% Setup of the initial configuration
p=0.2;
A0 = zeros(n,n); % A stands for the current adjacency matrix
for i = 1:n-1
  for j = i+1:n
     A0(i,j) = rand < p;
     A0(j,i) = A0(i,j);
  end
end

%% Micro model simulation: Gillespie path
A = A0;
microval(1) = sum(sum(A))*0.5/N;
T = 5000;
c_3_hat = c_3/(n-2);
t=0; % time counter
count=0;
Rmax = 5*10^2;   % stop after this many reactions
C = ones(n,n) - eye(n,n);    %represents a full graph
timeval(1)=0;   %time of each reaction

% Gillespie code
while t<T   
  
    Asq = A*A;   %number of common neighbours
  
    % First we will decide which of the three types
    % of reaction is going to fire next.
    % Evaluate the propensity for each of the 
    % three types.
    % Spontaneous death: omega for every edge   
    a(1) = sum(sum(c_2*A))/2; 
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Spontaneous birth: delta for every non-edge
    a(2) = sum(sum(c_1*(C-A)))/2;
     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Triadic closure-driven birth:  
    %  epsilon * # common neighbs for every non-edge
    a(3) = sum(sum(c_3_hat*Asq.*(C-A)))/2;
     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    asum = sum(a);
    
    j = min(find(rand<cumsum(a/asum)));
    tau = log(1/rand)/asum;
    %tau=tau*2;
    if j==1 %spontaneous death
        [x1,x2]=find(A==1);     % locate existing edges
        zz=find((x1<x2)==1);  % restrict to upper triangle 
        xx1=x1(zz);                 % row indices
        xx2=x2(zz);                 % column indices
        z=length(xx1);             % number to choose from
        index=randi(z);            % choose one at random
        % We now have an edge to remove
        A(xx1(index),xx2(index))=A(xx1(index),xx2(index))-1;
        A(xx2(index),xx1(index))=A(xx2(index),xx1(index))-1;
    elseif j==2 %spontaneous birth
        [x1,x2]=find(A==0);    % locate existing non-edges
        zz=find((x1<x2)==1); % restrict to upper triangle
        xx1=x1(zz);                % row indices
        xx2=x2(zz);                % column indices
        z=length(xx1);            % number to choose from
        index=randi(z);            % choose one at random
        % We now have a non-edge that will become an edge
        A(xx1(index),xx2(index))=A(xx1(index),xx2(index))+1;
        A(xx2(index),xx1(index))=A(xx2(index),xx1(index))+1;
    else %birth by triadic closure 
        [x1,x2]=find(A==0&Asq~=0); % locate candidate non-edges;
                                                    % these are non-edges
                                                    % where the two nodes
                                                    % that have a least one
                                                    % common neighbour.
        zz=find((x1<x2)==1);            % restrict to upper triangle
        xx1=x1(zz);                           % row indices
        xx2=x2(zz);                           % column indices
        % These candidate non-edges do not have a uniform 
        % probability of being chosen. Probability is 
        % proportional to number of common neigjbours.
        % So we need to check how many common neighbours 
        % exist for each pair of nodes forming a non-edge.
        weights=Asq(xx1,xx2);  % Number of common neighbours.
        s_w=sum(weights);        % Compute the sum, so that 
                                            % we can work out proportions.
        as_w=sum(s_w);            %
                                            % Each candidate non-edge 
                                            % has a probability of being
                                            % chosen that is portional to
                                            % its contribition to the overall 
                                            % sum:
        index=min(find(rand<cumsum(s_w/as_w)));
        % We now have a non-edge that will be created
        % according to the triadic closure principle.
        A(xx1(index),xx2(index))=A(xx1(index),xx2(index))+1;
        A(xx2(index),xx1(index))=A(xx2(index),xx1(index))+1;
        
     end
   
   count=count+1;    %update number of reactions
   microval(count+1) = sum(sum(A))*0.5/N; %update medge density 
   
 
   t=t+tau; % upodate the current t
   
   timeval(count+1)=t;    
      
end
timeval(end)=T;

%% Macro model simulation
V=[-1 1 1];
c_3_tilde = c_3/(N^2);
X(1) = sum(sum(A0))*0.5;
%intensity functions
tvals(1) = 0; % time for macro simulation
t=0;
count=1;


while t < T
     a(1) = c_2*X(count);
     a(2) = c_1*(N-X(count));
     a(3) = c_3_tilde*X(count)*(X(count)-1)*(N-X(count));
     asum = sum(a);
     j = min(find(rand<cumsum(a/asum)));
     tau = log(1/rand)/asum;
     X(count+1) = X(count) + V(:,j);
      
     count = count + 1;
     t = t + tau;
     tvals(count) = t;
end
tvals(end) = T;
macroval = X/N;



%% Mean-field approzimation
Dt = 1;
L=T/Dt;
mf = zeros(L+1,1);
mf(1) = macroval(1);

%% loop for time
for i=1:L
    mf(i+1) = (1-c_2)*mf(i)+(1-mf(i))*(c_1+c_3*(mf(i)^2));
end


%% Plot 
plot(timeval,microval,'red','Linewidth',2)
hold on
plot(tvals,macroval,'Linewidth',2)
plot([0:Dt:T],mf,'black','Linewidth',3)