%VitalSim: 
% This is a program to do a general stochastic PVA simulation
% including within-year, auto-, and cross-correlations
% and a choice of different distributions for the vital rates.
clear all;

%*****************Simulation Parameters***********************
% the data shown here are for H. montana (best Kendall estimates):
% the parameters are defined in a set of arrays
% with elements corresponding to these vital rates or matrix elements:
%1,	 	2,		3,	 	4,		5,		6		7		8		9		10	11	12	13	
%g>4,3 g>4,4	g>5,4	g>6,4	g>4,5	g>5,5	g>6,5	g>5,6	g>6,6	S3	S4	S5	S6
% and these:
%14	 15	  16  17	18	 19  20  21	 22	 23   24
%a11 a13  a14 a15   a16  a21 a23 a24 a25 a26 a32

% vrtypes identifies the distribution for each rate: 
% 1 = beta, 2 = stretched beta, 3 = lognormal
vrtypes= [ones(1,13),3,3,3,3,3,1,1,1,1,1,1];
% now means and variances of these rates or elements
vrmeans= [0.24	0.79	0.30	0.024	0.92	0.81	0.25	0.93 ...
	0.75	0.78798799	0.96452906	0.99519038	0.98076152...
	0.4995	4.5782	12.1425	22.3167	50.1895	0.0004	0.0033 ...
 	0.0088	0.0162	0.0364 0.4773];
vrvars=[0.00018836	0.00016164	0.00026851	0.00002805 ...
	0.00005477 	0.00184484	0.00937915	0.00002805 ...
	0.01333333	0.01790190	0.00000963	0.00000028 ...
 0.00000322 0	0.012	0.0843	0.2848	1.4403	0	0	0	0	0 0.0003];
% minimum and maximum values for each vital rate: 
% put in zeros for rates that are not stretched betas 
vrmins=zeros(1,24); vrmaxs=zeros(1,24);
np = length(vrmeans); % how many parameters are there? 
% you must load a file with a pair of correlation matrices: 
% corrin is the within yr correlation part of the matrix, 
% corrout is the between yr correlations. See Box 8.9 for a 
% description of the right format 
load hudcorrs
% NOTE: In this example, we use correlations between the first 13
% parameters, but not the other 11 reproductive parameters, which
% vary little. 
yrspan = 20; % the number of years of correlations to
% build into the M12 matrix

% you must also have a separate m-file that uses a vector of vital 
% rates (vrs) of the same form as vrmeans, above, to define the 
% elements of the populaton matrix (mx) (see Table 8.1 for example)
makemx='hudmxdef'; % hudmxdef.m defines the elements of mx

n0=[4264; 3; 30; 16; 25; 5];% initial population vector
Nx = 500; 	% quasi-extinction threshold
tmax = 50; 	% number of years to simulate;
np  = 13;  	% number vital rates in correlation matrix;
np2 = 11;	% number of uncorrelated rates. 
dims = 6; 	% dimensions of the population matrix;
runs = 1000; % how many trajectories to do
%*************************************************************

randn('state',sum(100*clock)); % seeds random numbers
Nstart = sum(n0);					 % starting population number
vrs=vrmeans;	% set vital rates to their mean values
eval(makemx);	% use matrix definition file to make mean matrix
lam0=max(eig(mx)); %find the deterministic population growth rate

%------------------------------------------------------------
% this section makes sets of beta or str. beta values to choose
% from during the simulations; it makes 99 values for 1% 
% increments of Fx for each parameter -- if you already have made 
% a set of beta values for this life history, you can recall
%  these, and save the time of recalculating them. 
yesno = ... 
input('type 0 to calculate betas, or 1 to get a stored set');
if yesno == 1
  betafile = ...
  input('type filename with stored betas; put in single quotes');
   load(betafile)% this line is corrected from that in Morris and Doak
else    % make a set of values for each beta or stretched beta
  parabetas=zeros(99,np+np2);
  for ii = 1:(np+np2)
	if vrtypes(ii) ~= 3
 	 for fx99 = 1:99
  		if vrtypes(ii) ==1; 
			parabetas(fx99,ii) = ...  
			 betaval(vrmeans(ii),sqrt(vrvars(ii)),fx99/100); end;
		if vrtypes(ii) ==2; 
			parabetas(fx99,ii) = ...
			 stretchbetaval(vrmeans(ii),sqrt(vrvars(ii)),...
			 vrmins(ii), vrmaxs(ii), fx99/100); end;
 	 end; % fx99
	end; % if vrtypes(ii)
  end; % ii loop
  yesno = input('type 1 to store the betas, or 0 if not');
  if yesno ==1 
    betafile = ...
    input('type filename to store betas, put in single quotes');
	save(betafile, 'parabetas'); % this line is corrected from that in Morris and Doak
  end; %if yesno 
 end; %else

% the next section is identical to that in the program 
% BetweenYrCorr in Box 8.9 - simply substitute in the lines 
% from that program that are delineated by same comment 
% lines as the two that follow here:
%--------creating and using the big correlation matrix, M--------
%--------creating and using the big correlation matrix, M----
% this set of loops makes the big correlation matrix (M)
% with multi-year correlations: the  if statements are used to 
% estimate the correct correlations with increasing time lags, 
% always assuming that all long-time-lag correlations are only 
% caused by within-year and one-time-step correlations
for ii = 1:yrspan;   
 for jj=1:yrspan;
  		if ii==jj, litmx = corrin;
  			else litmx =corrout; end;
		expo=1;
		if ii > jj, expo = ii-jj; litmx = (litmx)^expo; end;
		if ii < jj expo = jj-ii;  litmx = (litmx')^expo; end;
 		for ip=1:np;
 		  for jp = 1:np;
		  M(ip+np*(ii-1),jp+np*(jj-1)) = litmx(ip,jp);
		  end; % jp
		end; % ip
 end; % jj
end; % ii

% get the eigenvalues for calculating the M12 matrix 
% (here called zfull):
[W,D] = eig(M); % getting the eigenvalues and vectors

% now, a check for negative eigenvalues -- if you have them,
% it sets negatives and small positive values = 0 
check = min(min(D));	% are the smallest eigenvalues negative?
if check < 0
	disp('The min eigenvalue is < 0. Eigenvalues are:')
	disp(diag(D))
	disp('hit enter to continue with approximation')
	pause
	maxneg = max(max(abs(D(find(D<=0))))); 
	% maxneg is the largest negative eigenvalue
	D(find(abs(D<=maxneg))) = 0; % set the negatives = 0
	disp('Corrected eigenvalues are:')
	disp(diag(D))
	newfullmx = W*newD*W'; % make a corrected matrix
	for ii=1:np % change from covariances to correlations  
	 for jj=1:np
	if newfullmx(ii,ii) ==0 | newfullmx(jj,jj)==0; 
			newfullmx(ii,jj) =0;
		else
				newfullmx(ii,jj) = ...
newfullmx(ii,jj)/((newfullmx(ii,ii)*newfullmx(jj,jj))^0.5);
		end;
	 end;
	end;
	[W,D] = eig(newfullmx);  
end; %check < 0

M12 = W*abs(D.^0.5)*W';	% the M^(1/2) matrix  
sz = length(M12); % the total number of lines in M12
	
% get the lines from the middle of M12 
% to use to generate correlations
startcase = (round(yrspan/2)*np +1);  			
zvalold=real(M12(startcase:(startcase+np-1),:));
zvalnew=real(M12((startcase+np):(startcase+2*np-1),:));
clear M12 W D; % clearing memory  
% zvalold and zvalnew are each one year of rows in M12 
% to start the whole thing off, calculate a first set of
% normals, then multiply with zvalold to get correlated normals:
newns = randn(sz,1);	
oldxy = zvalold*newns;
%-----end of: creating and using the big correlation matrix------

%-----end of: creating and using the big correlation matrix------

% finally, do sets of runs to get growth rate and extinction risk
results = []; normresults = [];
PrExt = zeros(tmax,1);    % the extinction time tracker
logLam = zeros(runs,1);   % the tracker of log-lambda values
stochLam = zeros(runs,1); % tracker of stochastic lambda values
for xx = 1:runs;
 if round(xx/10) == xx/10 disp(xx); end; % displays progress
 nt = n0; % start at initial population vector
 extinct = 0;
 for tt = 1:tmax
	% make random normals
	newns = [newns((np+1):sz,1); randn(np,1)]; 
	newxy= zvalnew*newns;	% make new set of corr'ed normals
	normresults = [normresults; oldxy', newxy']; % these lines 
	oldxy = newxy;	      % save normals to check correlations 

	yrxy=[newxy;randn(np2,1)]; 
				%adds in randoms for uncorrelated vital rates. 
    for yy = 1:(np+np2)  % loop finds vital rate values
	 	if vrtypes(yy) ~= 3 % if not a lognormal rate
			index = round(100*stnormfx(yrxy(yy)));
			if index == 0 index = 1; end; % round at extremes
			if index ==100 index = 99; end;
      		vrs(yy) = parabetas(index,yy); % find stored value
		% else, calculated a lognormal value:
		else vrs(yy) = lnorms(vrmeans(yy),vrvars(yy),yrxy(yy));
		end;% if vrtypes(yy) ~= 3
    end; % yy loop
  	eval(makemx); 	% make a matrix with the new vrs values. 
  	nt = mx*nt;	 	% multiply by the population vector

	if extinct == 0	% check for extinction 
		Ntot = sum(nt);
		if Ntot <= Nx 
				PrExt(tt) = PrExt(tt) +1;
				extinct = 1;
		end; % if Ntot
	end; % if extinct

 end % time (tt) loop
logLam(xx) = (1/tmax)*log(sum(nt)/Nstart); % calculate loglambda 
stochLam(xx) = (sum(nt)/Nstart)^(1/tmax);	% and stoch. lambda
end; % runs (xx) loop

CDFExt = cumsum(PrExt./runs); % make the extinction CDF function
disp('This is the deterministic lambda value'); disp(lam0);
disp('And this is the mean stochastic lambda'); disp(exp(mean(logLam)));
disp('Below are mean and standard deviation of log lambda');
disp(mean(logLam)); disp(std(logLam));
disp('Next is a histogram of logLams'); Hist(logLam);
disp('And now, the extintion time CDF'); figure; plot(CDFExt);