function [averages, bincentres] = variogram1dvaribins(xdata, values, perbin)
% VARIOGRAM1DVARIBINS calculate and display experimental and modeled variogram
%     [averages, domains] = variogram1dvaribins(xdata, values, perbin)
% 
%     Calculate and display the experimental variogram of a 1d dataset, given
%     x-coordinates in xdata and in values the values at those
%     coordinates. The number of bins is variable; perbin determins 
%     the number of pairs per bin. The result is given as the
%     average covariances for each bin. 
% 
%     Yannick Kremer 2012



if size(xdata, 1) == 1, xdata = xdata'; end
if size(values, 1) == 1, values = values'; end


% Size above which a waitbar is displayed, to show progress.
slowlength = 1000;
showwaitbar = length(xdata) > slowlength;



if showwaitbar, 
    w=waitbar(0, 'this can take a while');
end

% Fill a square matrix with the values of xdata, and subtract its transpose, this
% is a fast (vectorized) way to calculate the differences in the X direction between all
% the data points. 





% Calculate the real distances using the x and y distances, using
% pythagoras. Distances contains the distances between datapoints.
% distances(1,6) for example contains the distance between the first and
% the sixth element. It is symmetric so distancess(6,1) contains the same
% distance

XCUBE = repmat(xdata, 1, length(xdata));
DISTANCES = abs( XCUBE - XCUBE' );

ifwa(0.6);

% Lowerhalf is used to filter the distances and vdiff2 matrix
LOWERHALF = logical( tril( ones(length(xdata)), -1) );
% lowerhalf now looks like this:
% 1 0 0
% 1 1 0
% 1 1 1


% The distances matrix is symmetric, the upper half above the diagonal is
% the same as the half below it. All values on the diagonal are zero, since
% this is the distance between an element and itself, these values need to
% be removed from the dataset.
DISTANCES = DISTANCES( LOWERHALF );

% Next calculate the variation between points in a similar fashion and
% store the result in vdiff2.
VCUBE = repmat(values, 1, length(values) );
ifwa(0.8);

VDIFF2 = (VCUBE - VCUBE').^2;
% figure(3), imagesc(vdiff2 .* lowerhalf)


VDIFF2 =  VDIFF2( LOWERHALF ) ;
 
% Calculate the sizes of the histogram based on the size of the data and
% the numberofdomains
% mindist = min( DISTANCES );
% maxdist = 0.8*max( DISTANCES );
% distrange = maxdist - mindist;


% Define the domains for the histogram
% domains =( 0 :  binsize : lags*binsize)';

ifwa(0.8);

% pre-allocate memory for the averages vector



maxbins = floor( 0.6*length(DISTANCES) / perbin );

averages = zeros(maxbins, 1);
bincentres = zeros(maxbins, 1);

[DISTANCES, i] = sort( DISTANCES);
VDIFF2 = VDIFF2( i );

for idx = 1 : maxbins,
    averages( idx ) = 0.5 * mean( VDIFF2( (1+(idx-1) * perbin):(idx * perbin) ) ) ;
    bincentres( idx ) = mean( DISTANCES( (1+(idx-1) * perbin):(idx * perbin) ) ) ;
end

% Save (half) the average variance for each histogram class in averages.
% for idx = 1: length(domains)-1,
%     %if length(VDIFF2( DISTANCES> domains(idx) & DISTANCES < domains(idx + 1) ) ) >0 ,
%         averages( idx ) = mean( VDIFF2( DISTANCES >= (domains(idx) - overlap) & DISTANCES < (domains(idx + 1) + overlap) ) )/2;
%     %else
% %         extend = 1.1;
% %         while length(VDIFF2( DISTANCES> domains(idx)/extend & DISTANCES < domains(idx + 1)*extend ) ) <5
% %             extend = extend + 0.1;
% %         end
% %         averages( idx ) = mean( VDIFF2( DISTANCES> domains(idx)/extend & DISTANCES < domains(idx + 1)*extend ) )/2;
% %     end
% end

if showwaitbar, close(w), end;


% histogram = domains(1:end-1) + binsize/2;
% fx = histogram( isfinite( averages ));
% fy = averages( isfinite( averages ));

% Generate data to plot the fitted polynomial (=modeled variogram)


% figure,

variogramplot( bincentres, averages, [], 'b+', false )

function ifwa(position)
if showwaitbar, 
    waitbar(position, w);
end
end

end