function nearest_neighbor_test ( )

%*****************************************************************************80
%
%% nearest_neighbor_test() tests nearest_neighbor().
%
%  Discussion:
%
%    nearest_neighbor() computes nearest neighbors to a set of points from a
%    set of candidate points by Euclidean distance.  Points are considered to 
%    be matrix columns.
% 
%    Points can be specified either as a separate matrix, or as indices 
%    into the candidate matrix.
%
%    Multiple neighbors can be requested.
%
%    Neighbors within a certain radius can be requested.
%
%    The Delaunay triangulation is used if it would be faster.
%
%  Licensing:
%
%    Copyright 2006 Richard Brown. This code may be freely used and
%    distributed, so long as it maintains this copyright line.
%
%  Modified:
%
%    17 August 2024
%
%  Author:
%
%    Richard Brown
%
  addpath ( '../nearest_neighbor' );

  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'nearest_neighbor_test():\n' );
  fprintf ( 1, '  MATLAB/Octave version %s\n', version() );
  fprintf ( 1, '  Test nearest_neighbor().\n' );

  test01 ( );
  test02 ( );
  test03 ( );
  test04 ( );
  test05 ( );
  test06 ( );
  test07 ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'nearest_neighbor_test():\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  timestamp ( );

  rmpath ( '../nearest_neighbor' );

  return
end
function test01 ( )

%*****************************************************************************80
%
%% test01() does a basic nearest neighbor search.
%
%  Licensing:
%
%    Copyright 2006 Richard Brown. This code may be freely used and
%    distributed, so long as it maintains this copyright line.
%
%  Modified:
%
%    17 August 2024
%
%  Author:
%
%    Richard Brown
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'test01():\n' );
  fprintf ( 1, '  A basic nearest neighbor search.\n' );
%
%  2D points of interest
%
  P = rand ( 2, 5 );
%
%  Candidate point set
%
  X = rand ( 2, 50 );
%
%  Index the nearest X to each P.
%
  I = nearest_neighbor ( P, X );

  disp ( '  P: candidate points:' );
  disp ( P );
  disp ( '  X(I): Nearest neighbors' );
  disp ( X(:,I) );
%
%  Plot the points,
%
  hold ( 'on' );
  plot ( X(1,:), X(2,:), 'b.', ...
         P(1,:), P(2,:), 'r.', 'MarkerSize', 15 );
  quiver ( P(1, :), P(2, :), X(1,I) - P(1, :), X(2, I) - P(2, :), 0, 'k' );
  hold ( 'off' );
  filename = 'test01.png';
  print ( '-dpng', filename );
  fprintf ( 1, '  Graphics saved as "%s"\n', filename );

  return
end
function test02 ( )

%*****************************************************************************80
%
%% test02() find nearest neighbors from a subset.
%
%  Licensing:
%
%    Copyright 2006 Richard Brown. This code may be freely used and
%    distributed, so long as it maintains this copyright line.
%
%  Modified:
%
%    18 August 2024
%
%  Author:
%
%    Richard Brown
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'test02():\n' );
  fprintf ( 1, '  Find nearest neighbors from a subset of full list.\n' );
%
%  Candidate 4D feature set.
%
  X = rand(4, 20);
%
%  Find nearest neighbor to the 2nd, 4th, 10th features
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Nearest X to X[2,4,10].\n' );

  I = nearest_neighbor([2 4 10], X)
%
% Find nearest neighbors to all columns of X
%
  fprintf ( 1, '  Nearest X to any X.\n' );
  I = nearest_neighbor(X)

  return
end
function test03 ( )

%*****************************************************************************80
%
%% test03() finds 4 nearest neighbors.
%
%  Licensing:
%
%    Copyright 2006 Richard Brown. This code may be freely used and
%    distributed, so long as it maintains this copyright line.
%
%  Modified:
%
%    18 August 2024
%
%  Author:
%
%    Richard Brown
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'test03():\n' );
  fprintf ( 1, '  Find four nearest neighbors.\n' );
%
%  three points of interest
%
  P = rand(2, 3);
%
%  Candidate feature set.
%
  X = rand(2, 50);
%
%  Calculate the 4 nearest neighbors to each point: I(:, 1) will be the 4
%  nearest neighbors to P(:, 1), etc.
%
  fprintf ( 1, '\n' )
  fprintf ( 1, '  Index of four nearest neighbors to 3 points\n' );

  I = nearest_neighbor ( P, X, 'NumberOfNeighbors', 4 )
%
%  Plot the points - show the neighbors for one of them
%
  clf ( );
  hold ( 'on' );
  plot(P(1,:), P(2, :), 'r.', X(1,:), X(2,:), 'b.', 'MarkerSize', 15)
  p1 = repmat(P(1,1), 1, 4); p2 = repmat(P(2,1), 1, 4);
  quiver(p1, p2, X(1, I(:, 1)) - p1, X(2, I(:, 1)) - p2, 0, 'k')
  hold ( 'off' );
  title ( 'Show 4 nearest neighbors for just one point.' );
  filename = 'test03.png';
  print ( '-dpng', filename );
  fprintf ( 1, '  Graphics saved as "%s"\n', filename );

  return
end
function test04 ( )

%*****************************************************************************80
%
%% test04() finds neighbors within a specific radius.
%
%  Licensing:
%
%    Copyright 2006 Richard Brown. This code may be freely used and
%    distributed, so long as it maintains this copyright line.
%
%  Modified:
%
%    17 August 2024
%
%  Author:
%
%    Richard Brown
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'test04():\n' );
  fprintf ( 1, '  Find all neighbors within a specific radius.\n' );
%
%  Calculate the nearest neighbors to three points within a radius of 0.2
%
%  Three points of interest
%
  P = rand(2, 3);
%
%  Candidate feature set.
%
  X = rand(2, 50);
%
%  Calculate the nearest neighbors within a radius of 0.2 to each point:
%  I(:, 1) will be the neighbors to P(:, 1), etc. Note that if different
%  numbers of neighbors are found, the columns of I are padded with zeros
%  at the end to make I rectangular
%
  I = nearest_neighbor ( P, X, 'Radius', 0.2 )
%
%  Use idx to select neighbors of point 1;
%  Suppress any dummy neighbors.
%
  idx = I(I(:, 1) ~= 0, 1);
%
%  Plot the points - show the neighbors for one of them
%
  clf ( );
  hold ( 'on' );
  plot(P(1,:), P(2, :), 'r.', X(1,:), X(2,:), 'b.', 'MarkerSize', 15)
  p1 = repmat ( P(1,1), 1, length(idx) );
  p2 = repmat ( P(2,1), 1, length(idx) );
  quiver ( p1, p2, X(1, idx) - p1, X(2, idx) - p2, 0, 'k' );
  hold ( 'off' );
  title ( 'Neighbors of point 1 within a radius of 0.2' );
  filename = 'test04.png';
  print ( '-dpng', filename );
  fprintf ( 1, '  Graphics saved as "%s"\n', filename );

  return
end
function test05 ( )

%*****************************************************************************80
%
%% test05() finds up to 5 neighbors within a specific radius.
%
%  Licensing:
%
%    Copyright 2006 Richard Brown. This code may be freely used and
%    distributed, so long as it maintains this copyright line.
%
%  Modified:
%
%    18 August 2024
%
%  Author:
%
%    Richard Brown
%
  fprintf ( 1, '\n' );

  fprintf ( 1, 'test05():\n' );
  fprintf ( 1, '  Find up to 5 neighbors within a specific radius.\n' );
%
% Point of interest.
%
  P = rand(2, 1);
%
% Candidate feature set.
%
  X = rand(2, 30);
%
% Calculate the 5 nearest neighbors to P, specifying that they must be
% at most 0.2 away by Euclidean distance.
%
  I = nearest_neighbor(P, X, 'n', 5, 'r', 0.2);
  n = length(I);
%
% Plot the points - show the neighbors for one of them
%
  clf ( );
  hold ( 'on' );
  plot(P(1,:), P(2, :), 'r.', X(1,:), X(2,:), 'b.', 'MarkerSize', 15 );
  quiver(repmat ( P(1), 1, n), repmat(P(2), 1, n), ...
    X(1, I) - P(1), X(2, I) - P(2), 0, 'k');
  hold ( 'off' );
  title ( '  Up to 5 neighbors within 0.2 radius' );
  filename = 'test05.png';
  print ( '-dpng', filename );
  fprintf ( 1, '  Graphics saved as "%s"\n', filename );

  return
end
function test06 ( )

%*****************************************************************************80
%
%% test06() looks at Delaunay mode for large data, small pointset.
%
%  Licensing:
%
%    Copyright 2006 Richard Brown. This code may be freely used and
%    distributed, so long as it maintains this copyright line.
%
%  Modified:
%
%    17 August 2024
%
%  Author:
%
%    Richard Brown
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'test06():\n' );
  fprintf ( 1, '  Delaunay decision for large data, small pointset.\n' );
%
%  We'll run a couple of examples and show how |nearest_neighbor|
%  automatically works out when to use the Delaunay triangulation. 
%  We'll try a large data set with few points to lookup - it shouldn't use 
%  Delaunay Triangulation as the overhead required to compute this should 
%  definitely offset any gains.
%
  d = 3;
  np = 10;
  nx = 5000;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Use a large data set, small number of points.\n' );
  fprintf ( 1, '  Spatial dimension d = %d\n', d );
  fprintf ( 1, '  Number of candidate points np = %d\n', np );
  fprintf ( 1, '  Number of data points nx = %d\n', nx );

  P = rand(d, np);
  X = rand(d, nx);
%
%  DelaunayMode off
%
  t = timing_test ( P, X, 'off' );
  fprintf ( 1, '  (DelaunayMode off) : %.2g seconds\n', t );
%
%  DelaunayMode on
%
  t = timing_test ( P, X, 'on' );
  fprintf ( 1, '  (DelaunayMode on)  : %.2g seconds\n', t );
%
%  DelaunayMode automatic.
%
  t = timing_test ( P, X, 'auto' );
  fprintf ( 1, '  (DelaunayMode auto): %.2g seconds\n\n', t );

  return
end
function test07 ( )

%*****************************************************************************80
%
%% test07() looks at the decision to use the Delaunay approach.
%
%  Licensing:
%
%    Copyright 2006 Richard Brown. This code may be freely used and
%    distributed, so long as it maintains this copyright line.
%
%  Modified:
%
%    17 August 2024
%
%  Author:
%
%    Richard Brown
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'test07():\n' );
  fprintf ( 1, '  Delaunay decision for small data, large pointset.\n' );
%
%  Small Data set, lots of candidate points.
%
  d = 3;
  np = 10000;
  nx = 500;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Use a small data set, lots of candidate points.\n' );
  fprintf ( 1, '  Spatial dimension d = %d\n', d );
  fprintf ( 1, '  Number of candidate points np = %d\n', np );
  fprintf ( 1, '  Number of data points nx = %d\n', nx );

  P = rand ( d, np );
  X = rand ( d, nx );
%
%  Run with DelaunayMode off
%
  t = timing_test ( P, X, 'off' );
  fprintf ( 1, '  (Delaunay mode off): %.2g seconds\n', t )
%
%  Run with DelaunayMode on
%
  t = timing_test ( P, X, 'on' );
  fprintf ( 1, '  (DelaunayMode on)  : %.2g seconds\n', t )
%
%  Run with DelaunayMode automatic
%
  t = timing_test ( P, X, 'auto' );
  fprintf ( 1, '  (DelaunayMode auto): %.2g seconds\n', t )

  return
end
function tAv = timing_test ( P, X, mode )

%*****************************************************************************80
%
%% timing_test() estimates average cputime required to run nearest_neighbor().
%
%  Example:
%
%    % A set of points
%    X = rand(2, 1000)
%    % A set of points to find the neighbors of
%    P = rand(2, 100)
%    % Compute the time to run with DelaunayMode 'on'
%    timing_test(P, X, 'on')
%
%  Licensing:
%
%    Copyright 2006 Richard Brown. This code may be freely used and
%    distributed, so long as it maintains this copyright line.
%
%  Modified:
%
%    17 August 2024
%
%  Author:
%
%    Richard Brown
%
  tRun = 1;

  t = 0;
  n = 0;
  while ( t < tRun )
    t0 = cputime ( );
    nearest_neighbor ( P, X, 'DelaunayMode', mode );
    t = t + cputime ( ) - t0;
    n = n + 1;
  end

  tAv = t / n;

  return
end
function timestamp ( )

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

  return
end

