function PVAL = tau_perm(L1,L2,N,poly,sims,D,s)


%TAU_PERM compares two point populations L1,L2 by doing cell counts for each with respect
%to N random reference points in a given area POLY containing the points. 

%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

% INPUTS:
%     (i)     L1   = locations of pop1
%     (ii)    L2   = locations of pop2
%     (iii)    N   = number of random reference points
%     (iv)  poly   = ([x11,x12],...,[xk-1,1,xk-1,2],[x11,x12])
%                  = (k:2)-matrix of points defining polygon
%                    (Note: last point must be same as first!)
%     (iv)  sims   = number of random permutations
%     (v)      D   = vector of distance values at which to measure k12-values.
%     (vi)     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

% Assertions for compiling

mbrealscalar(Dmax);
mbintscalar(N) ;




%Initial Calculations

n1 = size(L1,1);

n2 = size(L2,1);

k = length(D);


% Generate Points

      % Define containing rectangle
   
   	Mn = min(poly) ; 
   	Mx = max(poly) ;
   
   	x0 = Mn(1) ;
   	y0 = Mn(2) ;
   	x1 = Mx(1) ;
   	y1 = Mx(2) ;
   
   % Generate points
   
	disp(' ');
	disp('GENERATE POINTS');
	disp(' ');

   
   Pts = zeros(N,2);
   
   for i = 1:N
   
   	kk = 0 ;
   
   	counter = 1;
   
   	while kk == 0 ;
      
         r = rand(1);
         
         x = x0 + r*(x1 - x0) ;
            
         r = rand(1);            
            
         y = y0 + r*(y1 - y0) ;
      
      	if pt_in_poly([x,y],poly) == 1
         
         	Pts(i,:) = [x,y] ;
         
         	kk = 1;
         
         end %end if
         
         % Check to avoid infinite cycles
         
         counter = counter + 1;
         
         if (kk == 0)&(counter == 100)
            
           	error('looping error'); 
            
         end	         
         
      end %end point generation
      
   end %for
   

%Plot Points

xp = poly(:,1);

yp = poly(:,2);

%plot(xp,yp); hold on;

x = Pts(:,1);

y = Pts(:,2);

%plot(x,y,'.r');

%Compute Distances and Counts

disp(' ');
disp('COMPUTE DISTANCES AND COUNTS');
disp(' ');

C1 = zeros(N,k);

C2 = zeros(N,k);

for i = 1:N
   
   Diff1 = [L1 - ones(n1,1)*Pts(i,:)]';
   
   D1 = sqrt(sum(Diff1.^2))';
   
   Diff2 = [L2 - ones(n2,1)*Pts(i,:)]';
   
   D2 = sqrt(sum(Diff2.^2))';  
    
   C1(i,:) = tau_count(D1,D)';
   
   C2(i,:) = tau_count(D2,D)';
   
end

%Compute observed tau values

tau0 = zeros(k,1);

for j = 1:k
   
   x1 = C1(:,j);    
   x2 = C2(:,j);
   
   if (max(x1) == min(x1))|(max(x2) == min(x2))
      
      tau0(j) = NaN; %Kendall's Tau is not defined
      
   else
      
      tau0(j) = kendall_tau(x1,x2);
      
   end
      
end

%Now compute Tau Distributions

%First Construct Permutations for Testing

disp(' ');
disp('START PERMUTATIONS');
disp(' ');

PERM = zeros(N,sims);


for s = 1:sims
   
   PERM(:,s) = randperm(N)';
   
end

%Now construct TAU and P-Values

disp(' ');
disp('CONSTRUCT TAU');
disp(' ');

TAU = zeros(sims,k);

for j = 1:k
   
   if ~isnan(tau0(j))     
           
      for s = 1:sims
         
         x1 = C1(:,j);
         
         I = PERM(:,s);       
                 
         x2 = C2(I,j);
         
         TAU(s,j) = kendall_tau(x1,x2);
         
      end
      
   end %End if
   
end %End TAU Construction

disp(' ');
disp('CONSTRUCT P-VALUES');
disp(' ');

%Next construct PVAL

PVAL = zeros(k,1);

for j = 1:k
   
   v = - sort(-TAU(:,j)); %Sorts in descending order
   
   del = 1 ;
   
   row = 1 ;
   
   prob = 0 ;
   
   while del > 0
      
      if row == N+1

         
         prob = 1;
         
         del = 0;
         
      elseif (tau0(j) < v(row))
         
         row = row+1;
         
      elseif (tau0(j) == v(row)) %Treat Ties here
         
         I0 = row;
         
         index = 1;
         
         while index > 0 %Determine Rank
            
            row = row+1 ;
            
            if row == N+1
               
               I1 = N ;
               
               index = 0 ;
               
            elseif (tau0(j) == v(row))
               
               I1 = row;
               
            else %Must now have tau0(j) > v(row) 
               
               I1 = row;
               
               index = 0;
               
            end            
            
         end %end while
         
         rank = (I0 + I1)/2; %Define Rank
         
         prob = rank/(N+1);
         
         del = 0;
         
      else 
         
         prob = row/(N+1);  %No Ties
         
         del = 0;
         
      end %end rank procedure 
      
      PVAL(j) = prob;
      
   end %End while
   
   %disp(['PERCENT_DONE = ',num2str(100*j/k)]); 
   
end  %End PVAL construction