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
%
% Input:
%
% integer M, the spatial dimension.
%
% integer N, the number of cell generators.
%
% real A(M), B(M), the lower and upper limits of the region.
%
% real X(M), the point to be checked.
%
% real GENERATOR(M,N), the cell generators.
%
% logical MODULAR, is TRUE if modular arithmetic is to be used.
%
% Output:
%
% integer NEAREST, the index of the nearest cell generators.
%
% 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