function []=Simulation_Lotka_Volterra(V0,P0,MaxT)
%
% Simulation_Lotka_Volterra(N, MaxTime)
%       where V and P are the initial population sizes
%       MaxTime is the duration fo the dynamics
%

%
%  We assume that that consumption and birth terms are given by the
%  geometeric mean of the two parameters A and C.
%
%

if nargin<3
    MaxT=500;
    if nargin<2
        V0=5000; P0=500;
    end
end
 
% set the parameters
TimeFactor=1/5000;

P=ones(ceil(P0),1)*[1 V0]*TimeFactor;  % A and D 
V=ones(ceil(V0),1)*[P0 1]*TimeFactor;  % B and C

%Use tau-leap updating. 
step=0.05/max(max([P; V]));

for T=1:MaxT
    
    for t=step:step:1
        
        birthsV=random('binomial',ones(size(V,1),1),V(:,1)*step);  % only allow one birth per step.
        deathsV=random('binomial',ones(size(V,1),1),sqrt(V(:,2)*mean(P(:,1)))*size(P,1)*step);
        
        birthsP=random('binomial',ones(size(P,1),1),sqrt(P(:,1)*mean(V(:,2)))*size(V,1)*step); % only allow one birth per step.
        deathsP=random('binomial',ones(size(P,1),1),P(:,2)*step);
        
        %fprintf(1,'%g P=%d+%g-%g V=%d+%g-%g\n',T+t,size(P,1),mean(sqrt(P(:,1)*mean(V(:,2)))*size(V,1)*step),mean(P(:,2)*step),size(V,1),mean(V(:,1)*step),mean(sqrt(V(:,2)*mean(P(:,1)))*size(P,1)*step));
        
        bV=V(birthsV>0,:); V=V(deathsV==0,:);
        bP=P(birthsP>0,:); P=P(deathsP==0,:);
        
        %decide what parameters to mutate and throw-away invalid ones.
        bP(:,1)=bP(:,1)+0.1*TimeFactor*randn(size(bP,1),1); bP=bP(bP(:,1)>0,:);
        %bP(:,2)=bP(:,2)+0.1*TimeFactor*V0*randn(size(bP,1),1); bP=bP(bP(:,2)>0,:);  %proportional to initial value
        %bV(:,1)=bV(:,1)+0.1*TimeFactor*P0*randn(size(bV,1),1); bV=bV(bV(:,1)>0,:);  %proportional to initial value
        %bV(:,2)=bV(:,2)+0.1*TimeFactor*randn(size(bV,1),1); bV=bV(bV(:,2)>0,:);
        
        V=[V; bV]; P=[P; bP];  % this is very inefficient in terms of memory allocation.
        
    end
    
    figure(1);
    subplot(2,2,1);
    plot(V(:,1),V(:,2),'.g');
    xlabel 'Parameter B'
    ylabel 'Parameter C'
    
    subplot(2,2,2);
    plot(P(:,1),P(:,2),'.r');
    xlabel 'Parameter A'
    ylabel 'Parameter D'
    
    subplot(8,1,5);
    hist(P(:,1),100); xlabel 'Parameter A';
    subplot(8,1,6);
    hist(V(:,1),100); xlabel 'Parameter B';
    subplot(8,1,7);
    hist(V(:,2),100); xlabel 'Parameter C';
    subplot(8,1,8);
    hist(P(:,2),100); xlabel 'Parameter D';
    
    figure(2);
    KeepNumbers(T,1:2)=[size(V,1) size(P,1)];
    KeepParameters(T,1:4)=[mean(P(:,1)) mean(V(:,1)) mean(V(:,2)) mean(P(:,2))];
    subplot(2,1,1);
    h=plot([1:T],KeepNumbers(1:T,1),'-g',[1:T],KeepNumbers(1:T,2),'-r');
    legend(h,'Prey','Predators');
    xlabel 'Time';
    subplot(4,2,5);
    plot([1:T],KeepParameters(1:T,1),'-k'); ylabel 'Parameter A';
    subplot(4,2,6);
    plot([1:T],KeepParameters(1:T,2),'-k'); ylabel 'Parameter B';
    subplot(4,2,7);
    plot([1:T],KeepParameters(1:T,3),'-k'); ylabel 'Parameter C';
    subplot(4,2,8);
    plot([1:T],KeepParameters(1:T,4),'-k'); ylabel 'Parameter D';
    drawnow;
    
    if size(V,1)==0 || size(P,1)==0
        error('One of the population sizes hits zero');
    end
end
end