N = 10000;
% Declare network size (N.B. networks larger than 10^6 may require more
% memory than is available to Matlab...
runs = 10;
% Declare number of realisations

m0 = 5;
m = 5;
k0 = 4;
% Declare parameters
% k_0 == inital 'score' of a node
% m_0 == initial number of nodes
% m == number of connections each new node makes

if m0 < m
    m0 = m;
end
% the network must begin fully connected, so we cannot make more
% connections per step than there are nodes, hence the tacky solution:
% increase the number of initial nodes!

d = zeros(runs*N,1);
dk = zeros(runs*N,1);
dknni = zeros(runs*N,1);
% vector to hold degree sequences

p_coeff=zeros(runs,1);

for z=1:runs
    
    e = zeros(1000000,1);
    % edge list: clever method for creating a preferential attachement
    % network, each node appears k+k0 times in this edge list, we can then
    % pick a new node to
    length = 0;
    % actual length of edge list
    
    for i=1:m0
        for j=1:(k0+m0-1)
            e(j+length) = i;
        end
        length = length+(k0+m0-1);
    end
    % loop to populate the edge list with the initial condition
    
    A = sparse(N,N);
    % A = zeros(N);
    % adjacency matrix, stored as a sparse matrix (type 'help sparse' for
    % more information
    
    for i=1:m0
        for j=1:m0
            if i < j
                A(i,j) = 1;
            end
        end
    end
    % update the adjacency matrix to reflect the initial condition of the
    % network, which should be a fully connected network of m_0 nodes
    
    for i=m0+1:N
        
        R = e(randsample((length),m));
        % randomly pick m integers between 1 and (i-1) that will recieve
        % new connections from node i
        while size(unique(R),1)<m
            R = e(randsample((length),m));
        end
        % if there are fewer than m unique values, pick again
        % N.B. this is a TERRIBLE way of picking 4 random unique numbers
        % from e, in general one should NEVER use this approach...
        
        e(length+1:length+m) = R;
        % add new elements to the edge list corresponding to the selected
        
        e(length+m+1:length+m+m+k0) = i;
        % add the new node to the edge list
        
        length = length+m+k0+m;
        % update length of edge list
        
        A(R,i) = 1;
        % add new elements to the adjacency matrix
        
        
    end
    % loop to add new nodes and edges until the network size is N
    
    H = waitbar(z/runs);
    % cool loading bar :)
    
    A = triu(A) + triu(A)';
    % copy upper triangle into lower triangle
    
    deg = sum(full(A));
    
    B=full(A);
    
    degbis=B*deg';
    knni=degbis./deg';
    dk(z*N-(N-1):z*N) = deg;
    dknni(z*N-(N-1):z*N) = knni;
    
    p_coeff(z,1)=pearson(A);
    d(z*N-(N-1):z*N) = sum(A);
    % add degree sequence of current matrix to d
    
end

close(H);
% close loading bar :(

[F,X] = ecdf(d);
Y = 50*X.^(-2-(k0/m));
figure
loglog(X,1-F,'bo',X,Y,'r-')
figure
plot(dk,dknni,'bx');
figure
plot([1:runs],p_coeff,'bx');