function W = dist_wts(L,info)

%DIST_WTS computes a variety of weight matrices as functions of 
%         pairwise distances between centroids.All matrices are 
%         row normalized.

%Written by: TONY E. SMITH, 3/28/98 [modified 6/1/09]

%INPUTS: L = (Nx2) coordinate matrix (Xi,Yi),i=1:N
%
%        info = structure of weight matrix specifications.
%
%             info.type: specifies the TYPE of weight matrix desired. 
%                        
%             info.type = 1  yields the nearest-neighbor matrix, W, with 
%                       Wij > 0 denoting that j is a nearest neighbor
%                       of i. 
%
%             info.type = 2  yields a SYMMETRIC nearest-neighbor matrix, W,
%                        Wij > 0 iff either i or j is a nearest neighbor 
%                        of the other.
%
%             info.type = [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.
%                              
%             info.type = [4,a,s]  yields an exponential gravity-model form
%
%                              Wij = exp(-a d(ij)) 
%
%                              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.
%
%             info.type = [5,a,s]   yields an power gravity-model form
%
%                              Wij = d(ij)^(-a) 
%
%                              where dij is FRACTIONAL distance with dij = 1 
%                              for the maximum pairwise distance. 
%                              [NOTE: You will get a WARNING if there is at 
%                              least more than one point has the same location.
%                              To include the value Wii write s=1, and otherwise 
%                              write s=0.]
%   
%              info.type = [6,k,b,s]  yields a double-power form
%
%                              Wij = [1 - (dij/b)^k]^k , dij <= b
%                                  = 0  ,  dij > b
%
%                              where b is a bandwidth for the function, and
%                              where k is a positive integer power (usually
%                              k = 2 or higher).To include the value Wii 
%                              write s=1, and otherwise write s=0
%
%
%          info.norm: (optional) normalization. Default is NO normalization
%
%              info.norm = 1  if row normalization is desired
%
%              info.norm = 2  if max-eigenvalue normalization is desired
%
%
% OUTPUT: W = (NxN) weight matrix (in sparse form) for the option 
%                   specified.


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

D = distance_mat(L);  %Make distance matrix

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

W = sparse(N,N) ;

if info.type == 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). 
   
   W = sparse(u*u' - D) ; %Leaves ones in exactly the n-neighbor positions.        
    
elseif info.type == 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 = sparse(u*u' - D) ; %Leaves ones in exactly the n-neighbor positions.
   
   W = max(D,D') ; %SYMMETRIZES the matrix D so Dij = 1 iff either i or j
                   %is a nearest neighbor of the other. 

elseif info.type(1) == 3
    
   if length(info.type) ~= 3
    
    error('MUST HAVE FULL INFORMATION FOR TYPE 3');
    
    end
   
   d = info.type(2) ;
   
   s = info.type(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.
   
   W = sparse(u*u' - D); %Leaves ones in exactly the d-neighbor positions.   
    
   
elseif info.type(1) == 4  
    
    if length(info.type) ~= 3
    
    error('MUST HAVE FULL INFORMATION FOR TYPE 4');
    
    end
   
   a = info.type(2) ;
   
   s = info.type(3) ;
   
   d = max(max(D)) ;
   
   W = exp(-(a/d)*D) ;
   
   if s == 0
      
      W = W - diag(diag(W)) ;
      
   end   
   
     
elseif info.type(1) == 5
    
    if length(info.type) ~= 3
    
    error('MUST HAVE FULL INFORMATION FOR TYPE 5');
    
    end
    
    a = info.type(2) ;
    
    s = info.type(3);

    n = size(D,1);

    D = D/max(max(D)) + eye(n);

    D = D.^(-a) ;  

    if isinf(max(max(D)))

       error('ONE OR MORE DISTANCES IS ZERO');

    end

   if s == 0
      
      W = D - diag(diag(D)) ; % Remove diagonal elements
      
   end   
    
else % assumes info.type(1) == 6
    
    if length(info.type) ~= 4
    
    error('MUST HAVE FULL INFORMATION FOR TYPE 6');
    
    end
    
    k = info.type(2);
    
    b = info.type(3);
    
    s = info.type(4);
    
    d = max(max(D));

    d_pos = min(min(D + d*eye(N)));  %calculates min positive distance  

    if b < d_pos        

        error(['BANDWIDTH MUST EXCEED SMALLEST POSITIVE DISTANCE = ',num2str(d_pos)]);

    end

    U = ones(N); %Unit column vector

    W = (U - (D/b).^k).^k;

    Z = (D <= b);

    W = sparse(W.*Z); %makes all wij's with dij > b equal to zero.
    
    if s == 0
      
      W = W - diag(diag(W)) ; % Remove diagonal elements
      
   end
     
end


% Finally, normalize if desired

if isfield(info,'norm');
       
       if info.norm == 1
           
           W = row_norm(W) ; %Row normalizes D
           
       elseif info.norm == 2 
           
           W = sparse(W);
           
           if isinf(max(max(W)))
               
               error('W has one or more infinite elements');
               
           else % The following version of eigs SILENTLY computes
                % the largest eigenvalue of W
               
               options.disp = 0; 
               
               max_eig = abs(max(eigs(W,1,'lm',options)));
              
               W = W/max_eig;
              
           end
           
       end
       
   end