% Generalized, discrete Lotka-Voltera model
% Written for CS 477; Computer Simulation and Modeling
% 3/1/08 JVJ

% Prepare the environment
clear;
clf;


%=========================
% Model parameters
N = 100;                   % Number of species
Nat = 15;                  % Number of autotrophs
conversion = .1;           % Fraction of prey that is usable.
deathRate = .01;           % Fraction of heterotrophs dying each time
r = 2.5;                   % Basic reproductive number for autotrophs.
interactionCutoff = 1e-4;  % Weakest interaction
extenctionDensity = 1e-5;  % Populations less than this are set to zero.
evolutionStrength = 0.8;   % How fast the interaction matrix 'evolves'
timeSteps = 150;           % Number of time steps
%=========================

% Create interaction matrix M
M = -triu(rand(N,N),1); % Produce upper triangular, negatives (losses)

% Randomly asign preditors (positive) and prey (negative)
tmp = rand(size(M(Nat+1:end,:)));
tmp1 = M(Nat+1:end,:);
tmp1(tmp>.5)=tmp1(tmp>.5)*-1;
M(Nat+1:end,:) = tmp1;

% Complete the matrix with negative symmetric values.
M = M - M';

% Only set conversion factor of the things eaten
M(M>0)= M(M>0)*conversion;

% Create the initial population (this is important, but 'tuned')
x = rand(N,1)/N;
x(1:Nat)=N*x(1:Nat);

% Arrays to store data for plotting
plotArray=[]; % This array stores populations
ev=[];        % This array stores eigen values

for i = 1:timeSteps
    % This turns off very weak interactions
    Meff=M;
    tmp = Meff./max(Meff(:));
    Meff(abs(tmp)<interactionCutoff) = 0;

    % Interaction vector. This is sum over j  of M_{ij} x_i x_j
    interaction = x.*(M*x);
    
    % Handle updates to each population seperately. 
    % Autotrophs
    dNatdt = r * x(1:Nat) .* (1 - x(1:Nat)) + interaction(1:Nat);
    x(1:Nat) =  x(1:Nat) + dNatdt;

    % Heterotrophs
    dNhdt =  - deathRate * x(Nat+1:end) + interaction(Nat+1:end);
    x(Nat+1:end) = x(Nat + 1:end) + dNhdt;
    
    % All rates of change, for the evolution of links calculation
    dxdt = [dNatdt;dNhdt];
    
    % Can't be deader than extinct, no negative populations.
    x(x<1e-5)=0;
    
    % Evolution of links
    M = M + evolutionStrength * (repmat(dxdt',N,1) - repmat(dxdt,1,N));
    
    % Compute and store the set of eigen values
    % Note that this form of the matirx is NOT the whole story!!!
    ev=[ev eig(M.*repmat(x,1,N))];
    
    % Store data for plots
    plotArray = [plotArray x];
end

% A summary plot

% Populations on a semilog scale
subplot(2,3,1);
semilogy(plotArray');
title('Populations vs. Time')
xlabel('Time (steps)');
ylabel('logl0 Population');

% log log plot of interaction strengths
subplot(2,3,2);
[b c]=hist(M(M>0),30);
loglog(c,b,'ko');
title('loglog of interaction strength (matrix M entries)');
xlabel('log interaction strength');
ylabel('log number of species');

% color plot of the interaction Matrix
subplot(2,3,3);
imagesc(M);
colorbar;
title('Interaction Matrix Values');

% Total Number of species
subplot(2,3,4);
plot(mean(plotArray,2),'k');
title('Population of all Species');
xlabel('Species Number');
ylabel('Population');

% The eigen values
subplot(2,3,5);
pev=zeros(size(ev));
pev(real(ev)<0) =1; % Stable points will be 1 (red)
imagesc(pev);
title({'Eigen values','Red where real part of eigen values is negative (stable).',...
    'Blue where real part of eigen values is positive (un-stable).'});

subplot(2,3,6);
title('Final Relative Species Abundance Curve');
s=plotArray(:,end);
s=s(s>0);
s=s*1e10;
[Nrsa xrsa]=hist(log10(s),10);
plot(xrsa,Nrsa,'ko-');
xlabel('Abundance');
ylabel('Count');
fprintf('The number of species remaining is: %i of %i starting.\n',sum(x>0),N);