function [averages, bincentres] = variogram1dmanualbins(xdata, values, bins, showpaircount, style)
% VARIOGRAM calculate and display experimental and modeled variogram
%     [averages, domains, p] = variogram1dmanualbins(xdata, values, bins, showpaircount)
% 
%     Calculate and display the experimental variogram of a 2d dataset, given
%     x-coordinates in xdata, y-coordinates in ydata, and the values at those
%     coordinates. Numberofclasses defines the number of classes for which the
%     variogram (which is a histogram) is displayed. The result is given as the
%     average covariances for each domain in bins. p contains the
%     polynomials fitted to the variogram. If the optional argumant
%     showVariogram is false, no variogram is displayed, default = true.
% 
%     Yannick Kremer 2012

if nargin == 3,
    showpaircount = false;
    style = 'b+';
end

if nargin == 4,
    style = 'b+';
end

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

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 ) ;
 
ifwa(0.8);

% pre-allocate memory for the  vectors

averages = zeros(length(bins)-1, 1);
paircount = zeros(length(bins)-1, 1);

% Save (half) the average variance for each histogram class in averages.
for idx = 1: length(bins)-1,
    averages( idx )   = mean( VDIFF2( DISTANCES >= (bins(idx)) & DISTANCES < (bins(idx + 1)) ) )/2;
    paircount( idx )  = sum( DISTANCES >= (bins(idx)) & DISTANCES < (bins(idx + 1)) );
end

if showwaitbar, close(w), end;

bincentres = ( bins(1:end-1)+bins(2:end) ) .* 0.5;

%  plot results
variogramplot( bincentres, averages, paircount, style, showpaircount )

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

end