%Age from Stage Workshop
% Age determination workshop: Oct 2007, Rostock, code from Carol Horvitz 
%  
%Age from stage in a constant environment
% To prepare:
%  1)separate production of new individuals from fates of surviving individuals
%    A=Q+R
%  2)get rid of immortal stages: make column sums of Q be  < 1

clear all;
close all;

%get the Q matrix

load ardesc_q1.txt  %input data from text file
q=ardesc_q1;

%pinus_platt;       %read another  matlab script with the data in it
%q=tt1;

nstages=size(q,1);%find the number of stages

%let's see how this works:

maxage=100;% set the number of ages
%maxage=350;
ages=1:maxage;          %a vector from one to max age
Nzero=zeros(nstages,1); %column vector with nstages rows
Nzero(1,1)=1;           %all individuals start out in stage 1

Nx=zeros(nstages,maxage); %
Qx=zeros(nstages,nstages,maxage);                            
 for x=1:maxage            %for each age
    Nx(:,x)=q^x*Nzero;     %the population vector in stages at each age
    Qx(:,:,x)=q^x;         %the powers of Q: gives you the same information
 end
   Nx2=squeeze(Qx(:,1,:)); %all rows, column 1 of each age
   lx=sum(Nx);

 %find the stage structure at each age as a proportion
 Ux=zeros(nstages,maxage)
 for x=1:maxage
    Ux(:,x)=Nx(:,x)/lx(x)
 end
 
%get mortality trajectory
lxplus1=lx(2:maxage);
mu = log(lx(1:maxage-1)) - log(lxplus1);
plateau =-log(max(eig(q)));

% get distribution of age at death
deaths = lx(1:maxage-1)- lxplus1;
deaths= deaths/sum(deaths);

%get oneperiod survival
oneperiodsurvival=sum(q); %stage-specific survival

%get the Fundamental matrix and the Life Expectancies
I=eye(size(q)); %identity matrix of dimension q
N=(I-q)^-1; % see page 118, Caswell, 2001,
    %N is the "fundamental matrix"...how long an organism in each stage will spend in each stage
% IF newborns are in stage 1
% (otherwise this is NOT correct) then the first element of lifeexp will be expected lifespan of newborns
% lifespan=sum(N(:,i)); %expected lifespan of individuals in stage i
lifeexp=sum(N); %Expected lifespan of individuals in each stage
lifeexp_var=sum(2*N^2 - N) - sum(N).*sum(N); % see page 119, Caswell, 2001
lifeexp_sd=(lifeexp_var).^(0.5);
lifeexp_cv=lifeexp_sd./lifeexp;


figure(1);
subplot(2,1,1)
semilogy(ages,lx,'k','linewidth',2)
ylabel('Survivorship, {\it l(x)}','FontSize',14)
xlabel('Age','FontSize',14)

subplot(2,1,2)
bar(deaths)
ylabel('Proportion of deaths','FontSize',14)
xlabel('Age','FontSize',14)

figure(2)
subplot(3,1,3)
plot(ages(1:maxage-1),mu,'k','linewidth',2)
ylabel('Mortality, {\it \mu(x)}','FontSize',14)
xlabel('Age, yrs','Fontsize',14)
hold on;
plot(ages(1:maxage),plateau,'k:','linewidth',2)
text(10,1.2*plateau,['-log(\lambda)']) 
hold off;

subplot(3,1,2)
bar(oneperiodsurvival)
xlabel('Stage','FontSize',14)
ylabel('One period survival','FontSize',14)
axis([0 9 0 1.1])

subplot(3,1,1)
area(Ux')
legend('1','2','3','4','5','6','7','8')
xlabel('Age','FontSize',14)
ylabel ('Proportion in stage','FontSize',14)
axis([0 maxage-1 0 1])

figure(3)
plot(lifeexp)
ylabel('Life expectancy','FontSize',14)
xlabel('Stage','FontSize',14)

%now go back and try it for Pinus palustris