# include <math.h>
# include <stdio.h>
# include <stdlib.h>
# include <time.h>

int main ( );
double cpu_time ( );
int *find_closest1 ( int m, int nr, double r[], int ns, double s[] );
void timestamp ( );

/******************************************************************************/

int main ( )

/******************************************************************************/
/*
  Purpose:

    test_nearest() compares the performance of nearest neighbor routines.

  Discussion:

    We are given R, a set of NR points in M dimensions.

    We are given S, a set of NS points in M dimensions.

    For each S(I) in S, we seek the index J of the point R(J)
    which is nearest to S(I) over all points in R.

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    26 December 2012

  Author:

    John Burkardt

  Local:

    int M, the spatial dimension.

    int NR, the number of data points.

    int NS, the number of sample points.

    double R[M*NR], the data points.

    double S[M*NS], the sample points. 
*/
{
  int i;
  int j;
  int m;
  int m_test;
  int m_test_data[3] = { 2, 4, 8 };
  int m_test_num = 3;
  int *nearest;
  int nr;
  int nr_test_data[6] = {
    1000000, 100000, 10000,  1000,    100,      10 };
  int n_test_num = 6;
  int ns;
  int ns_test_data[6] = {
      10,    100,  1000, 10000, 100000, 1000000 };
  double *r;
  double *s;
  double t1;
  double t2;
  int test;

  timestamp ( );
  printf ( "\n" );
  printf ( "test_nearest():\n" );
  printf ( "  C version\n" );
  printf ( "  Consider various nearest neighbor algorithms.\n" );

  for ( m_test = 0; m_test < m_test_num; m_test++ )
  {
    m = m_test_data[m_test];

    for ( test = 0; test < n_test_num; test++ )
    {
      nr = nr_test_data[test];
      r = ( double * ) malloc ( m * nr * sizeof ( double ) );

      for ( i = 0; i < m; i++ )
      {
        for ( j = 0; j < nr; j++ )
        {
          r[i+j*m] = drand48 ( );
        }
      }

      ns = ns_test_data[test];
      s = ( double * ) malloc ( m * ns * sizeof ( double ) );

      for ( i = 0; i < m; i++ )
      {
        for ( j = 0; j < ns; j++ )
        {
          s[i+j*m] = drand48 ( );
        }
      }
      printf ( "\n" );
      printf ( "  M = %d, NR = %d,  NS = %d\n", m, nr, ns );

      t1 = cpu_time ( );
      nearest = find_closest1 ( m, nr, r, ns, s );
      t2 = cpu_time ( );
      printf ( "  #1 time: %g,  size = %d,  i[0] = %d\n", t2 - t1, ns, nearest[0] );

      free ( nearest );
      free ( r );
      free ( s );
    }
  }
/*
  Terminate.
*/
  printf ( "\n" );
  printf ( "test_nearest():\n" );
  printf ( "  Normal end of execution.\n" );
  printf ( "\n" );
  timestamp ( );

  return 0;
}
/******************************************************************************/

double cpu_time ( )

/******************************************************************************/
/*
  Purpose:

    cpu_time() returns the current reading on the CPU clock.

  Discussion:

    The CPU time measurements available through this routine are often
    not very accurate.  In some cases, the accuracy is no better than
    a hundredth of a second.  

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    06 June 2005

  Author:

    John Burkardt

  Output:

    double CPU_TIME, the current reading of the CPU clock, in seconds.
*/
{
  double value;

  value = ( double ) clock ( ) 
        / ( double ) CLOCKS_PER_SEC;

  return value;
}
/******************************************************************************/

int *find_closest1 ( int m, int nr, double r[], int ns, double s[] )

/******************************************************************************/
/*
  Purpose:

    find_closest1() finds the nearest R point to each S point.

  Discussion:

    We are given R, a set of NR points in M dimensions.

    We are given S, a set of NS points in M dimensions.

    For each S(I) in S, we seek the index J of the point R(J)
    which is nearest to S(I) over all points in R.

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    26 December 2012

  Author:

    John Burkardt

  Input:

    int M, the spatial dimension.

    int NR, the number of data points.

    double R[M*NR], the data points.

    int NS, the number of sample points.

    double S[M*NS], the sample points.

  Output:

    int FIND_CLOSEST1[NS], the index of the nearest data point.
*/
{
  double dist_sq;
  double dist_sq_min;
  int i;
  int jr;
  int js;
  int *nearest;

  nearest = ( int * ) malloc ( ns * sizeof ( int ) );

  for ( js = 0; js < ns; js++ )
  {
    dist_sq_min = HUGE_VAL;
    nearest[js] = -1;
    for ( jr = 0; jr < nr; jr++ )
    {
      dist_sq = 0.0;
      for ( i = 0; i < m; i++ )
      {
        dist_sq = dist_sq + pow ( r[i+jr*m] - s[i+js*m], 2 );
      }
      if ( dist_sq < dist_sq_min )
      {
        dist_sq_min = dist_sq;
        nearest[js] = jr;
      }
    }
  }

  return nearest;
}
/******************************************************************************/

void timestamp ( )

/******************************************************************************/
/*
  Purpose:

    timestamp() prints the current YMDHMS date as a time stamp.

  Example:

    31 May 2001 09:45:54 AM

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    24 September 2003

  Author:

    John Burkardt

  Parameters:

    None
*/
{
# define TIME_SIZE 40

  static char time_buffer[TIME_SIZE];
  const struct tm *tm;
  time_t now;

  now = time ( NULL );
  tm = localtime ( &now );

  strftime ( time_buffer, TIME_SIZE, "%d %B %Y %I:%M:%S %p", tm );

  fprintf ( stdout, "%s\n", time_buffer );

  return;
# undef TIME_SIZE
}