% DemoMetaSim. A program to do a multi-site stochastic
% demography simulation with correlations and a choice of 
% different distributions for the vital rates. 
% The code is now tailored to data for Coryphantha robbinsorum
% see Schmalzel et al. 1995
clear all;
%*****************Simulation Parameters***********************
% the order of the variables is survival and growth rates from 
% youngest to oldest, for each of the three sites (a,b,c), 
% and then the fecundity used in all the matrices:
% s1a	s2a	s3a	s4a	g2a	g3a	s1b	s2b	s3b	s4b
% g2b	g3b	s1c	s2c	s3c	s4c	g2c	g3c 	ff
load coryrates.txt; % contains the values for each rate
% (columns) for each year of data (rows)
vrates = coryrates; % assign data to the generic rate variable 
vrtypes= [ones(1,18),3]; % see VitalSim.m for explanation
vrmins=zeros(1,19); vrmaxs=zeros(1,19);
% Initial distribution across stages and sites:
n0 = [720, 26, 13, 24, 920, 6, 10, 23, 980, 12, 20, 31]'; 
makemx='makecorymx'; % assign matrix definition file: this
%  makes a separate matrix for each population (mx1 mx2 mx3)
%  and a matrix for all pop'ns, mx, same as G in Equ. 11.14
Ncap = [100, 100, 100]; % cap on juveniles and adults in pop'ns
Nx = [20 20 20]; % quasi-extinction threshold for each pop'n 
 % [this variable was called Ne in Morris and Doak -- changed to
 % be consistent with quasi-exinction variables in other programs]
tmax = 100; 	  	% number of years to simulate;
np = 19;  		% number of vital rates
dims = 12; 		% dimensions of total multi-site matrix;
runs = 10;  	% how many trajectories
popnum = 3; 		% number of separate populations
dimpop = 4; 		% number of classes per population
%*************************************************************
vrmeans = mean(vrates); 	% means of vital rates
vrvars=diag(cov(vrates)); % variances of vital rates
corrmatrix = corrcoef(vrates);   % correl. matrix for rates 
corrmatrix(:,19) = 0; % specific to this matrix: the last
corrmatrix(19,:) = 0; % parameter has no variance
corrmatrix(19,19) = 1;
for pp = 1:popnum % total each pop'n size w/o seeds
  popNstart(pp) = sum(n0((2+dimpop*(pp-1)):dimpop*pp)); end; 

randn('state',sum(100*clock)); rand('state',sum(150*clock)); 
[W,D] = eig(corrmatrix);  M12 = W*(sqrt(abs(D)))*W';

% INSERTED HERE: the section in VitalSim.m to make or retrieve
% sets of beta and stretched beta values to use in simulations: 
% replace np+np2 in this section with np. 
yesno = ... 
input('type 0 to calculate betas, or 1 to get a stored set');
if yesno == 1
  betafile = ...
input('type filename with stored betas; put in single quotes');
  load(betafile)% this line is corrected from that in Morris and Doak
else    % make a set of values for each beta or stretched beta
  parabetas=zeros(99,np);
  for ii = 1:(np)
	if vrtypes(ii) ~= 3
 	 for fx99 = 1:99
  		if vrtypes(ii) ==1; 
			parabetas(fx99,ii) = ...  
			 betaval(vrmeans(ii),sqrt(vrvars(ii)),fx99/100); end;
		if vrtypes(ii) ==2; 
			parabetas(fx99,ii) = ...
			 stretchbetaval(vrmeans(ii),sqrt(vrvars(ii)),...
			 vrmins(ii), vrmaxs(ii), fx99/100); end;
 	 end; % fx99
	end; % if vrtypes(ii)
  end; % ii loop
  yesno = input('type 1 to store the betas, or 0 if not');
  if yesno ==1 
    betafile = ...
        input('type filename to store betas, put in single quotes');
	save(betafile, 'parabetas'); % this line is corrected from that in Morris and Doak
  end; %if yesno 
 end; %else


% find the eigenvalues for each of the three populations:
vrs=vrmeans;			eval(makemx); 
[uu,lam1] = eig(mx1); 	lam1a = max(max(lam1));
[uu,lam1] = eig(mx2); 	lam1b = max(max(lam1));		
[uu,lam1] = eig(mx3); 	lam1c = max(max(lam1));
% initialize results trackers:
PrExt = zeros(tmax,1); 		logLam = zeros(runs,popnum); 
stochLam = zeros(runs,popnum); 

for xx = 1:runs;
 disp(xx); 
 nt = n0; 	% initialize population sizes
 extinct = zeros(1,popnum); % vector of extinction recorders
 for tt = 1:tmax
m = randn(1,np);
   	yrxy = (M12*(m'))';
% INSERTED HERE: the lines contained in the loop VitalSim.m 
%(for yy = 1:(np+np2)) to assign values for each rate, 
% each year. In this loop, replace np+np2 with np
 for yy = 1:(np)  % loops finds vital rate values
	 	if vrtypes(yy) ~= 3 % if not a lognormal rate
			index = round(100*stnormfx(yrxy(yy)));
			if index == 0 index = 1; end; % round at extremes
			if index ==100 index = 99; end;
      		vrs(yy) = parabetas(index,yy); % find stored value
		% else, calculated a lognormal value:
		else vrs(yy) = lnorms(vrmeans(yy),vrvars(yy),yrxy(yy));
		end;% if vrtypes(yy) ~= 3
    end; % yy loop
  

  	eval(makemx); % make a matrix with new vrs values. 
  	nt = mx*nt;	 % multiple by population vector
    for pp = 1:popnum % enforce population cap
        popNs(pp) = sum(nt((2+dimpop*(pp-1)):dimpop*pp)); 
        if popNs(pp) > Ncap(pp) 
            nt((2+dimpop*(pp-1)):dimpop*pp) = ... 
                nt((2+dimpop*(pp-1)):dimpop*pp)/Ncap(pp);
        end;
    end; % for pp
    if sum(extinct) < popnum % check for extinction   
	  for nn = 1:popnum
	      if popNs(nn) <= Nx(nn) extinct(nn) = 1;  end; 
      end;     
      if sum(extinct) == popnum PrExt(tt) = PrExt(tt) +1;  
       end; 
    end;     
 end % tt
logLam(xx,:) = (1/tmax)*log(popNs./popNstart);
stochLam(xx,:) = (popNs./popNstart).^(1/tmax);
end; % xx

CDFExt = cumsum(PrExt./runs); % make the extinction CDF
disp('deterministic lambdas for the three populations:');
display([lam1a, lam1b,lam1c]);  
% INSERT HERE: the final 6 lines in VitalSim.m to show results
disp('and this is the mean stochastic lambda'); disp(exp(mean(logLam)));
disp('below are mean and standard deviation of log lambda');
disp(mean(logLam)); disp(std(logLam));
disp('next is a histogram of logLams'); Hist(logLam);
disp('and then, the extinction time CDF'); figure; plot(CDFExt)