function PVal = k2_diff(loc,n1,sims,D,s)


%K2_DIFF does the calculations for K2_DIFF_PLOT.

%NOTE 1: By convention, P-Values at each distance h are heredefined to be  
%        the probability of getting a frequency count greater than the given
%        value IF pop1 were a random draw from pop1 + pop2. So low P-Values 
%        correspond to significant CLUSTERING of pop1 vs pop2.

%NOTE 2: Ties are treated neutrally, adding .5 rather 1 to rankings.

%Written by: TONY E. SMITH, 6/18/99

%Programs Called: k2ct_glo, rand_perm
%
% INPUTS:
%     (i)   loc  = combined location file [loc(i)=(Xi, Yi),i=1..n]
%                  where n = n1 + n2. It is assumed that points for
%                  pop1 are in the first n1 rows.]
%     (ii)  n1  = size of pop1.
%     (iii) sims = number of relabelings
%     (iv)  D    = vector of distance values at which to measure clustering.
%     (v)   s    = seed for random number generator

% OUTPUTS:  PVal = Vector of P-Values at each distance.


%compiling initializations

mbintscalar(n1);

mbintscalar(sims);

mbintscalar(s);


[n,junk] = size(loc) ;

n2 = n - n1;

mbintscalar(n2);

B = length(D);


%Compute counts for Pop1

Dist = dist_list(loc);

S1 = [1:n1];
S2 = [n1+1:n];

N1 = n1*(n1-1)/2;
N2 = n2*(n2-1)/2;

D1 = zeros(N1,3);
D2 = zeros(N2,3);

%Construct initial distance vector D1

k = 1;

for i = 1:n1-1
   
   I = S1(i);
      
   for j = i+1:n1
         
      J = S1(j);
         
      K = n*(I-1) - (I*(I-1)/2) + (J-I); %Counting rule for finding d(I,J) in Dist
         
      D1(k,:) = [Dist(K,1),i,j];
      
      k = k + 1;
      
   end
   
end

%Construct initial distance vector D2


k = 1;

for i = 1:n2-1
   
   I = S2(i);
      
   for j = i+1:n2
         
      J = S2(j);
         
      K = n*(I-1) - (I*(I-1)/2) + (J-I); %Counting rule for finding d(I,J) in Dist
                 
      D2(k,:) = [Dist(K,1),i,j];
      
      k = k + 1;
      
   end
   
end


C1o = k2ct_diff(n1,D1,D) ;
C2o = k2ct_diff(n2,D2,D) ;

C0 = (1/n1^2)*C1o - (1/n2^2)*C2o; %Density corrected


count = zeros(1,B);

%Start Simulations

i = 1 ;

index = 10 ;

ctr = 1;

while ctr <= sims
   
   %Form random relabeled sample
   
   [list,s] = rand_perm(n,s) ;
   
   S1 = list(1:n1) ;
   S2 = list(n1+1:n);
   
   %Construct distance vector D1

	k = 1;

	for i = 1:n1-1
   
  	 	I = S1(i);
      
   	for j = i+1:n1
         
      	J = S1(j);
         
      	K = n*(I-1) - (I*(I-1)/2) + (J-I); %Counting rule for finding d(I,J) in Dist
         
      	D1(k,:) = [Dist(K,1),i,j];
      
      	k = k + 1;
      
   	end
   
	end

	%Construct distance vector D2

	k = 1;

	for i = 1:n2-1
   
   	I = S2(i);
      
   	for j = i+1:n2
         
      	J = S2(j);
         
      	K = n*(I-1) - (I*(I-1)/2) + (J-I); %Counting rule for finding d(I,J) in Dist
         
      	D2(k,:) = [Dist(K,1),i,j];
      
      	k = k + 1;
      
      end
      
   end
   
     
   %Add to order counts
   
   
	C1 = k2ct_diff(n1,D1,D) ;
	C2 = k2ct_diff(n2,D2,D) ;
  
   C = (1/n1^2)*C1 - (1/n2^2)*C2; %Density Corrected
   
   inc = zeros(1,B);
   
   for j = 1:B
      
     	if C(j) > C0(j)
         
        	inc(j) = 1;
         
     	elseif C(j) == C0(j)
         
        	inc(j) = 0.5;
         
     	else
         
        	inc(j) = 0;
         
     	end
      
   end
   
   
   count = count + inc; 
   
   
   if 100*(ctr/sims)>= index
      
     	percent_done = index
      
      index = index + 10;
     
   end

   
   ctr = ctr + 1 ;
   
end %end while


%Compute p-values

u = ones(1,B);

PVal = count ./ ((sims+1)*u) ;