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