function results = space_time_pvals(LT,sims,s,t)

%SPACE_TIME_PVALS computes space-time K-function for a given space-time pattern,
%and in addition, computes the (cluster) P-values for a set of random time 
%relabelings of that pattern. 

%Written by: TONY E. SMITH, 2/7/05

%Functions called: space_time.m
%
% INPUTS:
%     (i)   LT   = location-time file [loc(i)=(Xi,Yi,Ti),i=1..n]
%     (ii)  sims = number of relabelings
%     (iii)    s = number of point values to use in space range.
%     (iv)     t = number of point values to use in time range.

% OUTPUTS: 1x4)-structure, results
%     (i)    results.cdf   = (s x t) matrix of CDF values for original pattern.
%     (ii)   results.pvals = (s x t) matrix of P-values.
%     (iii)  results.dist  = (s x 1) array of distance values.
%     (iv)   results.time  = (t x 1) array of time values.
%           [Here the max pairwise distance, Dmax, and time interval,
%            Tmax,are computed and (d0,t0) = (Dmax/2,Tmax/2)then 
%            defines the relevant ranges.]

m = size(LT,1) ;

%Compute pairwise distances and Dmax

DIST = zeros(m*(m-1)/2, 3) ; %Matrix of loc-dist pairs

i = 1 ;

k = 1 ;  %counter for DIST.

while i < m
   
   j = i + 1 ;
   
   while j <= m
      
      DIST(k,1) = norm(LT(i,1:2) - LT(j,1:2));

      DIST(k,2) = i ;
      
      DIST(k,3) = j ;  
                
      j = j + 1 ;
      
      k = k + 1 ;
      
   end
   
   i = i + 1 ;
   
end

Dmax = max(DIST(:,1)) ;

Tmax = max(LT(:,3)) - min(LT(:,3)) ;

%Compute CDF for LT

CDF = space_time_ct(LT,s,t,DIST) ;

COUNT = ones(s,t); %include LT as a sample

%Start Simulations

i = 1 ;

index = 10;

while i <= sims
   
   %Form randomly relabeled sample
   
   sample = randperm(m) ;   
        
   LT(:,3) = LT(sample,3) ;
   
   %Construct Envelopes
   
   X = space_time_ct(LT,s,t,DIST) ;
   
   diff = X - CDF;
   
   INC = (diff > 0) + .5*(diff == 0); %include half weights for ties.
   
   COUNT = COUNT + INC;
   
   percent = 100*(i/sims);
     
     if sims <= 20
         
         
         percent >= index
        
        disp(['percent_completed = ',num2str(index)]);
        
        index = index + 10;
        
     end
   
      
   i = i + 1 ;
   
end

%Compute P-values and Significance

PVal = COUNT/(sims + 1);

SIG = 1 - PVal;

%Compute space-time mesh

d0 = Dmax/2 ;

inc = d0/(s - 1) ;

D = 0:inc:d0 ;

t0 = Tmax/2 ;

inc = t0/(t - 1) ;

T = 0:inc:t0 ;

%Assign output

results.cdf = CDF;

results.pvals = PVal;

results.dist = D;

results.time = T;