%stochseq_codeforcolleagues is a copy with slightly modified comments
% of stochseq_Kalisz_Aug 21, 2004: combined code of Carol and Tulja,
%(line in module 1 that gets the nstate corrected relative to Aug 15 version
%  varmat, sdmat, cvmat and cov_ij added to the output of module 2, uses n weighting not n-1 )
% Stanford University

%authorship rights
% the stochastic programs and the stochastic simulation reported here is
% described in Tuljapurkar, Horvitz and Pascarella, 2003 AmNat
% and in       Horvitz, Tuljapurkar and Pascarella, in press (as of Dec 2005) Ecology 

%to run this program you need
  %a mat-file with the demograhy matrices
  %a mat-file with the habitat transition rules
        %c(i,j) = Prob{state at time t+1 is i GIVEN that state at time t is j}

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

%the program is organized in modules

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%(0) before all modules: load the data
%(1) module 1: 
    
%you choose
   
    %tr  the number of time steps that are considered transient dynamics
    %ts  the number of time steps that you will actually use

%program generates

    %t2     the number of time steps in the sequence generated by the
                %simulation
    %nstage  the number of life history stages
    %nstate  the number of different states of the habitat in the system
 
    %sets up matrices and vectors to receive
    
    %states  the sequence of environmental states generated by the
                %probability rules
    %growth  the sequence of one time-step growth rates generated by the
                %simulation
    %uvecs   the sequence population structures generated by the
                %simulation
    %vvecs   the sequence of reproductive value vectors generated by the
                %simulation
    %cstar  the equilibrium distribution of states
    %avmat  the average matrix, Abar
    %varmat   the variance of each matrix element
    %sdmat    the standard deviation of each matrix element
    %cvmat    the coefficient of variability of each matrix element
    %cov_ij   the covariance of matrix elements 
 
%(2) module 2: uses the habitat transition rules to generate
      
     %states   sequence of state
     %cstar    equilibrium distribution of habitats
     %avmat    the mean matrix, Abar
     %varmat   the variance of each matrix element
     %sdmat    the standard deviation of each matrix element
     %cvmat    the coefficient of variability of each matrix element
     %cov_ij   the covariance of matrix elements 
     
%(3) module 3: uses the sequence of states to generate
 
     %growth    sequence of one-time growth rates
     %uvecs     sequence of population structure
     %vvecs     sequence of reproductive value vectors
  
 %uses the sequence of growth rates to calculate   
   
      %a         stochastic growth rate in "r" currency
      %lam_S     stochastic growth rate in "\lambda" currency

%(4) module 4: uses the sequences of vectors and growth rates to generate

      %sensE    E^S          sensitivity of lam_S to proportional perturbations 
                                %of ij^th rate in all states 
      %sensA    E^{S,\mu}    sensitivity of lam_S to proportional perturbations 
                                %of the mean of the ij^th rate  
      %sensV    E^{S,\sigma} sensitivity of lam_S to proportional perturbations 
                                %of variability of ij^th rate 
      %eijSbs   E^S_\beta    sensitivity of lam_S to proportional perturbations 
                                %of ij^th rate in a given, particular state 
      %ebS      summed habitat elasticity   
                              %sensitivity of lam_S to proportional perturbations 
                                %of all rates in a given, particular state 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%pre module activities for Kalisz data... is shown
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% others substitute your "ms" (the matrices of the states of the environment in array form
%%%% and your "habs" the matrix that describes the transitions among states of the environment 
%%%% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear
%load vitalrates ms; %loads the matrices in an array from the mat file called vitalrates
load sites ms_bt; %loads the matrices in an array from the mat file called sites
ms=ms_bt; %renames the array 

load kalisz_habs; %loads four versions of habitat transition matrices 
         
cmat = hab3; %choose this one for this run
                    %hab4 uses all 12 matrices with equal probability 
                    %independent of current habitat state
                    %hab3 is periodic for 3 states         
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%MODULE 1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tr=2000;            %set the transients at 2000
ts=100000;   %set the number you actually want at 100,000
t2=ts + 2*tr; %simulate 100,000 plus 2000 transients at beginning and end, respectively
nstage=length(ms(:,:,1)); 
nstate=length(cmat);
states=zeros(1,t2+1); %to receive the sequence of states
growth=zeros(1,ts);   %to receive the sequence of one-time growth rates
uvecs=zeros(nstage,ts); %to receive the sequence of population structures
vvecs=zeros(nstage,ts); %to receive the sequence of reproductive value vectors
avmat=zeros(nstage);    %to receive the average matrix
varmat=zeros(nstage);   %the variance of each matrix element
sdmat=zeros(nstage);    %the standard deviation of each matrix element
cvmat=zeros(nstage);    %the coefficient of variability of each matrix element
cov_ij=zeros(nstage^2); %the covariance matrix of matrix elements 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%MODULE 2  
%guts of  mchain2.m, modified such that %p(i,j) defined as Prob (j->i );
%                                       and calculates Abar, 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%outputs: 
         %states   sequence of states
         %cstar    equilibrium distribution
         %avmat    the mean matrix, Abar
         %varmat   the variance of each matrix element
         %sdmat    the standard deviation of each matrix element
         %cvmat    the coefficient of variability of each matrix element
         %cov_ij   the covariance matrix of matrix elements 

         
%inputs: %cmat     habitat transition probability matrix,
                      %p(i,j) defined as Prob (j->i );
           %t2      length of run
     
% determine equilibrium dist of hab states
[vec,c] = eig(cmat);
[cmax,loc] =max(abs(diag(c)));
cstar = vec(:,loc);
cstar = cstar/sum(cstar);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%calculate Abar  
%calculate covariances
    %adapted from morris and doak, p.235, 
    %but uses n rather than n-1 weighting
   
% cov(x,y) = E(x*y) - E(x)*E(y),
                                 % you need the means of each parameter 
                                 %  & the mean of the products of parameters
Abar=zeros(nstage^2,1);
Exy=zeros(nstage^4,1);

for i=1:nstate
    extractmat=ms(:,:,i);
    A=extractmat(:); %chooses the matrix for state i, represent it as a column vector that is nstages^2 long
    Exy=Exy+cstar(i)*kron(A,A); %kron(A,A) gives a vector that is nstages^4 long, since A is nstages^2 long
                                % the pairwise products of parameters,
                                % weight by hab freq, accumulate to get mean
                                    %try a vector that is [1 2 3 4]' to see how
                                    %this works, in this simple example the
                                    %first 4 are 1 times [1 2 3 4]'
                                    %next  4 are 2 times [1 2 3 4]', 
                                    %next  4 are 3 times [1 2 3 4]', 
                                    %next  4 are 4 times [1 2 3 4]'
    Abar=Abar+cstar(i)*A;     %weight by hab freq, accumulate to get mean matrix
end
%C=(Exy-kron(Abar,Abar))*nstate/(nstate-1); % Carol isn't convinced it should be weighted this way
C=(Exy-kron(Abar,Abar)); % Carol thinks it should be weighted this way
cov_ij=reshape(C,nstage^2,nstage^2);  %puts it into the form of the variance-covariance matrix
avmat=reshape(Abar,nstage,nstage);%puts it into the form of the population projection matrix
varmat=diag(cov_ij);
varmat=reshape(varmat,nstage,nstage);
sdmat= sqrt(varmat);
cvmat = sqrt(varmat); %intermediate step 
cvmat(avmat>0) = cvmat(avmat>0)./avmat(avmat>0);%the coefficients of variability for the nonzero Abars
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%make rand nos
rx = rand(1,t2+1);

%%%%%%submodule combines "cumvec", "randcomp" %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
probvec = cstar;           %choose initial state at random from equilibrium vector
for i = 1:t2+1
   psum = cumsum(probvec); %find the cumulative probability for a given vector
   j = min(find(psum > rx(i))); %get random number at time i from the rx vector
                                % and decide which state it corresponds to
   states(i) = j;          %put the state into the sequence of states 
   probvec = cmat(:,j);    %choose new column from the matrix as the next probability vector
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%MODULE 3
%generate sequences of growth rates, population and reproductive value vectors,
%  calculate the stochastic growth rate
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% the spirit of projection... as we move through time
%  we project the population, but each time we change the matrix,
%  according to the sequence of states.
% Each time step we observe the size of the population and its shape.
% We can calculate the N(t+1)/N(t) for each step...producing a sequence of
%    one-time growth rates...from these we calculcate a 
%    stochastic population growth rate for the dynamic environment..

vec1=ones(1,nstage);
vec1=vec1/sum(vec1);
vec1=vec1'; %the 3 lines create a vector to use as an initial population vector
vec2 = vec1;
    
trun=ts+tr; %time to run is desired length of sample path plus transients

%set up a new cumulative growth for all trun
	mgrowth = 0;
    %mgrowth2 = 0;
       for i = 1:trun
          i2=states(i);
		  mat1 = ms(:,:,i2);
          vec1 = mat1* vec1;
          growth1 = sum(vec1);
          vec1= vec1/growth1;
          mgrowth = mgrowth + log(growth1);
          %mgrowth2 = mgrowth2 + log(growth1)*log(growth1);
		if( i > tr)      %only store the ones after the transients
            i1 = i - tr;
           uvecs(:,i1) = vec1;
           growth(i1) = growth1;
		end
		
		%in reverse order for the rv's
		j = (t2 - i+1);
		%to make v(3000) use X(3000) but this comes from environment 3001
         i2 = (states(j+1));
          mat1 = ms(:,:,i2);
		%the next line Carol rewrote to make 'left' vectors
          vec2 = (vec2'*mat1)';
          vec2= vec2/(sum(vec2));
		if (i > tr)        %only store the ones after the transients
            vvecs(:, j-tr) = vec2;
          end
		end


a=mean(log(growth(1:ts)));
lam_S=exp(a);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%MODULE 4 
%  generate elasticities: E^S,  E^{S,\mu}, E^{S,\sigma}, E^S_\beta, 
%                        and summed habitat elasticity
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%set up matrices to receive the output
sensE =zeros(nstage); %E^S, perturb the rate in every state by same percent 
sensA =zeros(nstage); %E^{S\mu}, perturb the mean only 
sensV =zeros(nstage); %E^{S\sigma}, pertub the variability only
eijabS=zeros(nstate*nstage); %make hab-stage matrices for e's, intermediate step
eijSbs=zeros(nstage,nstage,nstate); %E^S\beta

%start loop
for i = (tr+1):(tr+ts-1) %choose index for states, starting one time step after transients
                         %go as far as the state just before the (tr+ts)^th state 

    %find state -- Keep in mind that: phenology of disturbance determines
    %whether you have picked state(i) or state(i+1) for the matrix

    %In either case that matrix becomes X(t)...Tulja .. u(t)=X(t)*u(t-1)
    %  growth(t) is also named by the end of the interval...
    %  growth(t) is the amount of growth from t-1 to t

i2 = (states(i+1));
mat1 = ms(:,:,i2);    
itime = i+1-tr; % if this is: t
i1 = i-tr;      % then this is: t-1

mat2S = mat1;
mat2A = avmat;
mat2V = mat1 - avmat;

mat2S = diag(vvecs(:,itime))*mat2S*diag(uvecs(:,i1));
mat2A = diag(vvecs(:,itime))*mat2A*diag(uvecs(:,i1));
mat2V = diag(vvecs(:,itime))*mat2V*diag(uvecs(:,i1));
scale1=( (vvecs(:,itime))' )*(uvecs(:,itime));

mat2S = (1/(scale1*growth(itime)))*mat2S; %separate it out for hab-stage                                            %use the growth rate at position 3 (t+1)
sensE = sensE + mat2S; %accumulate these for E_{ij}

mat2A = (1/(scale1*growth(itime)))*mat2A; %separate it out                                          %use the growth rate at position 3 (t+1)
sensA = sensA + mat2A; %accumulate these for E_{ij}

mat2V = (1/(scale1*growth(itime)))*mat2V; %separate it out                                           %use the growth rate at position 3 (t+1)
sensV = sensV + mat2V; %accumulate these for E{ij}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% submodule: hab-stage blocking
%lines are from "sensimod2.m" 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
blockrow = states(i-1); blockcol = states(i-2);
%so we are in block (blockrow, blockcol)
rowdex = nstage*(blockrow-1) + (1:nstage); 
coldex = nstage*(blockcol-1) + (1:nstage);
eijabS(rowdex, coldex) = eijabS(rowdex, coldex) + mat2S;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

sensE = sensE/(ts-1);
sensA = sensA/(ts-1);
sensV = sensV/(ts-1);
eijabS = eijabS/(ts-1);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%% calculate E^S_{beta}, and summed habitat elasticity
eijS = zeros(nstage,nstage);
ebS = zeros(nstate); 
   for b = 1:nstate
      colx = nstage*(b-1);
      matrix1 = zeros(nstage,nstage);
     for a = 1:nstate
      rowx = nstage*(a-1);
      matrix1 = matrix1 +  eijabS( (rowx + (1:nstage)), (colx + (1:nstage)) );
     eijS = eijS + eijabS( (rowx + (1:nstage)), (colx + (1:nstage)) );
     end
  eijSbs(:,:,b) = matrix1; %puts the matrix into the array
  ebS(b) = ebS(b) + sum(sum(matrix1));
end

ebS=(ebS(:,1));

%%%%%%%%%%%%%%%%%%%end set for eijS_\beta***************
 
%%%%%%%%%%%%%a few figures%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(1);bar3(sensE);axis ij;title('E_{ij}^S')
figure(2);bar3(sensA);axis ij;title('E_{ij}^{S\mu}')
figure(3);bar3(sensV);axis ij;title('E_{ij}^{S\sigma}')


for i=1:nstate; %make figures of habitat-stage elasticities
figure(3+i);
bar3(eijSbs(:,:,i));title(['E_{ij,\beta=' num2str(i) '}^S']); axis ij;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%