function [DAT,BW] = var_robust(M,varargin)

%VARIOGRAM computes the ROBUST empirical (semi)variogram for a given
%set of spatial data (as in Cressie, p.75). Bins are determined internally.
%[If you want %to increase bin size, go to line 120 and reset 'binMax'].

%Written by: TONY E. SMITH, 2/17/98
%
%INPUTS:  M = (nx3)-matrix of locations and values,
%              M(i) = (Xi,Yi,Zi), i = 1:n
%   varargin = optional value of bandwidth 
%
%OUTPUTS: DAT = (bX3)-matrix with OUT(i)=(Di,Vi,Ni),i=1:b
%                Di = distance value for bin i
%                Vi = variogram value at distance Di
%                Ni = number of distinct pairs in bin i
%          BW = Bandwidth used

%See also VARIOGRAM_PLOT

n = length(M(:,1)) ; %number of locations

%Compute matrix, D, of pairwise distances

X = M(:,1) ; %column vector

Y = M(:,2) ; %column vector

U = ones(n,1) ; %column vector

XX = X*U' - U*X' ;

YY = Y*U' - U*Y' ;

D = sqrt(XX.^2 + YY.^2) ;

clear('X','Y','XX','YY') 

D = tril(D) ;  %take lower triangular part


%Now reduce to nonzero distances and determine range


[SS(:,1),SS(:,2),SS(:,3)] = find(D) ;

clear('D');

	%The matrix SS has rows (i,j,d) containing the row,
   %column, and values of all positive cells of D 


Dmin = min(min(SS(:,3))) ;

Dmax = max(max(SS(:,3))) ;

%Set BandWidth if not already specified

if (nargin < 2)
   
   if (Dmin < Dmax/2)
   
   d = Dmax/2 ; %maximum relevant distance
   
	else
   
      error('Dmin >= Dmax/2') ; 
      
   end
   
else
         
   d = varargin{:};
   
   if (Dmax < d)  %check to see if d is too large
      
      d = Dmax ; 
      
   end
   
end

BW = d;  %Bandwidth
   
%Now compute (Zi - Zj)^2

Z = M(:,3) ;

S(:,1) = (abs(Z(SS(:,1)) - Z(SS(:,2)))).^(1/2) ; %sqrt of absolute diffs 

S(:,2) = SS(:,3) ; %distances

clear('SS') ;

S = sortrows(S,2);

	%The matrix S now contains (Zi - Zj)^2 values in column 1
	%and distance values Dij in column 2 sorted from smallest
   %to largest
   
%Now eliminate all pairs with distance > d

Rows = find(S(:,2) <= d);

S = S(Rows,:);

%Now relabel

V = S(:,1) ; %sqrt of absolute diffs

D = S(:,2) ; %distances

clear('S') ;

	   
%Now aggregate into bins
   
i = 1 ; %point counter

b = 1 ; %bin counter

c = 0 ; %bin size counter

VV = 0 ; %cumulative values 

d1 = D(1) ; %left end of bin 1

%Choose bin size to be at least 30, and also large
%enough so that there are no more than 50 bins 

Length = length(D) ;

binMax = 100 ;

NUM = max([ceil(Length/binMax),30]);

%Grouping into bins


for i = 1:Length
   
   c = c + 1 ;
   
   VV = VV + V(i) ;
   
   if (c>=NUM)&((i == Length)|(D(i)<D(i+1)))
      
      DAT(b,1) = (D(i) + d1)/2 ; %assign bin distance
      
      DAT(b,2) = (1/2)*((VV/c)^4)/(.457 +(.494/c)) ; % assign bin value
                                                     % Cressie,p.75, exp.(2.4.12)      
      DAT(b,3) = c ; %assign bin count
          
      b = b+1; 
      
      d1 = D(i) ; %Reinitialize left-end bin distance
      
      c = 0 ; 
      
      VV = 0 ;
      
   end
   
      
end

%Make last bin

if (c > 0) %if last bin is not completed
   
   if (c < NUM)&(b > 1)  %if last bin is not filled, merge with previous one 
   
   DAT(b-1,1) = (D(i) + d1)/2 ; %reassign bin distance
   
   DAT(b-1,2) = (VV + DAT(b-1,2)*(2*DAT(b-1,3)))/(2*(c + DAT(b-1,3))) ;    
     
   DAT(b-1,3) = DAT(b-1,3) + c ; %reassign bin count
   
   else   
   
   DAT(b,1) = (D(i) + d1)/2 ; %assign bin distance
      
   DAT(b,2) = VV/(2*c) ; %assign bin value
      
   DAT(b,3) = c ; %assign bin count
   
   end
    
end