function [ nearest, xp ] = cvtp_find_closest ( m, n, a, b, x, generator, modular )

%*****************************************************************************80
%
%% CVTP_FIND_CLOSEST finds the Voronoi cell generator closest to a point X.
%
%  Discussion:
%
%    This routine finds the closest Voronoi cell generator by checking every
%    one.  For problems with many cells, this process can take the bulk
%    of the CPU time.  Other approaches, which group the cell generators into
%    bins, can run faster by a large factor.
%
%    For this routine, if MODULAR is TRUE, then distance is done in a modular 
%    sense, as though the points were on a generalized torus.  It's simple, 
%    really, we just need, in each coordinate, to consider 
%
%     X(I)-WIDTH(I), X(I), and X(I)+WIDTH(I).
%
%    The bad part is, to keep our sanity, we want to replace X on output
%    by the actual coordinates that got closest to some generator G,
%    even though some of these coordinates may lie outside the unit
%    hypercube.  This is the right thing to do, so that the averaging
%    process works correctly.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license. 
%
%  Modified:
%
%    31 July 2016
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, the spatial dimension.
%
%    Input, integer N, the number of cell generators.
%
%    Input, real A(M), B(M), the lower and  upper limits of the region.
%
%    Input, real X(M), the point to be checked.
%
%    Input, real GENERATOR(M,N), the cell generators.
%
%    Input, logical MODULAR, is TRUE if modular arithmetic is to be used.
%
%    Output, integer NEAREST, the index of the nearest cell
%    generators.
%
%    Output, real XP(M), the periodic image of X which achieves the
%    minimum distance to a generator in the region.  
%
  nearest = -1;

  dist_sq_min = Inf;
  xp = x;

  for i = 1 : n

    dist_sq = 0.0;

    for j = 1 : m

      if ( modular )

        side1 = abs ( generator(j,i) - x(j) );
        side2 = abs ( generator(j,i) - x(j) + b(j) - a(j) );
        side3 = abs ( generator(j,i) - x(j) - b(j) + a(j) );
      
        if ( side2 < side1 && side2 < side3 )
          side = side2;
          y(j) = x(j) - b(j) + a(j);
        elseif ( side3 < side1 && side3 < side2 )
          side = side3;
          y(j) = x(j) + b(j) - a(j);
        else
          side = side1;
          y(j) = x(j);
        end

      else

        side = abs ( generator(j,i) - x(j) );
        y(j) = x(j);

      end

      dist_sq = dist_sq + side ^ 2;

    end

    if ( dist_sq < dist_sq_min )
      dist_sq_min = dist_sq;
      nearest = i;
      xp(1:m) = y(1:m);
    end

  end

  return
end

