% Simulation of preferential attachment model (Chapter 2, Box C)
% Written by David Sumpter

%Number of positions (e.g. total number of papers, or sites for animals)
m=2000;

%Attachment constant
c=1;

%Exponent, always set to one in example in book.
ek=1

%Set up initial positions
x=zeros(m,1);
x(1:m)=1;

%Number of time steps
maxt=10000;

%Main loop
for t=1:maxt

    %Take a random position 
    pos=floor(rand*m)+1;
    %Stores how many were taken from that position and sets population at that position
    add(t)=x(pos);
    x(pos)=0;
    
    %Now add the removed individuals one by one to other positions.
    for rep=1:add(t)
        r=rand;
        i=1;
        p=x(i)+c;
        sx=sum((x+c).^ek);
        %Loop through until we find the position it should go in.
        while ((r>p/sx) & (i<m))
            i=i+1;
            p=p+(x(i)+c)^ek;
        end
        x(i)=x(i)+1;
    end

    mx(t)=max(x);

    %This just prints out how many time steps have gone.
    if (mod(t,100)==0)
        t 
    end

end


%Plots the frequency distribution
figure(1)
range=[1:1:4000]
h=histc(x,range)
hold on
x2=add(1000:maxt);
h=histc(x2,range)
hpl=plot(log(range),log(h./sum(h)),'kx')
set(hpl,'MarkerSize',10);
hold on
hpl=plot(log(range),-2.5*log(range),'k')
set(hpl,'MarkerSize',10);
hold off
axis([0 9 -11 0])
hlx=xlabel('log(Group size)');
hly=ylabel('log(Frequency)');
set(hlx,'FontSize',14);
set(hly,'FontSize',14);


%Plots the time series.
figure(2)
plot(add,'k')

hlx=xlabel('Time step')
hly=ylabel('Observed group size');
set(hlx,'FontSize',14);
set(hly,'FontSize',14);
axis([0 2000 0 200])