%StochSensSim: this program performs simulations to estimate 
% stochastic growth rate and extinction risk sensitivity 
% analyses for vital rate means, variances, and correlations. 
% Notes in the Box 8.9 program explain the stochastic 
% simulation methods used and the way to create some of the
% files that must be pre-made. 
%uses the function betaset.m

clear all;

%*****************Simulation Parameters***********************
%the data shown here are for the desert tortoise:
% order of rates: survival2-7, growth2-6, fertility5-7
vrtypes= [ones(1,11),3,3,3]; % types of distributions to use (see VitalSim.m)
vrmeans= [0.77 0.95		0.79		0.93		0.91		0.87 ...
0.3		0.31		0.19		0.22		0.052	0.42	0.69	0.69]; 
vrvars=[0.0006		0.0006		0.053		0.0006		0.049 ...
0.02		0.0096 0.017		0.048		0.0006		0.0006	0	0	0];
vrmins=zeros(1,14);
vrmaxs=zeros(1,14);
load tortcorrmx; % a .mat file that contains a premade,
% corrected correlation matrix (corrmx) for the vital rates. 
% You must also have a file with a MATLAB function to create 
% the matrix elements from the vital rates.
% The program assumes that this function uses a vector 'vrs' to 
% build the matrix. For tortoises, this file is maktortmx.m
makemx = 'maketortmx'; % (Box 8.10 has a prodecure to make these)
% Also, you must have a set of premade beta values to load:
load tortbetas
Nstart = 500;  	% total initial population size;
Nx = 50; 			% quasi-extinction threshold. 
tmax = 75; 		% number of years to simulate;
np = 14;  		% number vital rates;
dims = 8; 		% dimensions of the population matrix;
runs = 500;  	% how many trajectories to do
change = 0.05; 	% the proportional magnitude of the 
					% perturbations used to estimate sensitivities
%*************************************************************
rand('state',sum(100*clock)); randn('state',sum(150*clock)); 
[W,D] = eig(corrmx); % making the M12 matrix 
M12 = W*(sqrt(abs(D)))*W';
vrs = vrmeans;
eval(makemx) 			% make a matrix of mean values
[uu,lam1] = eig(mx); 	% find the eigenvalues and vectors
lam1s = max(lam1);		
[lam0,iilam] = max(lam1s); 	% find dominant eigenvalue
uvec = uu(:,iilam)/sum(uu(:,iilam));% dominant right e-vector
n0 = Nstart*sum(uvec)*uvec; % make the initial pop'n vector
ExtResults = [];
Lresults = [];
crow = 1; ccol = 2; % used to pick correlation mx elements
pertnum=(2*np + 0.5*(np^2 -np)); 
% initialize results variables
Extsens = zeros(tmax,pertnum); ExtElast = zeros(tmax,pertnum);
Lsens = zeros(pertnum,1); Lelast = zeros(pertnum,1);
    
% loop one-by-one through perturbations of each mean, var, 
% correlation. For each type of rate, a small change is made,
% the appropriate set of parameters are redone, and then the 
% new model is stochastically simulated.
for pert = 0:(2*np + 0.5*(np^2 -np)); 
  % set temporary parameters to estimated values
  Tvrmeans =vrmeans; Tvrvars=vrvars; 
  Tparabetas=parabetas;  TM12 = M12;
  %now, modify values, depending on what is changed:
  if pert > 0 & pert < np+1 % change a mean rate
	Tvrmeans(pert) = (1+change)*vrmeans(pert);
	% to guard against too high a mean:
	if Tvrmeans(pert)*(1-Tvrmeans(pert))<= vrvars(pert) 
		Tvrmeans(pert) = (1-change)*vrmeans(pert); end;
	vrdiff = Tvrmeans(pert) - vrmeans(pert);
    % now, use the function betaset.m to find new beta or stretched betas
    % values for the changed parameter (if not lognormal):
	Tparabetas(:,pert) =betaset(vrtypes(pert),Tvrmeans(pert),...
	Tvrvars(pert),vrmins(pert),vrmaxs(pert));
  elseif pert > np & pert < 2*np+1 % change a variance
	Tvrvars(pert-np) = (1+change)*vrvars(pert-np);
	vrdiff = Tvrvars(pert-np) -vrvars(pert-np);
	Tparabetas(:,pert-np) =betaset(vrtypes(pert-np), ...
 	Tvrmeans(pert-np), Tvrvars(pert-np),vrmins(pert-np), ...
	vrmaxs(pert-np));
  elseif pert > 2*np % change a correlation
	Tcorrmx = corrmx;
	Tcorrmx(crow,ccol) = corrmx(crow,ccol)*(1+change);
	Tcorr(ccol,crow) = Tcorrmx(crow,ccol);
	vrdiff = Tcorrmx(crow,ccol) - corrmx(crow,ccol);
	[W,D] = eig(Tcorrmx);		
	TM12 = W*(sqrt(abs(D)))*W';
	% advancing counters thru indices for the upper 1/2
	% of the correlation matrix
	if ccol == crow+1;  		  % advancing counters thru
 		ccol = ccol+1; crow = 1;    % indices for upper 1/2
		else crow = crow+1;		  % of the correlation matrix
	end
end; % if pert 

% now do a set of random runs for the perturbation
PrExt = zeros(tmax,1); % extinction time tracker
logLam = [0]; 			% log-stochastic lambda tracker
for xx = 1:runs;
 if round(xx/10) == xx/10 disp([pert,xx]); end;
 nt = n0; extinct = 0;
 for tt = 1:tmax
	m = randn(1,np); rawelems = (TM12*(m'))';	
   for yy = 1:np  % find vital rates for this year
	 	if vrtypes(yy) ~= 3
			index = round(100*stnormfx(rawelems(yy)));
			if index == 0 index = 1; end;
			if index ==100 index = 99; end;
      	vrs(yy) = Tparabetas(index,yy);
		else 
       vrs(yy) = lnorms(Tvrmeans(yy),Tvrvars(yy),rawelems(yy));
		end;% if vrtypes(yy)
   end; %yy
   eval(makemx); %make a matrix with the new vrs values. 
   nt = mx*nt;	 %multiply by the population vector
   if extinct == 0	
		Ntot = sum(nt);
		if Ntot <= Nx 
				PrExt(tt) = PrExt(tt) +1;
				extinct = 1;
		end; % if Ntot
	end; %if extinct
 end %time (xx)
 logLam = logLam +(1/tmax)*log(sum(nt)/Nstart);
 end; % xx

if pert == 0 % establish the base-line results
	baseext= cumsum(PrExt./runs); 
	baselam = exp(logLam/runs);
	origbaselam = baselam;
	% smooth the extinction CDF function:
	baseext=([baseext]+[baseext(1);baseext(1:tmax-1)] ...
		 + [baseext(2:tmax);baseext(tmax)])/3;
else % summarize results for each perturbed model
	cumext=cumsum(PrExt./runs);
	origcumext = cumext;
	cumext=([cumext]+[cumext(1);cumext(1:tmax-1)] ...
		 + [cumext(2:tmax);cumext(tmax)])/3;
	ExtResults(1:tmax,pert)= cumext;
	Lresults(pert,1) = exp(logLam/runs);

   if vrdiff ~=0
	 % calculate extinction time sensitivities and elasticities:
	 ExtSens(1:tmax,pert)=(cumext - baseext)/vrdiff;
     % next line will generate 'divide by zero' errors: its OK
     ExtElast(1:tmax,pert)=((cumext - baseext)./baseext)/ ...
		(change*vrdiff/abs(vrdiff)); 
	 % calculate stoch-lambda sensitivities and elasticities:
	 Lsens(pert,1) =(Lresults(pert,1)-baselam)/vrdiff;
	 Lelast(pert,1)=((Lresults(pert,1)-baselam)/ ...
 		(baselam))/(change*vrdiff/abs(vrdiff));
   end; %vrdiff
	%smooth extinction elasticities:
	sExtElast(1:tmax,pert)=([ExtElast(1:tmax,pert)]+ ...
		[ExtElast(1,pert);ExtElast(1:tmax-1,pert)]...
	 	+ [ExtElast(2:tmax,pert);ExtElast(tmax,pert)])/3;
end; %if, else
end; % pert loop
% ExtSens and sExtElast hold extinction sensitivities and 
% elasticities for each year (rows) for each rate (columns, 
% going through the vital rate means, variances, and then
% correlations in order).
% Lsens and Lelast hold the stochastic growth rate
% sensitivities and elasticities in the same order.
% Given their size, it is best to pull out pieces of these
% results for display. For example:
disp('Below are the stochastic growth elasticities for')
disp('survival rates from classes 2 to 7')
disp([Lelast(1:7)])
disp('The figure shows extinction elasticities for')
disp('survival rates from years 25 to tmax');
figure
 plot([25:tmax], sExtElast(25:tmax,1:7));	 
 xlabel('Time');
 ylabel('Extinction Elasticities');