IMPLEMENTATION MODULE OptimLib3;
(*------------------------------------------------------------------------*)
(* Minimirungsroutinen fuer Funktionen mit bekannter 1. Ableitung. *)
(* Minimization of functions with derivative *)
(*------------------------------------------------------------------------*)
(* Letzte Veraenderung *)
(* *)
(* 11.04.18, MRi: Erstellen der ersten Version von MinInDer *)
(* 12.04.18, MRi: Kosmetische Aenderungen *)
(* 23.04.18, MRi: Erstellen der ersten Versionen von VMMin und CGMin *)
(* in enger Anlehnung an die Pascal-Quellen *)
(* 23.04.18, MRi: Erstellen der ersten Versionen von VMMin und CGMin *)
(* mit uebergebenden Funktionen *)
(* 25.04.18, MRi: Umstellung von von VMMin und CGMin auf offene Felder *)
(* und dynamische Speicherverwaltung *)
(* 04.04.18, MRi: Pruefen auf NAN und INF in VMMin und CGMin eingefuegt *)
(*------------------------------------------------------------------------*)
(* Testumgebung in TstMinDer fuer MinInDer *)
(* Testumgebung in TstOptimG fuer VMMin,CGMin *)
(*------------------------------------------------------------------------*)
(* Offene Punkte *)
(* *)
(* - Testen gegen andere Funktionen *)
(* - Weitere Schnittstellenanpassungen *)
(* - Vollstaendiges umstellen auf offene Felder *)
(* - Doku und Referenzen *)
(* - CGMin mit Beale Sorenson gibt bei Testfall 12 NAN zurueck - wieso ? *)
(*------------------------------------------------------------------------*)
(* Diese Modul enthaelt Quellcode der aus der NUMAL-Bibiliothek *)
(* (Numerical Algol library) Stichting CWI abgeleitet wurde. *)
(* *)
(* This module contains code derived from NUMAL, Numerical Algol library *)
(* (Stichting CWI) *)
(* *)
(* [1] Fletcher, R.; "A new approach to variable metric *)
(* algorithms", Computer Journal, 13/3 pp. 317-322 (1970) *)
(* [2] Nash, J.C.; "Compact Numerical Methods for Computers", *)
(* Second Edition, Adam Hilger, Bristol, UK (1990) *)
(* [3] Fletcher,R.; "Practical Methods of Optimization", Second Edition, *)
(* Wiley, Hoboken, US (1987) ISBN 0471915475 *)
(* [4] Werner, J.; "��ber die globale Konvergenz von Variable-Metric *)
(* Verfahren mit nichtexakter Schrittweitenbestimmung", *)
(* Numer. Math., 31 pp. 321���334 (1978) *)
(* [5] Kelley, C.T.; "Iterative methods for optimization", Frontiers in *)
(* Appl. Math. , 18 , SIAM (1999) *)
(* [6] Dennis, J.E.; Schnabel, R.B.; "Numerical Methods for Nonlinear *)
(* Equations and Unconstrained Optimization", Classics in Applied *)
(* Math. 16, SIAM (1996) *)
(* see www.encyclopediaofmath.org (Broyden-Fletcher-Goldfarb-Shanno) *)
(*------------------------------------------------------------------------*)
(* Implementation : Michael Riedl *)
(* Licence : GNU Lesser General Public License (LGPL) *)
(*------------------------------------------------------------------------*)
(* $Id: OptimLib3.mod,v 1.3 2018/04/26 09:54:56 mriedl Exp mriedl $ *)
FROM SYSTEM IMPORT TSIZE;
FROM Storage IMPORT ALLOCATE,DEALLOCATE;
FROM Errors IMPORT FatalError;
FROM TestReal IMPORT IsReal8NaN,IsReal8INF;
FROM DynMat IMPORT AllocMat,DeAllocMat;
FROM Deklera IMPORT PMATRIX;
FROM LongMath IMPORT sqrt;
FROM LMathLib IMPORT MachEps,MachEpsR4,GoldSec,CGold,Funktion1,FunktionN,sqr;
IMPORT TIO;
FROM MatLib IMPORT ENorm;
FROM LibDBlasM2 IMPORT dcopy,ddot;
CONST reltest = 10.0; (* A relative size used to check equality of *)
(* numbers. Numbers x and y are considered equal if *)
(* the floating-point representation of reltest+x *)
(* equals that of reltest+y. *)
debug = FALSE;
TermOut = TRUE;
PROCEDURE MinInDer(VAR a,b : LONGREAL;
f : Funktion1;
df : Funktion1;
ftol : Funktion1;
VAR x,y : LONGREAL;
VAR xmin : LONGREAL;
VAR cnt : CARDINAL);
VAR c,e,d,ba,z,p,q,s,sgn : LONGREAL;
fa,fb,fu,dfa,dfb,dfu,tol : LONGREAL;
BEGIN
cnt := 0;
IF (b <= a) THEN c:=a; a:=b; b:=c; END;
x:=a; y:=b;
fa := f(a); dfa := df(a); INC(cnt);
fb := f(b); dfb := df(b); INC(cnt);
c := (3.0 - sqrt(5.0)) / 2.0; (* Goldener Schnitt *)
d := b - a; e := 2.0*d; z := 2.0*e;
LOOP
ba := b - a; tol := ftol(x);
IF (ba < 3.0*tol) THEN EXIT; END;
IF ABS(dfa) <= ABS(dfb) THEN
x := a; sgn := 1.0;
ELSE
x := b; sgn := -1.0;
END;
IF (dfa <= 0.0) AND (dfb >= 0.0) THEN
z := 3.0*(fa - fb)/ba + dfa + dfb;
s := sqrt(sqr(z) - dfa * dfb);
IF sgn = 1.0 THEN p := dfa - s - z ELSE p := dfb + s - z; END;
p := p * ba;
q := dfb - dfa + 2.0*s; z := e; e := d;
IF ABS(p) <= ABS(q) * tol THEN d := tol*sgn ELSE d := -p/q; END;
ELSE
d:= ba;
END;
IF (ABS(d) >= ABS(0.5*z)) OR (ABS(d) > 0.5*ba) THEN
e := ba; d := c*ba*sgn;
END;
x:=x + d;
fu := f(x); dfu := df(x); INC(cnt);
IF (dfu >= 0) OR ((fu >= fa) AND (dfa <= 0)) THEN
b := x; fb := fu; dfb := dfu;
ELSE
a := x; fa := fu; dfa := dfu;
END;
END; (* LOOP *)
IF fa <= fb THEN
x := a; y := b; xmin := fa;
ELSE
x := b; y := a; xmin := fb;
END;
END MinInDer;
PROCEDURE VMMin( Func : FunktionN;
n : CARDINAL;
Gradient : TGradient;
VAR X : ARRAY OF LONGREAL;
VAR Bvec : ARRAY OF LONGREAL;
VAR Fmin : LONGREAL;
VAR funcount : CARDINAL;
VAR gradcount : CARDINAL;
VAR iFehl : INTEGER);
(*----------------------------------------------------------------*)
(* Variable metric function minimiser, see Ref [1] and [2] *)
(*----------------------------------------------------------------*)
CONST stepredn = 0.2; (* factor to reduce stepsize in line search *)
acctol = 0.0001; (* acceptable point tolerance *)
MaxIter = 30000;
VAR B : PMATRIX;
c,g,t : POINTER TO ARRAY [0..MAX(INTEGER)-1] OF LONGREAL;
PROCEDURE FreeMem();
BEGIN
DEALLOCATE(c,(n+1)*TSIZE(LONGREAL));
DEALLOCATE(g,(n+1)*TSIZE(LONGREAL));
DEALLOCATE(t,(n+1)*TSIZE(LONGREAL));
DeAllocMat(B,n,n);
END FreeMem;
VAR accpoint : BOOLEAN;
count : CARDINAL;
D1,D2,f,s : LONGREAL;
gradproj : LONGREAL;
i,j,ilast : CARDINAL;
steplength : LONGREAL;
notcomp : BOOLEAN;
Iter : CARDINAL;
Fehler : BOOLEAN;
BEGIN
iFehl:=0;
f := Func(Bvec,n);
IF (f = MAX(LONGREAL)) THEN
TIO.WrStr(
'**** Function cannot be evaluated at initial parameters ****');
TIO.WrLn;
iFehl := -1;
RETURN;
END;
ALLOCATE(c,(n+1)*TSIZE(LONGREAL));
ALLOCATE(g,(n+1)*TSIZE(LONGREAL));
ALLOCATE(t,(n+1)*TSIZE(LONGREAL));
AllocMat(B,n,n);
IF (c = NIL) OR (g = NIL) OR (t = NIL) OR (B = NIL) THEN
FatalError("Kein Freispeicher vorhanden (CGMin)");
END;
Fmin := f;
funcount := 1;
Gradient(Bvec,n,g^); gradcount:=1;
ilast := gradcount;
Iter:=0;
REPEAT
INC(Iter);
IF (ilast = gradcount) THEN
FOR i:=1 TO n DO
FOR j:=1 TO n DO B^[i]^[j]:=0.0; END; B^[i]^[i]:=1.0;
END;
END;
IF debug THEN
TIO.WrCard(gradcount,1); TIO.WrChar(' ');
TIO.WrCard(funcount,1);
TIO.WrLngReal(Fmin,20,-9); TIO.WrLn; TIO.WrStr('parameters ');
FOR i:=0 TO n-1 DO TIO.WrLngReal(Bvec[i],10,5); TIO.WrChar(' '); END;
END;
FOR i:=0 TO n-1 DO
X[i] := Bvec[i];
c^[i] := g^[i];
END;
gradproj:=0.0;
FOR i:=0 TO n-1 DO
s:=0.0;
FOR j:=0 TO n-1 DO s:=s - B^[i+1]^[j+1]*g^[j]; END;
t^[i] := s; gradproj:=gradproj + s*g^[i];
END;
IF (gradproj < 0.0) THEN(* !! note change to floating point *)
steplength:=1.0;
accpoint:=FALSE;
REPEAT
count:=0;
FOR i:=0 TO n-1 DO
Bvec[i] := X[i] + steplength*t^[i];
IF ((reltest + X[i])=(reltest + Bvec[i])) THEN INC(count); END;
END;
IF (count < n) THEN
f := Func(Bvec,n); INC(funcount);
notcomp := (f = MAX(LONGREAL));
accpoint:=(NOT notcomp) AND
(f <= Fmin+gradproj*steplength*acctol);
IF NOT accpoint THEN
steplength := steplength*stepredn;
IF debug THEN TIO.WrChar('*'); END;
ELSE
IF debug THEN TIO.WrChar(' '); END;
END;
END;
UNTIL (count=n) OR accpoint;
IF (count < n) THEN
Fmin := f;
Gradient(Bvec,n,g^); INC(gradcount);
D1:=0.0;
FOR i:=0 TO n-1 DO
t^[i] := steplength*t^[i]; c^[i] := g^[i] - c^[i];
D1 := D1 + t^[i]*c^[i];
END;
IF (D1 > 0.0) THEN
D2:=0.0;
FOR i:=0 TO n-1 DO
s:=0.0;
FOR j:=0 TO n-1 DO s:=s + B^[i+1]^[j+1]*c^[j];END;
X[i] := s; D2:=D2 + s*c^[i];
END;
D2 := 1.0 + D2/D1;
FOR i:=0 TO n-1 DO
FOR j:=0 TO n-1 DO
B^[i+1]^[j+1] := B^[i+1]^[j+1] -
(t^[i]*X[j] + X[i]*t^[j] - D2*t^[i]*t^[j])/D1;
END;
END;
ELSE
TIO.WrStr(' update not possible'); TIO.WrLn;
ilast := gradcount;
(*
FreeMem();
iFehl:=-2; RETURN;
(****************)
*)
END;
ELSE
IF (ilast < gradcount) THEN
count:=0;
ilast := gradcount;
END;
END;
ELSE
TIO.WrStr('UPHILL SEARCH DIRECTION'); TIO.WrLn;
count:=0; (* !! order of statements *)
IF (ilast=gradcount) THEN count := n ELSE ilast := gradcount; END;
(* Resets Hessian inverse if it has not just *)
(* been set, otherwise forces a convergence. *)
END;
UNTIL ((count=n) AND (ilast=gradcount)) OR (Iter > MaxIter);
IF (Iter > MaxIter) THEN iFehl:=-3; END;
FreeMem();
(* Check if result is a valid number *)
i:=0; Fehler:=FALSE;
WHILE (i < n) AND NOT Fehler DO
Fehler := (IsReal8NaN(X[i]) OR IsReal8INF(X[i])); INC(i)
END;
IF Fehler THEN iFehl := 99; END;
END VMMin;
PROCEDURE CGMin( Func : FunktionN;
n : CARDINAL;
Gradient : TGradient;
VAR X : ARRAY OF LONGREAL;
VAR Bvec : ARRAY OF LONGREAL;
VAR Fmin : LONGREAL;
VAR funcount : CARDINAL;
VAR gradcount : CARDINAL;
VAR intol : LONGREAL;
setstep : LONGREAL;
methode : CARDINAL;
VAR iFehl : INTEGER);
(*----------------------------------------------------------------*)
(* Conjugate gradients function minimiser see Ref [2] *)
(*----------------------------------------------------------------*)
TYPE methodtype= (FletcherReeves, PolakRibiere, BealeSorenson);
CONST stepredn = 0.2; (* factor to reduce stepsize in line search *)
acctol = 0.0001; (* acceptable point tolerance *)
VAR i,count : CARDINAL;
cycle,cyclimit : CARDINAL;
f,tol : LONGREAL;
G1,G2,G3,gradproj : LONGREAL;
newstep,oldstep : LONGREAL;
steplength : LONGREAL;
accpoint,notcomp : BOOLEAN;
method : methodtype;
c,g,t : POINTER TO
ARRAY [0..MAX(INTEGER)-1] OF LONGREAL;
Fehler : BOOLEAN;
BEGIN
CASE methode OF
| 1: method := FletcherReeves;
| 2: method := PolakRibiere;
| 3: method := BealeSorenson;
ELSE
iFehl:=-99; RETURN;
END;
iFehl:=0;
ALLOCATE(c,(n+1)*TSIZE(LONGREAL));
ALLOCATE(g,(n+1)*TSIZE(LONGREAL));
ALLOCATE(t,(n+1)*TSIZE(LONGREAL));
IF (c = NIL) OR (g = NIL) OR (t = NIL) THEN
FatalError("Kein Freispeicher vorhanden (CGMin)");
END;
cyclimit := n;
IF (intol < 0.0) THEN intol := MachEps; END;
tol := intol*LFLOAT(n)*sqrt(intol);
steplength := 0.01; (* Wg. Compilerwarnung *)
TIO.WrStr('tolerance used in gradient test=');
TIO.WrLngReal(tol,12,-2); TIO.WrLn;
f := Func(Bvec,n);
IF (f = MAX(LONGREAL)) THEN
TIO.WrStr('**** Function cannot be evaluated at initial parameters ****');
TIO.WrLn;
iFehl := -1;
ELSE
Fmin := f;
funcount := 1;
gradcount := 0;
REPEAT
FOR i:=0 TO n-1 DO t^[i]:=0.0; c^[i]:=0.0; END;
cycle:=0;
oldstep:=1.0;
count:=0;
REPEAT
INC(cycle);
IF debug THEN
TIO.WrCard(gradcount,1); TIO.WrChar(' '); TIO.WrCard(funcount,1);
TIO.WrChar( ' '); TIO.WrLngReal(Fmin,12,-3); TIO.WrLn;
TIO.WrStr('parameters ');
FOR i:=1 TO n DO
TIO.WrLngReal(Bvec[i-1],10,-5); TIO.WrChar(' ');
IF (7 * (i DIV 7) = i) AND (i < n) THEN
TIO.WrLn;
END;
END;
TIO.WrLn;
END;
Gradient(Bvec,n,g^); INC(gradcount);
G1:=0.0; G2:=0.0;
FOR i:=0 TO n-1 DO
X[i] := Bvec[i];
CASE method OF
|FletcherReeves: G1:=G1 + sqr(g^[i]);
G2:=G2 + sqr(c^[i]);
|PolakRibiere : G1:=G1 + g^[i]*(g^[i] - c^[i]);
G2:=G2 + sqr(c^[i]);
|BealeSorenson : G1:=G1 + g^[i]*(g^[i] - c^[i]);
G2:=G2 + t^[i]*(g^[i] - c^[i]);
END;
c^[i] := g^[i];
END; (* FOR i *)
IF (G1 > tol) THEN
IF (G2 > 0.0) THEN G3 := G1/G2 ELSE G3 := 1.0; END;
gradproj:=0.0;
FOR i:=0 TO n-1 DO
t^[i]:=t^[i]*G3 - g^[i]; gradproj:=gradproj + t^[i]*g^[i];
END;
steplength := oldstep;
accpoint := FALSE;
REPEAT
count:=0;
FOR i:=0 TO n-1 DO
Bvec[i] := X[i] + steplength*t^[i];
IF (reltest + X[i]) = (reltest + Bvec[i]) THEN INC(count); END;
END;
IF (count < n) THEN
f := Func(Bvec,n); INC(funcount);
notcomp := (f = MAX(LONGREAL));
accpoint:=(NOT notcomp)
AND (f <= Fmin+gradproj*steplength*acctol);
IF NOT accpoint THEN
steplength:=steplength*stepredn;
IF debug THEN TIO.WrChar('*'); END;
ELSE
IF debug THEN TIO.WrChar(' '); END;
END;
END;
UNTIL (count = n) OR accpoint;
IF (count < n) THEN
newstep := 2.0*((f - Fmin) - gradproj*steplength);
IF (newstep > 0) THEN
newstep := -gradproj*sqr(steplength)/newstep;
FOR i:=0 TO n-1 DO
Bvec[i] := X[i] + newstep*t^[i];
END;
Fmin := f;
f := Func(Bvec,n); INC(funcount);
IF (f < Fmin) THEN
Fmin := f;
IF debug THEN TIO.WrStr(' i< '); END;
ELSE
IF debug THEN TIO.WrStr(' i> '); END;
FOR i:=0 TO n-1 DO Bvec[i] := X[i] + steplength*t^[i]; END;
END;
END;
END;
END;
oldstep := setstep*steplength;
IF (oldstep > 1.0) THEN oldstep:=1.0;END;
UNTIL (count = n) OR (G1 <= tol) OR (cycle = cyclimit);
UNTIL (cycle = 1) AND ((count=n) OR (G1 <= tol));
END;
DEALLOCATE(c,(n+1)*TSIZE(LONGREAL));
DEALLOCATE(g,(n+1)*TSIZE(LONGREAL));
DEALLOCATE(t,(n+1)*TSIZE(LONGREAL));
(* Check if result is a valid number *)
i:=0; Fehler:=FALSE;
WHILE (i < n) AND NOT Fehler DO
Fehler := (IsReal8NaN(X[i]) OR IsReal8INF(X[i])); INC(i)
END;
IF Fehler THEN iFehl := 99; END;
END CGMin;
END OptimLib3.