function OUT = geo_krige(y0,X0,L0,X1,L1,h)

%GEO_KRIGE uses the geoststistical regression procedure in GEO_REGR
%to estimate beta and covariogram parameters, and then uses these to
%construct kriging predictions at a set of locations, L.

%Written by: TONY E. SMITH, 2/30/98

%Functions called: distance_mat.m, var_spher.m, spher_mat.m, spher.m
%
%INPUTS: (i)   y0  = (n,1) vector of dependent variable values
%                   at n sample locations.
%        (ii)  X0 = (n,k) data matrix (Xi, i = 1:k) at sample
%                   locations [such as polynomial functions of 
%                   coordinates].
%        (iii) L0 = (n,2) matrix of sample location coordinates.
%        (iv)  X1  = (N,k) matrix of data {same as X0) for N 
%                   prediction locations.
%        (v)   L1  = (N,2) matrix of prediction location coordinates.
%        (vi)  h  = bandwidth for Kriging prediction.
%
% 					NOTE: To ensure that local covariance matrices are 
%              nonsingular, it is necessary to have at least k + 1 points 
%              in every kriging neighborhood. For added stability, it is
%              appropriate to require that there be at least k + 2 points.
%              So if the neighborhood, Ni, of a point i does not contain 
%              at least k + 2 points, then Ni is defined to be the k + 2 
%              nearest neighbors of i.
%
%OUTPUTS: OUT = (4x1) cell with
%         OUT{1} = [b,P] = (k+1) vectors of regression betas and P-values.
%                  Here b(1) = intercept, and b(2),..,b(k+1) are the slope
%                  coefficients in the same order as the columns of X0.
%         OUT{2} = (3,1) vector of estimated variogram parameters
%                  [r,s,a]
%         OUT{3} = N-vector of kriged values at locations in L
%         OUT{4} = N-vector of standard errors at locations in L

%Initial Settings

[n,k] = size(X0) ;

U = ones(n,1) ;

XX(:,1) = U ;

XX(:,2:k+1) = X0 ;


%Now rescale data to avoid possible instabilities

S = inv(diag(max(abs(XX)))) ;

Z = XX*S ;

%Define Locations and compute distances

M(:,1:2) = L0 ;

D = distance_mat(L0) ;

%Initial OLS estimates and residuals

bb = Z\y0 ;

b = S*bb ; %This factors out the rescaling effect

%EXPLANATION: Z = (XX)*S => bb = inv(Z'*Z)*Z'*y
%                              = inv[S'*(XX)'*(XX)*S]*S'*(XX)'*y
%                              = inv(S)*inv[(XX)'*(XX)]*inv(S')*S'*(XX)*y
%                              = inv(S)*{inv[(XX)'(XX)]*(XX)'*y}
%                              = inv(S)*b
%             => b = S*bb.



e = y0 - XX*b ;

%Construction of initial variogram estimates

M(:,3) = e ;

disp(' ');
disp('Beginning initial covariogram estimation');
disp(' ');

res = var_spher(M) ;

p = res.param ; %parameter vector for variogram

VARIOGRAM_BANDWIDTH = res.band %Print Bandwidth

%Begin Estimation Loop

disp(' ');
disp('Beginning estimation loop');
disp(' ');

del = 1.0 ;

while del > .001
   
   %Save current beta values
   
   b0 = b ;

	%Construction of Covariance Matrix
	
	C = p(2)*U*U' - spher_mat(p,D) + p(3)*eye(n) ;   
     
   %New b estimates and residuals
   
   A = chol(C)';
        
   AA = A\eye(n) ;   
           
   bb = (AA*Z)\(AA*y0) ;
   
   b = S*bb ; %This factors out the rescaling effect   
   
   M(:,3) = y0 - XX*b ;
   
   res = var_spher(M) ;
   
   p = res.param ;   
        
   del = max(abs((b - b0)./b0))       
      
end

% Update Covariance Matrix

C = p(2)*U*U' - spher_mat(p,D) + p(3)*eye(n) ; 


%se = sqrt(diag(inv(XX'*AA'*AA*XX)));

%Rescale to compute standard errors

cov_bb = inv(Z'*AA'*AA*Z);

cov_b = S*cov_bb*S';  %Recall that b = S*bb

se = sqrt(diag(cov_b));

z = abs(b./se);

P = 2*normcdf(-z);

clear('M','A','AA') ;


%Note: b and p are now defined for use in Kriging

%Compute Kriged Values and Standard Errors

disp(' ');
disp('Begin Kriging');
disp(' ');

N = length(L1(:,1)) ;

x(:,1) = ones(N,1) ;

x(:,2:k+1) = X1 ;

KRIGE = zeros(N,1) ;

STD_ERR = zeros(N,1) ;

RES = y0 - XX*b ;

for i = 1:N   
     
   %Identify points within distance h of location i
   
   x1 = L1(i,1)* U - L0(:,1) ;
   
   x2 = L1(i,2)* U - L0(:,2) ;
   
   if L1(i,:) == [460,340]
       
       i = i;
       
   end
   
   DD = sqrt(x1.^2 + x2.^2) ;
   
   pos = max(zeros(n,1),(h*U - DD)) ; 
   
   K = find(pos) ;  %List of nonzero positions in DD
   
   if length(K) < k+2  %Must have at least k+1 pts to ensure full rank.
      
      %find the k+2 closest points
      
      list(:,1) = [1:n]' ;
      
      list(:,2) = DD(:) ;
      
      ngh = sortrows(list,2) ; 
      
      K = ngh(1:k+2,1) ;
      
   end
                
      
   %Compute local covariances  
   
   CC = inv(C(K,K)) ; %inverse of local cov matrix
   
   d = DD(K) ; %vector of distances to i
   
   c = p(2)*U(K) - spher(p,d); %local cov vector
   
   %Correct c at any zero-distance points   
   
   KK = length(K) ;
   
   for j = 1:KK
      
      if d(j) == 0 
         
         c(j) = c(j) + p(3) ;
         
      end
      
   end
   
   KRIGE(i) = x(i,:)*b + c'*CC*RES(K) ;
   
   w = x(i,:)' ;
   
   W = XX(K,:) ;
   
   %Compute inverse, B = inv(W'*CC*W)
   
   ZZ = Z(K,:) ; %same as XX(K,:)*S 
   
   %Note: In the following calculation we employ the Pseudo 
   %      Inverse function, pinv(.), if ZZ'*CC*ZZ is numerically 
   %      less than full rank.
        
   if rank(ZZ'*CC*ZZ) < k+1
      
      B = S*pinv(ZZ'*CC*ZZ)*S ;      
           
   else
            
      B = S*inv(ZZ'*CC*ZZ)*S ;
      
   end   
  
      
   ERR_VAR = p(2) - c'*CC*c + (w - W'*CC*c)'*B*(w - W'*CC*c) ;
   
   STD_ERR(i) = sqrt(ERR_VAR) ;   
        
end


%Prepare Output


OUT{1} = [b,P] ;
OUT{2} = p ;
OUT{3} = KRIGE ;
OUT{4} = STD_ERR;