'***********************************************************
'* 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