% Program logregextsim.m. This does stochastic simulations to 
% estimate the extinction CDF for a patch-based metapopulation 
% model. All the data file requirements, parameter definitions,
% and simulation methods are the same as for logregA and 
% logregB in Boxes 11.2 and 11.1, so comments are kept to a 
% minimum here. 
clear all;
%*****************User-defined Parameters*********************
load delugA.txt; % get basic information file for the 
siteinfo=delugA; % metapopulation to simulate
simreps = 500; % how many simulations to do
maxtime = 200; % number of years to replicate 
% Values for parameters (same order and meaning as in logregA)
eeparams = [0.233908  0.000000  -2.862026 5.849953  0.000000];    
ccparams = [-1.107957 0.000000   0.652810]; 
%*************************************************************
rand('state',sum(100*clock)); 
XX = siteinfo(:,1); 	% latitude in km
YY = siteinfo(:,2); 	% longitude in km
FF1 = siteinfo(:,3);	% occupancy state in year 1 (0,1)
FF2 = siteinfo(:,4); 	% occupancy state in year 2 (0,1)
AA = siteinfo(:,5); 	% size of site
sitenum = length(XX);	% number of sites
x1 = XX(:,[ones(1,sitenum)]);  x2 = x1';
y1 = YY(:,[ones(1,sitenum)]);  y2 = y1';
Dists = sqrt((x1-x2).^2 + (y1-y2).^2); 
AAmx = AA(:,[ones(1,sitenum)])'; % size matrix 
ee = eeparams(1:5)'; 	cc = [ccparams,eeparams(4:5)]; 
simcount = 0;		PrExt = zeros(1,maxtime);

for reps = 1:simreps  % loop of simulations 
 FFnow = FF2; % set initial occupancy as it was at last census
 disp(reps);                        
 for tt=1:maxtime % time loop

 % Next, to run, insert the lines in the program in Box 11.2 
 % that are between the same two comment lines as follow here:
 %----------Start of basic calculations---------------------
 %----------Start of basic calculations---------------------
 % matrix of occupancies for all neighbors:
 FFnowmx = FFnow(:,[ones(1,sitenum)])'; 
 % Calculation of neighbor effects for extinction probability:
 bit1=-ee(4)*Dists; % first part of extinction equation
 bit1(find(bit1>100))=100;	% check for values too big or small
 bit1(find(bit1<-100))=-100;
 connmxe = exp(bit1).*(AAmx.^(ee(5))).*FFnowmx;  
 connmxe = connmxe-diag(diag(connmxe)); % remove site's self    
 connmxe(find(connmxe>20))=20; % check for too big or small
 connmxe(find(connmxe<(-20)))=-20;
 % summed, distance-discounted areas of neighbors for ext calc:
 conne = sum(connmxe')'; % 

 % Calculate neighbor effects for colonization probability:
 bit1=-cc(4)*Dists;
 bit1(find(bit1>100))=100;
 bit1(find(bit1<-100))=-100;
 connmxc = exp(bit1).*(AAmx.^(cc(5))).*FFnowmx;
 connmxc = connmxc-diag(diag(connmxc));  
 connmxc(find(connmxc>20))=20;
 connmxc(find(connmxc<(-20)))=-20;
 % summed, distance-and- areas of neighbors for col. calc:
 connc = sum(connmxc')'; 

 % make the linear functions ue and uc to go into logistics:
 uue = ee(1) + ee(2)*AA + ee(3)*conne;
 uuc = cc(1) + cc(2)*AA + cc(3)*connc;
 % initialize logistic equations:
 loge=zeros(sitenum,1); logc=zeros(sitenum,1);

 % find the predicted probabilities of ext and col:
 uue(find(uue >10))=10; uue(find(uue<-10))=-10;
 loge = exp(uue)./(exp(uue)+1);
 uuc(find(uuc >10))=10; uuc(find(uuc<-10))=-10;
 logc = exp(uuc)./(exp(uuc)+1);

 % calculate Negative Log-Likelihood for the observed
 % extinctions and colonizations:
 if ((tt == 1) & (reps ==1))
    LLec=FFnow.*( (1-FF2).*log(loge) +FF2.*log(1-loge) ) ...
    + (1-FFnow).*( (FF2).*log(logc) +(1-FF2).*log(1-logc) );
    NLLec = -sum(LLec); % neg LL of extinctions and colon.'s
 end; % if 

 % randomly make next year's simulated occupancy:
 transprobs=rand(sitenum,1); 
 FFnew = FFnow.*(ceil(transprobs-loge)) + ...
(1-FFnow).*(ceil(logc-transprobs));
%----------End of basic calculations---------------------
 %----------End of basic calculations-----------------------

 % check for extinction
 if isempty(find(FFnew>0))
    PrExt(tt) = PrExt(tt) +1;  end; 
 extantsites(tt,reps) = sum(FFnew); % stores number occupied 
 FFnow = FFnew;

 end; % tt loop 
end; % reps loop

CDFExt = cumsum(PrExt/simreps);
disp('Below is mean and standard deviation')
disp('of number of pop"s each year');
disp([mean(extantsites')', std(extantsites')']);
disp('And now, the extinction time CDF');
figure;
plot(CDFExt);