function cvtp_1d ( g_num, it_num, s_num )

%*****************************************************************************80
%
%% cvtp_1d() estimates a CVT in a 1D periodic region.
%
%  Discussion:
%
%    This program approximates a Centroidal Voronoi Tessellation, Periodic
%    (CVTP) on the periodic interval [0,1], using sampling rather than exact
%    geometric calculations.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    20 July 2016
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    integer G_NUM, the number of generators.
%
%    integer IT_NUM, the number of iterations.
%
%    integer S_NUM, the number of samples in [0,1] to use
%    when estimating the Voronoi regions.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'cvtp_1d():\n' );
  fprintf ( 1, '  Estimate a periodic centroidal Voronoi tessellation (CVTP)\n' );
  fprintf ( 1, '  in the 1D periodic unit interval [0,1].\n' );
  fprintf ( 1, '  The CVTP will have G_NUM generators.\n' );
  fprintf ( 1, '  There will be IT_NUM iterations taken.\n' );
  fprintf ( 1, '  Each iteration will consider S_NUM sample points.\n' );

  if ( nargin < 1 )
    g_num = input ( '  Enter number of generators (such as 10): ' );
  elseif ( ischar ( g_num ) )
    g_num = str2num ( g_num );
  end

  if ( nargin < 2 ) 
    it_num = input ( '  Enter number of iterations (such as 10): ' );
  elseif ( ischar ( it_num ) )
    it_num = str2num ( it_num );
  end

  if ( nargin < 3 ) 
    s_num = input ( '  Enter number of sample points (such as 1000): ' );
  elseif ( ischar ( s_num ) )
    s_num = str2num ( s_num );
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of generators is %d\n', g_num );
  fprintf ( 1, '  Number of iterations is %d\n', it_num );
  fprintf ( 1, '  Number of samples is %d\n', s_num );
%
%  Initialize the generators.
%  For convenience, make sure they are sorted.
%
  if ( true )
    fprintf ( 1, '\n' );
    fprintf ( 1, '  Using random initial generators\n' );
    g(1:g_num,1) = rand ( g_num, 1 );
    g = sort ( g );
  else
    fprintf ( 1, '\n' );
    fprintf ( 1, '  Using generators in tight range [0.01,0.02]\n' );
    g(1:g_num,1) = linspace ( 0.01, 0.02, g_num );
  end
%
%  Print the initial generators.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Initial generators:\n' );
  fprintf ( 1, '\n' );
  for k = 1 : g_num
    fprintf ( 1, '  %4d  %g\n', k - 1, g(k,1) );
  end
%
%  Initialize the plotting arrays.
%
  g_plot = zeros ( g_num, it_num + 1 );
  e_plot = nan ( it_num, 1 );
  m_plot = nan ( it_num, 1 );
%
%  Carry out the iterations.
%
  for it = 1 : it_num

    g_plot ( 1:g_num, it ) = g(1:g_num,1);

    s = rand ( s_num, 1 );
%
%  For each sample point, determine the nearest generator.
%  Here, "nearest" is in the periodic sense.  
%  In fact, we essentially start with a sample point S, and have
%  to consider which of three points S, S+1 and S-1 is nearest
%  to some generator G.
%
    g_new = zeros ( g_num, 1 );
    e_new = zeros ( g_num, 1 );
    w_new = zeros ( g_num, 1 );

    for i = 1 : s_num

      d = Inf;
      s2 = 0;
      k = 0;
      for j = 1 : g_num
        d1 = abs ( s(i) - g(j) );
        if ( d1 < d )
          d = d1;
          k = j;
          s2 = s(i);
        end
        d2 = abs ( s(i) + 1.0 - g(j) );
        if ( d2 < d )
          d = d2;
          k = j;
          s2 = s(i) + 1.0;
        end
        d3 = abs ( s(i) - 1.0 - g(j) );
        if ( d3 < d )
          d = d3;
          k = j;
          s2 = s(i) - 1.0;
        end
      end

      g_new(k) = g_new(k) + s2;
      w_new(k) = w_new(k) + 1;
      e_new(k) = e_new(k) + d^2;

    end
%
%  Most new G's will be the average of their sample points.
%  Some G's might have no hits, so they keep the old values.
%
    i1 = find ( w_new(1:g_num) == 0 );
    g_new(i1) = g(i1);
    i2 = find ( w_new(1:g_num) ~= 0 );
    g_new(i2) = g_new(i2) ./ w_new(i2);
%
%  Have to do a MOD so G's don't go out of the region.
%
    g_new = mod ( g_new, 1.0 );
%
%  Sum the total energy for this iteration.
%
    e_plot(it) = sum ( e_new(1:g_num) ) / s_num;
%
%  Compute the generator motion for this iteration.
%
    g_new = g_new(:);
    g = g(:);

    t = 0.0;
    t2 = 0.0;
    for k = 1 : g_num
      t =          abs ( g_new(k) - g(k) );
      t = min ( t, abs ( g_new(k) - g(k) + 1.0 ) );
      t = min ( t, abs ( g_new(k) - g(k) - 1.0 ) );
      t2 = t2 + t * t;
    end
    m_plot(it,1) = t2 / g_num;
%
%  Update the generators.
%
    g = g_new;

  end

  g_plot ( 1:g_num, it_num + 1 ) = g_new(1:g_num,1);
%
%  Print the current generators.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Final generators:\n' );
  fprintf ( 1, '\n' );
  for k = 1 : g_num
    fprintf ( 1, '  %4d  %g\n', k - 1, g(k,1) );
  end
%
%  Plot the energy over iterations.
% 
  figure ( 1 );
  plot ( 1:it_num, log ( e_plot ), 'm-*' );
  title ( 'Log (Energy)' );
  xlabel ( 'Step' );
  ylabel ( 'Energy' );
  grid ( 'on' );
  print ( '-dpng', 'cvtp_1d_energy.png' );
%
%  Plot the generator motion over iterations.
%
  figure ( 2 );
  plot ( 1:it_num, log ( m_plot ), 'm-*' );
  title ( 'Log (Average generator motion)' );
  xlabel ( 'Step' );
  ylabel ( 'Motion' );
  grid ('on' );
  print ( '-dpng', 'cvtp_1d_motion.png' );
%
%  Plot the generator locations over iterations.
%
  figure ( 3 );
  y = ( 0 : it_num );
  hold ( 'on' );
  for k = 1 : g_num
    plot ( g_plot(k,1:it_num+1), y );
  end
  grid ( 'on' );
  hold ( 'off' );
  title ( 'Generator evolution' );
  xlabel ( 'Generator positions' );
  ylabel ( 'Iterations' ); 
  print ( '-dpng', 'cvtp_1d_generators.png' );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'cvtp_1d():\n' );
  fprintf ( 1, '  Normal end of execution.\n' );

  return
end

