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

# include "ppc_array.h"
# include "ppc_nelder_mead.h"
# include "ppc_xmalloc.h"

# define N 5

int main ( );
void demo2d ( );
double demo2d_f ( double *x, int n, void *param );
void get_centroid_test ( );
int incr ( double *a, int n, int k );
int increment ( double *a, int n );
void print_row ( double *a, int n );
void print_simplex ( double **S, int n );
void rank_vertices_test ( );
void replace_row_test ( );
void shrink_test ( );
void test_result ( double *a, int n, int ia, int iy, int iz );
void timestamp ( );
void transform_test ( );

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

int main ( ) 

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

    ppc_nelder_mead_test() tests ppc_nelder_mead().

  Discussion:

    The code generates all vectors of length n with entries coming from
    the set {0,1,2,...,n-1}.

    Then it lets n vary from 2 to N, where N is #defined below.

    For each generated vector, it runs the given rank_vertices() function.

    If the result of rank_vertices() is incorrect, it prints a diagnostic,
    otherwise no output is generated.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    13 May 2024

  Author:

    John Burkardt.

  Reference:

    Rouben Rostamian,
    Programming Projects in C 
    for Students of Engineering, Science, and Mathematics,
    SIAM, 2014,
    ISBN: 978-1-611973-49-5
*/
{
  timestamp ( );

  printf ( "\n" );
  printf ( "ppc_nelder_mead_test():\n" );
  printf ( "  C version\n" );
  printf ( "  Test ppc_nelder_mead():\n" );
/*
  Test utilities.
*/
  get_centroid_test ( );
  rank_vertices_test ( );
  replace_row_test ( );
  shrink_test ( );
  transform_test ( );
/*
  Test nelder_mead()
*/
  demo2d ( );
/*
  Terminate.
*/
  printf ( "\n" );
  printf ( "ppc_nelder_mead_test():\n" );
  printf ( "  Normal end of execution.\n" );

  timestamp ( );

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

void demo2d ( ) 

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

    demo2d() is a simple demonstration of nelder_mead().

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    13 May 2024

  Author:

    John Burkardt.

  Reference:

    Rouben Rostamian,
    Programming Projects in C 
    for Students of Engineering, Science, and Mathematics,
    SIAM, 2014,
    ISBN: 978-1-611973-49-5
*/
{
  int evalcount;
  int n = 2;
  double value;
  double x[] = { 0.0, 0.0 };

  struct nelder_mead NM =
  {
    .f = demo2d_f,
    .n = n,
    .s = NULL,
    .x = x,
    .h = 0.1,
    .tol = 0.0001,
    .maxevals = 1000,
    .params = NULL
  };

  printf ( "\n" );
  printf ( "demo2d():\n" );
  printf ( "  A simple demonstration of nelder_mead().\n" );

  value = demo2d_f ( x, n, NULL );
  printf ( "  For initial x = (%g,%g)\n", x[0], x[1] );
  printf ( "  f(x) = %g\n", value );

  evalcount = nelder_mead ( &NM );

  if ( NM.maxevals < evalcount )
  {
    printf ( "\n" );
    printf ( "  nelder_mead() did not converge.\n" );
  }
  else
  {
    printf ( "\n" );
    printf ( "  nelder_mead() converged after %d function evaluations\n",
      evalcount );
    printf ( "  Computed solution:\n" );
    printf ( "  For x = (%g,%g)\n", NM.x[0], NM.x[1] );
    printf ( "  f(x) = %g\n", NM.minval );
  }
  return;
}
/******************************************************************************/

double demo2d_f ( double *x, int n, void *params ) 

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

    demo2d_f() evaluates the objective function for demo2d().
*/
{
  double value;

  value = pow ( x[0], 2 ) + pow ( x[1], 2 ) - 4.0 * x[0] - 2.0 * x[1] + 8.0;

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

void get_centroid_test ( ) 

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

    get_centroid_test() tests get_centroid().

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    13 May 2024

  Author:

    John Burkardt.

  Reference:

    Rouben Rostamian,
    Programming Projects in C 
    for Students of Engineering, Science, and Mathematics,
    SIAM, 2014,
    ISBN: 978-1-611973-49-5
*/
{
  double *C;
  int iz;
  int n = 2;
  double **S;

  printf ( "\n" );
  printf ( "get_centroid_test():\n" );
  printf ( "  get_centroid() returns the centroid of the simplex face\n" );
  printf ( "  opposite to vertex x(iz).\n" );

  make_matrix ( S, n + 1, n );

  S[0][0] = 0.0;
  S[0][1] = 0.0;
  S[1][0] = 1.0;
  S[1][1] = 0.0;
  S[2][0] = 0.0;
  S[2][1] = 3.0;

  printf ( "\n" );
  printf ( "  Simplex coordinates:\n" );
  print_simplex ( S, n );

  make_vector ( C, n );

  printf ( "\n" );
  printf ( "  For node iz, centroid of opposite face:\n" );
  printf ( "\n" );
  for ( iz = 0; iz < n + 1; iz++ )
  {
    get_centroid ( S, n, iz, C );
    printf ( "  %d  %10.4f  %10.4f\n", iz, C[0], C[1] );
  }
  return;
}
/******************************************************************************/

int incr ( double *a, int n, int k )

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

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    13 May 2024

  Author:

    Original C version by Rougen Rostamian.
    This version by John Burkardt.

  Reference:

    Rouben Rostamian,
    Programming Projects in C 
    for Students of Engineering, Science, and Mathematics,
    SIAM, 2014,
    ISBN: 978-1-611973-49-5
*/
{
  if ( k == n )
  {
    return 1;
  }

  if ( a[k] < n - 1 )
  {
    a[k]++;
    return 0;
  } 
  else
  {
    a[k] = 0;
    return incr ( a, n, k+1 );
  }
}
/******************************************************************************/

int increment ( double *a, int n )

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

    increment() increments a vector.

  Discussion:

    The code receives a vector a[] of length n containing entries from the set 0..(n-1).
    It "increments" the vector, as if the vector entries were the digits of
    an number in base n.  Returns 0 on success or 1 on failure, which happens
    when all vector entries have hit their max.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    13 May 2024

  Author:

    Original C version by Rouben Rostamian.
    This version by John Burkardt.

  Reference:

    Rouben Rostamian,
    Programming Projects in C 
    for Students of Engineering, Science, and Mathematics,
    SIAM, 2014,
    ISBN: 978-1-611973-49-5
*/
{
  return incr ( a, n, 0 );
}
/******************************************************************************/

void print_row ( double *a, int n )

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

    print_row() prints a row of values.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    13 May 2024

  Author:

    John Burkardt.

  Reference:

    Rouben Rostamian,
    Programming Projects in C 
    for Students of Engineering, Science, and Mathematics,
    SIAM, 2014,
    ISBN: 978-1-611973-49-5
*/
{
  int i;

  for ( i = 0; i < n; i++ )
  {
    printf( "%g ", a[i] );
  }
  return;
}
/******************************************************************************/

void print_simplex ( double **S, int n )

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

    print_simplex() prints the rows of the simplex matrix.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    13 May 2024

  Author:

    John Burkardt.

  Reference:

    Rouben Rostamian,
    Programming Projects in C 
    for Students of Engineering, Science, and Mathematics,
    SIAM, 2014,
    ISBN: 978-1-611973-49-5
*/
{
  int i;
  int j;

  printf ( "\n" );
  for ( i = 0; i < n + 1; i++ )
  {
    for ( j = 0; j < n; j++ )
    {
      printf ( "  %10.4f", S[i][j] );
    }
    printf ( "\n" );
  }
  return;
}
/******************************************************************************/

void rank_vertices_test ( )

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

    rank_vertices_test() tests rank_vertices().

  Discussion:

    The code generates all vectors of length n with entries coming from
    the set {0,1,2,...,n-1}.

    Then it lets n vary from 2 to N, where N is #defined below.

    For each generated vector, it runs the given rank_vertices() function.

    If the result of rank_vertices() is incorrect, it prints a diagnostic,
    otherwise no output is generated.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    13 May 2024

  Author:

    John Burkardt.

  Reference:

    Rouben Rostamian,
    Programming Projects in C 
    for Students of Engineering, Science, and Mathematics,
    SIAM, 2014,
    ISBN: 978-1-611973-49-5
*/
{
  double a[N+1];
  int i;
  int ia;
  int iy;
  int iz;
  int n;

  printf ( "\n" );
  printf ( "rank_vertices_test():\n" );
  printf ( "  Test rank_vertices(), which returns the indices of the\n" );
  printf ( "  lowest, second highest, and highest values in a vector.\n" );
 
  for ( n = 2; n <= N; n++ )
  {
    for ( i = 0; i < n; i++ )
    {
      a[i] = 0;
    }
    do 
    {
      rank_vertices ( a, n, &ia, &iy, &iz );
      test_result ( a, n, ia, iy, iz );
    } while ( increment ( a, n ) != 1 );
  }

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

void replace_row_test ( ) 

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

    replace_row_test() tests replace_row().

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    13 May 2024

  Author:

    John Burkardt.

  Reference:

    Rouben Rostamian,
    Programming Projects in C 
    for Students of Engineering, Science, and Mathematics,
    SIAM, 2014,
    ISBN: 978-1-611973-49-5
*/
{
  int iz;
  int n = 2;
  double *R;
  double **S;

  printf ( "\n" );
  printf ( "replace_row_test():\n" );
  printf ( "  replace_row() replaces a row in the simplex matrix.\n" );

  make_matrix ( S, n + 1, n );

  S[0][0] = 0.0;
  S[0][1] = 0.0;
  S[1][0] = 1.0;
  S[1][1] = 0.0;
  S[2][0] = 0.0;
  S[2][1] = 3.0;

  printf ( "\n" );
  printf ( "  Simplex coordinates:\n" );
  print_simplex ( S, n );

  make_vector ( R, n );
  R[0] = 1.0;
  R[1] = 2.0;
  printf ( "\n" );
  printf ( "  Replacement row R: " );
  print_row ( R, n );
  printf ( "\n" );

  iz = 2;
  replace_row ( S, iz, &R );
  printf ( "\n" );
  printf ( "  Simplex after replacing row %d by R:\n", iz );
  print_simplex ( S, n );

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

void shrink_test ( ) 

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

    shrink_test() tests shrink().

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    13 May 2024

  Author:

    John Burkardt.

  Reference:

    Rouben Rostamian,
    Programming Projects in C 
    for Students of Engineering, Science, and Mathematics,
    SIAM, 2014,
    ISBN: 978-1-611973-49-5
*/
{
  int ia;
  int n = 2;
  double **S;

  printf ( "\n" );
  printf ( "shrink_test():\n" );
  printf ( "  shrink() shrinks a simplex towards a given vertex ia.\n" );

  make_matrix ( S, n + 1, n );

  S[0][0] = 0.0;
  S[0][1] = 0.0;
  S[1][0] = 1.0;
  S[1][1] = 0.0;
  S[2][0] = 0.0;
  S[2][1] = 3.0;

  printf ( "\n" );
  printf ( "  Simplex coordinates:\n" );
  print_simplex ( S, n );

  ia = 1;
  shrink ( S, n, ia );

  printf ( "\n" );
  printf ( "  Simplex coordinates after shrinking towards node 1:\n" );
  print_simplex ( S, n );

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

void test_result ( double *a, int n, int ia, int iy, int iz )

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

    test_result() tests the user's "rank_vertices()" function on a test vector.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    13 May 2024

  Author:

    Original C version by Rouben Rostamian.
    This version by John Burkardt.

  Reference:

    Rouben Rostamian,
    Programming Projects in C 
    for Students of Engineering, Science, and Mathematics,
    SIAM, 2014,
    ISBN: 978-1-611973-49-5
*/
{
  int i;
/* 
  Are the indices distinct? 
*/
  if ( ( n == 2 && iy == iz ) || (n > 2 && (ia == iy || iy == iz || iz == ia ) ) )
  {
    printf ( "vector: " );
    print_row ( a, n );
    printf ( " => " );
    printf ( "ia = %u, iy = %u, iz = %u\n", ia, iy, iz );
    printf ( "  The indices are not distinct!\n" );
    return;
  }
/*
  Are the values ordered? 
*/
  if ( ! ( a[ia] <= a[iy] && a[iy] <= a[iz] ) )
  {
    printf ( "vector: " );
    print_row ( a, n );
    printf ( "ia = %u, iy = %u, iz = %u\n", ia, iy, iz );
    printf ( "a[%u] <= a[%u] && a[%u] <= a[%u] failed!\n", ia, iy, iy, iz );
    printf ( "  The values are not properly ordered!\n" );
    return;
  }
/* 
  is ia a true min? 
*/
  for ( i = 0; i < n; i++ )
  {
    if ( i != ia && a[i] < a[ia] )
    {
      printf ( "vector: " );
      print_row ( a, n );
      printf ( "ia = %u, iy = %u, iz = %u\n", ia, iy, iz );
      printf ( "ia = %u but a[%u] < a[%u]!\n", ia, i, ia );
      printf ( "  ia does not index a true minimum!\n" );
      return;
    }
  }
/*
  is iz a true max? 
*/
  for ( i = 0; i < n; i++ )
  {
    if ( i != iz && a[i] > a[iz] )
    {
      printf ( "vector: " );
      print_row ( a, n );
      printf ( "ia = %u, iy = %u, iz = %u\n", ia, iy, iz );
      printf ( "iz = %u but a[%u] > a[%u]!\n", iz, i, iz );
      printf ( "  iz does not index a true maximum!\n" );
      return;
    }
  }
/* 
  is iy correct? 
*/
  for ( i = 0; i < n; i++ )
  {
    if ( i != iz && i != iy && a[i] > a[iy] )
    {
      printf ( "iy = %u but a[%u] > a[%u]!\n", iy, i, iy );
      printf ( "  iy is not correct!\n" );
      return;
    }
  }
/*
  This vector was handled correctly!
*/
  printf ( "OK!: vector: " );
  print_row ( a, n );
  printf ( " => " );
  printf ( "ia = %u, iy = %u, iz = %u\n", ia, iy, iz );

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

void timestamp ( )

/******************************************************************************/
/*
  Purpose:
 
    timestamp() prints the current YMDHMS date as a time stamp.

  Example:

    17 June 2014 09:45:54 AM

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    01 May 2021

  Author:

    John Burkardt
*/
{
# 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 );

  printf ( "%s\n", time_buffer );

  return;
# undef TIME_SIZE
}
/******************************************************************************/

void transform_test ( ) 

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

    transform_test() tests transform().

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    13 May 2024

  Author:

    John Burkardt.

  Reference:

    Rouben Rostamian,
    Programming Projects in C 
    for Students of Engineering, Science, and Mathematics,
    SIAM, 2014,
    ISBN: 978-1-611973-49-5
*/
{
  double beta;
  int i;
  double *P;
  double *Q;
  double *R;
  int n = 2;
 
  printf ( "\n" );
  printf ( "transform_test():\n" );
  printf ( "  transform() carries out the transform R=(1-beta)*P+beta*Q.\n" );

  make_vector ( P, n );
  make_vector ( Q, n );
  make_vector ( R, n );

  P[0] = 0.0;
  P[1] = 1.0;

  Q[0] = 2.0;
  Q[1] = 3.0;

  printf ( "\n" );
  printf ( "  P:  " );
  print_row ( P, n );
  printf ( "\n" );

  printf ( "\n" );
  printf ( "  Q:  " );
  print_row ( Q, n );
  printf ( "\n" );

  printf ( "\n" );
  printf ( "  beta   R\n" );
  printf ( "\n" );

  for ( i = -2; i < 6; i++ )
  {
    beta = i / 4.0;
    transform ( P, Q, n, beta, R );
    printf ( "  %6.2f  ( %10.4f,  %10.4f )\n", beta, R[0], R[1] );
  }
  return;
}