function OUT = dist_wts(L,varargin)

%DIST_WTS computes a variety of weight matrices as functions of 
%         pairwise distances between centroids.

%Written by: TONY E. SMITH, 3/28/98

%INPUTS: L = (Nx2) coordinate matrix (Xi,Yi),i=1:N
%        'varargin'specifies the weight matrix desired. In the 
%        following options, the equality 'varargin = a,b,c' means
%        means that the arguments in dist_wts are (L,a,b,c).
%   
%         varargin = 1  yields the nearest-neighbor matrix, W, with 
%                       Wij > 0 denoting that j is a nearest neighbor
%                       of i. Rows are normalized, so that there are
%                       k nearest neighbors, each will have value 1/k)
%
%         varargin = 2  yields a SYMMETRIC nearest-neighbor matrix, W,
%                       Wij > 0 iff either i or j is a nearest neighbor 
%                       of the other.
%
%         varargin = 3,d,s  yields a matrix indicating all neighbors, j, 
%                           within distance d of i. To include i in the 
%                           count, write s=1 and otherwise, write s=0.
%                           Here, Wij = 1/k if there are k such neighbors.
%
%         varargin = 4,a,s  yields an exponential gravity-model form
%
%                           Wij = exp(-a d(ij))/SUMk(exp(-a d(ik)) 
%
%                    where dij is FRACTIONAL distance with dij = 1 for the
%                    maximum pairwise distance. To include the value Wii 
%                    write s=1, and otherwise write s=0.
%
%
%OUTPUTS: OUT{1} = (NxN) weight matrix (in sparse form) for the option 
%                   specified.
%         OUT{2} = (Nx2) matrix of (ordered)nonzero point pairs
%                   [useful for certain types of analyses]

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

D = distance(L);  %Make distance matrix

u = ones(N,1); %Unit column vector

W = zeros(N) ;

if varargin{1} == 1
   
   D = D + diag(D*u); %Eliminate diagonal elements as n-neighbors.   
    
   D = D - min(D')'*u' ; %Substract Row min from each element(nearest
                         %neighbors are now exactly the zero elements).
   
   D = spones(sparse(D + eye(N))); %Converts all nonzero elements to ones
   										  %(adding Identity matrix preserves the
                                   %dimension of the sparse matrix). 
   
   D = u*u' - D ; %Leaves ones in exactly the n-neighbor positions.
        
   D = inv(diag(D*u)) * D ; %Row normalizes D
   
elseif varargin{1} == 2
   
   D = D + diag(D*u); %Eliminate diagonal elements as n-neighbors.   
    
   D = D - min(D')'*u' ; %Substract Row min from each element(nearest
                         %neighbors are now exactly the zero elements).
   
   D = spones(sparse(D + eye(N))); %Converts all nonzero elements to ones
   										  %(adding Identity matrix preserves the
                                   %dimension of the sparse matrix). 
   
   D = u*u' - D ; %Leaves ones in exactly the n-neighbor positions.
   
   D = max(D,D') ; %SYMMETRIZES the matrix D so Dij = 1 iff either i or j
                   %is a nearest neighbor of the other.
        
   D = inv(diag(D*u)) * D ; %Row normalizes D
   

elseif varargin{1} == 3
   
   d = varargin{2} ;
   
   s = varargin{3} ;
   
   if s == 0
      
      D = D + diag(D*u); %Eliminate diagonal elements as d-neighbors.
      
   end
   
   D = max(D - d*u*u',zeros(N)) ; %sets all desired points to zero.
   
   D = spones(sparse(D)); %sets all undesired points to one.
   
   D = u*u' - D ; %Leaves ones in exactly the d-neighbor positions.
   
   for i = 1:N
      
      if max(D(i,:)) > 0
         
         D(i,:)  = D(i,:)/(D(i,:)*u) ; %Normalize all nonzero rows
         
      end
      
   end
   
else % Assumed here that varargin{1}= 3
   
   a = varargin{2} ;
   
   s = varargin{3} ;
   
   d = max(max(D)) ;
   
   D = exp(-(a/d)*D) ;
   
   if s == 0
      
      D = D - diag(diag(D)) ;
      
   end
   
   D = inv(diag(D*u)) * D ; %Row normalizes D
   
end

if varargin{1} == 4
   
   OUT{1} = D ;
   
   OUT{2} = [] ;
   
else
   
   OUT{1} = sparse(D) ;
   
   [x,y,z] = find(D) ;
   
   DD(:,1) = x ;
   DD(:,2) = y ;
   
   OUT{2} = sortrows(DD);
   
end