% Simulation of Niwa's model (Chapter 2, Box B)
% Written by David Sumpter, corrections by Qi Ma

% Number of sites 
s=2000;
k=1;

%Initialise sites, one individual per site.
x=zeros(1,s);
x(1,1:s)=1;

%Split probability
p=0.3;

%Number of time steps
maxt=2000;
%Keeps track of how many individuals are at site 1.
track=zeros(maxt,1);


for t=1:maxt

    y=zeros(1,s);

    for site=1:s

        r=rand;

        if ((x(site)>0) & (r<p))
            %A split has occurred at this site.
            %Determine the size of the two splitting groups
            n1=floor(rand*x(site))+1;
            n2=x(site)-n1;  
            %Put the groups at two new sites chosen at random
            site1=floor(rand*s)+1;
            site2=floor(rand*s)+1;
            y(site1)=y(site1)+n1;
            y(site2)=y(site2)+n2;
        else
            %Whole group moves to a new randomly chosen site.
            site1=floor(rand*s)+1;
            y(site1)=y(site1)+x(site);
        end
    end 
    
    %Update sites.
    x=y;
    track(t)=x(1);
    
    %Prints out the time step.
    if (mod(t,1000)==0)
        t 
    end

end

%Calculate Niwa's predicted distribution.
range=[1:1:2000];
h=histc(track(1:maxt),range);
P1=h'.*range/sum(h'.*range);
NP=sum(P1.*range)
fitw=exp(-range/NP.*(1-exp(-range/NP)/2))./range;
fitw=fitw/sum(fitw);

%Plot simulation vs analytic prediction.
figure(1)
hpl=plot(log(range),log(h/sum(h)),'kx')
set(hpl,'MarkerSize',10);
hold on
hpl=plot(log(range),log(fitw),'k')
set(hpl,'MarkerSize',10);
hold off
axis([0 4.5 -7 0])
hlx=xlabel('log(Group size)');
hly=ylabel('log(Frequency)');
set(hlx,'FontSize',14);
set(hly,'FontSize',14);

%Plot time series of group sizes at site 1. (since the individuals move each time step this is a 'randomly chosen' site.
figure(2)
plot(track(1:maxt),'k')
hlx=xlabel('Time step')
hly=ylabel('Observed group size');
set(hlx,'FontSize',14);
set(hly,'FontSize',14);
axis([0 maxt 0 max(track)+10])