/**********************************************************
*   Multidimensional minimization of a function FUNC(X)   *
*  where X is an NDIM-dimensional vector, by the downhill *
*  simplex method of Nelder and Mead.                     *
* ------------------------------------------------------- *
* SAMPLE RUN: Find a minimum of function F(x,y):          *
*             F=Sin(R)/R, where R = Sqrt(x*x+y*y).        *
*                                                         *
*  Number of iterations: 22                               *
*                                                         *
*  Best NDIM+1 points:                                    *
*   4.122686  1.786915                                    *
*   4.166477  1.682698                                    *
*   4.142454  1.741176                                    *
*                                                         *
*  Best NDIM+1 mimimum values:                            *
*   -0.2172336265                                         *
*   -0.2172336281                                         *
*   -0.2172336271                                         *
*                                                         *
* ------------------------------------------------------- *
* REFERENCE: "Numerical Recipes, The Art of Scientific    *
*             Computing By W.H. Press, B.P. Flannery,     *
*             S.A. Teukolsky and W.T. Vetterling,         *
*             Cambridge University Press, 1986"           *
*             [BIBLI 08].                                 *
*                                                         *
*                      C++ Release By J-P Moreau, Paris.  *
*                            (www.jpmoreau.fr)            *
**********************************************************/
#include <stdio.h>
#include <math.h>

#define   MP  22
#define   NP  21    //Maximum value for NDIM=20 

typedef   double MAT[MP][NP];

MAT P;
double Y[MP], PT[MP];
int    I,ITER,J,NDIM;
double FTOL;

// user defined function to minimize
double FUNC(double *P) {
  double R;
  R=sqrt(P[1]*P[1]+P[2]*P[2]);
  if (fabs(R) < 1e-12)
    return 1.0;
  else
    return sin(R)/R;
}


void AMOEBA(MAT P, double *Y, int NDIM, double FTOL, int *ITER) {
/*-------------------------------------------------------------------
! Multidimensional minimization of the function FUNC(X) where X is
! an NDIM-dimensional vector, by the downhill simplex method of
! Nelder and Mead. Input is a matrix P whose NDIM+1 rows are NDIM-
! dimensional vectors which are the vertices of the starting simplex
! (Logical dimensions of P are P(NDIM+1,NDIM); physical dimensions
! are input as P(NP,NP)). Also input is the vector Y of length NDIM
! +1, whose components must be pre-initialized to the values of FUNC
! evaluated at the NDIM+1 vertices (rows) of P; and FTOL the fractio-
! nal convergence tolerance to be achieved in the function value. On
! output, P and Y will have been reset to NDIM+1 new points all within
! FTOL of a minimum function value, and ITER gives the number of ite-
! rations taken.
!-------------------------------------------------------------------*/
// Label:  e1
const NMAX=20, ITMAX=500;
//Expected maximum number of dimensions, three parameters which define
// the expansions and contractions, and maximum allowed number of
//iterations.
  double PR[MP], PRR[MP], PBAR[MP];
  double ALPHA=1.0, BETA=0.5, GAMMA=2.0;
  int I,IHI,ILO,INHI,J,MPTS;
  double RTOL,YPR,YPRR;
  MPTS=NDIM+1;
  *ITER=0;
e1:ILO=1;
  if (Y[1] > Y[2]) {
    IHI=1;
    INHI=2;
  }
  else {
    IHI=2;
    INHI=1;
  }
  for (I=1; I<=MPTS; I++) {
    if (Y[I] < Y[ILO])  ILO=I;
    if (Y[I] > Y[IHI]) {
      INHI=IHI;
      IHI=I;
    }
    else if (Y[I] > Y[INHI])
      if (I != IHI)  INHI=I;
  }
//Compute the fractional range from highest to lowest and return if
//satisfactory.
  RTOL=2.0*fabs(Y[IHI]-Y[ILO])/(fabs(Y[IHI])+fabs(Y[ILO]));
  if (RTOL < FTOL)  return;  //normal exit
  if (*ITER == ITMAX) {
    printf(" Amoeba exceeding maximum iterations.\n");
    return;
  }
  *ITER= (*ITER) + 1;
  for (J=1; J<=NDIM; J++)  PBAR[J]=0.0;
  for (I=1; I<=MPTS; I++)
    if (I != IHI)
      for (J=1; J<=NDIM; J++)
        PBAR[J] += P[I][J];
  for (J=1; J<=NDIM; J++) {
    PBAR[J] /= 1.0*NDIM;
    PR[J]=(1.0+ALPHA)*PBAR[J] - ALPHA*P[IHI][J];
  }
  YPR=FUNC(PR);
  if (YPR <= Y[ILO]) {
    for (J=1; J<=NDIM; J++)
      PRR[J]=GAMMA*PR[J] + (1.0-GAMMA)*PBAR[J];
    YPRR=FUNC(PRR);
    if (YPRR < Y[ILO]) {
      for (J=1; J<=NDIM; J++) P[IHI][J]=PRR[J];
      Y[IHI]=YPRR;
    }
    else {
      for (J=1; J<=NDIM; J++) P[IHI][J]=PR[J];
      Y[IHI]=YPR;
	}
  }
  else if (YPR >= Y[INHI]) {
	if (YPR < Y[IHI]) {
      for (J=1; J<=NDIM; J++)  P[IHI][J]=PR[J];
      Y[IHI]=YPR;
    }
    for (J=1; J<=NDIM; J++)  PRR[J]=BETA*P[IHI][J] + (1.0-BETA)*PBAR[J];
    YPRR=FUNC(PRR);
    if (YPRR < Y[IHI]) {
      for (J=1; J<=NDIM; J++)  P[IHI][J]=PRR[J];
      Y[IHI]=YPRR;
    }
    else
      for (I=1; I<=MPTS; I++)
		if (I != ILO) {
		  for (J=1; J<=NDIM; J++) {
            PR[J]=0.5*(P[I][J] + P[ILO][J]);
	        P[I][J]=PR[J];
		  }
          Y[I]=FUNC(PR);
		}
  }
  else {
    for (J=1; J<=NDIM; J++)  P[IHI][J]=PR[J];
    Y[IHI]=YPR;
  }
  goto e1;
}


void main()  {

  NDIM=2;      // 2 variables
  FTOL=1e-8;   // User given tolerance

  //define NDIM+1 initial vertices (one by row)
  P[1][1]= 1.0; P[1][2]=2.0;
  P[2][1]=-2.0; P[2][2]=-3.0;
  P[3][1]= 4.0; P[3][2]=2.0;

//Initialize Y to the values of FUNC evaluated 
//at the NDIM+1 vertices (rows] of P
  for (I=1; I<=NDIM+1; I++) {
    PT[1]=P[I][1]; PT[2]=P[I][2];
    Y[I]=FUNC(PT);
  }

  //call main function
  AMOEBA(P,Y,NDIM,FTOL,&ITER);

  //print results
  printf("\n Number of iterations: %d\n\n", ITER);
  printf(" Best NDIM+1 points:\n");
  for (I=1 ; I<=NDIM+1; I++) {
    for (J=1; J<=NDIM; J++)  printf("  %f", P[I][J]);
    printf("\n");
  }
  printf("\n Best NDIM+1 mimimum values:\n");
  for (I=1; I<=NDIM+1; I++) printf(" %14.10f\n", Y[I]);
  printf("\n");

}
    
// end of file tamoeba.cpp