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

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

# define REFLECT  1.0
# define EXPAND   2.0
# define CONTRACT 0.5
# define SHRINK   0.5

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

int done ( double **s, int n, double *y, int ia, int iz, double err2 )

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

    done() determines whether the algorithm has completed.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    11 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

  Input:

    double **s: the (n+1)*n simplex matrix.

    int n: the spatial dimension.

    double y[n+1]: the objective function values at the vertices.

    int ia, iz: the indices of the vertices with lowest and highest values.

    double err2: the value ( h * tol )^2;

  Output:

    int done: 1 if done, 0 otherwise.
*/
{
  int j;
  double t;

  if ( err2 < fabs ( y[iz] - y[ia] ) )
  {
    return 0;
  }

  t = 0.0;
  for ( j = 0; j < n; j++ )
  {
    t = t + pow ( s[iz][j] - s[ia][j], 2 );
  }

  if ( err2 < t )
  {
    return 0;
  }
  return 1;
}
/******************************************************************************/

void get_centroid ( double **s, int n, int iz, double *C )

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

    get_centroid() computes the centroid of the face opposite vertex x(iz).

  Discussion:

    For our intended 2-dimensional case, this is simply computing the
    midpoint of the edge opposite to vertex x(iz).

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    11 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

  Input:

    double **s: the (n+1)*n simplex matrix.

    int n: the spatial dimension.

    int iz: the index of a vertex of the simplex.

  Output:

    double *C: the coordinates of the face of the simplex opposite vertex iz.
*/
{
  int i;
  int j;

  for ( j = 0; j < n; j++ )
  {
    C[j] = 0.0;
  }

  for ( i = 0; i < n + 1; i++ )
  {
    for ( j = 0; j < n; j++ )
    {
      if ( i != iz )
      {
        C[j] = C[j] + s[i][j];
      }
    }
  }

  for ( j = 0; j < n; j++ )
  {
    C[j] = C[j] / n;
  }
  return;
}
/******************************************************************************/

int nelder_mead ( struct nelder_mead *nm )

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

    nelder_mead() seeks to minimize a function using the Nelder Mead method.

  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

  Input:

    struct nelder_mead *nm:

  Output:

    struct nelder_mead *nm:
*/
{
  double *C;
  int fevalcount = 0;
  double h = nm->h;
  int i;
  int ia;
  int iy;
  int iz;
  int j;
  int n = nm->n;
  double *Pc;
  double *Pe;
  double *Pr;
  double **s = nm->s;
  int simplex_to_be_freed = 0;
  double tol = nm->tol;
  double *y;
  double yc;
  double ye;
  double yr;

  double err2 = ( h * tol ) * ( h * tol );

  make_vector ( y, n + 1 );
  make_vector ( Pr, n );
  make_vector ( Pe, n );
  make_vector ( Pc, n );
  make_vector ( C, n );

  if ( s == NULL )
  {
    make_matrix ( s, n + 1, n );
    simplex_to_be_freed = 1;

    for ( i = 0; i < n + 1; i++ )
    {
      for ( j = 0; j < n; j++ )
      {
        if ( i == j + 1 )
        {
          s[i][j] = nm->x[j] + h;
        }
        else
        {
          s[i][j] = nm->x[j];
        }
      }
    }
  }

  for ( i = 0; i < n + 1; i++ )
  {
    y[i] = nm->f ( s[i], n, nm->params );
  }
  fevalcount = fevalcount + 1;
/*
  Loop.
*/
  while ( fevalcount <= nm->maxevals )
  {
    rank_vertices ( y, n + 1, &ia, &iy, &iz );

    if ( done ( s, n, y, ia, iz, err2 ) )
    {
      nm->minval = y[ia];
/*
?????
*/
/*
      nm->x = s[ia];
*/
      for ( j = 0; j < n; j++ )
      {
        nm->x[j] = s[ia][j];
      }
/*
  COPY BEST VERTEX INTO NM->x
*/
      break;
    }

    get_centroid ( s, n, iz, C );
    transform ( C, s[iz], n, -REFLECT, Pr );
    yr = nm->f ( Pr, n, nm->params );
    fevalcount = fevalcount + 1;
/*
  Case 1:
  Replace vertex x[iz] with the better of x[r] and x[e] (gotten by expansion).
*/
    if ( yr < y[ia] )
    {
      transform ( C, Pr, n, EXPAND, Pe );
      ye = nm->f ( Pe, n, nm->params );
      fevalcount = fevalcount + 1;
      if ( ye < yr )
      {
        replace_row ( s, iz, &Pe );
        y[iz] = ye;
      }
      else
      {
        replace_row ( s, iz, &Pr );
        y[iz] = yr;
      }
    }
/*
  Case 2:
  Replace the vertex x[iz] with x[r].
*/
    else if ( yr < y[iy] )
    {
      replace_row ( s, iz, &Pr );
      y[iz] = yr;
    }
/*
  Case 3 and 4:
  Push x[r] toward the centroid by a factor of beta, getting x[c].
  Apply shrink
*/
    else
    {
      transform ( C, Pr, n, CONTRACT, Pc );
      yc = nm->f ( Pc, n, nm->params );

      if ( yr < y[iz] )
      {
        if ( yc < yr )
        {
          replace_row ( s, iz, &Pc );
          y[iz] = yc;
        }
        else
        {
          shrink ( s, n, ia );
        }
      }
      else if ( y[iz] <= yr )
      {
        if ( yc < y[iz] )
        {
          replace_row ( s, iz, &Pc );
          y[iz] = yc;
        }
        else
        {
          shrink ( s, n, ia );
        }
      }
    }  
  }
/*
  Free memory.
*/
  free_vector ( Pc );
  free_vector ( Pe );
  free_vector ( Pr );
  free_vector ( y );

  if ( simplex_to_be_freed )
  {
    free_matrix ( s );
  }

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

void rank_vertices ( double *y, int m, int *ia, int *iy, int *iz )

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

    rank_vertices() identifies the best, middle, and worst vertices.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    06 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

  Input:

    int m: the length of the vector y.

    double y[m]: the function values at the simplex's vertices.

  Output:

    int *ia, *iy, *iz: the indices in y of the lowest, second highest,
    and highest values.  
*/
{
  int i;
/*
  Initial sort.
*/
  *ia = 0;

  if ( y[*ia] < y[1] )
  {
    *iy = 1;
  } 
  else
  {
    *ia = 1;
    *iy = 0;
  }

  if ( 2 < m )
  {
    if ( y[*iy] < y[2] )
    {
      *iz = 2;
    } 
    else if ( y[*ia] < y[2] )
    {
      *iz = *iy;
      *iy = 2;
    }
    else
    {
      *iz = *iy;
      *iy = *ia;
      *ia = 2;
    }
  }
  else
  {
    *iz = *iy;
    *iy = *ia;
  }

  for ( i = 3; i < m; i++ )
  {
/*
  New value is lower than previous lowest?
*/
    if ( y[i] < y[*ia] )
    {
      *ia = i;
    }
/*
  New value is higher than previous highest?
*/
    else if ( y[*iz] < y[i] )
    {
      *iy = *iz;
      *iz = i;
    }
/*
  New value is higher than previous next highest?
*/
    else if ( y[*iy] < y[i] )
    {
      *iy = i;
    }

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

void replace_row ( double **s, int i, double **r )

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

    replace_row() replaces row i in 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

  Input:

    double **s: the (n+1)*n simplex matrix.

    int i: the index of the row to be replaced.

    double **r: a pointer to the replacement row.

  Output:

    double **s: the simplex matrix after the row was replaced.
*/
{
  double *tmp;

  tmp = s[i];
  s[i] = *r;
  *r = tmp;

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

void shrink ( double **s, int n, int ia )

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

    shrink() shrinks the simplex towards the vertex ia.

  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

  Input:

    double **s: the (n+1)*n simplex matrix.

    int n: the spatial dimension.

    int ia: the index of the simplex vertex with lowest function value

  Output:

    double **s: the simplex matrix after shrinking.
*/
{
  double beta;
  int i;

  beta = SHRINK;

  for ( i = 0; i < n + 1; i++ )
  {
    if ( i != ia )
    {
      transform ( s[ia], s[i], n, beta, s[i] );
    }
  }
  return;
}
/******************************************************************************/

void transform ( double *P, double *Q, int n, double beta, double *R )

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

    transform() computes a transformed point R = (1-beta)*P + beta*Q.

  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

  Input:

    double P[n], Q[n]: two points.

    int n: the spatial dimension.

    double beta: the combination factor. 

  Output:

    double R[n]: the transformed point.  
*/
{
  int i;

  for ( i = 0; i < n; i++ )
  {
    R[i] = ( 1.0 - beta ) * P[i] + beta * Q[i];
  }
  return;
}