function PVAL = exact_clustering(loc_0,loc_1,D);

%EXACT_CLUSTERING considers the null hypothesis that loc_0 is a random
%sample from a larger population of N points, and for each
%point x in loc and distance radius d in D calculates the exact
%probability of getting as many points from loc_0 within distance d as the
%number observed. These serve as P-values for defining local clustering --
%a low P-value indicates strong clustering at x
% 
% INPUTS:   loc_0 = target pattern
%           loc_1 = pattern of all other points
%             D = selected radii {d1,..,dk}
%              
% OUTPUT:  PVAL = matrix of P-values, where Pij = prob of getting at least 
%                 as many points from loc_0 within distance dj of xi as the 
%                 number observed.       
         

n0 = length(loc_0); 

loc = [loc_0;loc_1]; %so loc_0 is the first n0 points in loc

n = length(loc);

%Compute pairwise distances


b = length(D); %number of distance values specified

C = zeros(n,b) ; %Matrix of counts

ub = ones(b,1);

Dmax = D(b) ; %Maximum distance value

PVAL = ones(n0,b);

% Start

for i = 1:n0
   
   % Compute distances from i to loc points
   
   x = zeros(1,b) ; %Vector of frequency counts in loc
   
   x0 = zeros(1,b) ; %Vector of frequency counts in loc_0
   
   for j = 1:n
      
      Dist = sqrt((loc_0(i,1) - loc(j,1))^2 + ...
         (loc_0(i,2) - loc(j,2))^2); 
      
      %Add to point frequencies      
            
      if Dist <= Dmax  %ignores distances larger than max in D
         
         for k = 1:b        
                                   
            if Dist <= D(k)
               
               x(k:end) = x(k:end) + 1 ; %adds to frequency counts to x
               
               if j <= n0
                   
                   x0(k:end) = x0(k:end) + 1 ; %adds to frequency counts to x0
                   
               end
                              
               break
               
            end
            
         end         
         
      end %End if
               
   end %end frequency count for point i  
   
   for k = 1:b
       
       PVAL(i,k) = 1 - hygecdf(x0(k) - 1,n,x(k),n0);  
       
%        NOTE: hygecdf(x0,n,x,n0) = prob of getting no more than x0 members 
%        of loc_0 in a sample of size x from loc.
       
   end
   
end %End procedure