function energy = cvt_energy ( ndim, n, batch, sample, initialize, ...
  sample_num, 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 MIT license.
%
%  Modified:
%
%    04 December 2004
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    integer NDIM, the spatial dimension.
%
%    integer N, the number of generators.
%
%    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.
%
%    integer SAMPLE, specifies how the sampling is done.
%    -1, 'RAND', using MATLAB's RAND function;
%     0, 'UNIFORM', using a simple uniform RNG;
%     2, 'GRID', points from a grid;
%
%    logical INITIALIZE, is TRUE if the pseudorandom process may need to be
%    reinitialized.
%
%    integer SAMPLE_NUM, the number of sample points.
%
%    real R(NDIM,N), the Voronoi cell generators.
%
%  Output:
%
%    real ENERGY, the estimated CVT energy.
%

%
%  Generate the sampling points S in batches.
%
  have = 0;
  energy = 0.0;
  
  while ( have < sample_num )

    get = min ( sample_num - have, batch );

    s = cvt_sample ( ndim, sample_num, get, sample, initialize );
  
    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
