%% CODE DESCRIPTION
% This script plots the spy of the adjacency matrix A(t) in a bistable regime at four
% different time levels: at the initial time, at the time of high level, at
% the time of low level, at the final time.

% The script impolements the Gillespie algorithm to obtain A(t)

% The roots are: 0.1718, 0.3066, 0.5216

% List of parameters
% n: number of nodes
% N: number of possible edges
% c_1: birth rate
% c_2: death rate
% c_3: triadic closure effect
%c_3_hat: size dependent rate
% T: final time of simulations
% p: initial probability for edge to exist

% Initial conditions: E-R condiotions



clear all
clc

rng(100) % setting of the random number generator

%% Parameter setting
n=30; % number of nodes
N=n*(n-1)*0.5; % number of total edges
c_1 = 0.025; %birth constant
c_2 = 0.25; % death constant
c_3 = 0.91; % traidic closur constant
c_3_hat = c_3/(n-2); % size-dependent rate
T = 10000; % final time of the simulation


%%Initializations
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

%% Initial configurations
p=0.3;
A = zeros(n,n); % A stands for the current adjacency matrix
for i = 1:n-1
  for j = i+1:n
     A(i,j) = rand < p;
     A(j,i) = A(i,j);
  end
end

microval(1)=sum(sum(A))*0.5/N; % initial density

i=1; % counter for the plots
ban_low = 0; % ban for checking the spy of the matrix
ban_up =  0;
spy_time(1) = 0; % time of the spy


%% Plot of the initial density of A
 figure(i)
 spy(A)
 i=i+1;


%% Gillespie simulation
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;    
   
   %% spy of the mnatrix
   if microval(count+1)>0.17 && microval(count+1)<0.18 && ban_low==0
       figure(i)
       spy(A)
       ban_low = 1;
       spy_time(i) = t;
       i=i+1;
   end
   
    if microval(count+1)>0.5 && microval(count+1)<0.6 && ban_up==0
       figure(i)
       spy(A)
       ban_up=1;
       spy_time(i) = t;
       i=i+1;
   end
   
    
end
timeval(end)=T;
figure(4)
spy(A)
spy_time(4) = T;






