function cvtm_1d ( g_num, it_num, s_num )

%*****************************************************************************80
%
%% CVTM_1D estimates a CVT in a 1D mirror-periodic region.
%
%  Discussion:
%
%    This program approximates a Centroidal Voronoi Tessellation, 
%    Mirror-Periodic (CVTM) on the interval [0,1], using sampling 
%    rather than exact geometric calculations.
%
%    Although we claim to use S_NUM sample points, we actually generate
%    3*S_NUM sample points because of mirror images.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 August 2016
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer G_NUM, the number of generators.
%
%    Input, integer IT_NUM, the number of iterations.
%
%    Input, integer S_NUM, the number of samples in [0,1] to use
%    when estimating the Voronoi regions.
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'CVTM_1D\n' );
  fprintf ( 1, '  MATLAB version\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 );
    sa = 0.0 - s;
    sb = 2.0 - s;
%
%  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
%
%  For each "real" sample point, find the closest generator.
%
      d = Inf;
      k = -1;

      for j = 1 : g_num
        dj = abs ( s(i) - g(j) );
        if ( dj < d )
          d = dj;
          k = j;
        end
      end

      da = Inf;
      ka = -1;

      for j = 1 : g_num
        daj = abs ( sa(i) - g(j) );
        if ( daj < da )
          da = daj;
          ka = j;
        end
      end

      db = Inf;
      kb = -1;

      for j = 1 : g_num
        dbj = abs ( sb(i) - g(j) );
        if ( dbj < db )
          db = dbj;
          kb = j;
        end
      end
%
%  If KA = K, then replace S by 0.0.
%  If KB = K, then replace S by 1.0.
%
      if ( ka == k )
        s(i) = 0.0;
      elseif ( kb == k )
        s(i) = 1.0;
      end

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

    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);
%
%  Compute generator motion here, before we apply the MOD function!
%
    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;
%
%  Have to do a MOD so G's don't go out of the region.
%
    g_new = max ( g_new, 0.0 );
    g_new = min ( g_new, 1.0 );
%
%  Sum the total energy for this iteration.
%
    e_plot(it) = sum ( e_new(1:g_num) ) / s_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-*', 'Linewidth', 2 )
  title ( 'Log (Energy)' )
  xlabel ( 'Step' )
  ylabel ( 'Energy' )
  grid on
  print ( '-dpng', 'cvtm_1d_energy.png' );
%
%  Plot the generator motion over iterations.
%
  figure ( 2 )
  plot ( 1:it_num, log ( m_plot ), 'm-*', 'Linewidth', 2 )
  title ( 'Log (Average generator motion)' )
  xlabel ( 'Step' )
  ylabel ( 'Motion' )
  grid on
  print ( '-dpng', 'cvtm_1d_motion.png' );
%
%  Plot the generator locations over iterations.
%
  figure ( 3 )
  y = ( 0 : it_num );
  for k = 1 : g_num
    plot ( g_plot(k,1:it_num+1), y, 'Linewidth', 2 )
    hold on;
  end
  grid on
  hold off;
  title ( 'Generator evolution' );
  xlabel ( 'Generator positions' );
  ylabel ( 'Iterations' ); 
  print ( '-dpng', 'cvtm_1d_generators.png' );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'CVTM_1D\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end
function timestamp ( )

%*****************************************************************************80
%
%% TIMESTAMP prints the current YMDHMS date as a timestamp.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  t = now;
  c = datevec ( t );
  s = datestr ( c, 0 );
  fprintf ( 1, '%s\n', s );

  return
end
