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

%*****************************************************************************80
%
%% CVT_3D_LUMPING carries out the Lloyd algorithm with lumping in 3D.
%
%  Discussion:
%
%    This program is an extension of the CVT_1D_LUMPING program to the 3D case.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 October 2016
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of generators.
%    N must be at least 4.
%
%    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.  This number is cubed to
%    get the total number of 3D sample points.
%
%    Input, function value = MU_FUN ( X, Y, Z ), the name of the M file 
%    which accepts points (X,Y,Z) and returns a real vector
%    VALUE of density function values.
%
%    Input, integer INIT, the initialization option.
%    *1, random sorted values;
%
%    Output, real G(N,3), the locations of the generators after the
%    iterations have been carried out.
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'CVT_3D_LUMPING\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Use lumping to implement a version of Lloyd''s algorithm\n' );
  fprintf ( 1, '  in the 3D square [-1,1]x[-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 < 4.
%
  if ( n < 4 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'CVT_3D_LUMPING - Fatal error!\n' );
    fprintf ( 1, '  N must be at least 4, but input N = %d\n', n );
    error ( 'CVT_3D_LUMPING - Fatal error!' );
  end
%
%  Initialize the generators.
%
  g(1:n,1:3) = 2.0 * rand ( n, 3 ) - 1.0;
%
%  Print the initial generators.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Initial generators (uniform random):\n' );
  fprintf ( 1, '\n' );

  for k = 1 : n
    fprintf ( 1, '  %4d  %8.4f  %8.4f  %8.4f\n', k, g(k,1), g(k,2), g(k,3) );
  end
%
%  Initialize some arrays.
%
  g_new = zeros ( n, 3 );
  step = ( 1 : it_num )';
  e = nan ( it_num, 1 );
  gm = nan ( it_num, 1 );
%
%  SX, SY, SZ contains the XYZ coordinates of the sample points.  
%  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 );
%
%  Make separate vectors of X, Y and Z coordinates for the 3D mesh,
%  and a combined array S of dimension [SIZE Y * SIZE X * SIZE Z]
%  for some godforsaken reason which doesn't matter to us at the
%  moment since all three coordinates have the same range and same
%  number of points...but any variation from this will cause a
%  disaster.
%
  [ sx_mat, sy_mat, sz_mat ] = meshgrid ( s_1d, s_1d, s_1d );
  sx_vec = reshape ( sx_mat, s_num * s_num * s_num, 1 );
  sy_vec = reshape ( sy_mat, s_num * s_num * s_num, 1 );
  sz_vec = reshape ( sz_mat, s_num * s_num * s_num, 1 );
  s = [ sx_vec, sy_vec, sz_vec ];
%
%  R is the density function RHO(X,Y,Z).
%
%  We will NOT assume it is normalized.
%
%  We will assume it is the cube of some other density function MU(X,Y,Z).
%
%  In this way, we expect the asymptotic distribution of CVT generators
%  to follow the distribution MU(X,Y,Z).
%
  mu_mat = mu_fun ( sx_mat, sy_mat, sz_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 3d this becomes a fifth power.
%
  r_mat = mu_mat .^ 5;
  r = reshape ( r_mat, s_num * s_num * s_num, 1 );
%
%  Now carry out the iteration.
%
  for it = 1 : it_num
%
%  Compute the Delaunay triangle information T for the current nodes.
%  G must be an MxN array of M points in N dimensional space.
%
    t = delaunayn ( g );
%
%  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:3) 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);
    g_new(:,3) = accumarray ( k, r(:) .* s(:,3) ) ./ 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 ...
                  + ( r(:) .* s(:,3) - g(k(:,1),3) ) .^ 2 ) / s_num;
%
%  Compute the generator motion.
%
    gm(it,1) = sum ( ( g_new(:,1) - g(:,1) ) .^ 2 ...
                   + ( g_new(:,2) - g(:,2) ) .^ 2 ...
                   + ( g_new(:,3) - g(:,3) ) .^ 2 ) / n;
%
%  Update the generators.
%
    g = g_new;

  end
%
%  Print the Final generators.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Final generators:\n' );
  fprintf ( 1, '\n' );

  for k = 1 : n
    fprintf ( 1, '  %4d  %8.4f  %8.4f  %8.4f\n', k, g(k,1), g(k,2), g(k,3) );
  end
%
%  Display the energy.
%
  figure ( 3 )
  plot ( step, log ( e ), 'm-*' )
  title ( 'Log (Energy)' )
  xlabel ( 'Iteration' )
  ylabel ( 'Energy' )
  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, '  Saved graphics file "%s"\n', filename );
%
%  Display the generator motion.
%
  figure ( 4 )
  plot ( step, log ( gm ), 'm-*' )
  title ( 'Log (Average generator motion)' )
  xlabel ( 'Iteration' )
  ylabel ( 'Motion' )
  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, '  Saved graphics file "%s"\n', filename );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'CVT_3D_LUMPING\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end
