% mat3prosretr_08 
%   May, 2008
%     uses function file  "senstuff.m"
%     uses mat-file "Caretta.mat"
%
%Demographic lessons
%  prospective analysis   (_sensu_ Horvitz, Schmeske and Caswell 1997)
%      sensitiviy
%      elasticity
%  retrospective analyses (_sensu_ Horvitz, Schmeske and Caswell 1997)
%      variance   by sensitivity-squared
%      covariance by sensitivity-product
%      LTRE contributions analysis: deviation by sensitivity

%MATLAB lessons
%  refering to matrix elements
%  naming matrices in a numerical series
%  max function  
%  num2str 
%    used in a do-loop with a numbered series
%  subplot   for puttting more than one figure on a page
%  matrixname(:)  for making a single vector from a matrix
%  means, standard deviation and variances
%  rearranging the elements of a vector
%  element by element multiplication

%this example uses  Human demography (Caswell 2001, Example 4.3, p.78)
% for some analyses 

%start by clearing the workspace
clear all
close all

% Human population, USA 1966, 
%5 yr projection interval, Keyfitz and Flieger 1971
Fi=[0
   0.00102
   0.08515
   0.30574
   0.40002
   0.28061
   0.15260
   0.06420
   0.01483
   0.00089]

Pi=[0.99670
   0.99837
   0.99780
   0.99672
   0.99607
   0.99472
   0.99240
   0.98867
   0.98274]

humans=zeros(10)%set up a 10 by 10 matrix filled with zeros
humans(1,:)=Fi';  %change the first row to fertilities

%write a do loop to construct the matrix
for i=1:9
   humans(i+1,i)=Pi(i);
end

humans

save humans humans    
%saves the variable 'humans' in the file 'humans.mat'    

%use senstuff get lambda, etc
% prospective analyses
% aymptotic properties
% sensitivity and elasticity matrices for the humans example

[lhumans ssdhumans rvhumans senhumans elahumans]=senstuff(humans);

nstages=size(humans,1);

figure(1)
subplot(2,1,1)
bar(ssdhumans)
ylabel('Proportion of population')
title('Stable age-class distribution')
subplot(2,1,2)
plot(rvhumans)
ylabel('Relative reproductive value')
title('Reproductive value vector')
xlabel('Age interval (for humans each is 5 yrs long)')


%pick out the sensitivities for the survival parameters 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  CAUTION: THE NEXT SECTION REFERS  TO an age-structured matrix only
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
P_structure=zeros(nstages,nstages); %set up matrix
for i=1:nstages-1;
    P_structure(i+1,i)=1
end
P_sen=P_structure.*senhumans;%get sensitivities for survival only
P_ela=P_structure.*elahumans;%get elasticities 
P_sen=max(P_sen); %there is only one per column that is nonzero
P_ela=max(P_ela); % similarly for elasticities

figure(2)
subplot(2,1,1) % make two figures, this is the first
semilogy(P_sen,':v','LineWidth',2);
hold on
semilogy(senhumans(1,:),'-o','LineWidth',2);% the F sensitivities
ylabel('Sensitivity of \lambda')
legend('P_i','F_i')
title('Sensitivity and elasticity of \lambda for humans')
hold off
subplot(2,1,2) % make two figures, this is the second
plot(P_ela,':v','LineWidth',2);
hold on
plot(elahumans(1,:),'-o','LineWidth',2);% the F elasticities
ylabel('Elasticity of \lambda')
xlabel('Age interval (for humans each is 5 yrs long)')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  END AGE-STRUCTURED EXAMPLE FOR HUMANS ONLY
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  THE FOLLOWING IS MORE GENERAL AND WILL HELP VISUALIZE 
%  SENSITIVITIES AND ELASTICITIES FOR ANY MATRIX
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(3)
subplot(1,2,1)
bar3(senhumans)
xlabel('Stage at time t')
ylabel('Stage at time t+1')
zlabel('Sensitivity of \lambda')
subplot(1,2,2)
bar3(elahumans)
xlabel('Stage at time t')
ylabel('Stage at time t+1')
zlabel('Elasticity of \lambda')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

elamat=elahumans;
%sum the elasticities in various ways
elafec=sum(elamat(1,:)); %row1 sums
elastage=sum(elamat); %column sums
elastasis=sum(diag(elamat));


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Sensitivity analysis for Caretta caretta
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all 

load Caretta
[lam_turtle ssd_turtle rv_turtle sen_turtle ela_turtle] = senstuff(A);

rv_turtle=rv_turtle/rv_turtle(1)%scales as in Crowder et al.

figure(4)
subplot(2,1,1)
bar(ssd_turtle)
ylabel('Proportion of population')
title('Stable age-class distribution for Caretta')
subplot(2,1,2)
plot(rv_turtle)
ylabel('Relative reproductive value')
title('Reproductive value vector for Caretta')
xlabel('Stage')



figure(5)
subplot(1,2,1)
bar3(sen_turtle)
xlabel('Stage at time t')
ylabel('Stage at time t+1')
zlabel('Sensitivity of \lambda')
subplot(1,2,2)
bar3(ela_turtle)
xlabel('Stage at time t')
ylabel('Stage at time t+1')
zlabel('Elasticity of \lambda') 

figure(6)
plot(diag(ela_turtle),'-rs')
hold on
plot(ela_turtle(1,:),'-gv')
plot(diag(ela_turtle,-1),'-b^')
legend('P','F','G')
ylabel('Elasticity')
xlabel('Stage')
title('Caretta (Figure 1, Crowder et al. 1994)')

%sum the elasticities in various ways
elafec_turtle=sum(ela_turtle(1,:)); %row1 sums
elastage_turtle=sum(ela_turtle); %column sums
elastasis_turtle=sum(diag(ela_turtle));




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Retrospective analyses

%imagine that there are three different habitats for the humans
% and that we have studied the demography in each
% just for the sake of saving time in writing these
%  we start with equal ones to later make alterations

load humans

humansenv1 = humans;
humansenv2 = humans;
humansenv3 = humans;

%change the top row and first column, ie reproduction and survival of young
humansenv2(1,:) = humans(1,:)* 10;
humansenv2(:,1) = humans(:,1)* .5

humansenv3(1,:) = humans(1,:)* 100;
humansenv3(:,1) = humans(:,1)* 1.5

%save them for tomorrow in the matfile 'humansenvs.mat'
save humans_envs humansenv1 humansenv2 humansenv3

%how do they differ?  how do their lambdas differ?
% how do their differences contribute to their different lanbdas

         
lhumansenv=[0 0 0]%these will be reset inside the loop

 for i = 1:3
i2 = num2str(i);
mat = eval(['humansenv' i2]);
lhumansenv(i)=senstuff(mat);
end

%the lambdas of all three environments
figure(7)
bar(lhumansenv)
xlabel('environments')
ylabel('\lambda')




%make three columns for all the entries to get means and variances,etc
humansall=[humansenv1(:) humansenv2(:) humansenv3(:)];
%take the transpose 
humansall=humansall';
%then take the mean and std
humansbar=mean(humansall);
humansstd=std(humansall);


%reassemble as a matrix the humansbar
humansbar=humansbar';
humansbar= reshape(humansbar,10,10)
%reassemble as a matrix the humansstd
humansstd=humansstd';
humansstd= reshape(humansstd,10,10)

%find the lambda and sensitivity of the mean matrix
[lhumansbar ssdhumansbar rvhumansbar senhumansbar elabar]... 
      =senstuff(humansbar);

 
 %show the two components: sensitivity and std
 figure(8)
 subplot(3,1,1), plot(humansstd(:),'r')
 ylabel('standard deviation')
 subplot(3,1,2), plot(senhumansbar(:),'b')
 ylabel('sensitivity')
 
 %so you can see these components differ alot
 %what happens when you look at their product?
 senbystd= senhumansbar.* humansstd;
 subplot(3,1,3), plot(senbystd(:),'c')
 ylabel('contribution')
 xlabel('transition parameters,from (1,1) to (n,n)')
 
 figure(9)
 bar3(senhumansbar)
 title('Sensitivity for mean matrix') 

 figure(10)
 bar3(humansstd)
 title('Standard deviation of matrix elements')
 
 figure(11)
 bar3(senbystd)
 title('Sensitivity weighted by deviation')
 
 % similarly we can look at sensitivity squared, weighted by variance
 sensq=senhumansbar.^2;
 humansvar=humansstd.^2;
 
 figure(12)
 subplot(3,1,1), plot(humansvar(:),'r')
 ylabel('variance')
 subplot(3,1,2), plot(sensq(:),'b')
 ylabel('sensitivity squared')
  %so you can see these components differ alot
 %what happens when you look at their product?
  sensqbyvar= sensq.* humansvar;
 subplot(3,1,3), plot(sensqbyvar(:),'c')
 ylabel('contribution')
 xlabel('transition parameters,from (1,1) to (n,n)')
 
  figure(13)
 bar3(sensq)
 title('Sensitivity of mean matrix squared') 

 figure(14)
 bar3(humansvar)
 title('Variance matrix elements')
 
 figure(15)
 bar3(sensqbyvar)
 title('Sensitivity squared weighted by variance')
 
 
%similarly one could look at the variance-covariance matrix
% and products of the sensitivities (ij,kl)...
%the dimension would be 100 by 100...