% mat4mega_08 
% May, 2008
%     uses mat-file "jpdata"
%     uses function file "senstuff"

%Disturbance mosaics (treefall gaps, hurricanes, etc.): 
%  Average over many patches or over many populations experiencing the same probabilities
%   of environments

% Demographic lessons
%		Megamatrix.
%		Megamatrix population growth rate and what influences it.
%		Invent a disturbance mosaic for your invented organism.
%	
% MATLAB lessons
%   eye  = identity matrix
%   zeros = matrix of zeros
%   kron = kronecker product 
%   num2str
%   function m-files in contrast to m-scripts
%   3-d figures
%   some graphics commands


% The steps are:
% 1. Load a set of population matrices representing different environmental states 
%    or patch types, say
%       m1, m2, m3...

% 2. Write a matrix for the probabilities of transitions among environments
%      the probabilities of going from m1 to m2, from m1 to m3 
%                                 from m2 to m1, from m2 to m3...etc

%     Let's call it 
%                   c 

% 3.  Write a time line for the phenology of your census and 
%         for the phenology of enviromental change... to decide when to apply
%         the environmental dynamic process 

% 4. Make the megamatrix that  includes the dynamics of environments and the dynamics 
%         of populations within each environmental state.
%         The megamatrix describes the probability that an individual will go from,
%         one environment-by-stage class to another in one time unit.
%  5. Analyze the asymptotic dynamics of the megamatrix, eigenvectors, eigenvalues,
%         sensitivity and elasticity.


clear all
close all

load jpdata
%loads the matfile that contains John Pascarella's dissertation data
% reference:
% Pascarella, J.B. and C.C. Horvitz. 1998. Hurricane disturbance 
% and the population dynamics of a tropical understory shrub:
% megamatrix elasticity analysis: Ecology 79: 547-563
%   some of these parameters were altered in subsequent papers
%   to make some reducible
%   matrices irreducible, etc, but the version you have here is 
%  the one published in this reference

m1
m2
m3
m4
m5
m6
m7

cmat

 %m1,m2, m3... are the demographic matrices for the 1, 2, 3... environments 
 %cmat is the environmental transitions

nstate=size(cmat,1)
nstage=size(m1,1)
%nstate is the number of environmental states
%nstage is the number of life history stages
% Note that the total number of patch-stages is given by the product nstate*nstage

%Which is the highest quality state of the habitat for the plant? 
% this will be answered in Figure 1

%  Use the command "eval" together with a string that contains another
%  command within which is an embedded "num2str" to get a variable coded with ending
%  numbers like m1, m2, m3...
for i=1:nstate
    eval([ 'lams(i)= senstuff(m' num2str(i) ')']); %the command we are evaluating is
         % "lams(i)=senstuff(mi)"
end
figure(1)
plot(lams,'o','LineWidth',3)
xlabel('Habitat')
ylabel('\lambda_{each habitat}')
title('Habitat quality')

%Which is the most frequent state of the habitat?
% This will be answered in figure 2
[lam, ssd, rv, sens, elas] = senstuff(cmat); % find the eigenvectors of cmat 
         %not these outputs all are relevant here...BUT it is convenient
fstar=ssd; %rename the right eigenvector as the equilibrium frequency vector
figure(2)
bar(fstar)
xlabel('Habitat')
ylabel('Frequency')
title('Equilibrium habitat frequency')

%NEXT: build the megamatrix
ckron = kron(cmat, eye(nstage));
%eye(nstage) is an identity matrix, nstage X nstage
% eye is a MATLAB function!, type 'help eye' 

%ckron is a matrix we have defined that has dimension environment-stage X environment-stage

% it results from taking 
% the Kronecker Product of the c matrix with an nXn identity matrix.
%     kron is the MATLAB function for Kronecker Product, type 'help kron'
% ckron looks like: 
%      it is subdivided into nstate*nstate blocks
%       each of which is the same size as the demography matrices (nstageXnstage); 
%       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 if you could think, in a case of 7 environmental state and 8 stages 
%     define a matrix that is zeros everywhere except the diagonal
%     which is entirely 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 C15 C16 C17
%         C21 C22 C23 C24 C25 C26 C27 
%         C31 C32 C33 C34 C35 C36 C37 
%         C41 C42 C43 C44 C45 C46 C47 
%         C51 C52 C53 C54 C55 C56 C57 
%         C61 C62 C63 C64 C65 C66 C67 
%         C71 C72 C73 C74 C75 C76 C77 ]
% where each of these is an 8x8 matrix, as is C11

%the following loop creates a large matrix that is also of the full patch-stage dimension

%if you want to look at ckron try highlighting and looking at these ONE at A TIME
ckron(1:8,1:56) %the top set of block matrices =    [C11 C12 C13 C14 C15 C16 C17]

ckron(9:16,1:56) %the second set of block matrices =[C21 C22 C23 C24 C25 C26 C27]

ckron(17:24,1:56)%the third set of block matrices = [C31 C32 C33 C34 C35 C36 C37]

ckron(25:32,1:56)%the fourth set of block matrices =[C41 C42 C43 C44 C45 C46 C47]
%etc

%the following loop creates a large matrix that is also of the full patch-stage dimension
for i = 1:nstate
i2 = num2str(i);
mat1 = eval(['m' i2]); %another way to use the "eval" command
rowindex1=nstage*(i-1)+1; rowindex2=nstage*i;
matrix2(rowindex1:rowindex2,rowindex1:rowindex2) = mat1;
end

%this matrix2 contains within it all the demographic matrices;
%     is also subdivided into nstate*nstate blocks, all of these blocks are
%      zeros everywhere except the blocks along the diagonal..and these are
%      the demographic matrices for each patch type
% an alternative way to get this matrix for our example
%    follows...it might be more intuitive to you
% of course you would have to enter it by hand if you changed dimensions of either
% demography or environment, whereas the above loop is more general

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

megL = ckron * matrix2;
megE = matrix2 * ckron;

%use "senstuff" to analyze the megamatrix
[lmat wmat vmat senmat elamat]=senstuff(megL);


%use the "text" command to place a label on a figure
%    the format is "text(x,y,['text of interest'])"
%        where x and y are the coordinates of where the string should go
% use the "num2str" (number to string) command to incorporate 
%      the value of a variable into your text label
figure(3)
plot(lmat,'ro')
title('Growth rate averaged over many patches: the Megamatrix \lambda')
text(1.1,lmat,['\lambda_M=' num2str(lmat)])

figure(4)
subplot(2,1,1)
bar(wmat)
title('stable environment-by-stage distribution, megL')
subplot(2,1,2)
bar(vmat)
title('reproductive values of each environment-by-stage class, megL')

%Now a figure of the elaticities
figure(5)

waterfall(elamat)
%this is one kind of 3-d graph
%try "help waterfall" 
%to compare to other 3-d graphs "help mesh", "help surface"

axis ij, grid off;
%the "axis ij" controls whether the arrangement is like a matrix 
% contrast with "axis xy" which puts the arrangement in cartesian mode
% try "help axis" 

AXIS([1 56 1 56 0 0.2]);
% manually sets  the x, the y and the z axis
% the x goes from 1 to 56, because there are 56 habitat-stage classes
% the y goes from 1 to 56, "
% the z goes from 0 0.2, set this high to be able to compare to another

shading faceted;
% shades the planes of the figure

%NOW for the axis labels
ylabel('Habitat-stages at t+1','FontName','Helvetica','FontSize',10,...
'FontWeight','bold')
xlabel('Habitat-stages at t','FontName','Helvetica','FontSize',10,...
   'FontWeight','bold')
zlabel('Elasticity','FontName','Helvetica','FontSize',16,...
   'FontWeight','bold')

%the title
title('Megamatrix, Late','FontName','Helvetica','FontSize',18,...
   'FontWeight','bold')

%NOW set up how many ticks there are and how to label each tick

%the following is sets some special specific text for labelling the ticks
% for this case in particular
xvalues=['1';'2';'3';'4';'5';'6';'7']; %labels the Habitat states

set(gca,'XTick',[1 9 17 25 33 41 49],'FontName','Helvetica','FontSize',10)
set(gca,'XTickLabel',{xvalues},'FontName','Helvetica','FontSize',10)
set(gca,'YTick',[1 9 17 25 33 41 49],'FontName','Helvetica','FontSize',10)
set(gca,'YTickLabel',{xvalues},'FontName','Helvetica','FontSize',10)


% OR GENERICALLY:
%  Following is a loop for  Megamatrix figure labels based only on
%  number of habitat states and number of life history stages

for i=1:nstate
    my_xvalues(i,:)=num2str(i);
    my_xticks(i)=nstage*(i-1)+1; % 8*(i-1) + 1 for JP data, try it for yours!
end

%then you can use
set(gca,'XTick',my_xticks,'FontName','Helvetica','FontSize',10)
set(gca,'XTickLabel',{my_xvalues},'FontName','Helvetica','FontSize',10)
set(gca,'YTick',my_xticks,'FontName','Helvetica','FontSize',10)
set(gca,'YTickLabel',{my_xvalues},'FontName','Helvetica','FontSize',10)

%QUESTIONS TO THINK OVER AS YOU SEE THESE RESULTS:

%Where is the "action" in this figure? 
% What is the megamatrix lambda most sensitive to?
% Is it for events in the most frequent habitat state?
% Is if for events associated with the highest quality habitat state?

%See Tuljapurkar, Horvitz and Pascarella, AmNat 2003 for further discussion
%of the kind of average population growth that is obtained here.