function [ energy, seed ] = cvt_energy ( ndim, n, batch, sample, initialize, ...
sample_num, seed, r )
%*****************************************************************************80
%
%% CVT_ENERGY computes the CVT energy of a dataset.
%
% Discussion:
%
% For a given number of generators, a CVT is a minimizer (or at least
% a local minimizer) of the CVT energy. During a CVT iteration,
% it should generally be the case that the CVT energy decreases from
% step to step, and that perturbations or adjustments of an
% approximate CVT will almost always have higher CVT energy.
%
% Licensing:
%
% This code is distributed under the GNU LGPL license.
%
% Modified:
%
% 04 December 2004
%
% Author:
%
% John Burkardt
%
% Parameters:
%
% Input, integer NDIM, the spatial dimension.
%
% Input, integer N, the number of generators.
%
% Input, integer BATCH, sets the maximum number of sample points
% generated at one time. It is inefficient to generate the sample
% points 1 at a time, but memory intensive to generate them all
% at once. You might set BATCH to min ( SAMPLE_NUM, 10000 ), for instance.
% BATCH must be at least 1.
%
% Input, integer SAMPLE, specifies how the sampling is done.
% -1, 'RAND', using MATLAB's RAND function;
% 0, 'UNIFORM', using a simple uniform RNG;
% 1, 'HALTON', from a Halton sequence;
% 2, 'GRID', points from a grid;
%
% Input, logical INITIALIZE, is TRUE if the SEED must be reset to SEED_INIT
% before computation. Also, the pseudorandom process may need to be
% reinitialized.
%
% Input, integer SAMPLE_NUM, the number of sample points.
%
% Input, integer SEED, the random number seed.
%
% Input, real R(NDIM,N), the Voronoi cell generators.
%
% Output, real ( kind = 8 ) ENERGY, the estimated CVT energy.
%
% Output, integer SEED, an updated seed for the random number generator.
%
%
% Generate the sampling points S in batches.
%
have = 0;
energy = 0.0;
while ( have < sample_num )
get = min ( sample_num - have, batch );
[ s, seed ] = cvt_sample ( ndim, sample_num, get, sample, initialize, ...
seed );
initialize = 0;
have = have + get;
%
% Find the index N of the nearest cell generator to each sample point S.
%
nearest = find_closest ( ndim, n, get, s, r );
%
% Add S to the centroid associated with generator N.
%
for j = 1 : get
energy = energy + sum ( ( r(1:ndim,nearest(j)) - s(1:ndim,j) ).^2 );
end
end
energy = energy / sample_num;
return
end