% This program calculates the near-term sensitivities presented in the
% paper 
% J. Yearsley, "Transient population dynamics and short-term sensitivity analysis of matrix
% population models", Ecological Modelling.
%
% The paper uses the cactus data (site A) from Fox and Gurevitch "Population numbers count:
% Tools for near-term demographic analysis.", Am. Nat., 156, p242--256
% 
% The program calculates the decomposition of the near-term sensitivity (S and DeltaS) and
% shows how the near-term sensitivity can be constructed from the
% decomposition.
%
% This script requires the files near_term.m, S.m and deltaS.m
%
% A number of additional Leslie matrices from a couple of sources are also
% included.
%
% Jon Yearsley 18/2/2004
% j.yearsley@macaulay.ac.uk
%
clear all

%% The cactus Leslie Matrices 

% Fox & Gurevitch cactus model  Site A
A = [[0.67  0 0.56]; [0.02 0.84 0]; [0 0.14 0.97]];
% Fox & Gurevitch cactus model  Site B
%A = [[0.49  0 0.56]; [0.01 0.73 0]; [0 0.23 0.99]];
% Fox & Gurevitch cactus model  Site C
%A = [[0.43  0 0.56]; [0.33 0.61 0]; [0 0.30 0.96]];

%% Additional Leslie Matrices for other classic examples

% Leslie matrix from Caswel p 166, lambda=1.0291
% A = [[0 0 10];[0.5 0.675 0];[0 0.075 0]];

% Teasel data from Caswell p60
% A = [[0 0 0 0 0 322.38];
%     [0.966 0 0 0 0 0];
%     [0.013 00.010 0.125 0 0 3.448];
%     [0.007 0 0.125 0.238 0 30.170];
%     [0.001 0 0 0.245 0.167 0.862];
%     [0 0 0 0.023 0.750 0]];

% A = [[0 0 0 0 0 0 431];
%      [0.749 0 0 0 0 0 0];
%      [0 0.966 0 0 0 0 0];
%      [0.008 0.013 0.01 0.125 0 0 0];
%      [0.07 0.007 0 0.125 0.238 0 0];
%      [0.002 0.001 0 0.036  0.245 0.167 0];
%      [0 0 0 0 0.023 0.750 0]];

%% Input parameters
% Element to vary = A(i,j)
i = 2;
j = 2;

% Time after perturbation
t = 5;
% And the initial population structure
n0 = [10, 5, 2]';

%% Perform short-term sensitivity analysis

% Calculate right eigenvectors and eigenvalues
[w,lambda] = eig(A);
[~, ind]=sort(real(diag(lambda)));
lambda = diag(lambda);
lambda = lambda(ind);

% Sort the right eigenvectors into order
w = w(:,ind);
v = inv(w)';
N = length(lambda);

% Prinicipal eigenvectors
w1 = w(:,N);
v1 = v(:,N);

% The equilibrium stage distribution
n_infty = sum(n0)/sum(w1);


% Deviation of initial population structure from asymptotic stable stage
% structure
delta_n = n0 - sum(n0)/sum(w1)*w1;

% Near-term sensitivity of population size for non-equilibrium population
S_near = near_term(A,i,j,t,n0,2);

% Near-term sensitivity for total population growth rate
S_lambda = near_term(A,i,j,t,n0,1);

% Near-term sensitivity for population growth rate stage by stage
S_lambda_i = near_term(A,i,j,t,n0,0)';

% Near-term sensitivity assuming initial population at equilibrium
S_equilib = near_term(A,i,j,t,sum(n0)*w1/sum(w1),2);


% Calculate the terms in the decomposition of the near-term sensitivity
s = S(A,i,j,t);

% s * sum(n0) / sum(w1) = S_equilib

Delta_s = delta_S(A,i,j,t);
Delta_s_t2 = delta_S(A,i,j,t+1);
nt = A^t*n0;
nt2 = A*nt;
GR = sum(nt2)./sum(nt);

% Check sesnitivity of population growth rate with a second calculation
S_lambda_check = sum(Delta_s_t2 - GR*Delta_s)*n0 / sum(nt);
S_lambda_i_check = (Delta_s_t2 - diag(nt2./nt)*Delta_s)*n0 ./ nt;

% Delta_s * delta_n = S_near - S_equilib

% Here's the decomposition explicitly calculated
S_near2 = s * sum(n0)/sum(w1) + Delta_s * delta_n;

%% Run a simulation where one element changes after t=t_perturb;

t_perturb = 20;
t_end = 50;
n = zeros(length(n0),t_end);
n(:,1) = n0;
for time = 1:t_end,
    if time==t_perturb,
        A(i,j) = 0.8*A(i,j);
    end
   n(:,time+1) = A*n(:,time); 
end

plot((n(:,2:end)./n(:,1:end-1))')
axis tight
xlabel('Time')
ylabel('n(t)/n(t-1)')