function PVal = k12_perm(L1,L2,sims,D,s)

%K12_PERM compares two point populations L1,L2 by randlomly permuting their labels
%and computing K12-values for a range D of distance values. P-values are then computed, 
%reflecting significant attraction and/or repulsion between these patterns.

%NOTE 1: By convention, low P-Values here correspond to significant ATTRACTION between 
%        patterns.
%NOTE 2: Ties are treated neutrally, adding .5 rather 1 to rankings.

%Written by: TONY E. SMITH, 12/2/01 [modified 1/4/07]

% INPUTS:
%     (i)   L1   = locations of pop1
%     (ii)  L2   = locations of pop2
%     (iii) sims = number of random permutations
%     (iv)  D    = vector of distance values at which to measure k12-values.
%     (v)   s    = random seed

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

% SCREEN OUTPUT:  A plot of P-Values against distances in D.

% PROGRAMS CALLED: dist_vec

% Initial Calculations

rand('seed',s);

B = length(D);
n1 = length(L1(:,1));
n2 = length(L2(:,1));
n = n1 + n2;

loc = [L1;L2];

Dist = dist_vec(loc);

%mbrealvector(Dist);

Dmax = D(B);

Dmin = D(1);

C = zeros(n1,B) ;

for i = 1:n1  
   
   for j = n1+1:n     
            
       K = n*(i-1) - (i*(i-1)/2) + (j-i); %Counting rule for finding d(i,j) in Dist
   
   	   dist = Dist(K); 
         
      if dist <= Dmax
         
         if dist <= Dmin %j in all rings
            
            C(i,:) = C(i,:) + 1;
            
         else
            
            b = B;
            
            while D(b) >= dist
               
               C(i,b) = C(i,b) + 1;
               
               b = b - 1; %b = 1 already done
               
            end
            
         end %if <= Dmin
         
      end %if <= Dmax 
      
   end % for j
   
end %for i


%Compute Output

C0 = sum(C) ; %sum columns of C to get total counts

count = ones(1,B); %include C0 >= C0 as first count [1/4/07]

%Start Random Permutations

index = 10 ;

i = 1;

while i <= sims   
      
   C = k12_perm_count(loc,n1,Dist,D);
        
   inc = zeros(1,B);
   
   for j = 1:B
      
      if C(j) > C0(j)
         
         inc(j) = 1;
         
      elseif C(j) == C0(j)
            
         inc(j) = .5;
            
      else  
         
         inc(j) = 0;
         
      end
      
   end
   
   
   count = count + inc; 
   
   
    if 100*(i/sims)>= index
      
      disp(['percent_done = ', num2str(index)]);
      
      index = index + 10;
     
    end

   
   i = i + 1 ;
   
end %while


%Compute p-values 


PVal = count / (sims+1) ; %Small if attraction is strong