% Limitsens.m  Calculates sensitivities over a range of values 
% for each vital rate, accounting either for variation or 
% uncertainty in rate values.  
% This example assumes a uniform distribution of values for
% each rate between their high and low bounds

%IMPORTANT NOTE: you can't run this program without having 
% the Symbolic toolbox for Matlab

clear all;

%*****************Simulation Parameters***********************
% Emperor goose data from Schumtz et al. 1997
meanvr = [0.1357    0.8926    0.6388  0.8943]; % best estimates 
vrhi = [0.38 0.90	0.64	1.16]; % highest estimates 
vrlo = [0.018 0.50 0.57 0.80];   % lowest estimates
syms  Ss0 Ss1 Sf2 Sf3 % symbolic variable definitions
Svr = [Ss0 Ss1 Sf2 Sf3]; % vector of symbolic variables
% symbolic definition of matrix:
mx = 	[	0, 		0,		Sf2*Ss1, 	Sf3*Ss1;
       	Ss0, 	0, 		0,			0;
	       0, 	Ss1, 	0,			0;
          0,		0,		Ss1,		Ss1];       
reps = 500; %number of replicate matrices to make and analyze
%************************************************************

rand('state',sum(100*clock));	% random seed
numvrs = length(meanvr); 		% how many vital rates?

% make uniform randoms and convert to random vital rates
allvrs = rand(reps,numvrs);
allvrs =allvrs.*repmat(vrhi-vrlo,reps,1) + repmat(vrlo,reps,1); 

for rr = 1:reps  % do calculations for each random matrix
 disp(rr); 		 % display progress
 realmx = subs(mx,Svr,allvrs(rr,:)); % Make a random matrix
 [lambdas,lambda1,W,w,V,v]= eigenall(realmx); 
 sensmx = v*w'/(v'*w);  	% sensitivities of elements
 elastmx = (sensmx.*realmx)/lambda1;  % element elasticities    
 alllams(rr,1) = lambda1; 			    % save lambda value
 for xx=1:numvrs % loop to make sensitivities for vital rates
		diffofvr = double(subs(diff(mx,Svr(xx)),Svr,allvrs(rr,:)));%changed 4/22/04	
		vrsens(xx) = sum(sum(sensmx.*diffofvr));
 end; % xx 
 allelasts(rr,:) = ((vrsens.*allvrs(rr,:))/lambda1); 
end; % rr

% now get sensitivities and elasticities for the best estimates
realmx = subs(mx,Svr,meanvr); 
[lambdas,lambda1,W,w,V,v]= eigenall(realmx); 
meanlam1 = lambda1;  	sensmx = v*w'/(v'*w);  
elastmx = (sensmx.*realmx)/lambda1; 
meansens = zeros(1,numvrs);
for xx = 1:numvrs
	diffofvr = double(subs(diff(mx,Svr(xx)),Svr,meanvr)); %changed 4/22/04
	meansens(xx) = sum(sum(sensmx.*diffofvr));
end; 
meanelast = ((meansens.*meanvr)/lambda1);	

% ask in turn how much the maximum value for each rate would 
% change lambda, if all other rates are at best estimate
maxlams=zeros(1,numvrs);
for rate = 1:numvrs  
	vrates = meanvr; 	vrates(rate) = vrhi(rate);
	realmx = subs(mx,Svr,vrates);
	[lambdas,lambda1,W,w,V,v]= eigenall(realmx); 
	maxlams(rate) = lambda1;  % save resulting lambdas 
end ;

% Now display the results
disp('Below are the maximum lambdas and max proportional')
disp('change in lambdas from changing each vital rate')
disp(Svr);	disp(maxlams)
disp((maxlams-ones(1,numvrs)*meanlam1)/meanlam1)


disp('Below are the r^2 values for lambda and each')
disp('vital rate, a measure of influence on population')
disp('growth for the random simulated matrices')
correls = corrcoef([allvrs,alllams]);
disp(Svr); disp((correls(numvrs+1,1:numvrs)).^2);

disp('Below are the elasticities for the mean vital rates,')
disp('and then the min, max, mean, and st. deviation of the')
disp('elasticity values from the random matrices')
disp(Svr); disp(meanelast);
disp(min(allelasts)); disp(max(allelasts))
disp(mean(allelasts)); disp(std(allelasts))

disp('Next, upper and lower 95% confidence limits for the')
disp('elasticities calculated from the random matrices')
disp(Svr);  CLup=zeros(1,numvrs); CLlo=zeros(1,numvrs);
for vr=1:numvrs
	x=sort(allelasts(:,vr));
	CLup(vr)=x(1+round((reps-1)*0.975));
	CLlo(vr)=x(1+round((reps-1)*0.025));
end;
disp(CLup); disp(CLlo);

disp('Finally, squared correlations between the elasticity')
disp('of each rate and the values of other vital rates,') 
disp('a measure of the influence of each vital rate on the') 
disp('elasticity results')
for rate = 1:numvrs
	disp('For elasticity of rate:');
	disp(Svr(rate));
	correls = corrcoef([allvrs,allelasts(:,rate)]);
	disp(Svr);
	disp((correls(numvrs+1,1:numvrs)).^2);
end;