%% CODE DESCRIPTION
% This script plots the probability that the link i<->j exists at time t, i.e., 
%P_ij(t) = Prob(A_{ij}(t) = 1), for all i,j=1,...,n, in a bistable regime.

% We compute M paths of Gillespie simulation and copmpute the probability
% as a the frequence

% 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
% tot_edges: number of existing edges at time 0
% T: final time of simulations
% Dt: stepsize of the simulation
% M: number of Gillespie paths


% Initial config: we choose at random tot_edges links

clear all
clc

rng('shuffle')
%% 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
c_3_hat = c_3/(n-2);
T = 250; % final time of the simulation
Dt = 1;
tem = [Dt:Dt:T];
L=T/Dt;
M=250; %& numebr of paths

%% Initial configuration
A_in = zeros(n,n);
tot_edges = 217;
nedges = 0;
while nedges < tot_edges
    i = floor(rand*n)+1;
    j = floor(rand*n)+1;
    if i ~= j && A_in(i,j) == 0
            nedges = nedges+1;
            A_in(i,j) = 1;
            A_in(j,i) = 1;
    end
end

prob_vec = zeros(L+1,N); % probability matrix, times x prob 
in = 1;
for i=1:n-1
    for j=i+1:n
        prob_vec(1,in) = A_in(i,j);
        in = in+1;
    end
end

steps = length(tem); % number of time steps
summation = zeros(steps,N); % sum of ones for every edges for different paths
temp = zeros(steps,N);

%% Loop for M paths
for ind=1:M
    
    ind
    t = 0;
    A = A_in;
    %% Loop for Gillespie simulations
    for m=1:steps 
        
        
        Rmax = 5*10^2;   % stop after this many reactions
        C = ones(n,n) - eye(n,n);    %represents a full graph
      

        %% Gillespie simulation
        while t<tem(m)

            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

           t=t+tau; % upodate the current t

        end % end of Gillepsie simulation
        
        in = 1;
        for i1=1:n-1
            for j1=i1+1:n
                summation(m,in) = A(i1,j1);
                in = in+1;
            end
        end

    end % end loop for differenti gillespie simulations
    
    
  temp = summation+temp;
    
    
    
end % end cycle for M paths
prob_vec(2:end,:) = temp/M;
plot([0,tem],prob_vec)

%% Superimposition of the E[P_{ij(t)]
hold on
plot([0,tem],mean(prob_vec'))





