'***********************************************************
'*   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.122685894370079   1.786915332078934                 *
'*   4.166477359831333   1.682698175311089                 *
'*   4.142454409040511   1.741175873205066                 *
'*                                                         *
'*  Best NDIM+1 mimimum values:                            *
'*  -.21723362651456                                       *
'*  -.2172336281071481                                     *
'*  -.2172336271378329                                     *
'*                                                         *
'* ------------------------------------------------------- *
'* 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].                                 *
'*                                                         *
'*                    Basic Release By J-P Moreau, Paris.  *
'*                            (www.jpmoreau.fr)            *
'***********************************************************
'PROGRAM TEST_AMOEBA

DEFDBL A-H, O-Z
DEFINT I-N

MP = 21
NP = 20  'Maximum value for NDIM=20

DIM P(MP, NP)
DIM Y(MP), PT(MP)

  NDIM = 2    '2 variables
  FTOL = .00000001#'User given tolerance

  'define NDIM+1 initial vertices (one by row)
  P(1, 1) = 1#: P(1, 2) = 2#
  P(2, 1) = -2#: P(2, 2) = -3#
  P(3, 1) = 4#: P(3, 2) = 2#

  'Initialize Y to the values of FUNC evaluated
  'at the NDIM+1 vertices (rows] of P
  FOR I = 1 TO NDIM + 1
    PT(1) = P(I, 1): PT(2) = P(I, 2)
    'Y(I)=FUNC(PT)
    GOSUB 500: Y(I) = FUNC
  NEXT I

  'call main subroutine
  GOSUB 1000 'call AMOEBA(P,Y,MP,NP,NDIM,FTOL,ITER)

  'print results
  CLS
  PRINT
  PRINT " Number of iterations: "; ITER
  PRINT
  PRINT " Best NDIM+1 points:"
  FOR I = 1 TO NDIM + 1
    FOR J = 1 TO NDIM
      PRINT " "; P(I, J);
    NEXT J
    PRINT
  NEXT I
  PRINT
  PRINT " Best NDIM+1 mimimum values:"
  FOR I = 1 TO NDIM + 1
    PRINT " "; Y(I)
  NEXT I
  PRINT

END 'of main program


'user defined function to minimize}
500 'FUNCTION FUNC(PT)
  R = SQR(PT(1) * PT(1) + PT(2) * PT(2))
  IF ABS(R) < 1E-12 THEN
    FUNC = 1#
  ELSE
    FUNC = SIN(R) / R
  END IF
RETURN


1000 'Subroutine AMOEBA(P, Y, NDIM, FTOL, 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:  1
  NMAX = 20: ALPHA = 1#: BETA = .5: GAMMA = 2#: ITMAX = 500
' Expected maximum number of dimensions, three parameters which define
' the expansions and contractions, and maximum allowed number of
' iterations.
DIM PR(MP), PRR(MP), PBAR(MP)
  MPTS = NDIM + 1
  ITER = 0
1 ILO = 1
  IF Y(1) > Y(2) THEN
    IHI = 1
    INHI = 2
  ELSE
    IHI = 2
    INHI = 1
  END IF
  FOR I = 1 TO MPTS
    IF Y(I) < Y(ILO) THEN ILO = I
    IF Y(I) > Y(IHI) THEN
      INHI = IHI
      IHI = I
    ELSEIF Y(I) > Y(INHI) THEN
      IF I <> IHI THEN INHI = I
    END IF
  NEXT I
' Compute the fractional range from highest to lowest and return if
' satisfactory.
  RTOL = 2# * ABS(Y(IHI) - Y(ILO)) / (ABS(Y(IHI)) + ABS(Y(ILO)))
  IF RTOL < FTOL THEN RETURN   'normal return
  IF ITER = ITMAX THEN
    PRINT " Amoeba exceeding maximum iterations."
    RETURN
  END IF
  ITER = ITER + 1
  FOR J = 1 TO NDIM
    PBAR(J) = 0#
  NEXT J
  FOR I = 1 TO MPTS
    IF I <> IHI THEN
      FOR J = 1 TO NDIM
	PBAR(J) = PBAR(J) + P(I, J)
      NEXT J
    END IF
  NEXT I
  FOR J = 1 TO NDIM
    PBAR(J) = PBAR(J) / (1# * NDIM)
    PR(J) = (1# + ALPHA) * PBAR(J) - ALPHA * P(IHI, J)
  NEXT J
  'YPR=FUNC(PR)
  PT(1) = PR(1): PT(2) = PR(2): GOSUB 500: YPR = FUNC
  IF YPR <= Y(ILO) THEN
    FOR J = 1 TO NDIM
      PRR(J) = GAMMA * PR(J) + (1# - GAMMA) * PBAR(J)
    NEXT J
    'YPRR=FUNC(PRR)
    PT(1) = PRR(1): PT(2) = PRR(2): GOSUB 500: YPRR = FUNC
    IF YPRR < Y(ILO) THEN
      FOR J = 1 TO NDIM
	P(IHI, J) = PRR(J)
      NEXT J
      Y(IHI) = YPRR
    ELSE
      FOR J = 1 TO NDIM
	P(IHI, J) = PR(J)
      NEXT J
      Y(IHI) = YPR
    END IF
  ELSEIF YPR >= Y(INHI) THEN
    IF YPR < Y(IHI) THEN
      FOR J = 1 TO NDIM
	P(IHI, J) = PR(J)
      NEXT J
      Y(IHI) = YPR
    END IF
    FOR J = 1 TO NDIM
      PRR(J) = BETA * P(IHI, J) + (1# - BETA) * PBAR(J)
    NEXT J
    'YPRR=FUNC(PRR)
    PT(1) = PRR(1): PT(2) = PRR(2): GOSUB 500: YPRR = FUNC
    IF YPRR < Y(IHI) THEN
      FOR J = 1 TO NDIM
	P(IHI, J) = PRR(J)
      NEXT J
      Y(IHI) = YPRR
    ELSE
      FOR I = 1 TO MPTS
	IF I <> ILO THEN
	  FOR J = 1 TO NDIM
	    PR(J) = .5 * (P(I, J) + P(ILO, J))
	    P(I, J) = PR(J)
	  NEXT J
	  'Y(I)=FUNC(PR)
	  PT(1) = PR(1): PT(2) = PR(2): GOSUB 500: Y(I) = FUNC
	END IF
      NEXT I
    END IF
  ELSE
    FOR J = 1 TO NDIM
      P(IHI, J) = PR(J)
    NEXT J
    Y(IHI) = YPR
  END IF
  GOTO 1
RETURN

'end of file tamoeba.bas