% simulate the Moran process

% 
clear all
set(0,                           ...
   'defaultaxesfontsize', 20,   ...
   'defaultaxeslinewidth', 1.0, ...
   'defaultlinelinewidth', 1.2, ...
   'defaultpatchlinewidth', 0.7);


for n = 1:6
   
N = 10*n;
nj(n) = N;
%the total number of cells must be an even number !!
M = 1000; % number of trials; 
X = N*ones(M,1)/2; % initial position is N/2
XS = X;% stores the positions at each time step
 
% Specify the movement vector
xm(1) = 1; %to the right
xm(2) = -1; %to the left
xm(3) = 0; % don't move 

nz = sum(X==0);
nN= sum(X==N);
nt = nz+nN;
j = 0;
while nt<M
    j = j+1;
    R = rand(M,1); % generate a random number for each trial
    for k = 1:M
        pij=X(k)*(N-X(k))/N^2;
        c = [pij,2*pij,1]; %this is the reaction vector to determine the actual move

       ndx= min(find(R(k)<=c)); % determine how to move
       X(k) = X(k) + xm(ndx); %update the position for particle k
    end
       XS(:,j+1) = X; % store the answer
     nz = sum(X==0);
nN= sum(X==N);
nt = nz+nN
  
end
 T=[0:j];
 figure(1)
plot(T,XS(1:10,:)','linewidth',2)
xlabel('time step','fontsize',20)
ylabel('position','fontsize',20)
axis([0 T(end) 0 N])

% now determine the time to fixation
for j = 1:M
    xend = X(j);
    xt(j) = min(find(XS(j,:)==xend));
end
hist(xt,N)
mn(n) =mean(xt);
end

alpha = sum(mn.*nj.^2)/sum(nj.^4)
loglog(nj,mn,'*',nj,log(2)*nj.^2,'--','linewidth',2 )
xlabel('N','fontsize',20)
ylabel('E(T)','fontsize',20)
  %compare alpha with log(2)
  log(2)