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

m0 = 3;
m = 4;
k0 = 2;
% Delcare 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);
% vector to hold degree sequences

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);
    % adjacency matrix, stored as a sparse matrix (type 'help sparse' for
    % more information
    
    for i=1:m0
        for j=1:m0
            if i ~= j && 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
        
        %         x = e(randi(length));
        %         y = e(randi(length));
        %         % pick two numbers from e
        %
        %         while x == y
        %             y = e(randi(length));
        %         end
        %         % make sure they are unique
        %
        %         e(length+1) = x;
        %         e(length+2) = y;
        %         e(length+3:length+3+m+k0) = i;
        %         % update edge list
        %
        %         length = length+3+m+k0;
        %         % update the length
        %
        %         A(x,i) = 1;
        %         A(y,i) = 1;
        %         % update the upper triangle of 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
    
    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 = 10*X.^(-2-(k0/m));
loglog(X,1-F,'bo',X,Y,'r-')
% plot the results