'**************************************************************
'* Minimization of a Function FUNC of N Variables By Powell's *
'* Method Discarding the Direction of Largest Decrease *
'* ---------------------------------------------------------- *
'* SAMPLE RUN: Find a minimum of function F(x,y): *
'* F=Sin(R)/R, where R = Sqr(x*x+y*y). *
'* *
'* Number of iterations: 2 *
'* *
'* Minimum Value: -.2172336282112217 *
'* *
'* at point: 3.177320295222616 3.177320295222616 *
'* *
'* ---------------------------------------------------------- *
'* 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 TPOWELL
DEFDBL A-H, O-Z
DEFINT I-N
NP = 20
DIM P(NP), XI(NP, NP), P1(NP)
DIM PCOM(NP), XICOM(NP) 'PCOM,XICOM,NCOM are common variables
DIM XT(NP) 'for LINMIN and F1DIM only
N = 2
P(1) = 2#: P(2) = 2#
XI(1, 1) = 1#: XI(1, 2) = 1#
XI(2, 1) = 1#: XI(2, 2) = 1#
FTOL = .000001#
GOSUB 1000 'call POWELL(P,XI,N,NP,FTOL,ITER,FRET)
CLS
PRINT
PRINT " Number of iterations: "; ITER
PRINT
PRINT " Minimum value: "; FRET
PRINT
PRINT " at point: "; P(1); " "; P(2)
PRINT
END 'of main program
'user defined function to minimize
500 'FUNCTION FUNC(P1)
R1 = SQR(P1(1) * P1(1) + P1(2) * P1(2))
IF ABS(R1) < .000000000001# THEN
FUNC = 1#
ELSE
FUNC = SIN(R1) / R1
END IF
RETURN
600 'FUNCTION XMAX(a,b)
IF A >= B THEN
XMAX = A
ELSE
XMAX = B
END IF
RETURN
610 'FUNCTION XMIN(a,b)
IF A <= B THEN
XMIN = A
ELSE
XMIN = B
END IF
RETURN
650 'FUNCTION Sign(a,b)
IF B >= 0 THEN
Sign = ABS(A)
ELSE
Sign = -ABS(A)
END IF
RETURN
1000 'Subroutine POWELL(P, XI, N, FTOL, ITER, FRET)
'-----------------------------------------------------------
' Minimization of a function FUNC of N variables (FUNC is
' not an argument, it is a fixed function name). Input con-
' sists of an initial starting point P that is a vector of
' length N: an initial matrix XI whose logical dimensions
' are N by N, physical dimensions NP by NP, and whose columns
' contain the initial set of directions (usually the N unit
' vectors): and FTOL, the fractional tolerance in the func-
' tion value such that failure to decrease by more than this
' amount on one iteration signals doneness. On output, P is
' set to the best point found, XI is the then-current direc-
' tion set, FRET is the returned function value at P, and
' ITER is the number of iterations taken. The routine LINMIN
' is used.
'-----------------------------------------------------------
' Label: 1
ITMAX = 200
DIM PT(NP), PTT(NP), XIT(NP)
'FRET=FUNC(P)
P1(1) = P(1): P1(2) = P(2): GOSUB 500: FRET = FUNC
FOR J = 1 TO N
PT(J) = P(J): 'Save initial point
NEXT J
ITER = 0
1 ITER = ITER + 1
FP = FRET
IBIG = 0
DEL = 0# 'Will be the biggest function decrease.
FOR I = 1 TO N 'In each iteration, loop over all directions in the set.
FOR J = 1 TO N 'Copy the direction.
XIT(J) = XI(J, I)
NEXT J
FPTT = FRET
GOSUB 1500 'call LINMIN(P,XIT,N,FRET) (Minimize along it)
IF ABS(FPTT - FRET) > DEL THEN
DEL = ABS(FPTT - FRET)
IBIG = I
END IF
NEXT I
IF 2# * ABS(FP - FRET) <= FTOL * (ABS(FP) + ABS(FRET)) THEN RETURN'Termination criterion
IF ITER = ITMAX THEN
PRINT " Powell exceeding maximum iterations."
RETURN
END IF
FOR J = 1 TO N
PTT(J) = 2# * P(J) - PT(J)'Construct the extrapolated point and the average
XIT(J) = P(J) - PT(J) 'direction moved. Save the old starting point
PT(J) = P(J)
NEXT J
'FPTT=FUNC(PTT) (Function value at extrapolated point)
P1(1) = PTT(1): P1(2) = PTT(2): GOSUB 500: FPTT = FUNC
IF FPTT >= FP THEN GOTO 1: 'One reason not to use new direction.
T = 2# * (FP - 2# * FRET + FPTT) * (FP - FRET - DEL) * (FP - FRET - DEL) - DEL * (FP - FPTT) * (FP - FPTT)
IF T >= 0! THEN GOTO 1: 'Other reason not to use new direction.
GOSUB 1500 'call LINMIN(P,XIT,N,FRET) 'Move to the minimum of the new direction.
FOR J = 1 TO N 'and save the new direction.
XI(J, IBIG) = XIT(J)
NEXT J
GOTO 1
RETURN
1500 'Subroutine LINMIN(P, XIT, N, FRET)
'----------------------------------------------------------
' Given an N dimensional point P and a N dimensional direc-
' tion XI, moves and resets P to where the function FUNC(P)
' takes on a minimum along the direction XIT from P, and
' replaces XIT by the actual vector displacement that P was
' moved. Also returns as FRET the value of FUNC at the
' returned location P. This is actually all accomplished by
' calling the routines MNBRAK and BRENT.
'----------------------------------------------------------
TOL = .0001#
NCOM = N
FOR J = 1 TO N
PCOM(J) = P(J)
XICOM(J) = XIT(J)
NEXT J
AX = 0#
BX = 1#
CX = 2#
GOSUB 2000 'call MNBRAK(AX,BX,CX,FA,FB,FC)
GOSUB 3000 'call BRENT(AX,BX,CX,TOL,XMIN)
FRET = BRENT
FOR J = 1 TO N
XIT(J) = XMIN * XIT(J)
P(J) = P(J) + XIT(J)
NEXT J
RETURN
1600 'FUNCTION F1DIM(X1)
FOR J = 1 TO NCOM
XT(J) = PCOM(J) + X1 * XICOM(J)
NEXT J
'F1DIM = FUNC(XT)
P1(1) = XT(1): P1(2) = XT(2): GOSUB 500
F1DIM = FUNC
RETURN
2000 'Subroutine MNBRAK(AX,BX,CX,FA,FB,FC)
'--------------------------------------------------------------------
'Given a Function F1DIM(X), and given distinct initial points AX and
'BX, this routine searches in the downhill direction (defined by the
'Function as evaluated at the initial points) and returns new points
'AX, BX, CX which bracket a minimum of the Function. Also returned
'are the Function values at the three points, FA, FB and FC.
'--------------------------------------------------------------------
' Label: 10
GOLD = 1.618034: GLIMIT = 100#: TINY = 9.999999999999999D-21
'The first parameter is the default ratio by which successive intervals
'are magnified: the second is the maximum magnification allowed for
'a parabolic-fit step.
'FA=F1DIM(AX)
X1 = AX: GOSUB 1600: FA = F1DIM
'FB=F1DIM(BX)
X1 = BX: GOSUB 1600: FB = F1DIM
IF FB > FA THEN
DUM = AX
AX = BX
BX = DUM
DUM = FB
FB = FA
FA = DUM
END IF
CX = BX + GOLD * (BX - AX)
'FC=F1DIM(CX)
X1 = CX: GOSUB 1600: FC = F1DIM
10 IF FB >= FC THEN
R = (BX - AX) * (FB - FC)
Q = (BX - CX) * (FB - FA)
A = Q - R: B = TINY: GOSUB 600
A = XMAX: B = Q - R: GOSUB 650
U = BX - ((BX - CX) * Q - (BX - AX) * R) / (2# * Sign)
ULIM = BX + GLIMIT * (CX - BX)
IF (BX - U) * (U - CX) > 0! THEN
'FU=F1DIM(U)
X1 = U: GOSUB 1600: FU = F1DIM
IF FU < FC THEN
AX = BX
FA = FB
BX = U
FB = FU
GOTO 10
ELSEIF FU > FB THEN
CX = U
FC = FU
GOTO 10
END IF
U = CX + GOLD * (CX - BX)
'FU=F1DIM(U)
X1 = U: GOSUB 1600: FU = F1DIM
ELSEIF (CX - U) * (U - ULIM) > 0! THEN
'FU=F1DIM(U)
X1 = U: GOSUB 1600: FU = F1DIM
IF FU < FC THEN
BX = CX
CX = U
U = CX + GOLD * (CX - BX)
FB = FC
FC = FU
'FU=F1DIM(U)
X1 = U: GOSUB 1600: FU = F1DIM
END IF
ELSEIF (U - ULIM) * (ULIM - CX) >= 0! THEN
U = ULIM
'FU=F1DIM(U)
X1 = U: GOSUB 1600: FU = F1DIM
ELSE
U = CX + GOLD * (CX - BX)
'FU=F1DIM(U)
X1 = U: GOSUB 1600: FU = F1DIM
END IF
AX = BX
BX = CX
CX = U
FA = FB
FB = FC
FC = FU
GOTO 10
END IF
RETURN
3000 'Function BRENTAX,BX,CX,TOL,XMIN)
'-------------------------------------------------------------------
'Given a function F1DIM, and a bracketing triplet of abscissas
'AX,BX,CX (such that BX is between AX and CX, and F(BX) is less
'than both F(AX) and F(CX)), this routine isolates the minimum
'to a fractional precision of about TOL using Brent's method.
'The abscissa of the minimum is returned in XMIN, and the minimum
'function value is returned as BRENT, the returned function value.
'-------------------------------------------------------------------
'Labels: 11,2,3
ITMAX = 100: CGOLD = .381966#: ZEPS = 1E-10
'Maximum allowed number of iterations: golden ratio: and a small
'number which protects against trying to achieve fractional accuracy
'for a minimum that happens to be exactly zero.
'A=MIN(AX,CX)
IF AX <= CX THEN
A = AX
ELSE
A = CX
END IF
'B=MAX(AX,CX)
IF AX >= CX THEN
B = AX
ELSE
B = CX
END IF
V = BX
W = V
X = V
E = 0#
'FX=F1DIM(X)
X1 = X: GOSUB 1600: FX = F1DIM
FV = FX
FW = FX
FOR ITER1 = 1 TO ITMAX 'main loop
XM = .5 * (A + B)
TOL1 = TOL * ABS(X) + ZEPS
TOL2 = 2# * TOL1
IF ABS(X - XM) <= (TOL2 - .5 * (B - A)) THEN GOTO 3'Test for done here
IF ABS(E) > TOL1 THEN 'Construct a trial parabolic fit
R = (X - W) * (FX - FV)
Q = (X - V) * (FX - FW)
P = (X - V) * Q - (X - W) * R
Q = .2 * (Q - R)
IF Q > 0# THEN P = -P
Q = ABS(Q)
ETEMP = E
E = D
IF (ABS(P) >= ABS(.5 * Q * ETEMP)) OR (P <= Q * (A - X)) OR (P >= Q * (B - X)) THEN GOTO 11
' The above conditions determine the acceptability of the
' parabolic fit. Here it is o.k.:
D = P / Q
U = X + D
IF XM - X >= 0 THEN
Sign = ABS(TOL1)
ELSE
Sign = -ABS(TOL1)
END IF
IF (U - A < TOL2) OR (B - U < TOL2) THEN D = Sign
GOTO 2
END IF
11 IF X >= XM THEN
E = A - X
ELSE
E = B - X
END IF
D = CGOLD * E
2 IF ABS(D) >= TOL1 THEN
U = X + D
ELSE
IF D >= 0 THEN
Sign = ABS(TOL1)
ELSE
Sign = -ABS(TOL1)
END IF
U = X + Sign
END IF
'FU=F1DIM(U) (This the one function evaluation per iteration)
X1 = U: GOSUB 1600: FU = F1DIM
IF FU <= FX THEN
IF U >= X THEN
A = X
ELSE
B = X
END IF
V = W
FV = FW
W = X
FW = FX
X = U
FX = FU
ELSE
IF U < X THEN
A = U
ELSE
B = U
END IF
IF (FU <= FW) OR (W = X) THEN
V = W
FV = FW
W = U
FW = FU
ELSEIF (FU <= FV) OR (V = X) OR (V = W) THEN
V = U
FV = FU
END IF
END IF
NEXT ITER1
PRINT " Brent exceed maximum iterations."
3 XMIN = X 'exit section
BRENT = FX
RETURN
'end of file tpowell.bas