function [PVal,C0] = k_count_loc_ref(loc,ref,sims,D,M,poly)

% K_COUNT_LOC_REF is the same as K_COUNT_LOC, but now tests are run with respect
% to a reference set of (grid) points in an an areal lattice over the given region.
% This approach has the advantage of allowing an interpolated p-value map
% to be constructed.
%
% (1) A vector of distance values, D, must be prespecified
% (2) A set of reference locations, ref, must be prespecified
% (3) The output is a matrix of P-values at each refence point and each distance
%    value under the 'Randomness Hypothesis' with respect to the reference measure.

%Written by: TONY E. SMITH, 5/11/99 [modified 12/12/07]

%Programs Called: count_loc.m, polyform.m, rand_loc.m
%
% INPUTS:
%     (i)   loc  = population location file [loc(i)=(Xi, Yi),i=1:N]
%	  (ii)  ref  = reference location file [ref(k) = (Xk,Yk),k=1:n]
%     (iii) sims = number of simulations
%     (iv)  D    = set of distance values (in ASCENDING order)
%     (v)   M    = k-vector of measure values for each of k polygons
%     (vi)  poly = matrix describing boundaries of k polygons in the form:
%
%          [1   n1
%           x11 y11 
%           x12 y12
%           .
%           .
%           x11 y11
%           2   n2
%           x21 x22
%           .
%           .
%           .
%           k   nk
%           xk1 yk1
%           .
%           .
%           xk1 yk1]

% DATA OUTPUTS: PVal = (nxb)-matrix of P-Values, where b is the number of
%                       distances specified, and where each PVal(i,k) is the
%                       the P-value under the null hypothesis at reference 
%                       point i and distance k.
%
%                 C0 = Matrix of point counts at each distance

%Compiling assertions

%mbintscalar(sims);

% VARIABLE LIST

% s,N,n,b,k,C0,Z,u,M,P,F,Num,
% out,x,Freq,LL,CDF,d,inc,D0

%Initialize values

%s = 23311 ; %random seed

N = length(loc(:,1)) ; % number of locations

n = length(ref(:,1)) ; % number of reference points

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

k = length(M) ; % number of polygons



stage1 = 'Construct Initial Counts for given points in loc'

C0 = count_loc(loc,ref,D);


stage2 = 'Construct Polygon Sampling Distribution'

Z = polyform(poly) ;

u = ones(1,k);

m = u*M ;

P = M/m ;

% Form cdf

F = zeros(k,1);

F(1) = P(1) ;

i = 2 ;

while i <= k
   
   F(i) = P(i) + F(i-1) ;
   
   i = i + 1;
   
end


%Compute Sample Counts plus P-Values

FREQ = zeros(n,b) ;

i = 1 ;

stage = 'Begin Simulations'

while i <= sims
   
   % Construct Polygon Sample Numbers

	Num = zeros(k,1) ;

	ctr = 1 ;

	while ctr <= N
   
   	x = rand(1,1) ;
        
   	j = 1 ;
   
   	% Locate Polygon Interval containing x
   
   	while j <= k
      
      	if x <= F(j)
         
         	Num(j) = Num(j) + 1 ;
         
         	break
         
      	end
      
      	j = j + 1 ;
      
   	end
   
   	ctr = ctr + 1 ;
   
	end
  
   
   %Form random point sample
   
   LL = rand_loc(poly,Z,Num) ;        
      
   %Construct Cumulative Counts
   
   C = count_loc(LL,ref,D) ;
   
   Freq = zeros(n,b);
   
   for j = 1:n
      
      for kk = 1:b
         
         if C(j,kk) > C0(j,kk)
         
            Freq(j,kk) = Freq(j,kk) + 1;
            
         elseif C(j,kk)==C0(j,kk)
            
             Freq(j,kk) = Freq(j,kk) + .5;
            
         end       
         
         
      end   
   
   end
   
   FREQ = FREQ + Freq;
   
   disp(['simulation = ',num2str(i)]);
   
   i = i + 1;
   
end

PVal = (1 + FREQ) ./(sims + 1) ;