% Program joincount.m. This program estimates the correlation
% in binary conditions or fates in different distance 
% categories, using a randomization procedure to test for
% significance. See Sokal and Oden 1978a for background
% The program requires information on a binary variable 
% (occupancy at one time, or colonization or extinction 
% across a time step) and the location of each site 
% in the sample. 

clear all;
%*****************User-defined Parameters*******************
distcats = [ 1 3 5 7 12]; % the upper limit of distance
% categories with which to partition the data (km)
% Load stored information on each site's location and fate or 
% state. Each row of data should be for one site, and the
% columns should be: X-coordinate of site, Y-coordinate, fate 
% or state (0 = no, 1 = yes). We use data on D. uliginosum 
load delug90s.txt;            % from Harrison et al (2000)
siteinfo = delug90s; % assign to the siteinfo variable       
randruns = 1000; % number of random rearrangements to 
% use for significance estimation
%*************************************************************
rand('state',sum(100*clock)); randn('state',sum(100*clock));
XX = siteinfo(:,1); % making variables to use in calculations:
YY = siteinfo(:,2); % XX,YY coordinates
FF = siteinfo(:,3); % FF fates
sitenum = length(XX); % how many sites
x1 = XX(:,[ones(1,sitenum)]); % making matrices of 
x2 = x1';			    % x and y values
y1 = YY(:,[ones(1,sitenum)]); % to calculate distances
y2 = y1';
Dists = sqrt((x1-x2).^2 + (y1-y2).^2); % pairwise distances 

for dd =1:length(distcats) % loop thru distance classes 
 disp('dd ='); disp(dd);
 updist = distcats(dd); % Upper and lower distance bounds 
 if dd == 1;   downdist = 0;
     else downdist = distcats(dd-1); end;
 ii=[]; % initialize the list of sites to include 
 for jj=1:sitenum   % find coordinates of sites to include
     ind = ...
      find((Dists(:,jj) <= updist)&(Dists(:,jj) >downdist));    
     if ~isempty(ind)  ii = [ii;jj];  end;
 end;
 XXd = XX(ii); % define the reduced sample that have joins 
 YYd = YY(ii); % of right length.
 FFd = FF(ii);
 sitenumd = length(XXd); % number used in this distance class
 
 x1 = XXd(:,[ones(1,sitenumd)]);	% making matrices of 
 x2 = x1';			   	% x and y values,
 y1 = YYd(:,[ones(1,sitenumd)]);	% then all the pairwise
 y2 = y1';			   	% distances between sites
 Distsd = sqrt((x1-x2).^2 + (y1-y2).^2);   
    
 f1 = FFd(:,[ones(1,sitenumd)]); % matrix of combined fates:
 f2 = f1';	% (0,0)=0, (0,1)and(1,0)=1, and (1,1)=2
 fatesmx = f1+f2; % making combined fates matrix
 fates=fatesmx(find((Distsd <= updist)&(Distsd >downdist))); 
  % a column of the combined fates with correct distances
 fatehist = hist(fates,3);  % find numbers of (0,0), (0,1), 
 realhist(dd,:) = fatehist; % and (1,1) combinations of fates.
       
 % Now, generate random rearrangements of fates to make 
 % distributions of expected joins with no correlations 
 for pp = 1:randruns % make all the random rearrangements
    FFr = FFd(randperm(sitenumd));  % randomize fates
    f1 = FFr(:,[ones(1,sitenumd)]);  
    f2 = f1';
    fatesmx = f1+f2; 		   % combined fates matrix
    fates=fatesmx(find((Distsd <= updist)&(Distsd >downdist)));      
    fatehist = hist(fates,3); 
    randhist(pp,:) = fatehist; % store results  
 end; % pp    
  % sort random results in ascending order for 00 pairs
  randnum00 = sort(randhist(:,1)); 
  if realhist(dd,1) <= min(randnum00)  % find ranking of real
        probs(dd,1) = 0	;		% result
    elseif realhist(dd,1) >= max(randnum00)
        probs(dd,1)=1;
    else rank = find(randnum00 < realhist(dd,1)); 
         probs(dd,1) = max(rank)/randruns; 
  end; % if realhist
        
  randnum01 = sort(randhist(:,2)); % repeat for 01 pairs
  if realhist(dd,2) <= min(randnum01)
        probs(dd,2) = 0;
    elseif realhist(dd,2) >= max(randnum01)
        probs(dd,2)=1;
    else rank = find(randnum01 < realhist(dd,2)); 
         probs(dd,2) = max(rank)/randruns; 
  end;
        
  randnum11 = sort(randhist(:,3)); % repeat for 11 pairs
  if realhist(dd,3) <= min(randnum11)
        probs(dd,3) = 0;
    elseif realhist(dd,3) >= max(randnum11)
        probs(dd,3)=1;
    else rank = find(randnum11 < realhist(dd,3)); 
         probs(dd,3) = max(rank)/randruns; 
  end;
  expected(dd,:) = [mean(randhist)]; % expected results       
end; % dd loop    

disp('Below are the observed numbers of 00 prs, the expected') 
disp('number, and the percentile rank for the observed');
disp('(Rows are for each distance class, short to long)')
disp([realhist(:,1),expected(:,1),probs(:,1)]);

disp('Below are the observed numbers of 01 prs, the expected') 
disp('number, and the percentile rank for the observed');
disp('(Rows are for each distance class, short to long)')
disp([realhist(:,2),expected(:,2),probs(:,2)]);

disp('Below are the observed numbers of 11 prs, the expected') 
disp('number, and the percentile rank for the observed');
disp('(Rows are for each distance class, short to long)')
disp([realhist(:,3),expected(:,3),probs(:,3)]);

disp('Finally, the numbers of pairs of 00,01,11'); 
disp('and the total in each distance category');
disp([realhist,sum(realhist')']);