clear all;

% Supplement to Chapter 8 in book on Ranaviruses (Gray and Chinchar, eds.) 
% published in 2014 by Springer.

% Code is written by Julia E. Earl.  Feel free to contact her with questions- 
% https://sites.google.com/site/juliaeearl/

% This model is a stage structured population matrix model for exploring 
% the effects of Ranavirus exposure in year one on population dynamics.  
% The model requires input of demographic parameters (survival, fecundity)
% for uninfected populations, survival given ranavirus exposure, adult
% carrying capacity, and average frequency of ranavirus exposure.

% The current form of the model is based on wood frogs (Lithobates 
% sylvaticus), but a different stage structure could be implemented by
% altering the population matrix below.

% This is a female-only, post-breeding model, meaning that the numbers of
% individuals in each class for a year are that just after eggs are laid.

% Extinction is defined very conservatively as zero living individuals.


% POPULATION PARAMETERS
runs=1000; % each run is a simulation over a certain number of years.
years=1000; % number of years for each run.
a=4; % number of age classes.
xi=[300;300;50;50]; %initial age distribution vector of individuals in 
%population representing the number of individuals in each age class


%the mean survival probability for each age class.
p=[0.035 0.3775 0.249 0.12]; 
% p[1]= probability that an individual survives from egg to year one.
% p[2]= probability than an individual survives from year one to year two.
% p[3]= probability than an individual survives from year two to year three.
% P[4]= probability that an individual in the three year class continues to
% survive.

% standard deviation for survival probability for each age class.  If you want no
% variability in survival, set psd values to zero.
psd=[0.0245 0.102 0.1575 0.05]; 
f=[0 0.215 334.5 338.3]; % mean fecundity for each age class
%note that fecundities should be adjusted for female survival
fsd=[0 0.0548 31.94 31.62]; % sd for fecundity for each age class
% If you want no variability in fecundity, set fsd values to zero.
K=500; % Adult carrying capacity


% DISEASE PARAMETERS
% Set disprob to zero or dissurv to one to run the model with no ranavirus
%exposure.
disprob=0.1;  % This is the probability that ranavirus exposure occurs in a 
% given year.  
dissurv=0.17;  % This is the proportion that survive if exposed to 
% ranavirus assuming no mortality from any other source, such that the 
% model would run the same as no ranavirus exposure if dissurv is one.


%preallocating for speed
    N=zeros(1,years);
    Extinct=zeros(runs,1);
    finalAdult=zeros(runs,1);
    ExtinctYear=zeros(runs,1);
    finalAdult=zeros(runs,1);
 
    

% RUNNING SIMULATIONS
for i=1:runs % Number of times to run whole model

Xt=zeros(a,years); %Reserving memory where a is #age classes and 
% years is number of time steps
% X is the number of individuals in each age class at each time step
N=zeros(1,years); % The total number of individuals in all age classes 
% combined in each year.  Preallocating for space.


% Placing values for starting conditions at years=1
Xt(:,1)=xi; %placing the initial age distribution in the the first column 
%of the matrix X (X(:,1).... means "every row, first column")
N(1,1)=sum(xi); % placing the initial population size in the N vector.

  
%%Begin run of year 2 through final year%%
for t=2:years 
   
    % Drawing values of demographic parameters from distributions specified
    % above for each year of the model.  All are from a normal
    % distribution.
    P = abs(normrnd(p(1),psd(1)));
    F1 = abs(normrnd(f(1),fsd(1)));
    F2 = abs(normrnd(f(2),fsd(2)));
    F3 = abs(normrnd(f(3),fsd(3)));
    F4 = abs(normrnd(f(4),fsd(4)));
    P2 = abs(normrnd(p(2),psd(2)));
    P3 = abs(normrnd(p(3),psd(3)));
    P4 = abs(normrnd(p(4),psd(4)));
  
 
    
%  POPULATION MATRIX

if binornd(1,disprob)==1 % Determines whether it is a year with ranavirus exposure.
    LX=[F1 F2 F3 F4; 
    P*dissurv 0 0 0; 
    0 P2 0 0; 
    0 0 P3 P4]; % Matrix with ranavirus exposure
else LX=[F1 F2 F3 F4; 
    P 0 0 0; 
    0 P2 0 0; 
    0 0 P3 P4]; % Matrix without ranavirus exposure
end

%Putting your population Xt, through the Matrix for one time step.
Xt(:,t)=floor((LX)*Xt(:,t-1)); % floor means that there can't be fractions 
%of individuals


if Xt(3,t)+Xt(4,t)>K  
% Checks whether adult population size is greater than the carrying capacity.
    Xt(3,t)=floor(K/2); 
    Xt(4,t)=floor(K/2);
    % Sets adult population size equal to carrying capacity equally split
    % between two adult stage classes.
end

N(t)=sum(Xt(:,t));  % Calculates the total number of individuals each year

end

finalX(i,:) = Xt(:,years); % The number of individuals in each age class 
%in the final year.

finalN(i,:) = N(years);  % Total population size in the final year.

if finalN(i) < 1  % Determines whether the population is extinct.
    Extinct(i,1) = 1;  % Extinct = 0 is not extinct, Extinct = 1 is extinct.
end

if Extinct(i,1) == 1
    idxs = find(N==min(min(N)));
    ExtinctYear(i,1)=idxs(1,1);  % calculates the year of extinction
else ExtinctYear(i,1)=NaN;  % if not extinct the year is recorded as a missing value.
end

if Extinct(i,1)==1
     finalAdult(i,1)=NaN;
else finalAdult(i,1)=finalX(i,3)+finalX(i,4);  
    % calculates the number of adults in the final time step
 end

end


% PRINT OUTPUT IN MATLAB DISPLAY
Extpops=sum(Extinct);  % Number of model runs resulting in extinction
disp('Extinction probability')
ExtProb=sum(Extinct)/runs;
disp(ExtProb)
ExtTime=nanmean(ExtinctYear);
disp('Average time to extinction')
disp(ExtTime)
stderrExtTime=nanstd(ExtinctYear)/sqrt(sum(Extinct));
% Standard error of the time to extinction

disp('Average final population size')
avFinN=mean(finalN);  
disp(avFinN)
stderrFinN=std(finalN)/sqrt(runs);
%Standard error of the final population size across model runs
disp('Average final adult population size')
avFinAdult=nanmean(finalAdult); %Average of the final adult population size across model runs
disp(avFinAdult)
stderrAdult=nanstd(finalAdult)/sqrt(runs-Extpops);
%Standard error of the final population size across model runs


    
% OUTPUT TO CSV FILE
Output=zeros(runs,4);%creating matrix for output.  
%The number next to runs is equal to the total number of variables that you are tracking.

Output(:,1)=ExtinctYear;
Output(:,2)=finalN;
Output(:,3)=Extinct;
Output(:,4)=finalAdult;

Output;

csvwrite('RASYout.csv',Output); % Creates csv file.