%%%%%%%%%%%%%2 oct 06 %%%%% work in progress %%%%%%%%%%%%%%
%  (based on 14 sept 06 which was base on corrected code from calmex program)
%  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% 
%initial coding in july05 to implement analytical stuff from Tulja's notes 14june05
% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all
close all

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  MODULE 1
%
% in this module, we should 
%  
%  1) load the demographic and habitat trans matrices
%  2) specify the structures
%  3) correct for perfect one-period survival
%  4) choose the subsets of demographic matrices to be analyzed
%  5) choose the habitat dynamics matrix
%  6) get the dimensions of the set you have chosen
%demographic matrices must be in an array
%habitat transition must be in a matrx
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%name the species, set up
name=['Ardisia escallonioides, hurricane model'];
maxage=200; % pick a max age to which you will estimate survivorship
c=1;       %choose the demographic matrix set
chab=1;    %choose the habitat dynamics transition matrix
matfilename=['ares_demo' num2str(c) 'hab' num2str(chab)];% choose a name for 
                                                             % the output
                                                           % file
%load the data from a mat file
load cosmbase_mothnatural  %you will need to make the array

canopymats(:,:,1)=m1;
canopymats(:,:,2)=m2;
canopymats(:,:,3)=m3;
canopymats(:,:,4)=m4;
canopymats(:,:,5)=m5;
canopymats(:,:,6)=m6;
canopymats(:,:,7)=m7;
                   
%get the dimensions of your input system from one of the matrix arrays
nstages=size(canopymats,2);
nstates=size(canopymats,3);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%name the demo cases: how demography responds to environmental state change
cases= [', seven canopy states   '];                    

casemats=['canopymats']; %these are the actual array names

%chose the set of demo matrices to be analyzed and rename it mats
%c=1% choose case 1  %moved choice to top of program   

mats=eval(casemats(c,:));

%%%% correct for perfect survival %%%%%%%%%%%
for i=1:nstates;
mat=mats(:,:,i);
%stages are, eg:
%  seeds, sdlgs, juv prerep rep1 rep2
%  time step = 1 year (Aug to Aug)

matrix=mat;
perfects=find(sum(mat(2:nstages,:)) >= 1); % finds any columns that have perfect survival
%alter the perfect (or more!)survival columns ....make them .998 survival
m=length(perfects); %how many columns have perfect survival?
    for j=1:m         % for each column with perfect survival, designated, perfect(i), alter all entries
    matrix(2:nstages,perfects(j):nstages)=mat(2:nstages,perfects(j):nstages)*.998;
    end;

As(:,:,i)=matrix;
end;
% A's are like the original matrices, but no perfect survival

%%%%%%%%% define structures as in Cochran and Ellner %%%%%%%%%%

struct=zeros(nstages); %set the birth process matrices to zeros

%Put seeds in the B matrix by setting entries to 1's
B=struct;
B(1,2:nstages)=1; %specific to case

%Put small clonal offspring in the C matrix by setting entries to 1's
C=struct; %none in this case

%Put large clonal offspring in the V matrix by setting entries to 1's
V=struct; %none in this case

%Find the entries that are simply transitions among live states
T=ones(nstages); %set the transitions matrix to 1's
T(1,1)=0; % specific to the case at hand, no dormant seeds
T=T-B-C-V;       %subtract birth processes

%scale it up to the array
Ts=zeros(nstages,nstages,nstates); %
for i=1:nstates
    Ts(:,:,i)=T;
end

Ts= As.*Ts;  % the transition entries only 


%%% choose habitat transition matrix %%%%%%%%%%%%%%%%
%
%name the habitat dynamics cases: environmental state matrix 
%cmat will have been loaded with the above file
hab_cases=['historical hurricane frequency'];
               
                            
habmats= ['cmat '];  %in this case, cmat is also used as the name in the input file
                     %these are the names of the variables


%chose the matrix that describes the habitat dynamics and rename it cmat
%chab=1% choose case 1 %moved choice to top of program   
%chab=2% choose case 2
cmat=eval(habmats(chab,:)); %note the use of cmat as the output from this module

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% letters used correspond to Tulja's 14june05 notes
qs=Ts;                     %the q matrices in an array, one for each state, in the notes

%get the dimension of the sets you have actually chosen
nstates=size(cmat,1);        %S in the notes
nstages=size(qs(:,:,1),1);   %K in the notes

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% output variables of interest from module 1 
%         
%         maxage
%         name, c, cases, chab, hab_cases, matfilename
%         nstates, nstages, qs, cmat
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% MODULE 2 
%  
% ANALYZE habitat dynamics matrix
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%get the stationary habitat state distribution, pi in the notes and rho
[wmat,dmat]=eig(cmat);%asks for eigenvectors and eigenvalues
[lmat,j]=max(max(dmat));%gets the largest eigenvalue and its location,j
wmat=wmat(:,j);%gets the eigenvector in the jth location of wmat
pis=wmat./sum(wmat);%rescales the eigenvector to sum to 1

%get rho, the serial autocorrelation %%%%
d=diag(dmat); % the eigenvalues
absd=abs(d);   % gets their absolute magnitudes
[y,index]=sort(absd); %sort by absolute magnitude and get index for each 
rho=(d(index(nstates-1))); %finds the subdominant, which is in position nstates-1 = serial correlation
hab_eigs=d; %to save the eigenvalue spectrum of the habitat process

%%% note on 15 sep 06
%ASK T: can rho be a complex root...if so do we report its abs magnitude
%%% as the autocorrelation???? 
%ANSWER on 20 sep 06
%          yes, also take a look at the whole eigenvalue spectrum in the
%          complex plane
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% output variables of interest from module 2
%   
%         pis, rho, hab_eigs
%         
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%MODULE 3
%Deterministic analyses
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The individual environmental states as constant evironments
%%%%        constant environment survivorship, mortality 

ages=1:maxage; %a vector from one to max age
Nzero=zeros(nstages,1);
Nzero(1,1)=1;

qcon_surv=zeros(nstages,maxage,nstates);
for i=1:nstates
 for x=1:maxage %for individual environments
    qcon_surv(:,x,i)=sum(qs(:,:,i)^x)'; %all stages, each row is a stage
 end
end


%set up arrays to receive output
qcon_lx=zeros(nstates,maxage);
qcon_lxplus1=zeros(nstates,maxage-1);
qcon_mu=zeros(nstates,maxage-1);

for i=1:nstates
 qcon_lx(i,:)=qcon_surv(1,:,i); %take the stage 1 survivorships only
 qcon_lxplus1(i,:)=qcon_lx(i,2:maxage);
 qcon_mu(i,:) = log(qcon_lx(i,1:maxage-1)) - log(qcon_lxplus1(i,:));
 qcon_plateau(i)=-log(max(eig(qs(:,:,i))));
 qcon_oneperiodsurvival(i,:)=sum(qs(:,:,i)); %stage-specific survival
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%The Average q: the so-called "deterministic" case
% q's weighted by their frequency in habitat

qbar=zeros(nstages);
for i=1:nstates;
qbar=qbar+pis(i)*qs(:,:,i); 
end;      

%%%% "deterministic" survivorship, mortality, 
%%%%        and stage structure at each age %%%%%%%%%%%%%

ages=1:maxage; %a vector from one to max age
Nzero=zeros(nstages,1);
Nzero(1,1)=1;

for x=1:maxage 
    qbar_surv(x,:)= sum(qbar^x);
    qbar_trans(:,x)= qbar^x*Nzero/sum(qbar^x*Nzero); %stage structure at each age
end
qbar_lx=qbar_surv(:,1)'; %survivorship to age x, conditional on starting in stage 1,
                  % the other columns are conditional on starting in stages
                  % 2, 3, 4 etc
                  
qbar_lxplus1=qbar_lx(2:maxage);
qbar_mu = log(qbar_lx(1:maxage-1)) - log(qbar_lxplus1);

qbar_plateau=-log(max(eig(qbar)));

qbar_oneperiodsurvival=sum(qbar); %stage-specific survival

% distribution of age at death
qbar_deaths = qbar_lx(1:maxage-1)-qbar_lxplus1;
qbar_deaths=qbar_deaths/sum(qbar_deaths);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% output variables of interest from module 3
%    ages                 % a vector of ages from 1 to maxage
%    
%    qbar_lx
%    qbar_lxplus1
%    qbar_mu
%    qbar_plateau
%    qbar_deaths
%    qbar_oneperiodsurvival
%    qbar_trans           % the stage sructure at each age, transient dynamics
%
%    qcon_lx
%    qcon_lxplus1
%    qcon_mu
%    qcon_plateau
%    qcon_deaths
%    qcon_oneperiodsurvival
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%MODULE 4  
%  MEGAMATRIX FOR FUNDAMENTAL MATRIX FOR DOUBLY MARKOV PROCESS
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%borrowed kronecker product development of megamatrix 
% from mat4megapitt.m in the fivedaycourse folder
%  added notes for current application
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
ckron = kron(cmat, eye(nstages));
%eye(nstages) is an identity matrix, S X S
%ckron is a matrix we have defined that has dimension SK x SK
% the Kronecker Product of the c matrix with an S X S identity matrix.
% ckron looks like: 
%      it is subdivided into K X K blocks
%       each of which is the same size as the q matrices (S X S); 
%       in each block, all are zeroes, except the diagonal,
%       which has the same value throughout the block, 
%             one particular entry of the c-matrix
%  SO...  in a case of 4 environmental states and 8 stages 
%     define an 8 x8 matrix that is zeros everywhere except the diagonal
%         which is one c-matrix parameter
% 
%C11=[c11 0   0   0   0   0   0   0 
%       0   c11 0   0   0   0   0   0 
%       0   0   c11 0   0   0   0   0 
%       0   0   0   c11 0   0   0   0 
%       0   0   0   0   c11 0   0   0 
%       0   0   0   0   0   c11 0   0 
%       0   0   0   0   0   0   c11 0
%       0   0   0   0   0   0   0   c11 ]
%
% then 
%ckron = [C11 C12 C13 C14 
%         C21 C22 C23 C24  
%         C31 C32 C33 C34  
%         C41 C42 C43 C44] 
%        
% where each of these is an 8x8 matrix, as is C11

%the following loop creates a large matrix 
% that is also of the full SK x SK dimension

for i = 1:nstates
mat1 = qs(:,:,i); %choose the ith q matrix
rowindex1=nstages*(i-1)+1; rowindex2=nstages*i; %gets you into the block
matrix2(rowindex1:rowindex2,rowindex1:rowindex2) = mat1;
end

%this matrix2 contains within it all the q matrices;
%     is also subdivided into S x S blocks, all of these blocks are
%      zeros everywhere except the blocks along the diagonal..
%        and these are the q matrices, one for each state 

% an alternative way to get this matrix for our example
%    follows...it might be more intuitive to some
% of course one would have to enter it by hand 
%    if one changed dimensions of either demography or environment, 
%    whereas the above loop is more general

%z8=zeros(8)
%matrix2a=[m1 z8 z8 z8 
%          z8 m2 z8 z8  
%          z8 z8 m3 z8 
%          z8 z8 z8 m4] 
%    

m = matrix2 * ckron;           %the m matrix, of eqn 23 of T's 14june notes
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
chekm=(qs(:,:,2)*cmat(2,1))-m(nstages+1:2*nstages,1:nstages) %
%this should be exactly 0's, it is
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Itilda=eye(nstages*nstates)    %eqn 24

eatildas=zeros(nstates*nstages,nstages,nstates); % set up array to receive eatildas
for i=1:nstates
eatildas(:,:,i)=Itilda(:,(i-1)*nstages+1:nstages*i);
end;                           %array of eatilidas, eqn 25
                               % all rows, columns read from Itilda 1-8, 9-16,
                               %                                    17-24, etc

qatildas=zeros(nstates*nstages,nstages,nstates); % set up array to receive qatildas
for i=1:nstates
qatildas(:,:,i)=matrix2(:,(i-1)*nstages+1:nstages*i);
end;                            %array of qatildas, eqn 26

etilda=zeros(nstates*nstages,nstages); % set up matrix for etilda, only one not an array
for i=1:nstates
etilda((i-1)*nstages+1:nstages*i,:)=eye(nstages);
end;                            %etilda, eqn 27

I=eye(nstages);                 %identity matrix of dimension K x K

Ns_markov=zeros(nstages,nstages,nstates); %set up for array of Ns
for i=1:nstates;
Ns_markov(:,:,i)=I + etilda'*(Itilda-m)^-1*qatildas(:,:,i); %eqn 28, conditional on initial state 1
end;

%%%%average over all initial states %%%%%

Nbar_markov=zeros(nstages);
for i=1:nstates;
Nbar_markov=  Nbar_markov+pis(i)*Ns_markov(:,:,i); 
end;                         %eqn 29, weight each 
                                 % fundamental matrix by expected frequency 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%lifeexp, column sums
lifeexp_markov=zeros(nstates,nstages); %set up array
for i=1:nstates
lifeexp_markov(i,:)=sum(Ns_markov(:,:,i));
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% output variables of interest from module 4
%
%       m
%       Ns_markov
%       Nbar_markov
%       lifeexp_markov
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% MODULE 5  the variance  of the FUNDAMENTAL MATRIX
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Ntildas=zeros(nstages,nstages,nstates); %set up for Ntildas array
for i=1:nstates
Ntildas(:,:,i)=etilda'*m*(Itilda-m)^-1*eatildas(:,:,i); %eqn 30, one for each initial state
end;

for a=1:nstates
   for b=1:nstates
Nhat(:,:,b,a)=(eatildas(:,:,b))'*m*(Itilda-m)^-1*eatildas(:,:,a);%eqn 31, one for each a,b
   end;
end;

%now for the expectation of the nu_ij^2's given initial environment a
% eqn32
%       try to do it for one initial state, a=1....
%           start with the second term of eqn. 32
%           diag(diag(Ntildas(:,:,b)) *Nhat(:,:,b,a)
%   
%              set up arrays to calculate the variances conditional on each initial state

secondterm=zeros(nstages,nstages,nstates);
J=zeros(nstages,nstages,nstates);
Ysq_markov=zeros(nstages,nstages,nstates);
Yvar_markov=zeros(nstages,nstages,nstates);
Ycv_markov=zeros(nstages,nstages,nstates);

for i=1:nstates;
    
    sumterm=zeros(nstages,nstages);
    for b=1:nstates
    sumterm= diag(diag(Ntildas(:,:,b))) *Nhat(:,:,b,i) +sumterm;
    end    
    
  secondterm(:,:,i)=sumterm;
  %secondterm(:,:,i)=2*secondterm(:,:,i)  
  % I think this is an error because this is now multiplied by 2 as a
  % component of the sum J; check this with T
 

J(:,:,i)= secondterm(:,:,i)...  
         + diag(diag(Ns_markov(:,:,i)-I)) ...
         + diag(diag(Ntildas(:,:,i)))*qs(:,:,i);

Ysq_markov(:,:,i)=Ns_markov(:,:,i) + 2.*J(:,:,i); 
Yvar_markov(:,:,i)=Ysq_markov(:,:,i)-Ns_markov(:,:,i).*Ns_markov(:,:,i); 
Ycv_markov(:,:,i) = sqrt(Yvar_markov(:,:,i))./Ns_markov(:,:,i);
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%code for lifeexp var :equations developed by T  in 6july05 notes
lifeexpvar_markov=zeros(nstates,nstages);% set up matrix 
lifeexpcv_markov=zeros(nstates,nstages);
secondthing=zeros(nstages,nstages,nstates);
JJ=zeros(nstages,nstages,nstates);
KK=zeros(nstages,nstages,nstates);
Etasq_markov=zeros(nstates,nstages);

%expectation of squares of eta
for i=1:nstates;
    
    sumterm2=zeros(nstages,nstages);
    for b=1:nstates
    sumterm2= Ntildas(:,:,b) *Nhat(:,:,b,i) + sumterm2;
    end    
   
  secondthing(:,:,i)=sumterm2;

JJ(:,:,i)=                        secondthing(:,:,i) + Ns_markov(:,:,i) + Ntildas(:,:,i)*qs(:,:,i);
KK(:,:,i)= Ns_markov(:,:,i) + 2.* JJ(:,:,i);   

Etasq_markov(i,:)= sum(KK(:,:,i)); 
% I think the equation in the notes was wrong as etilda' is not a row of 1's...it is a
% block matrix of identity matrices...[I I I I];

%%%ok... now on the the variance and Cv
    lifeexpvar_markov(i,:) = Etasq_markov(i,:) - sum(Ns_markov(:,:,i)).*sum(Ns_markov(:,:,i));
    lifeexpcv_markov(i,:) = (sqrt(lifeexpvar_markov(i,:)))./sum(Ns_markov(:,:,i));
end;    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%used the iid case and substituted the markov N's

%finally equation 20, averaging over all initial conditions
thirdterm=0;
for i=1:nstates    
thirdterm=thirdterm + pis(i)*(Ns_markov(:,:,i).*Ns_markov(:,:,i));  %eqn 20
end 

Ybar_var_markov=2*diag(diag(Nbar_markov))*Nbar_markov - Nbar_markov - thirdterm;
Ybar_cv_markov = sqrt(Ybar_var_markov)./Nbar_markov;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% output variables of interest from module 5  (fundamental matrix variance)
%
% lifeexpvar_markov
% lifeexpcv_markov
% 
% Ybar_var_markov
% Ybar_cv_markov
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% MODULE 6
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%survivorship and hazard rates  7july05 notes %%%

%maxage=80; % maxage set at in module 1
ages=1:maxage; %a vector from one to max age
Thats=zeros(nstages,nstages,nstates,maxage);

for a=1:nstates;
    
for i=1:maxage;
       Thats(:,:,a,i)= etilda'*m^(i-1)*qatildas(:,:,a); %I think this should be i-1, Tulja had i-2
    end;

end;   


surv_allstages=sum(Thats,1);%survivorship of each stage to each age
chek=surv_allstages(:,:,2,10)-sum(Thats(:,:,2,10))
surv_1=surv_allstages(1,1,:,:); %get the element for those that start in stage 1
surv=squeeze(surv_1); %age-specific survivorships 
                        %rows are the states and the columns are the ages
        %so to reference surv for state 1...surv(1,:)
        %                     for state 2...surv(2,:)

shift_surv=surv(:,2:size(surv,2)); %lxplus1
mu = log(shift_surv) - log(surv(:,1:size(shift_surv,2))); % see next step, too!
mu_markov = -mu;
        %mu for state 1 ...mu(1,:)
        %mu for state 2 ...mu(2,:)
surv_markov = surv;      

surv_markov_early=surv_markov(:,1:10);
surv_markov_midearly=surv_markov(:,11:20);
surv_markov_midlate=surv_markov(:,21:30);
surv_markov_late=surv_markov(:,31:40);

%get better look at effects of initial state on early survival
surv_markov_early1=repmat(surv_markov_early(1,:),nstates,1); %repeats the row nstates times
rel_surv_markov_early=surv_markov_early./surv_markov_early1;

%get better look at effects of initial state on mid survival
surv_markov_midearly1=repmat(surv_markov_midearly(1,:),nstates,1);
rel_surv_markov_midearly=surv_markov_midearly./surv_markov_midearly1;

% distribution of age at death
deaths = surv(:,1:size(shift_surv,2))-shift_surv;
sumdeath=sum(deaths,2); %sum for each row
for i=1:nstates
    deaths_markov(i,:)=deaths(i,:)/sumdeath(i);
end;   %to divide each element by the corresponding row sum

 %%%% "killing rate" (sensu Steinsaltz) 
 %%%%     and "quasistationary distribtuion" 
lams_markov=zeros(1,1);
slope_markov=zeros(1,1);
us_markov=zeros(nstages,nstates);
vs_markov=zeros(nstages,nstates);
zs_markov=zeros(1,nstates); 

 
[wmat,dmat]=eig(m);%asks for eigenvectors and eigenvalues of megamatrix
[lmat,j]=max(max(dmat));%gets the largest eigenvalue and its location,j
wmat=wmat(:,j);%gets the eigenvector in the jth location of wmat
u=wmat./sum(wmat);%rescales the eigenvector to sum to 1
lams_markov=lmat;

[vmat,dmat]=eig(m');%asks for eigenvectors and eigenvalues
[lmat,j]=max(max(dmat));%gets the largest eigenvalue and its location,j
vmat=vmat(:,j);%gets the eigenvector in the jth location of vmat
v=vmat./sum(vmat);%rescales the eigenvector to sum to 1

us_markov=reshape(u,nstages,nstates);      
u_markov=sum(us_markov,2);%row sums
vs_markov=reshape(v,nstages,nstates);

slope_markov = -log(lams_markov); %%%%% - slope of log(survivorship) and zs
plateau_markov = slope_markov;

for i=1:nstates
    vq(i,:)=vs_markov(:,i)'*qs(:,:,i);
    zs_markov(i)=log(vq(i,1)) - log(v'*u);  %
end   


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% output variables of interest from module 6
%    
%    Thats
%    surv_markov
%      surv_markov_early
%      surv_markov_midearly
%      surv_markov_midlate
%      surv_markov_late
%    shift_surv
%    mu_markov
%    lams_markov
%    plateau_markov
%    deaths_markov
%    us_markov
%    vs_markov
%    zs_markov
%    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% MODULE 7
% SURVIVORSHIP VARIANCE
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%  18 july 05 %%%%%%%%%%%%%%%%%%%%%%%% 
Bi=eye(nstages^2);
B=zeros(nstages^2,nstages^2);

%A1=kron(cmat,kron(eye(nstages),eye(nstages)));
%alternatively and equivalent 
A1=kron(cmat,Bi); % this is parallel to ckron for the megamatrix


qkron=zeros(nstages^2,nstages^2,nstates);
qhat=zeros(nstages^2,nstages^2*nstates,nstates); %not yet transposed

ehat=Bi;
for a=2:nstates;
    ehat=[ehat Bi]; %not yet transposed
end;

for a=1:nstates
    qkron(:,:,a)=kron(qs(:,:,a),qs(:,:,a));
end

%changed to a loop...this is parallel to "matrix2" for the megamatrix 
% but the unit of the block size is nstages^2
  
%here was the code for a case with only 4 states of the habitat
%A4=[qkron(:,:,1) B B B];
%A4=[A4; B qkron(:,:,2) B B];
%A4=[A4; B B qkron(:,:,3) B ];
%A4=[A4; B B B qkron(:,:,4)];

%the following loop does the job by creating a large matrix 
 % you don't have to specify the size of the matrix ahead of time...
 % just increase the indices of the blocks and the rest is filled in with
 % zeros by MATLAB
 
for i = 1:nstates
mat5 = qkron(:,:,i); %choose the ith qkron matrix
rowindex1=nstages^2*(i-1)+1; rowindex2=nstages^2*i; %gets you into the block
matrix5(rowindex1:rowindex2,rowindex1:rowindex2) = mat5;
end

A4=matrix5;
%this matrix5 contains within it all the qkron matrices;
%     is subdivided into blocks, all of these blocks are
%     zeros everywhere except the blocks along the diagonal..
%     and these are the qkron matrices, one for each state 

R=A4*A1; %this is parallel to the megamatrix

for a=1:nstates
    qhat(:,((a-1)*nstages^2+1):a*nstages^2,a)=qkron(:,:,a)';
end    

Rx=eye(size(R));

 for x=1:maxage;
     for a=1:nstates;
         dummy=sum(ehat*Rx*qhat(:,:,a)');
      L(a,x)=dummy(1);
  end;
 Rx=R*Rx;
end;
  
%the variance at each age (equation 40)
for a=1:nstates
var_surv_markov(a,:)= L(a,:) - surv_markov(a,:).^2;
end;

%find negative numbers  
%   convert them to zero, in all test cases these negative numbers have
%   been tiny  and essentially zero... so no let's just set them to zero
%   but we can print them out just in case

[row,col]=find(var_surv_markov <0)

negs=[ ]
for i=1:length(row)
    negs(i)=[var_surv_markov(row(i),col(i))];
end

%print location and values of negatives
format short e
rowcolnegs=[row col negs']
format

%convert negatives to zeros
for i=1:length(row)
     var_surv_markov(row(i),col(i))=0;
end

cv_surv_markov=(var_surv_markov.^0.5)./surv_markov;
% 
%   
% beta as lambda of R
beta_markov=max(eig(R));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% output variables of interest from Module 7
%
%  R
%  var_surv_markov
%  cv_surv_markov
%  beta_markov
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

eval(['save ' matfilename])% saves a mat file of all the variables under the 
                           % name chosen in module 1



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%stop stop stop stop stop stop stop stop stop stop
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% NO MODULE 8 
% figures are now made in a separate file  
%   figure making loop
% this code needs to be standardized for any number of states, stages
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%STOP HERE...WE WILL DELETE THESE LINES EVENTUALLY
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(1)
colormap gray;
subplot(4,1,1);
bar(Ns_markov(:,:,1)');
%xlabel('groups=starting stage, stage j')
ylabel('N = E[\nu_{ij}]');
title(['Markov case ', cases(c,:),': expected yrs in stage i before death, \nu_{ij}']) 
text(0.1,7,'initial state = 1');
axis([0 9 0 15])
legend('i=1','i=2','i=3','i=4','i=5','i=6','i=7','i=8');

subplot(4,1,2);
bar(Ns_markov(:,:,2)');
%xlabel('initial stage, j')
ylabel('N = E[\nu_{ij}]');
text(0.1,7,'initial state = 2');
axis([0 9 0 15])
%legend('i=1','i=2','i=3','i=4','i=5','i=6','i=7','i=8');

subplot(4,1,3);
bar(Ns_markov(:,:,3)');
%xlabel('initial stage, j')
ylabel('N = E[\nu_{ij}]');
text(0.1,7,'initial state = 3');
axis([0 9 0 15])
%legend('i=1','i=2','i=3','i=4','i=5','i=6','i=7','i=8');

subplot(4,1,4);
bar(Ns_markov(:,:,4)');
xlabel('initial stage, j')
ylabel('N = E[\nu_{ij}]');
text(0.1,7,'initial state = 4');
axis([0 9 0 15])
%legend('i=1','i=2','i=3','i=4','i=5','i=6','i=7','i=8');


figure(2)
colormap gray;
subplot(3,1,1)
bar(lifeexp_markov');
ylabel('E[\eta_{j}]')
title(['Markov case ', cases(c,:),': life expectancy, E[\eta_{j}]'])
legend('initial state=1','initial state=2','initial state=3','initial state=4')
%xlabel('initial stage, j')
axis([0 9 0 25])

subplot(3,1,2)
bar(lifeexpvar_markov');
ylabel('Var[\eta_{j}]')
%title(['Markov case ', cases(c,:),': life expectancy, \eta_{j}'])
%legend('initial state=1','initial state=2','initial state=3','initial state=4')
%xlabel('initial stage, j')
%axis([0 9 0 25])

subplot(3,1,3)
bar(lifeexpcv_markov');
ylabel('CV[\eta_{j}]')
%title(['Markov case ', cases(c,:),': life expectancy, \eta_{j}'])
%legend('initial state=1','initial state=2','initial state=3','initial state=4')
xlabel('initial stage, j')
%axis([0 9 0 25])


%now one for the averaging over all initial states
figure(3)
colormap gray;
subplot(3,1,1)
bar(Nbar_markov');
%xlabel('initial stage, j')
ylabel('N = E[\nu_{ij}]');
title(['Markov case ', cases(c,:),': average over initial states' ])
legend('i=1','i=2','i=3','i=4','i=5','i=6','i=7','i=8');
axis([0 9 0 8])

subplot(3,1,2)
colormap gray;
bar(Ybar_cv_markov');
%xlabel('initial stage, j')
ylabel('CV[\nu_{ij}]');
%title(['markov case ', cases(c,:),': Average over initial states ' ]);
%legend('i=1','i=2','i=3','i=4','i=5','i=6','i=7','i=8');
axis([0 9 0 100])

subplot(3,1,3)
bar(sum(Nbar_markov));
ylabel('E[\eta_{j}]')
%title(['markov case ', cases(c,:),': Average over initial states, life expectancy'])
xlabel('initial stage, j')
axis([0 9 0 20])

%explore with figure
%figure(4)
%colormap gray;
%subplot(4,1,1)
%bar(Yvar_markov(:,:,1)')
%title(['Markov case ', cases(c,:),': Var(\nu_{ij})'])
%text(0.1,9,['initial state = 1' ]);
%legend('i=1','i=2','i=3','i=4','i=5','i=6','i=7','i=8')
%%xlabel('initial stage, j')
%ylabel('var[\nu_{ij}]')

%subplot(4,1,2)
%bar(Yvar_markov(:,:,2)')
%text(0.1,9,['initial state = 2' ]);
%%legend('i=1','i=2','i=3','i=4','i=5','i=6','i=7','i=8')
%%xlabel('initial stage, j')
%ylabel('var[\nu_{ij}]')

%subplot(4,1,3)
%bar(Yvar_markov(:,:,3)')
%text(0.1,9,['initial state = 3' ]);
%%legend('i=1','i=2','i=3','i=4','i=5','i=6','i=7','i=8')
%%xlabel('initial stage, j')
%ylabel('var[\nu_{ij}]')

%subplot(4,1,4)
%bar(Yvar_markov(:,:,4)')
%text(0.1,9,['initial state = 4' ]);
%%legend('i=1','i=2','i=3','i=4','i=5','i=6','i=7','i=8')
%xlabel('initial stage, j')
%ylabel('var[\nu_{ij}]')

figure(5)
colormap gray;
subplot(4,1,1)
bar(Ycv_markov(:,:,1)')
title(['Markov case ', cases(c,:),': CV(\nu_{ij})'])
text(0.1,75,['initial state = 1' ]);
legend('i=1','i=2','i=3','i=4','i=5','i=6','i=7','i=8')
%xlabel('initial stage, j')
ylabel('CV[\nu_{ij}]')
axis([0 9 0 150])

subplot(4,1,2)
bar(Ycv_markov(:,:,2)')
text(0.1,75,['initial state = 2' ]);
%legend('i=1','i=2','i=3','i=4','i=5','i=6','i=7','i=8')
%xlabel('initial stage, j')
ylabel('CV[\nu_{ij}]')
axis([0 9 0 150])

subplot(4,1,3)
bar(Ycv_markov(:,:,3)')
text(0.1,75,['initial state = 3' ]);
%legend('i=1','i=2','i=3','i=4','i=5','i=6','i=7','i=8')
%xlabel('initial stage, j')
ylabel('CV[\nu_{ij}]')
axis([0 9 0 150])

subplot(4,1,4)
bar(Ycv_markov(:,:,4)')
text(0.1,75,['initial state = 4' ]);
%legend('i=1','i=2','i=3','i=4','i=5','i=6','i=7','i=8')
xlabel('initial stage, j')
ylabel('CV[\nu_{ij}]')
axis([0 9 0 150])

%figure(6)
%colormap gray;
%bar(Ybar_var_markov');
%xlabel('initial stage, j')
%ylabel('Var[\nu_{ij}]');
%title(['Markov case, ', cases(c,:),': average over initial states ' ]);
%legend('i=1','i=2','i=3','i=4','i=5','i=6','i=7','i=8');

figure(7)
colormap white;
bar(pis)
xlabel('Environmental state')
ylabel('Frequency')
axis([0 5 0 1])
title(['Markov case ', cases(c,:),': stationary state distribution, \pi(i)'])
text(.5,.9,['\rho = ', num2str(rho)])

figure(8);
subplot(2,1,1)
semilogy(surv_markov');
ylabel('Survivorship (\it l_x)')
%xlabel('Age, yrs')
title(['Markov case ', cases(c,:)])
axis([0 85 10^-7 10^0])
legend('initial state=1','initial state=2','initial state=3','initial state=4')
text(10, 10^-1,['slope = ' num2str(slope_markov)])
text(10, 10^-1.5,['w      = ' num2str(ws_markov)])

subplot(2,1,2)
plot(mu_markov');
ylabel('Mortality rate (\it{\mu_x})')
xlabel('Age (\it{x})')

figure(9)
colormap gray;
subplot(2,1,1)
bar(us_markov')
title(['Markov case ', cases(c,:), ': quasi-stationary vectors'])
ylabel('Structure')
legend('stage=1','stage=2','stage=3','stage=4','stage=5','stage=6','stage=7','stage=8')
subplot(2,1,2)
bar(vs_markov')
ylabel('Reproductive value')
xlabel('Initial state')

%%%figures for examining early, mid and late survivorship
figure(10)
subplot(2,1,1)
bar(surv_markov_early')
ylabel('Survivorship (\it l_x)')
xlabel('x, yrs')
title(['Markov case: 1-10 yrs'])
legend('initial state=1','initial state=2','initial state=3','initial state=4')

subplot(2,1,2)
bar(surv_markov_midearly')
ylabel('Survivorship (\it l_x)')
xlabel('x+10, yrs')
title(['Markov case ', cases(c,:), ': 11-20 yrs'])


figure(11)
subplot(2,1,1)
bar(surv_markov_midlate')
ylabel('Survivorship (\it l_x)')
xlabel('x+20, yrs')
title(['Markov case ', cases(c,:), ': 21-30 yrs'])
legend('initial state=1','initial state=2','initial state=3','initial state=4')

subplot(2,1,2)
bar(surv_markov_late')
ylabel('Survivorship (\it l_x)')
xlabel('x+30, yrs')
title(['Markov case ', cases(c,:), ': 31-40 yrs'])

figure(12)
subplot(2,1,1)
%plot(rel_surv_markov_early')
bar(rel_surv_markov_early')
%ylabel('Survivorship ({\it l_x}) / Survivorship for initial state 1')
%xlabel('x, yrs')
title(['Markov case ', cases(c,:), ': 1-10 yrs'])
legend('initial state=1','initial state=2','initial state=3','initial state=4')
subplot(2,1,2)
plot(rel_surv_markov_early')
%bar(rel_surv_markov_early')
ylabel('Survivorship ({\it l_x}) / Survivorship for initial state 1')
xlabel('x, yrs')
%title(['Markov case ', cases(c,:), ': 1-10 yrs'])
legend('initial state=1','initial state=2','initial state=3','initial state=4')

figure(13)
subplot(2,1,1)
%plot(rel_surv_markov_midearly')
bar(rel_surv_markov_midearly')
%ylabel('Survivorship ({\it l_x}) / Survivorship for initial state 1')
%xlabel('x, yrs')
title(['Markov case ', cases(c,:), ': 11-20 yrs'])
legend('initial state=1','initial state=2','initial state=3','initial state=4')
subplot(2,1,2)
plot(rel_surv_markov_midearly')
%bar(rel_surv_markov_midearly')
ylabel('Survivorship ({\it l_x}) / Survivorship for initial state 1')
xlabel('x+10, yrs')

figure(14)
bar(zs_markov)
ylabel('{\it z_a}')
xlabel('initial state')
%legend('cycle=1234','cycle=2341','cycle=3412','cycle=4123')
title(['Markov case ', cases(c,:), ': intercepts'])
%legend('initial state=1','initial state=2','initial state=3','initial state=4')

figure(15)
subplot(2,1,1)
semilogy(var_surv_markov')
ylabel('Variance of survivorship, {\it l_x}')
title(['Markov case ', cases(c,:) ])
%xlabel('Age, yrs')
subplot(2,1,2)
plot(cv_surv_markov')
ylabel('CV of survivorship, {\it l_x}')
xlabel('Age, yrs')
text(1,4,['\beta = ' num2str(beta_markov) ])