function cvt_1d_lumping ( n, it_num, s_num, mu_fun, init )

%*****************************************************************************80
%
%% cvt_1d_lumping() carries out the Lloyd algorithm with lumping.
%
%  Discussion:
%
%    This program is a variation of the CVT_1D_LLOYD method
%    for a problem with an associated density rho(x).
%
%    Instead of sampling points uniformly with respect to the density,
%    it samples points uniformly with respect to geometry,
%    and weights them by rho(x).
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    22 August 2016
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    integer N, the number of generators.
%
%    integer IT_NUM, the number of CVT iterations.
%
%    integer S_NUM, the number of samples in [-1,1] to use
%    when estimating the Voronoi regions.
%
%    function value = MU_FUN ( X ), the name of the M file 
%    which accepts a real vector of points X and returns a real vector
%    VALUE of density function values.
%
%    integer INIT, the initialization option.
%    *1, random sorted values;
%    *2, Chebyshev Zero Nodes;
%    *3, uniform staggered values.
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'CVT_1D_LUMPING\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Use lumping to implement a version of Lloyd''s algorithm\n' );
  fprintf ( 1, '  in the 1D interval [-1,1].\n' );

  if ( nargin < 1 )
    n = input ( '  Enter number of generators: ' );
  elseif ( ischar ( n ) )
    n = str2num ( n );
  end

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

  if ( nargin < 3 ) 
    s_num = input ( '  Enter number of sample points: ' );
  elseif ( ischar ( s_num ) )
    s_num = str2num ( s_num );
  end
%
%  I don't know how to deal with a missing function argument.
%
  if ( nargin < 4 )
  end

  if ( nargin < 5 )
    init = input ( '  Enter initialization option 1/2/3: ' );
  elseif ( ischar ( init ) )
    init = str2num ( init );
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of generators is %d\n', n );
  fprintf ( 1, '  Number of iterations is %d\n', it_num );
  fprintf ( 1, '  Number of samples is %d\n', s_num );
%
%  Define the Chebyshev Zero nodes.
%
  gc = zeros ( n, 1 );
  gc(1:n,1) = cos ( pi * ( 2 * n - 2 * (1:n)' + 1 ) / ( 2 * n ) );
%
%  Initialize the generators.
%
  if ( init == 1 )
    fprintf ( 1, '\n' );
    fprintf ( 1, '  Initial generators (Random, sorted):\n' );
    g(1:n,1) = 2.0 * rand ( n, 1 ) - 1.0;
    g = sort ( g );
  elseif ( init == 2 )
    fprintf ( 1, '\n' );
    fprintf ( 1, '  Initial generators (Chebyshev Zero nodes):\n' );
    g(1:n,1) = gc(1:n,1);
  else
    fprintf ( 1, '\n' );
    fprintf ( 1, '  Initial generators (Uniform, staggered):\n' );
    g(1:n,1) = r8vec_linspace2 ( n, -1.0, +1.0 );
  end
%
%  Print the initial generators.
%
  fprintf ( 1, '\n' );
  for k = 1 : n
    fprintf ( 1, '  %4d  %f\n', k, g(k,1) );
  end
%
%  Initialize some arrays.
%
  g_new = zeros ( n, 1 );
  g_plot = zeros ( n, it_num + 1 );
  step = ( 1 : it_num )';
  e = nan ( it_num, 1 );
  gm = nan ( it_num, 1 );
%
%  S is the uniform sample points.  Gee, that's easy!
%
  s = linspace ( -1.0+eps, 1.0-eps, s_num );
  s = s';
%
%  R is the density function RHO(X).
%
%  We will NOT assume it is normalized.
%
%  We will assume it is the cube of some other density function MU(X).
%
%  In this way, we expect the asymptotic distribution of CVT generators
%  to follow the distribution MU(X).
%
  mu = mu_fun ( s );
  r = mu .^ 3;
%  
%  Plot MU(X) and RHO(X), but truncate at 10 so that the plot
%  is useful rather than correct.
%
  i1 = find ( r < 10 );
  i2 = find ( mu < 10 );
  figure ( 1 );
  plot ( s(i1), r(i1), 'r-', s(i2), mu(i2), 'b-', 'Linewidth', 3 );
  axis ( [ -1, +1, 0, 10 ] );
  title ( 'Densities \mu(x) = blue and \rho(x) = red:' );
  grid ( 'on' );
  filename = 'density.png';
  print ( '-dpng', filename );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Saved graphics file "%s"\n', filename );
%
%  Start the timer for the iteration.
%
  tic;
%
%  Now carry out the iteration.
%
  for it = 1 : it_num

    g_plot ( 1:n, it ) = g(1:n,1);
%
%  Bracket the generators.
%
    gb = zeros ( n + 1, 1 );
    gb(1) = -1.0;
    for j = 2 : n
      gb(j) = 0.5 * ( g(j-1) + g(j) );
    end
    gb(n+1) = +1.0;

    left = r8vec_bracket4 ( s_num, s, n + 1, gb );
%
%  Compute the new generators.
%
    e(it,1) = 0.0;

    k2 = 0;

    for j = 1 : n

      k1 = k2 + 1;
      k2 = left(j+1) - 1;

      if ( k2 <= k1 )
        g_new(j) = g(j);
      else
        r_new = sum ( r(k1+1:k2) );
        g_new(j) = r(k1+1:k2)' * s(k1+1:k2) / r_new;
        for k = k1 + 1 : k2
          e(it,1) = e(it,1) + r(k) * ( s(k) - g(j) ).^2;
        end
      end

    end
%
%  If N is odd, try to enforce symmetry.
%
    if ( mod ( n, 2 ) == 1 )
      g_new((n+1)/2) = 0.0;
    end

    e(it,1) = e(it,1) / s_num;
%
%  Compute the generator motion.
%
    gm(it,1) = sum ( ( g_new(:) - g(:) ).^2 ) / n;
%
%  Update the generators.
%
    g = g_new;

  end
%
%  Stop the timer for the iteration.
%
  duration = toc;
%
%  Save the final generators.
%
  g_plot(1:n,it_num+1) = g(1:n,1);
%
%  Display the energy.
%
  figure ( 2 );
  plot ( step, log ( e ), 'm-*' );
  title ( 'Log (Energy)' );
  xlabel ( 'Iteration' );
  ylabel ( 'Energy' );
  grid ( 'on' );
  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 ( 3 );
  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, '  Graphics saved as "%s"\n', filename );
%
%  Plot the evolution of the locations of the generators.
%
  figure ( 4 );
  hold ( 'on' );
  y = ( 0 : it_num );
  for k = 1 : n
    plot ( g_plot(k,1:it_num+1), y, 'Linewidth', 2 );
  end
  grid ( 'on' );
  hold ( 'off' );

  title ( 'Generator evolution' );
  xlabel ( 'Generator positions' );
  ylabel ( 'Iterations' );
 
  if ( init == 1 )
    filename = sprintf ( 'evolution_%03d_r.png', n );
  elseif ( init == 2 )
    filename = sprintf ( 'evolution_%03d_c.png', n );
  elseif ( init == 3 )
    filename = sprintf ( 'evolution_%03d_u.png', n );
  end
  print ( '-dpng', filename );
  fprintf ( 1, '  Saved graphics file "%s"\n', filename );
%
%  Print the current generators and energy distribution.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  CVT data:\n' );
  fprintf ( 1, '      K          G             H          rho H^3          mu H\n' );
  fprintf ( 1, '\n' );
  for k = 1 : n
    h = gb(k+1) - gb(k);
    rh3 = mu_fun ( g(k,1) )^3 * h^3;
    muh = mu_fun ( g(k,1) )   * h;
    fprintf ( 1, '  %4d  %12.4f  %12.4g  %12.4g  %12.4g\n', k, g(k,1), h, rh3, muh );
  end
%
%  Compare with Chebyshev zero points:
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Chebyshev Zero Node data for comparison\n' );
  fprintf ( 1, '     K          G             H          rho H^3          mu H\n' );
  fprintf ( 1, '\n' );
  gbc = zeros ( n + 1, 1 );
  gbc(1) = -1.0;
  for k = 2 : n
    gbc(k) = 0.5 * ( gc(k-1) + gc(k) );
  end
  gbc(n+1) = +1.0;
  for k = 1 : n
    h = gbc(k+1) - gbc(k);
    rh3 = mu_fun ( gc(k) )^3 * h^3;
    muh = mu_fun ( gc(k) ) * h;
    fprintf ( 1, '  %4d  %12.4f  %12.4g  %12.4g  %12.4g\n', k, gc(k), h, rh3, muh );
  end
%
%  Compute norm of difference between generators and Chebyshev Zero nodes.
%
  t = norm ( g(1:n,1) - gc(1:n,1) );
  fprintf ( '\n' );
  fprintf ( '  |G-ChebyshevZero| = %g\n', t );
%
%  Write CVT generators to file.
%
  if ( init == 1 )
    filename = sprintf ( 'cvt_%03d_r.txt', n );
  elseif ( init == 2 )
    filename = sprintf ( 'cvt_%03d_c.txt', n );
  elseif ( init == 3 )
    filename = sprintf ( 'cvt_%03d_u.txt', n );
  end

  r8vec_write ( filename, n, g );

  fprintf ( 1, '  Saved CVT coordinates in file "%s"\n', filename );
%
%  Write Chebyshev Zero nodes to file.
%
  filename = sprintf ( 'chebyzero_%03d.txt', n );

  r8vec_write ( filename, n, gc );

  fprintf ( 1, '  Saved Chebyshev Zero node coordinates in file "%s"\n', filename );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Iteration required %g seconds.\n', duration );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'CVT_1D_LUMPING\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end
function left = r8vec_bracket4 ( nt, t, ns, s )

%*****************************************************************************80
%
%% r8vec_bracket4() finds the interval to each of a vector of values.
%
%  Discussion:
%
%    An R8VEC is a vector of R8's.
%
%    The routine always returns the index LEFT of the sorted array
%    T with the property that either
%    *  T is contained in the interval [ T(LEFT), T(LEFT+1) ], or
%    *  T < T(LEFT) = T(1), or
%    *  T > T(LEFT+1) = T(N).
%
%    The routine is useful for interpolation problems, where
%    the abscissa must be located within an interval of data
%    abscissas for interpolation, or the "nearest" interval
%    to the (extreme) abscissa must be found so that extrapolation
%    can be carried out.
%
%  Licensing:
%
%    This code is distributed under the MIT license. 
%
%  Modified:
%
%    30 April 2009
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    integer NT, length of the input array.
%
%    real T(NT), an array that has been sorted
%    into ascending order.
%
%    integer NS, the number of points to be bracketed.
%
%    real S(NS), values to be bracketed by entries of T.
%
%  Output:
%
%    integer LEFT(NS).
%    LEFT(I) is set so that the interval [ T(LEFT(I)), T(LEFT(I)+1) ]
%    is the closest to S(I); it either contains S(I), or else S(I)
%    lies outside the interval [ T(1), T(NT) ].
%

%
%  Check the input data.
%
  if ( nt < 2 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8VEC_BRACKET4 - Fatal error!\n' );
    fprintf ( 1, '  NT must be at least 2.\n' );
    error ( 'R8VEC_BRACKET4 - Fatal error!' );
  end

  for i = 1 : ns

    left(i) = floor ( ( nt + 1 ) / 2 );
%
%  CASE 1: S < T(LEFT):
%  Search for S in [T(I), T(I+1)] for intervals I = 1 to LEFT-1.
%
    if ( s(i) < t(left(i)) )

      if ( left(i) == 1 )
        continue
      elseif ( left(i) == 2 )
        left(i) = 1;
        continue
      elseif ( t(left(i)-1) <= s(i) )
        left(i) = left(i) - 1;
        continue
      elseif ( s(i) <= t(2) )
        left(i) = 1;
        continue
      end
%
%  ...Binary search for S in [T(I), T(I+1)] for intervals I = 2 to LEFT-2.
%
      low = 2;
      high = left(i) - 2;

      while ( true )
  
        if ( low == high )
          left(i) = low;
          break
        end

        mid = floor ( ( low + high + 1 ) / 2 );

        if ( t(mid) <= s(i) )
          low = mid;
        else
          high = mid - 1;
        end

      end
%
%  CASE2: T(LEFT+1) < S:
%  Search for S in [T(I),T(I+1)] for intervals I = LEFT+1 to N-1.
%
    elseif ( t(left(i)+1) < s(i) )

      if ( left(i) == nt - 1 )
        continue
      elseif ( left(i) == nt - 2 )
        left(i) = left(i) + 1;
        continue
      elseif ( s(i) <= t(left(i)+2) )
        left(i) = left(i) + 1;
        continue
      elseif ( t(nt-1) <= s(i) )
        left(i) = nt - 1;
        continue
      end
%
%  ...Binary search for S in [T(I), T(I+1)] for intervals I = LEFT+2 to NT-2.
%
      low = left(i) + 2;
      high = nt - 2;

      while ( true )

        if ( low == high )
          left(i) = low;
          break
        end

        mid = floor ( ( low + high + 1 ) / 2 );

        if ( t(mid) <= s(i) )
          low = mid;
        else
          high = mid - 1;
        end

      end
%
%  CASE3: T(LEFT) <= S <= T(LEFT+1):
%  S is in [T(LEFT), T(LEFT+1)].
%
    else

    end

  end

  return
end
function x = r8vec_linspace2 ( n, a, b )

%*****************************************************************************80
%
%% r8vec_linspace2() creates a vector of linearly spaced values.
%
%  Discussion:
%
%    An R8VEC is a vector of R8's.
%
%    4 points evenly spaced between 0 and 12 will yield 2, 4, 6, 8, 10.
%
%    In other words, the interval is divided into N+1 even subintervals,
%    and the endpoints of internal intervals are used as the points.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    17 September 2012
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    integer N, the number of entries in the vector.
%
%    real A, B, the first and last entries.
%
%  Output:
%
%    real X(N,1), a vector of linearly spaced data.
%
  x = zeros ( n, 1 );

  for i = 1 : n
    x(i) = ( ( n  - i + 1 ) * a ...
           + (      i     ) * b ) ...
           / ( n      + 1 );
  end

  return
end
function r8vec_write ( output_filename, n, x )

%*****************************************************************************80
%
%% r8vec_write() writes an R8VEC file.
%
%  Discussion:
%
%    An R8VEC is a vector of R8's.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    10 July 2011
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    string OUTPUT_FILENAME, the output filename.
%
%    integer N, the number of points.
%
%    real X(N), the vector.
%

%
%  Open the file.
%
  output_unit = fopen ( output_filename, 'wt' );

  if ( output_unit < 0 ) 
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8VEC_WRITE - Error!\n' );
    fprintf ( 1, '  Could not open the output file.\n' );
    error ( 'R8VEC_WRITE - Error!' );
  end

  for j = 1 : n
    fprintf ( output_unit, '  %24.16g\n', x(j) );
  end
%
%  Close the file.
%
  fclose ( output_unit );

  return
end
function timestamp ( )

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

  return
end

