function [ generator_new, change_l2 ] = cvtp_iteration ( m, n, a, b, ...
  generator, modular, sample_num_cvt )

%*****************************************************************************80
%
%% CVTP_ITERATION takes one step of the CVT iteration.
%
%  Discussion:
%
%    The routine is given a set of points, called "generators", which
%    define a tessellation of the region into Voronoi cells.  Each point
%    defines a cell.  Each cell, in turn, has a centroid, but it is
%    unlikely that the centroid and the generator coincide.
%
%    Each time this CVT iteration is carried out, an attempt is made
%    to modify the generators in such a way that they are closer and
%    closer to being the centroids of the Voronoi cells they generate.
%
%    A large number of sample points are generated, and the nearest generator
%    is determined.  A count is kept of how many points were nearest to each
%    generator.  Once the sampling is completed, the location of all the
%    generators is adjusted.  This step should decrease the discrepancy
%    between the generators and the centroids.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license. 
%
%  Modified:
%
%    01 August 2016
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, the spatial dimension.
%
%    Input, integer N, the number of Voronoi cells.
%
%    Input, real A(M), B(M), the lower and upper region limits.
%
%    Input, real GENERATOR(M,N), the Voronoi cell generators.
%
%    Input, logical MODULAR, is TRUE if modular arithmetic is to be used.
%
%    Input, integer ( kind = 4 ) SAMPLE_NUM_CVT, the number of sample points.
%
%    Output, real GENERATOR_NEW(M,N), the updated generators.
%
%    Output, real ( kind = 8 ) CHANGE_L2, the sum of the L2 norms of the
%    change in each generator's position.
%
  generator_new = zeros ( m, n );
  count = zeros ( n, 1 );

  for j = 1 : sample_num_cvt
%
%  Generate a sampling point X.
%
    x = cvtp_region_sampler ( m, 1, a, b );
%
%  Find the cell generator G that is nearest to XP, some
%  periodic image of X.
%
    [ nearest, xp ] = cvtp_find_closest ( m, n, a, b, x, generator, modular );
%
%  Add XP to the averaging data for GENERATOR(*,NEAREST).
%
    generator_new(1:m,nearest) = generator_new(1:m,nearest) + xp(1:m);

    count(nearest) = count(nearest) + 1;

  end
%
%  Average the sample points to get new generators.
%
  for j = 1 : n
    if ( count(j) ~= 0 )
      generator_new(1:m,j) = generator_new(1:m,j) / count(j);
    end
  end
%
%  While some generators may have drifted out of the main region,
%  our "change" calculation doesn't mind that.
%  
  change_l2 = norm ( generator_new - generator, 'fro' );
%
%  It's possible that a generator would go outside the box.
%  Use modular arithmetic to fix that.
%
  if ( modular )

    for j = 1 : n
      for i = 1 : m

        if ( generator_new(i,j) < a(i) )
          generator_new(i,j) = generator_new(i,j) + b(i) - a(i);
        elseif ( b(i) < generator_new(i,j) )
          generator_new(i,j) = generator_new(i,j) - ( b(i) - a(i) );
        end

      end
    end

  end

  return
end

