% Simulation of Sibly's model (Chapter 2, Box A)
% Written by David Sumpter

%Number of time steps
time=3000;

%This is the parameter in the optimal group size curve. 
beta=0.1;

sites=2000;

%This is the range 
maxnum=4;

%This term allows an error to be introduced in individual's assement of how many others there are the site.
error=0;

%Initial group size configuration.
x=zeros(time,sites);
x(1,:)=[8+floor(rand(1,sites*1/2)*(maxnum+1)) zeros(sites*1/2,1)'];

for t=1:time
    x(t+1,:)=x(t,:);
    
    %On each time step 10 individuals move to a new site.
    for it=1:10
        %Evaluate optimal group size fitness function
        f=(x(t+1,:)+1).*exp(-beta*(x(t+1,:)+1)).*(1+rand(1,length(x(1,:)))*error);
        
        %Pick a random individual
        r=ceil(rand*sum(x(t+1,:)));
        i=1;
        count=x(t+1,i);
        while (count<r)
            i=i+1;
            count=count+x(t+1,i);
        end
        
        %Put the individual at the site with the maximum payoff.
        %If more than one site is maximal then choose one of those at random.
        funp=unique((f==max(f)).*[1:sites]);
        fp=funp(ceil(rand*(length(funp)-1))+1);
        fm=f(fp);
        if fm>f(i)
            x(t+1,i)=x(t+1,i)-1;
            x(t+1,fp)=x(t+1,fp)+1;
        end
    end
    
    %Make a histogram of the positions of the individuals
    h=[0:1:50];
    xh(t,:)=histc(x(t,:),h);
    
    %Output current time step.
    if (mod(t,100)==0)
        t 
    end
end

%Shows how the group size distribution changes with time.
figure(1)
colormap([1 [0.9:-1/63:0] ; 1 [0.9:-1/63:0] ; 1 [0.9:-1/63:0]]')
imagesc(h,[0:10:10*(size(xh,1)-1)],xh,[0 800] );
axis([0 25 0 3000])
hlx=ylabel('Time');
hly=xlabel('Group size');
set(hlx,'FontSize',14);
set(hly,'FontSize',14);

figure(2)
hold off
bar(h(2:length(h)),xh(t,2:length(h))./sum(xh(t,2:length(h))),'w')
hold on
h=[0:50];
fh=0.1*h.*exp(-beta*(h));
plot(h,fh,'kx')
plot([0 50],[fh(2) fh(2)],'k:')
axis([0 50 0 0.8])
hlx=ylabel('Frequency')
hly=xlabel('Group Size')
set(hlx,'FontSize',14);
set(hly,'FontSize',14);