function g = cvt_2d_lumping ( n, it_num, s_num, mu_fun, init )

%*****************************************************************************80
%
%% CVT_2D_LUMPING carries out the Lloyd algorithm with lumping in 2D.
%
%  Discussion:
%
%    This program is an extension of the CVT_1D_LUMPING program to the 2D case.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 May 2019
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of generators.
%    N must be at least 3.
%
%    Input, integer IT_NUM, the number of CVT iterations.
%
%    Input, integer S_NUM, the number of samples in [-1,1] to use
%    when estimating the Voronoi regions.
%
%    Input, function value = MU_FUN ( X, Y ), the name of the M file 
%    which accepts points X, Y and returns a real vector
%    VALUE of density function values.
%
%    Input, integer INIT, the initialization option.
%    *1, random sorted values;
%
%    Output, real G(N,2), the locations of the generators after the
%    iterations have been carried out.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'CVT_2D_LUMPING\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Use lumping to implement a version of Lloyd''s algorithm\n' );
  fprintf ( 1, '  in the 2D square [-1,1]x[-1,+1].\n' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of generators is %d\n', n );
  fprintf ( 1, '  Number of iterations is %d\n', it_num );
  fprintf ( 1, '  Sample rate in 1D    is %d\n', s_num );
%
%  Fail if N < 3 (because we need to use the Delaunay triangulation
%  to do searches, and we're too lazy to write special cases for N = 1, 2.
%
  if ( n < 3 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'CVT_2D_LUMPING - Fatal error!\n' );
    fprintf ( 1, '  N must be at least 3, but input N = %d\n', n );
    error ( 'CVT_2D_LUMPING - Fatal error!' );
  end
%
%  Initialize the generators in [-1,+1]x[-1,+1].
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Initial generators (uniform random):\n' );
  g(1:n,1:2) = 2.0 * rand ( n, 2 ) - 1.0;
%
%  Print the initial generators.
%
  fprintf ( 1, '\n' );
  for k = 1 : n
    fprintf ( 1, '  %4d  %8.4f  %8.4f\n', k, g(k,1), g(k,2) );
  end
%
%  Initialize some arrays.
%
  g_new = zeros ( n, 2 );
  step = ( 1 : it_num )';
  e = nan ( it_num, 1 );
  gm = nan ( it_num, 1 );
%
%  SX, SY contains the X and Y coordinates of the sample points in [-1,1]x[-1,+1].  
%  We choose to use a uniform mesh, and we need it to avoid the boundary
%  of the region.
%
  x_min = - 1.0 + eps;
  x_max = + 1.0 - eps;

  s_1d = linspace ( x_min, x_max, s_num );
  [ sx_mat, sy_mat ] = meshgrid ( s_1d, s_1d );
  sx_vec = reshape ( sx_mat, s_num * s_num, 1 );
  sy_vec = reshape ( sy_mat, s_num * s_num, 1 );
  s = [ sx_vec, sy_vec ];
%
%  R is the density function RHO(X,Y).
%
%  We will NOT assume it is normalized.
%
%  We will assume it is the cube of some other density function MU(X,Y).
%
%  In this way, we expect the asymptotic distribution of CVT generators
%  to follow the distribution MU(X,Y).
%
  mu_mat = mu_fun ( sx_mat, sy_mat );
%
%  MU blows up.  That's why we pulled the sample points in, away from
%  the boundary.  But it still blows up.  So try truncating it, say at 10.
%
  mu_mat = min ( mu_mat, 10 );
%
%  Rho was the cube of MU in 1D, but in 2d this becomes a fourth power.
%
  r_mat = mu_mat .^ 4;
  r = reshape ( r_mat, s_num * s_num, 1 );
%  
%  Contour plot of MU(X,Y).
%
  if ( true )
    figure ( 1 )
    surf ( sx_mat, sy_mat, mu_mat, 'EdgeColor', 'none' );
    title ( 'Density \mu(x)', 'Fontsize', 16 );
    filename = 'density_mu.png';
    print ( '-dpng', filename );
    fprintf ( 1, '\n' );
    fprintf ( 1, '  Graphics saved as "%s"\n', filename );
  end
%  
%  Contour plot of RHO(X,Y).
%
  if ( true )
    figure ( 2 )
    surf ( sx_mat, sy_mat, r_mat, 'EdgeColor', 'none' );
    title ( 'Density \rho(x):', 'Fontsize', 16 );
    filename = 'density_rho.png';
    print ( '-dpng', filename );
    fprintf ( 1, '  Graphics saved as "%s"\n', filename );
  end
%
%  Now carry out the iteration.
%
  for it = 1 : it_num
%
%  Compute the Delaunay triangle information T for the current nodes.
%
    t = delaunay ( g(:,1), g(:,2) );
%
%  For each sample point, find K, the index of the nearest generator.
%  
    k = dsearchn ( g, t, s );
%
%  M(I,1) contains the total mass of sample points closest to generator I.
%
    m(:,1) = accumarray ( k, r );
%
%  G_NEW(I,1:2) is the mass-weighted average of the sample points
%  closest to generator I.
%
    g_new(:,1) = accumarray ( k, r(:) .* s(:,1) ) ./ m(:,1);
    g_new(:,2) = accumarray ( k, r(:) .* s(:,2) ) ./ m(:,1);
%
%  The energy is the mass weighted distance.
%
    e(it,1) = sum ( ( r(:) .* s(:,1) - g(k(:,1),1) ) .^ 2 ...
                  + ( r(:) .* s(:,2) - g(k(:,1),2) ) .^ 2 ) / s_num;
%
%  Compute the generator motion.
%
    gm(it,1) = sum ( ( g_new(:,1) - g(:,1) ) .^ 2 ...
                   + ( g_new(:,2) - g(:,2) ) .^ 2 ) / n;
%
%  Update the generators.
%
    g = g_new;

  end
%
%  Display the energy.
%
  figure ( 3 )
  plot ( step, log ( e ), 'm-*' )
  title ( 'Log (Energy)', 'Fontsize', 16 )
  xlabel ( 'Iteration', 'Fontsize', 16 )
  ylabel ( 'Energy', 'Fontsize', 16 )
  grid
  if ( init == 1 )
    filename = sprintf ( 'energy_%03d_r.png', n );
  elseif ( init == 2 )
    filename = sprintf ( 'energy_%03d_c.png', n );
  elseif ( init == 3 )
    filename = sprintf ( 'energy_%03d_u.png', n );
  end
  print ( '-dpng', filename );
  fprintf ( 1, '  Graphics saved as "%s"\n', filename );
%
%  Display the generator motion.
%
  figure ( 4 )
  plot ( step, log ( gm ), 'm-*' )
  title ( 'Log (Average generator motion)', 'Fontsize', 16 )
  xlabel ( 'Iteration', 'Fontsize', 16 )
  ylabel ( 'Motion', 'Fontsize', 16 )
  grid
  if ( init == 1 )
    filename = sprintf ( 'motion_%03d_r.png', n );
  elseif ( init == 2 )
    filename = sprintf ( 'motion_%03d_c.png', n );
  elseif ( init == 3 )
    filename = sprintf ( 'motion_%03d_u.png', n );
  end
  print ( '-dpng', filename );
  fprintf ( 1, '  Graphics saved as "%s"\n', filename );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'CVT_2D_LUMPING\n' );
  fprintf ( 1, '  Normal end of execution.\n' );

  return
end
