% Vitalsens.m  Takes vital rate means and a matrix      
% definition to calculate deterministic sensitivities and 
% elasticities of lambda to vital rates, using symbolic
% functions to take derivatives. This example uses emperor
% goose demography: fall to fall, with the youngest class as 
% fledglings, using Schmutz et al. 1997 survey set 1
% The program calls the function eigenall.m (Box 7.1)

%IMPORTANT NOTE: you can't run this program without having 
% the Symbolic toolbox for Matlab
% this program was changed on 4/33/04 to conform to newer MATLAB conventions
clear all;
%**************** SIMULATION PARAMETERS **********************
vr = [0.1357    0.8926    0.6388    0.8943];  % vital rates
syms  Ss0 Ss1 Sf2 Sf3  % vital rates as symbolic variables
Svr = [Ss0 Ss1 Sf2 Sf3]; % vector of symbolic vital rates 

% Next, a symbolic definition of the matrix
mx = 	[	0, 		0,		Sf2*Ss1, 	Sf3*Ss1;
       	Ss0, 	0, 		0,			0;
          0, 	Ss1, 	0,			0;
          0,		0,		Ss1,		Ss1];       
%*************************************************************

% make a matrix of the mean numerical values using subs 
realmx = subs(mx,Svr,vr);  
% use eigenall.m to get eigenvalues
[lambdas,lambda1,W,w,V,v]= eigenall(realmx); 
sensmx = v*w'/(v'*w);   % sensitivities of matrix elements
elastmx = (sensmx.*realmx)/lambda1; % element elasticities
numvrs = length(vr); % how many vital rates?
vrsens = zeros(1,numvrs); % initialize vital rate sens.

% a loop to calculate sensitivity for each vital rate
for xx=1:numvrs 
	% derivatives of elements with respect to vital rates
	diffofvr = double(subs(diff(mx,Svr(xx)),Svr,vr)); % this was changed: 4/22/04
	% sum up to get vital rate sensitivities
    
	vrsens(xx) = double(sum(sum(sensmx.*diffofvr))); % this was changed: 2/22/2003
end; % xx

vrelast = ((vrsens.*vr)/lambda1); % calculate elasticities
disp('Matrix element sensitivities and elasticities:')
disp(sensmx); disp(elastmx);
disp('Below are the vital rate results');
disp('vital rates:');
disp(Svr)
disp('sensitivities');
disp(vrsens)
disp('elasticities');
disp(vrelast)