IMPLEMENTATION MODULE OptimLib1;
(*------------------------------------------------------------------------*)
(* Routinen fuer die Minimierung von ein- und mehrdimensionalen *)
(* Funktionen. *)
(* Module for finding the minimum of a function of one variable and *)
(* finding the minimum of a function of multible variables without *)
(* without derivatives. *)
(*------------------------------------------------------------------------*)
(* References: *)
(* *)
(* FMin : original Fortran version is due to Richard Brent and *)
(* George Forsythe (see desciption in definition section) *)
(* FMinBr : original C version is due to an unknown author, the *)
(* function is based on Ref. [5] *)
(* Zero : original Fortran versions are due to Richard Brent *)
(* GloMin : original Alogol 60 versions is due to Richard Brent *)
(* LocalMin : with modifications by John Burkardt *)
(* NelMin : original Fortran version of MelMin is due to Richard *)
(* O'Neill with modifications by John Burkardt *)
(* Praxis : original Pascal version of Praxis is due to Karl *)
(* Gegenfurtner based on a AlgolW version of Richard Brent *)
(* Fletcher- *)
(* Powell : routine is in tune based on ideas found in Gaussian 80 *)
(* Trudge : original Fortran version of Trudge is due to Michael *)
(* Depuis and coworkers of the Quantum Chemical package *)
(* Hondo *)
(*------------------------------------------------------------------------*)
(* Letzte Aenderung: *)
(* *)
(* FMinBr,FMin,GloMin,LocalMin,Zero *)
(* *)
(* 11.12.15, MRi: Erstellen der ersten Version von FMinBr *)
(* 19.06.16, MRi: Fehler in fminbir korrigiert. *)
(* 19.06.16, MRi: First version of FMin translated from Fortran source *)
(* 20.06.16, MRi: Erstellen der 1. Version von LocalMin *)
(* 21.06.16, MRi: Erstellen der 1. Version von Zero *)
(* 04.07.16, MRi: Erstellen der 1. Version von GloMin *)
(*------------------------------------------------------------------------*)
(* Nelder-Mead *)
(* *)
(* 13.06.16, MRi: Erstellen der ersten NelMin M2 Version (uebersetzbar) *)
(* 14.06.16, MRi: Korrekturen und erste Tests - noch nicht befriedigend *)
(* 15.06.16, MRi: Erste halbwegs brauchbaren Ergebnisse - aber noch nicht *)
(* vergleichbar mit F77 Original *)
(* 16.06.16, MRi: Fehler bei der initialisierung des Simplex gefunden *)
(* p^[nn,i] statt p^[n,i] in der ersten FOR Schleife, *)
(* Ergebisse scheinen nun der F77 Version zu entsprechen *)
(* Inneren LOOP durch REPEAT-Schleife ersetzt. *)
(*------------------------------------------------------------------------*)
(* Praxis *)
(* *)
(* 17.02.86, KGe: Initial Version (Turbo Pascal) *)
(* 20.06.16, MRi: "Umbau" der Pascal-Version um Konvertierung nach M2 *)
(* vorzubereiten. *)
(* 22.06.16, MRi: Erstellen der 1. Modula-2 Version *)
(* 23.06.16, MRi: Einfuehrung von MaxF,iFehl *)
(* Alle Testfunktionen werden durchlaufen, bei einigen *)
(* wird aber wegen Ueberschreiten der Maximalzahl der *)
(* Funktionsausrufe (MaxF) abgebrochen. Ergebisse zu- *)
(* friedenstellen - ... *)
(* 29.06.16, MRi: Verbesserte Kommentierung, Argumentenlisten fuer *)
(* interne Prozeduren *)
(*------------------------------------------------------------------------*)
(* Trudge *)
(* *)
(* June 95, MRi: Initial Version *)
(* 01.02.15, MRi: Error corrections, adopting to XDS,GM2 *)
(* 18.12.15, MRi: Falt in Trudge mit Func initialisiert *)
(* 04.06.16, MRi: Aufruf von LFLOAT durch VAL(LONGREAL,...) ersetzt *)
(* 08.06.16, MRi: Dokumentation im Definitionsmodul deutlich verbessert *)
(* Vorgabewerte fuer TolF,TolR,Noise eingefuegt *)
(*------------------------------------------------------------------------*)
(* Fletcher Powell *)
(* *)
(* Jan 95, MRi: Erste Version *)
(* 05.06.95, MRi: Modul INIT geloescht. Divisionen durch Null *)
(* abgefangen - alle Tests laufen durch, 6 Tests OK, *)
(* 4 Naeherungsweise, 4 ohne akzeptables Ergebniss. *)
(* Lokales MATaus geloescht. *)
(*------------------------------------------------------------------------*)
(* 19.06.16, MRi: Zusamennfassen der unter freier Lizenz stehenden *)
(* Module Trudge,FletcherPowell, NelderMead und FMinBr *)
(* 02.05.18, MRi: MinBrack aus dem Modul entfernt *)
(* 30.08.18, MRi: Anpassen von Formatanweisungen fuer WriteN *)
(*------------------------------------------------------------------------*)
(* Offene Punkte *)
(* *)
(* - In FletcherPowell wird auf die Matrix V in Fortran Speicherordnung *)
(* zugegriffen - koennte man korrigieren *)
(* - Bereinigen, IO durchsehen *)
(* - Testen *)
(* - Restart-Option in Praxis einfuegen (sichern von v,d und ???) *)
(* - Restart-Option in Trudge einfuegen *)
(* - Make routine output and desciption at least bi-lingual *)
(* - Literatur-Referenezen ausbauen *)
(*------------------------------------------------------------------------*)
(* Test routines for all 1 dimensional optimization routines are provided *)
(* in Module "TestMin". Test routines for all multiple demensional *)
(* optimization routines are given in Modul "TestOpt" together with *)
(* some test problems provided either in TestOpt or in module ObjFunkt *)
(*------------------------------------------------------------------------*)
(* References *)
(* *)
(* [5] Forsythe, George; Malcolm, Michael; Moler, Cleve; "Computer *)
(* methods for mathematical computations.", Prentice Hall, Upper *)
(* Saddle River, US (1977) *)
(*------------------------------------------------------------------------*)
(* Implementation : Michael Riedl *)
(* Licence : GNU Lesser General Public License (LGPL) *)
(*------------------------------------------------------------------------*)
(* $Id: OptimLib1.mod,v 1.3 2018/06/30 09:49:32 mriedl Exp mriedl $ *)
FROM SYSTEM IMPORT TSIZE;
FROM Storage IMPORT ALLOCATE,DEALLOCATE;
FROM Errors IMPORT Fehler,Fehlerflag,ErrOut;
IMPORT Deklera;
IMPORT Errors;
FROM LongMath IMPORT sqrt,power;
IMPORT LongMath;
FROM LMathLib IMPORT Small,CGold,MachEps,sqrt5,Pow10;
FROM MatLib IMPORT Stuerz;
FROM EigenLib1 IMPORT HhQL;
FROM RandomLib IMPORT Zufall;
(*
FROM EigenLib2 IMPORT Balance,BalBak;
*)
FROM Timer IMPORT GetUserTime;
FROM TFormAus IMPORT Write1,Write2,Write3,Write5;
FROM FileSystem IMPORT Output;
FROM FMatEA IMPORT MATaus;
IMPORT TIO;
VAR debug : BOOLEAN;
PROCEDURE FMin( A,B : LONGREAL;
Func : Funktion1;
Tol : LONGREAL;
VAR iter : CARDINAL;
VAR Min : LONGREAL);
PROCEDURE DSIGN(X,Y : LONGREAL) : LONGREAL;
BEGIN
IF (Y >= 0.0) THEN RETURN ABS(X); ELSE RETURN -ABS(X); END;
END DSIGN;
(* c is the squared inverse of the golden ratio *)
CONST C = 0.5*(3.0 - sqrt5);
MaxIter = 1000;
VAR a,b,d,e,p,q,r,u,v,w : LONGREAL;
eps,tol1,tol2 : LONGREAL;
fu,fv,fw,fx,x,xm : LONGREAL;
flag,goldsec : BOOLEAN;
BEGIN
eps:=sqrt(MachEps);
(* initialization *)
IF (Tol <= 0.0) THEN Tol:=1.0E-12; END;
a := A;
b := B;
IF (b < a) THEN
w:=a; a:=b; b:=w;
END;
v := a + C*(b - a);
w := v;
x := v;
e := 0.0;
fx := Func(x);
fv := fx;
fw := fx;
d := 0.0;
iter:=0;
LOOP
xm := 0.5*(a+b);
tol1 := eps*ABS(x) + Tol/3.0; (* Actual tolerance *)
tol2 := 2.0*tol1;
(* check stopping criterion *)
IF (ABS(x-xm) <= (tol2-0.5*(b-a))) THEN
EXIT; (* Acceptable approx. is found *)
END;
INC(iter);
IF (iter > MaxIter) THEN
Errors.Fehlerflag:="MaxIter (1000) uerberschritten (FMin)";
Errors.Fehler:=TRUE;
EXIT;
END;
(* is golden-section necessary *)
goldsec:=TRUE;
IF (ABS(e) > tol1) THEN
(* fit parabola *)
r := (x - w)*(fx - fv);
q := (x - v)*(fx - fw);
p := (x - v)*q - (x - w)*r;
q := 2.0*(q - r);
IF (q > 0.0) THEN
p := - p;
END;
q := ABS(q);
r := e;
e := d;
(* is parabola acceptable *)
IF (ABS(p) < ABS(0.5*q*r)) THEN
IF (p > q*(a-x)) THEN
IF (p < q*(b-x)) THEN
(* a parabolic interpolation step *)
d := p / q;
u := x + d;
(* f must not be evaluated too close to Xa or Xb *)
IF ((u-a) < tol2) THEN
d := DSIGN(tol1,xm-x)
END;
IF ((b-u) < tol2) THEN
d := DSIGN(tol1,xm-x)
END;
goldsec:=FALSE;
END; (* IF *)
END; (* IF *)
END; (* IF *)
END; (* IF *)
IF (goldsec) THEN (* a golden-section step *)
IF (x >= xm) THEN
e := a - x
END;
IF (x < xm) THEN
e := b - x
END;
d := C*e;
END; (* IF *)
(* F must not be evaluated too close to x *)
IF (ABS(d) >= tol1) THEN
u := x + d
ELSE
u := x + DSIGN(tol1,d);
END;
fu := Func(u); (* !!! *)
(* update a, b, v, w, and x *)
IF (fu > fx) THEN
IF (u < x) THEN
a := u
ELSE
b := u;
END;
flag:=TRUE;
IF (fu > fw) THEN
IF (w # x) THEN
IF ((fu <= fv) OR (v = x)) THEN
v := u;
fv := fu;
flag:=FALSE;
ELSIF (v # w) THEN
flag:=FALSE;
ELSE
v := u;
fv := fu;
flag:=FALSE;
END;
END;
END;
IF (flag) THEN
v := w;
fv := fw;
w := u;
fw := fu;
END;
ELSE
IF (u >= x) THEN
a := x;
END;
IF (u < x) THEN
b := x;
END;
v := w;
fv := fw;
w := x;
fw := fx;
x := u;
fx := fu;
END; (* IF *)
END; (* LOOP *)
Min := x;
END FMin;
PROCEDURE FMinBr( a,b : LONGREAL;
f : Funktion1;
tol : LONGREAL;
VAR iter : CARDINAL;
VAR Min : LONGREAL);
CONST MaxIter = 1000;
VAR x,fx,fv,fw,ft,atol : LONGREAL;
p,q,r,t,v,w : LONGREAL;
range,middlerange : LONGREAL;
newstep : LONGREAL;
BEGIN
IF (tol <= 0.0) THEN tol:=1.0E-12; END;
IF (b < a) THEN
t:=a; a:=b; b:=t;
END;
r := (3.0 - sqrt(5.0)) / 2.0; (* Gold section ratio *)
v := a + r*(b - a);
fv := f(v); (* First step - always gold section*)
x := v;
w := v;
fx:=fv;
fw:=fv;
iter:=0;
LOOP
range := b-a; (* Range over which the minimum is seeked for *)
middlerange := (a+b)/2.0;
atol := sqrt(MachEps)*ABS(x) + tol/3.0; (* Actual tolerance *)
IF (ABS(x-middlerange) + 0.5*range <= 2.0*atol) THEN
EXIT; (* Acceptable approx. is found *)
END;
INC(iter);
IF (iter > MaxIter) THEN EXIT; END;
IF (x < middlerange) THEN (* Obtain the gold section step *)
newstep := r*(b-x);
ELSE
newstep := r*(a-x);
END;
(* Decide IF the interpolation can be tried *)
IF (ABS(x-w) >= atol) THEN
(* If x and w are distinct interpolatiom may be tried *)
t := (x-w) * (fx-fv);
q := (x-v) * (fx-fw);
p := (x-v)*q - (x-w)*t;
q := 2.0*(q-t);
IF (q > 0.0) THEN (* q was calculated with the *)
p := -p; (* opposite sign; make q positive *)
ELSE (* and assign possible minus to p *)
q := -q;
END;
IF (ABS(p) < ABS(newstep*q)) AND (* If x+p/q falls in [a,b] *)
(p > q*(a-x+2.0*atol)) AND (* not too close to a and *)
(p < q*(b-x-2.0*atol)) (* b, and isn't too large *)
THEN
newstep := p/q; (* it is accepted *)
END;
(* If p/q is too large then the gold section procedure can *)
(* reduce [a,b] range to more extent *)
END;
IF (ABS(newstep) < atol) THEN (* Adjust the step to be *)
IF (newstep > 0.0) THEN (* not less than tolerance *)
newstep := atol;
ELSE
newstep := -atol;
END;
END;
(* Obtain the next approximation to min *)
(* and reduce the enveloping range *)
t := x + newstep; (* Tentative point for the min *)
ft := f(t);
IF (ft <= fx) THEN (* t is a better approximation *)
IF (t < x) THEN (* Reduce the range so that *)
b := x; (* t would fall within it *)
ELSE
a := x;
END;
v := w; w := x; x := t; (* Assign the best approx to x *)
fv:=fw; fw:=fx; fx:=ft;
ELSE (* x remains the better approx *)
IF (t < x) THEN (* Reduce the range enclosing x *)
a := t;
ELSE
b := t;
END;
IF (ft <= fw) OR (w = x) THEN
v := w; w := t;
fv:=fw; fw:=ft;
ELSIF (ft <= fv) OR (v = x) OR (v = w) THEN
v := t;
fv:=ft;
END;
END;
END; (* LOOP *)
Min := x; (* Acceptable approx. is found *)
END FMinBr;
PROCEDURE GloMin( a,b,c : LONGREAL;
m : LONGREAL;
e : LONGREAL;
t : LONGREAL;
f : Funktion1;
VAR x : LONGREAL) : LONGREAL;
VAR a0,a2,a3,d0,d1,d2 : LONGREAL;
h,m2,p,q,qs,r,s : LONGREAL;
y,y0,y1,y2,y3,yb : LONGREAL;
z0,z1,z2 : LONGREAL;
k : LONGINT;
retry,done : BOOLEAN;
BEGIN
(* Initialization *)
x :=b;
a0:=b;
a2:=a;
y0:=f(b);
yb:=y0;
y2:=f(a);
y :=y2;
IF (y0 < y) THEN y:=y0 ELSE x:=a; END;
IF (m > 0.0) AND (a < b) THEN
(* Nontrivial case (m > 0, a < b) *)
m2:=0.5*(1.0 + 16.0*MachEps)*m;
IF (c <= a) OR (c >= b) THEN c:=0.5*(a + b); END;
y1:=f(c); k:=3; d0:=a2 -c; h:= 9.0/11.0;
IF (y1 < y) THEN x:=c; y:=y1; END;
REPEAT (* Main loop *)
d1 := a2 - a0;
d2 := c - a0;
z2 := b - a2;
z0 := y2 - y1;
z1 := y2 - y0;
r := d1*d1*z0 - d0*d0*z1;
p := r;
qs := 2.0*(d0*z1 - d1*z0);
q := qs;
(* Try to find a lower value of f using quadratic interpolation *)
retry:=FALSE;
REPEAT
IF (NOT ((NOT retry) AND (((k > 1000000) AND (y < y2))))) THEN
IF (q*(r*(yb-y2)+z2*q*((y2-y)+t)) < (z2*m2*r*(z2*q-r))) THEN
a3:=a2 + r/q;
y3:=f(a3);
IF (y3 < y) THEN x:=a3; y:=y3; END;
END;
END;
(* With a probability of about 0.1 do a random search for *)
(* a lower value of f. Any reasonable random number generator *)
(* can be used in place of the one here (it need not to be *)
(* very good) *)
k:= ((1611*k) MOD 1048576);
r:=(b - a)*(VAL(LONGREAL,k) / 100000.0);
q:=1.0;
retry := (r < z2);
UNTIL (NOT retry);
(* Prepare to step as far as possible *)
r := m2*d0*d1*d2;
s := sqrt(((y2 - y) + t) / m2);
h := 0.5*(1.0 + h);
p := h*(p + 2.0*r*s);
q := r + 0.5*qs;
r := -0.5*(d0 + (z0 + 2.01*e) / (d0*m2));
IF (r >= s) AND (d0 >= 0.0) THEN r:=a2 + r ELSE r:=a2 + s; END;
(* It is safe to step to r, but we may try to step further *)
IF (p*q > 0.0) THEN a3 := a2 + (p / q) ELSE a3 := r; END;
REPEAT
IF (a3 < r) THEN a3:=r; END;
IF (a3 >= b) THEN
a3 := b;
y3 := yb;
ELSE
y3 := f(a3);
END;
IF (y3 < y) THEN x:=a3; y:=y3; END;
d0 := a3 - a2;
done:=TRUE;
IF (a3 > r) THEN
(* Inspect the parabolic lower bound on f in (a2,a3) *)
p:=2.0*(y2 - y3) / (m*d0);
IF ((ABS(p) < (1.0 + 9.0*MachEps)*d0) AND
((0.5*m2*(d0*d0 + p*p)) > ((y2 - y) + (y3 - y) + 2.0*t)))
THEN
(* Halve the step and try again *)
a3 := 0.5*(a2 + a3);
h := 0.9*h;
done:=FALSE;
END;
END;
UNTIL (done);
IF (a3 < b) THEN
(* Prapare for the next step *)
a0 := c;
c := a2;
a2 := a3;
y0 := y1;
y1 := y2;
y2 := y3;
END;
UNTIL (a3 >= b);
END;
RETURN y;
END GloMin;
PROCEDURE LocalMin( A,B : LONGREAL;
eps : LONGREAL;
T : LONGREAL;
F : Funktion1;
VAR x : LONGREAL) : LONGREAL;
CONST C = 0.5*(3.0 - sqrt5);
(* C is the square of the inverse of the golden ratio. *)
VAR d,e,u,v,w : LONGREAL;
fu,fv,fw,fx : LONGREAL;
m,p,q,r,t2,tol : LONGREAL;
BEGIN
x := A + C*(B - A);
w := x; v := w;
e := 0.0;
fx := F(x); fw := fx; fv := fw;
d := MAX(LONGREAL);
LOOP
m := 0.5*(A + B);
tol := eps*ABS (x) + T;
t2 := 2.0 * tol;
(* Check the stopping criterion. *)
IF (ABS(x - m) <= t2 - 0.5*(B - A)) THEN
EXIT;
END;
(* Fit a parabola. *)
r:=0.0; q:=0.0; p:=0.0;
IF (tol < ABS(e)) THEN
r := (x - w)*(fx - fv);
q := (x - v)*(fx - fw);
p := (x - v)*q - (x - w)*r;
q := 2.0*(q - r);
IF (0.0 < q) THEN
p := - p;
END;
q := ABS (q);
r := e; e := d; (* d nicht initialisiert *)
END;
IF (ABS(p) < ABS(0.5*q*r)) AND (q*(A-x) < p) AND (p < q*(B-x)) THEN
(* Take the parabolic interpolation step. *)
d := p / q;
u := x + d;
(* F must not be evaluated too close to A or B. *)
IF ((u-A) < t2) OR ((B-u) < t2) THEN
IF (x < m) THEN d := tol; ELSE d := - tol; END;
END;
ELSE
(* A golden-section step. *)
IF x < m THEN e := B - x; ELSE e := A - x; END;
d := C*e;
END;
(* F must not be evaluated too close to X. *)
IF (tol <= ABS(d)) THEN
u := x + d;
ELSIF (0.0 < d) THEN
u := x + tol;
ELSE
u := x - tol;
END;
fu := F(u);
(* Update A, B, V, W, and X. *)
IF (fu <= fx) THEN
IF (u < x) THEN
B := x;
ELSE
A := x;
END;
v := w; fv := fw;
w := x; fw := fx;
x := u; fx := fu;
ELSE
IF (u < x) THEN
A := u;
ELSE
B := u;
END;
IF (fu <= fw) OR (w = x) THEN
v := w; fv := fw;
w := u; fw := fu;
ELSIF (fu <= fv) OR (v = x) OR (v = w) THEN
v := u; fv := fu;
END;
END;
END; (* LOOP *)
RETURN fx;
END LocalMin;
PROCEDURE Zero(a,b : LONGREAL;
t : LONGREAL;
F : Funktion1) : LONGREAL;
(* Zero seeks the root of a function f(x) in an interval [A,B]. *)
VAR c,d,e : LONGREAL;
fa,fb,fc : LONGREAL;
m : LONGREAL;
tol : LONGREAL;
p,q,r,s : LONGREAL;
sa,sb : LONGREAL;
BEGIN
Fehler := FALSE;
(* Make local copies of A and B. *)
sa := a;
sb := b;
fa := F(sa);
fb := F(sb);
c := sa;
fc := fa;
e := sb - sa;
d := e;
LOOP
IF ABS(fc) < ABS(fb) THEN
sa := sb;
sb := c;
c := sa;
fa := fb;
fb := fc;
fc := fa;
END;
tol := 2.0*MachEps*ABS(sb) + t;
m := 0.5 * (c - sb);
IF (ABS(m) <= tol) OR (fb = 0.0) THEN
IF (ABS(fb) > 100.0*MachEps) THEN (* MRi, 01.05.18 *)
Fehler := TRUE;
Fehlerflag:="Kein Nullstelle gefunden (OptimLib1.Zero)";
END;
EXIT;
END;
IF (ABS(e) < tol) OR (ABS(fa) <= ABS(fb)) THEN
e := m;
d := e;
ELSE
s := fb / fa;
IF sa = c THEN
p := 2.0 * m * s;
q := 1.0 - s;
ELSE
q := fa / fc;
r := fb / fc;
p := s*(2.0*m*q*(q - r) - (sb - sa)*(r - 1.0));
q := (q - 1.0)*(r - 1.0)*(s - 1.0);
END;
IF 0.0 < p THEN
q := - q;
ELSE
p := - p;
END;
s := e;
e := d;
IF (2.0*p < 3.0*m*q - ABS(tol*q)) AND (p < ABS(0.5*s*q)) THEN
d := p / q;
ELSE
e := m;
d := e;
END;
END;
sa := sb;
fa := fb;
IF tol < ABS(d) THEN
sb := sb + d;
ELSIF 0.0 < m THEN
sb := sb + tol;
ELSE
sb := sb - tol;
END;
fb := F(sb);
IF (0.0 < fb) AND (0.0 < fc) OR (fb <= 0.0) AND (fc <= 0.0) THEN
c := sa;
fc := fa;
e := sb - sa;
d := e;
END;
END; (* LOOP *)
RETURN sb;
END Zero;
PROCEDURE NelMin( n : INTEGER;
VAR Start : ARRAY OF LONGREAL; (* <==> *)
VAR XMin : ARRAY OF LONGREAL; (* ==> *)
VAR YnewLo : LONGREAL;
FNProc : MinProc;
ReqMin : LONGREAL;
VAR Step : ARRAY OF LONGREAL; (* <== *)
Konvge : INTEGER;
KCount : INTEGER;
VAR ICount : INTEGER;
VAR NumRes : INTEGER;
VAR IFault : INTEGER);
PROCEDURE sqr(x : LONGREAL) : LONGREAL; BEGIN RETURN x*x; END sqr;
CONST ccoeff = 0.5; (* contraction *)
rcoeff = 1.0; (* reflection *)
ecoeff = 2.0; (* extension *)
eps = 1.0E-03;
VAR i,j,l,nm1 : INTEGER;
ihi,ilo,jcount : INTEGER;
del,dn,dnn,rq,x,z : LONGREAL;
ystar,y2star,ylo : LONGREAL;
pstar,p2star,pbar,y : POINTER TO
ARRAY [0..MAX(INTEGER)-1] OF LONGREAL;
p : POINTER TO
ARRAY [0..MaxVar] OF
ARRAY [0..MaxVar-1] OF LONGREAL;
done,check : BOOLEAN;
BEGIN
IF (ReqMin <= 0.0) THEN (* Check the input parameters *)
IFault:=1;
RETURN;
ELSIF (n < 1) THEN
IFault:=1;
RETURN;
ELSIF (n >= MaxVar) THEN
IFault:=1;
RETURN;
ELSIF (Konvge < 1) THEN
IFault:=1;
RETURN;
END;
NEW(p);
ALLOCATE(y, (n+1)*TSIZE(LONGREAL));
ALLOCATE(pstar ,n*TSIZE(LONGREAL));
ALLOCATE(p2star,n*TSIZE(LONGREAL));
ALLOCATE(pbar ,n*TSIZE(LONGREAL));
ICount:=0;
NumRes:=0;
jcount:=Konvge;
dn := VAL(LONGREAL,n);
nm1 := n - 1;
dnn := VAL(LONGREAL,n+1);
del := 1.0;
rq := ReqMin*dn;
(* Initial or restarted loop *)
LOOP
(* construction of initial simplex *)
FOR i:=0 TO nm1 DO
p^[n,i] := Start[i];
END;
FNProc(Start,n,y^[n]); INC(ICount);
FOR j:=0 TO nm1 DO
x := Start[j];
Start[j]:=Start[j] + Step[j]*del;
FOR i:=0 TO nm1 DO
p^[j,i] := Start[i];
END;
FNProc(Start,n,y^[j]); INC(ICount);
Start[j] := x;
END;
(* The simplex construction is complete. *)
(* Find highest and lowest y values YnewLo := y[ihi] *)
(* indicates the vertex of the simplex to be replaced *)
ylo := y^[0];
ilo := 0;
FOR i:=1 TO n DO
IF (y^[i] < ylo) THEN
ylo := y^[i];
ilo := i;
END;
END;
REPEAT
done:=TRUE;
YnewLo := y^[0];
ihi := 0;
FOR i:=1 TO n DO
IF (YnewLo < y^[i]) THEN
YnewLo := y^[i];
ihi := i;
END;
END;
(* Calculate pbar, the centroid of the simplex vertices *)
(* excepting the vertex with Y value YnewLo. *)
FOR i:=0 TO nm1 DO
z:=0.0;
FOR j:=0 TO n DO (* Hier wird auf p^[n,.] zugegriffen *)
z:=z + p^[j,i];
END;
z:=z - p^[ihi,i];
pbar^[i] := z / dn;
END;
FOR i:=0 TO nm1 DO (* Reflection through the centroid *)
pstar^[i] := pbar^[i] + rcoeff*(pbar^[i] - p^[ihi,i]);
END;
FNProc(pstar^,n,ystar); INC(ICount);
check:=TRUE;
IF (ystar < ylo) THEN (* Successful reflection, so extension *)
FOR i:=0 TO nm1 DO
p2star^[i] := pbar^[i] + ecoeff*(pstar^[i] - pbar^[i]);
END;
FNProc(p2star^,n,y2star); INC(ICount);
IF (ystar < y2star) THEN (* Check extension *)
FOR i:=0 TO nm1 DO
p^[ihi,i] := pstar^[i];
END;
y^[ihi] := ystar;
ELSE (* Retain extension or contraction *)
FOR i:=0 TO nm1 DO
p^[ihi,i] := p2star^[i];
END;
y^[ihi] := y2star;
END; (* IF *)
ELSE (* No extension *)
l := 0;
FOR i:=0 TO n DO
IF (ystar < y^[i]) THEN
INC(l);
END;
END;
IF (1 < l) THEN
FOR i:=0 TO nm1 DO
p^[ihi,i] := pstar^[i];
END;
y^[ihi] := ystar;
ELSIF (l = 0) THEN
(* Contraction on the Y(IHI) side of the centroid *)
FOR i:=0 TO nm1 DO
p2star^[i] := pbar^[i] + ccoeff*(p^[ihi,i] - pbar^[i]);
END;
FNProc(p2star^,n,y2star); INC(ICount);
IF (y^[ihi] < y2star) THEN
(* Contract the whole simplex *)
FOR j:=0 TO n DO
FOR i:=0 TO nm1 DO
p^[j,i] := (p^[j,i] + p^[ilo,i])*ccoeff;
XMin[i] := p^[j,i];
END;
FNProc(XMin,n,y^[j]); INC(ICount);
END;
IF (KCount < ICount) THEN
(* Max. Anzahl Funktionsauswertungen ueberschritten *)
done:=TRUE;
ELSE
ylo := y^[0];
ilo := 0;
FOR i:=1 TO n DO
IF (y^[i] < ylo) THEN
ylo := y^[i];
ilo := i;
END;
END;
done :=FALSE;
END;
check:=FALSE;
ELSE (* Retain contraction *)
FOR i:=0 TO nm1 DO
p^[ihi,i] := p2star^[i];
END;
y^[ihi] := y2star;
END;
ELSIF (l = 1) THEN
(* Contraction on the reflection side of the centroid *)
FOR i:=0 TO nm1 DO
p2star^[i] := pbar^[i] + ccoeff*(pstar^[i] - pbar^[i]);
END;
FNProc(p2star^,n,y2star); INC(ICount);
IF (y2star <= ystar) THEN (* Retain reflection ? *)
FOR i:=0 TO nm1 DO
p^[ihi,i] := p2star^[i];
END;
y^[ihi] := y2star;
ELSE
FOR i:=0 TO nm1 DO
p^[ihi,i] := pstar^[i];
END;
y^[ihi] := ystar;
END; (* IF *)
END; (* IF *)
END; (* IF *)
IF (check) THEN (* Check if ylo improved *)
IF (y^[ihi] < ylo) THEN
ylo := y^[ihi];
ilo := ihi;
END;
DEC(jcount);
IF (jcount # 0) THEN
done:=FALSE;
ELSE
(* Check to see if minimum reached *)
IF (ICount <= KCount) THEN
jcount := Konvge;
z:=0.0;
FOR i:=0 TO n DO
z:=z + y^[i];
END;
x := z / dnn;
z:=0.0;
FOR i:=0 TO n DO
z:=z + sqr(y^[i] - x); (* Hier geht es schief !!! *)
END;
IF (rq < z) THEN
done:=FALSE;
END;
END; (* IF *)
END; (* IF *)
END; (* IF check *)
UNTIL (done);
(* Factorial tests to check that YnewLo is a local minimum *)
FOR i:=0 TO nm1 DO
XMin[i] := p^[ilo,i];
END;
YnewLo := y^[ilo];
IF (KCount < ICount) THEN
IFault:=2;
EXIT;
END;
IFault:=0;
i:=0;
LOOP (* FOR i:=0 TO nm1 DO *)
IF (i >= n) THEN EXIT; END;
del := Step[i]*eps;
XMin[i]:=XMin[i] + del;
FNProc(XMin,n,z); INC(ICount);
IF (z < YnewLo) THEN
IFault:=2;
EXIT; (* local loop *)
END;
XMin[i]:=XMin[i] - del - del;
FNProc(XMin,n,z); INC(ICount);
IF (z < YnewLo) THEN
IFault:=2;
EXIT; (* local loop *)
END;
XMin[i]:=XMin[i] + del;
INC(i);
END; (* local LOOP *)
IF (IFault = 0) THEN
EXIT; (* outer loop *)
END;
(* Restart the procedure *)
FOR i:=0 TO nm1 DO
Start[i] := XMin[i];
END;
del := eps;
INC(NumRes);
END; (* LOOP *)
DISPOSE(p);
DEALLOCATE(y, (n+1)*TSIZE(LONGREAL));
DEALLOCATE(pstar ,n*TSIZE(LONGREAL));
DEALLOCATE(p2star,n*TSIZE(LONGREAL));
DEALLOCATE(pbar ,n*TSIZE(LONGREAL));
END NelMin;
(*=================== SECTION Praxis =====================================*)
CONST Epsilon = 1.0E+03*MachEps;
TYPE Vector = ARRAY [1..MaxVar] OF LONGREAL;
PROCEDURE sqr(x : LONGREAL) : LONGREAL; BEGIN RETURN x*x; END sqr;
PROCEDURE MinFit( N : INTEGER;
eps : LONGREAL;
Tol : LONGREAL;
VAR ab : SuchMatrix;
VAR q : Vector);
VAR i,j,k,l,kt : INTEGER;
c,f,g,h,s,x,y,z : LONGREAL;
e : Vector;
TestFconvergence : BOOLEAN;
Cancellation : BOOLEAN;
BEGIN
(* Householders reduction to bidiagonal form *)
l := 0;
x := 0.0;
g := 0.0;
FOR i:=1 TO N DO
e[i] := g;
s := 0.0;
l := i+1;
FOR j:=i TO N DO
s:=s + sqr(ab[j,i]);
END;
IF (s < Tol) THEN
g := 0.0;
ELSE
f := ab[i,i];
IF (f < 0.0) THEN
g := sqrt(s);
ELSE
g := - sqrt(s);
END;
h := f*g-s;
ab[i,i] := f - g;
FOR j:=l TO N DO
f := 0.0;
FOR k:= i TO N DO
f:=f + ab[k,i]*ab[k,j];
END;
f := f / h;
FOR k:=i TO N DO
ab[k,j] := ab[k,j] + f*ab[k,i];
END;
END; (* j *)
END; (* IF *)
q[i] := g;
s := 0.0;
FOR j:=l TO N DO
s:=s + sqr(ab[i,j]);
END;
IF (s < Tol) THEN
g := 0.0;
ELSE
f := ab[i,i+1];
IF (f < 0.0) THEN
g := sqrt(s);
ELSE
g := - sqrt(s);
END;
h := f*g-s;
ab[i,i+1] := f-g;
FOR j:=l TO N DO
e[j] := ab[i,j] / h;
END;
FOR j:=l TO N DO
s := 0.0;
FOR k:=l TO N DO
s:=s + ab[j,k]*ab[i,k];
END;
FOR k:=l TO N DO
ab[j,k] := ab[j,k] + s*e[k];
END;
END;
END; (* IF *)
y := ABS(q[i])+ABS(e[i]);
IF (y > x) THEN x := y; END;
END; (* i *)
(* accumulation of right hand transformations *)
FOR i:=N TO 1 BY -1 DO
IF (g <> 0.0) THEN
h := ab[i,i+1]*g;
FOR j:=l TO N DO
ab[j,i] := ab[i,j] / h;
END;
FOR j:=l TO N DO
s := 0.0;
FOR k:=l TO N DO
s := s + ab[i,k]*ab[k,j];
END;
FOR k:=l TO N DO
ab[k,j] := ab[k,j] + s*ab[k,i];
END;
END; (* j *)
END; (* IF *)
FOR j:=l TO N DO
ab[j,i] := 0.0;
ab[i,j] := 0.0;
END;
ab[i,i] := 1.0;
g := e[i];
l := i;
END; (* i *)
(* diagonalization to bidiagonal form *)
eps := eps*x;
FOR k:=N TO 1 BY -1 DO
kt := 0;
LOOP (* Test for splitting *)
kt := kt + 1;
IF (kt > 30) THEN
e[k] := 0.0;
TIO.WrStr('+++ QR failed'); TIO.WrLn;
END;
TestFconvergence := FALSE;
Cancellation := FALSE;
l := k+1;
REPEAT
l:=l - 1;
IF (ABS(e[l]) <= eps) THEN
TestFconvergence := TRUE;
ELSIF (ABS(q[l-1]) <= eps) THEN
Cancellation := TRUE;
END;
UNTIL (l = 1) OR TestFconvergence OR Cancellation;
IF Cancellation THEN
c := 0.0;
s := 1.0;
i:=l;
REPEAT
f := s*e[i];
e[i] := c*e[i];
IF (ABS(f) <= eps) THEN
TestFconvergence:=TRUE;
ELSE
g := q[i];
IF (ABS(f) < ABS(g)) THEN
h := ABS(g)*sqrt(1.0 + sqr(f / g));
ELSIF (f <> 0.0) THEN
h := ABS(f)*sqrt(1.0 + sqr(g / f));
ELSE
h := 0.0;
END;
q[i] := h;
IF (h = 0.0) THEN
h := 1.0;
g := 1.0;
END;
c := g / h;
s := -f / h;
END;
i:=i + 1;
UNTIL (i > k) OR TestFconvergence;
END;
(* TestFconvergence *)
z := q[k];
IF (l = k) THEN EXIT; END; (* convergence *)
(* shift from bottom 2*2 minor *)
x := q[l];
y := q[k-1];
g := e[k-1];
h := e[k];
f := ((y-z)*(y+z) + (g-h)*(g+h)) / (2.0*h*y);
g := sqrt(sqr(f) + 1.0);
IF (f <= 0.0) THEN
f := ((x-z)*(x+z) + h*(y / (f-g)-h)) / x;
ELSE
f := ((x-z)*(x+z) + h*(y / (f+g)-h)) / x;
END;
(* next QR transformation *)
s := 1.0;
c := 1.0;
FOR i:= l+1 TO k DO
g := e[i];
y := q[i];
h := s*g;
g := g*c;
IF (ABS(f)<ABS(h)) THEN
z := ABS(h)*sqrt(1.0 + sqr(f / h))
ELSIF (f <> 0.0) THEN
z := ABS(f)*sqrt(1.0 + sqr(h / f))
ELSE
z := 0.0;
END;
e[i-1] := z;
IF (z = 0.0) THEN
f := 1.0;
z := 1.0;
END;
c := f / z;
s := h / z;
f := x*c + g*s;
g := - x*s + g*c;
h := y*s;
y := y*c;
FOR j:=1 TO N DO
x := ab[j,i-1];
z := ab[j,i];
ab[j,i-1] := x*c + z*s;
ab[j,i] := - x*s + z*c;
END;
IF (ABS(f)<ABS(h)) THEN
z := ABS(h)*sqrt(1.0 + sqr(f / h));
ELSIF (f # 0.0) THEN
z := ABS(f)*sqrt(1.0 + sqr(h / f));
ELSE
z := 0.0;
END;
q[i-1] := z;
IF (z = 0.0) THEN
z := 1.0;
f := 1.0;
END;
c := f / z;
s := h / z;
f := c*g + s*y;
x := -s*g + c*y;
END; (* i *)
e[l] := 0.0;
e[k] := f;
q[k] := x;
END; (* LOOP *)
(* Convergence *)
IF (z < 0.0) THEN
q[k] := -z;
FOR j:=1 TO N DO
ab[j,k] := - ab[j,k];
END;
END;
END; (* k *)
END MinFit;
PROCEDURE Praxis(VAR X : ARRAY OF LONGREAL;
N : INTEGER;
F : MinProc;
VAR fx : LONGREAL;
t : LONGREAL;
h : LONGREAL;
Ktm : INTEGER;
Scbd : LONGREAL;
IllC : BOOLEAN;
Prin : CARDINAL;
MaxF : CARDINAL;
VAR nf : CARDINAL;
VAR iFehl : INTEGER);
VAR i,j,k,k2,nl,kl,kt : INTEGER;
s,sl,dn,dmin,f1 : LONGREAL;
lds,ldt,sf,df : LONGREAL;
qf1,qd0,qd1,qa,qb,qc : LONGREAL;
m2,m4,ldfac,t2 : LONGREAL;
small,vsmall : LONGREAL;
large,vlarge : LONGREAL;
d,y,z,q0,q1 : Vector;
v : SuchMatrix;
PROCEDURE Sort( N : INTEGER;
VAR d : Vector;
VAR v : SuchMatrix);
(* sort d and v in descending order *)
VAR k,i,j : INTEGER;
s : LONGREAL;
BEGIN
FOR i:=1 TO N-1 DO
k := i;
s := d[i];
FOR j:=i+1 TO N DO
IF (d[j] > s) THEN
k := j;
s := d[j];
END;
END;
IF (k > i) THEN
d[k] := d[i];
d[i] := s;
FOR j:=1 TO N DO
s := v[j,i];
v[j,i] := v[j,k];
v[j,k] := s;
END;
END; (* IF *)
END; (* FOR *)
END Sort;
PROCEDURE Print;
(* print a line of traces *)
VAR i : INTEGER;
BEGIN
TIO.WrStr('------------------------------------------------------');
TIO.WrLn;
TIO.WrStr('Chi Square reduced to ... '); TIO.WrLngReal(fx,20,-10);
TIO.WrLn;
TIO.WrStr(' ... after '); TIO.WrInt(nf,1);
TIO.WrStr(' function calls ...'); TIO.WrLn;
TIO.WrStr(' ... including '); TIO.WrInt(nl,1);
TIO.WrStr(' linear searches.'); TIO.WrLn;
TIO.WrStr('Current Values of X ...'); TIO.WrLn;
FOR i:=1 TO N DO
TIO.WrLngReal(X[i-1],20,-10);
IF ((i MOD 3) = 0) THEN TIO.WrLn; END;
END;
TIO.WrLn;
END Print;
PROCEDURE MatPrint(s : ARRAY OF CHAR;
v : SuchMatrix;
n,m : INTEGER);
VAR k,i : INTEGER;
BEGIN
TIO.WrLn; TIO.WrStr(s); TIO.WrLn; TIO.WrLn;
FOR k:=1 TO n DO
FOR i:=1 TO m DO
TIO.WrLngReal(v[k,i],20,-10);
END;
TIO.WrLn;
END;
TIO.WrLn;
END MatPrint;
PROCEDURE VecPrint( s : ARRAY OF CHAR;
VAR v : ARRAY OF LONGREAL;
n : INTEGER);
VAR i : INTEGER;
BEGIN
TIO.WrLn; TIO.WrStr(s); TIO.WrLn; TIO.WrLn;
FOR i:=0 TO n-1 DO
TIO.WrLngReal(v[i],20,-10);
IF (((i+1) MOD 3) = 0) THEN TIO.WrLn; END;
END;
TIO.WrLn;
END VecPrint;
PROCEDURE Min( j,Nits : INTEGER;
VAR d2,x1 : LONGREAL;
F : MinProc;
f1 : LONGREAL;
fk : BOOLEAN);
VAR k,i : INTEGER;
dz : BOOLEAN;
x2,xm,f0,f2,fm : LONGREAL;
d1,t2,s,sf1,sx1 : LONGREAL;
PROCEDURE Flin(l : LONGREAL) : LONGREAL;
VAR i : INTEGER;
t : Vector;
fx : LONGREAL;
BEGIN
IF (j > 0) THEN (* linear search *)
FOR i:=1 TO N DO
t[i] := X[i-1] + l*v[i,j];
END;
ELSE
(* search along parabolic space curve *)
qa := l*(l-qd1) / (qd0*(qd0+qd1));
qb := (l+qd0)*(qd1-l) / (qd0*qd1);
qc := l*(l+qd0) / (qd1*(qd0+qd1));
FOR i:=1 TO N DO
t[i] := qa*q0[i] + qb*X[i-1] + qc*q1[i];
END;
END;
F(t,N,fx); INC(nf);
RETURN fx;
END Flin;
VAR l0,l1 : BOOLEAN;
BEGIN (* Min *)
sf1 := f1;
sx1 := x1;
k := 0;
xm := 0.0;
fm := fx;
f0 := fx;
dz := (d2 < Epsilon);
(* Find step size *)
s:=0.0;
FOR i:=1 TO N DO
s:=s + sqr(X[i-1]);
END;
s := sqrt(s);
IF dz THEN
t2 := m4*sqrt(ABS(fx) / dmin + s*ldt) + m2*ldt
ELSE
t2 := m4*sqrt(ABS(fx) / d2 + s*ldt) + m2*ldt;
END;
s := m4*s + t;
IF (dz AND (t2 > s)) THEN t2 := s; END;
IF (t2 < small) THEN t2 := small; END;
IF (t2 > (0.01*h)) THEN t2 := 0.01*h; END;
IF (fk AND (f1 <= fm)) THEN
xm := x1;
fm := f1;
END;
IF (NOT fk OR (ABS(x1)<t2)) THEN
IF (x1 > 0.0) THEN
x1 := t2;
ELSE
x1 := - t2;
END;
f1 := Flin(x1);
END;
IF (f1 <= fm) THEN
xm := x1;
fm := f1;
END;
REPEAT
IF dz THEN
IF (f0 < f1) THEN
x2 := - x1;
ELSE
x2 := 2.0*x1;
END;
f2 := Flin(x2);
IF (f2 <= fm) THEN
xm := x2;
fm := f2;
END;
d2 := (x2*(f1-f0) - x1*(f2-f0)) / (x1*x2*(x1-x2));
END;
d1 := (f1-f0) / x1 - x1*d2;
dz := TRUE;
IF (d2 <= small) THEN
IF (d1 < 0.0) THEN
x2 := h;
ELSE
x2 := -h;
END;
ELSE
x2 := -0.5*d1 / d2;
END;
IF (ABS(x2)>h) THEN
IF (x2 > 0.0) THEN
x2 := h;
ELSE
x2 := -h;
END;
END;
l0:=FALSE;
REPEAT
f2 := Flin(x2);
l1:=TRUE;
IF (k < Nits) AND (f2 > f0) THEN
k:=k + 1;
IF (f0 < f1) AND ((x1*x2) > 0.0) THEN
l0:=TRUE;
ELSE
x2 := 0.5*x2;
END;
l1:=FALSE;
END;
UNTIL (l0 OR l1);
UNTIL (NOT l0);
nl:=nl + 1;
IF (f2 > fm) THEN
x2 := xm;
ELSE
fm := f2;
END;
IF (ABS(x2*(x2 - x1)) > small) THEN
d2 := (x2*(f1 - f0) - x1*(fm - f0)) / (x1*x2*(x1 - x2));
ELSIF (k > 0) THEN
d2 := 0.0;
END;
IF (d2 <= small) THEN d2 := small; END;
x1 := x2;
fx := fm;
IF (sf1 < fx) THEN
fx := sf1;
x1 := sx1;
END;
IF (j > 0) THEN
FOR i:=1 TO N DO
X[i-1]:=X[i-1] + x1*v[i,j];
END;
END;
END Min;
PROCEDURE Quad( N : INTEGER;
VAR X : ARRAY OF LONGREAL);
(* look for a minimum along the curve q0, q1, q2 *)
VAR i : INTEGER;
l,s : LONGREAL;
BEGIN
s := fx;
fx := qf1;
qf1 := s;
qd1 := 0.0;
FOR i:=1 TO N DO
s := X[i-1];
l := q1[i];
X[i-1] := l;
q1[i] := s;
qd1 := qd1 + sqr(s - l);
END;
s := 0.0;
qd1 := sqrt(qd1);
l := qd1;
IF (qd0 > 0.0) AND (qd1 > 0.0) AND (nl >= (3*N*N)) THEN
Min(0,2,s,l,F,qf1,TRUE);
qa := l*(l - qd1) / (qd0*(qd0 + qd1));
qb := (l + qd0) * (qd1 - l) / (qd0*qd1);
qc := l*(l + qd0) / (qd1*(qd0 + qd1));
ELSE
fx := qf1;
qa := 0.0;
qb := 0.0;
qc := 1.0;
END;
qd0 := qd1;
FOR i:=1 TO N DO
s := q0[i];
q0[i] := X[i-1];
X[i-1] := qa*s + qb*X[i-1] + qc*q1[i];
END;
END Quad;
VAR done : BOOLEAN;
BEGIN (* Praxis *)
(* initialization *)
iFehl:=0;
qd1:=0.0;
small := sqr(Epsilon);
vsmall := sqr(small);
large := 1.0 / small;
vlarge := 1.0 / vsmall;
m2 := sqrt(Epsilon);
m4 := sqrt(m2);
IF IllC THEN
ldfac := 0.1;
ELSE
ldfac := 0.01;
END;
nl := 0;
kt := 0;
nf := 0;
(* F(X,N,fx); INC(nf); wird uebergeben *)
qf1 := fx;
t2 := small + ABS(t);
t := t2;
dmin := small;
IF (h < (100.0*t)) THEN h := 100.0*t; END;
ldt := h;
FOR i:=1 TO N DO (* initial search directions is the identity matrix *)
FOR j:=1 TO N DO
v[i,j] := 0.0;
END;
v[i,i] := 1.0;
END;
d[1] := 0.0;
qd0 := 0.0;
FOR i:=1 TO N DO
q1[i] := X[i-1];
END;
IF (Prin > 1) THEN
TIO.WrLn;
TIO.WrStr('---------- enter function praxis ------------'); TIO.WrLn();
TIO.WrStr('Current paramter settings are:'); TIO.WrLn();
TIO.WrStr('... scaling ... '); TIO.WrLngReal(Scbd,20,-10); TIO.WrLn();
TIO.WrStr('... Epsilon ... '); TIO.WrLngReal(Epsilon,20,-10); TIO.WrLn();
TIO.WrStr('... Tol ... '); TIO.WrLngReal(t,20,-10); TIO.WrLn();
TIO.WrStr('... MaxStep ... '); TIO.WrLngReal(h,20,-10); TIO.WrLn();
TIO.WrStr('... IllC ... '); TIO.WrBool (IllC); TIO.WrLn();
TIO.WrStr('... Ktm ... '); TIO.WrInt (Ktm,1); TIO.WrLn();
END;
IF (Prin > 0) THEN Print; END;
LOOP
sf := d[1];
s := 0.0;
d[1] := 0.0;
(* minimize along first direction *)
Min(1,2,d[1],s,F,fx,FALSE);
IF (s <= 0.0) THEN
FOR i:=1 TO N DO
v[i,1] := - v[i,1];
END;
END;
IF (sf <= (0.9*d[1])) OR ((0.9*sf) >= d[1]) THEN
FOR i:=2 TO N DO
d[i] := 0.0;
END;
END;
FOR k:=2 TO N DO
FOR i:=1 TO N DO
y[i] := X[i-1];
END;
sf := fx;
IllC := IllC OR (kt>0);
done:=TRUE;
REPEAT
kl := k;
df := 0.0;
IF IllC THEN
(* random step to get off resolution valley *)
FOR i:=1 TO N DO
z[i] := (0.1*ldt + t2*Pow10(kt))*(Zufall() - 0.5);
s := z[i];
FOR j:=1 TO N DO
X[j-1]:=X[j-1] + s*v[j,i];
END;
END; (* I *)
F(X,N,fx); INC(nf);
IF (nf > MaxF) THEN
Errors.Fehlerflag:=
"Maximalzahl der Funktionsauswertungen ueberschritten (Praxis)";
Errors.ErrOut(Errors.Fehlerflag);
iFehl:=1;
EXIT;
END;
END; (* IF *)
FOR k2:=k TO N DO
(* minimize along non-conjugate directions *)
sl := fx;
s := 0.0;
Min(k2,2,d[k2],s,F,fx,FALSE);
IF IllC THEN
s := d[k2]*sqr(s + z[k2])
ELSE
s := sl - fx;
END;
IF (df < s) THEN
df := s;
kl := k2;
END;
END; (* K2 *)
IF NOT IllC AND (df < ABS(100.0*Epsilon*fx)) THEN
(* If there was not much improvement on the first try, *)
(* set IllC = true and start the inner loop again. *)
IllC := TRUE;
done := FALSE;
END;
UNTIL done;
IF (k = 2) AND (Prin > 1) THEN
VecPrint('New Direction ...', d, N);
END;
FOR k2:=1 TO k-1 DO
(* minimize along conjugate directions *)
s := 0.0;
Min(k2,2,d[k2],s,F,fx,FALSE);
END; (* K2 *)
f1 := fx;
fx := sf;
lds := 0.0;
FOR i:=1 TO N DO
sl := X[i-1];
X[i-1] := y[i];
y[i] := sl - y[i];
sl := y[i];
lds := lds + sqr(sl);
END;
lds := sqrt(lds);
IF (lds > small) THEN
(* Discard direction v[*,kl]. *)
(* If no random step was taken, v[*,kl] is the "non-conjugate" *)
(* direction along which the greatest improvement was made. *)
FOR i:=kl-1 TO k BY -1 DO
FOR j:=1 TO N DO
v[j,i+1] := v[j,i];
END;
d[i+1] := d[i];
END;
d[k] := 0.0;
FOR i:=1 TO N DO
v[i,k] := y[i] / lds;
END;
(* Minimize along the new "conjugate" direction v[*,k], *)
(* which is the normalized vector: (new x) - (old x). *)
Min(k,4,d[k],lds,F,f1,TRUE);
IF (lds <= 0.0) THEN
lds := -lds;
FOR i:=1 TO N DO
v[i,k] := -v[i,k];
END;
END;
END; (* IF *)
ldt := ldfac*ldt;
IF (ldt < lds) THEN ldt := lds; END;
IF (Prin > 1) THEN Print; END;
t2 := 0.0;
FOR i:=1 TO N DO
t2:=t2 + sqr(X[i-1]);
END;
t2 := m2*sqrt(t2) + t;
(* See whether the length of the step taken since *)
(* starting the inner loop exceeds half the tolerance. *)
IF (ldt > (0.5*t2)) THEN
kt := 0;
ELSE
kt := kt+1;
END;
IF (kt > Ktm) THEN EXIT; END;
END; (* K *)
(* try quadratic extrapolation in case *)
(* we are stuck in a curved valley *)
Quad(N,X);
dn := 0.0;
FOR i:=1 TO N DO
d[i] := 1.0 / sqrt(d[i]);
IF (dn < d[i]) THEN dn := d[i]; END;
END;
IF (Prin > 2) THEN
MatPrint('New Matrix of Directions ...', v,N,N);
END;
FOR j:=1 TO N DO
s := d[j] / dn;
FOR i:=1 TO N DO
v[i,j] := s*v[i,j];
END;
END;
IF (Scbd > 1.0) THEN
(* scale axis to reduce condition number *)
s := vlarge;
FOR i:=1 TO N DO
sl := 0.0;
FOR j:=1 TO N DO
sl := sl + sqr(v[i,j]);
END;
z[i] := sqrt(sl);
IF (z[i] < m4) THEN z[i] := m4; END;
IF (s > z[i]) THEN s := z[i]; END;
END; (* i *)
FOR i:=1 TO N DO
sl := s / z[i];
z[i] := 1.0 / sl;
IF (z[i] > Scbd) THEN
sl := 1.0 / Scbd; (* Wofuer *)
z[i] := Scbd;
END;
(*
FOR j:=1 TO N DO v[i,j]:=v[i,j]*sl; END; (* MRi, 01.05.18 *)
*)
END;
END; (* IF *)
(* Calculate a new set of orthogonal directions *)
(* before repeating the main loop. *)
FOR i:=2 TO N DO (* Transpose v for MinFit *)
FOR j:=1 TO i-1 DO
s := v[i,j];
v[i,j] := v[j,i];
v[j,i] := s;
END;
END;
(* Call MInFit to find the singular value decomposition of V. *)
(* This gives the principal values and principal directions of *)
(* the approximating quadratic form without squaring the *)
(* condition number. *)
MinFit(N,Epsilon,vsmall,v,d);
IF (Scbd > 1.0) THEN (* unscale the axes *)
FOR i:=1 TO N DO
s := z[i];
FOR j:=1 TO N DO
v[i,j] := v[i,j]*s;
END;
END;
FOR i:=1 TO N DO
s:=0.0;
FOR j:=1 TO N DO
s:=s + sqr(v[j,i]);
END;
s := sqrt(s);
d[i] := s*d[i];
s := 1.0 / s;
FOR j:=1 TO N DO
v[j,i] := s*v[j,i];
END;
END;
END;
FOR i:=1 TO N DO
IF (dn*d[i] > large) THEN
d[i] := vsmall;
ELSIF ((dn*d[i]) < small) THEN
d[i] := vlarge;
ELSE
d[i] := power(dn*d[i],-2.0);
END;
END;
Sort(N,d,v);
(* The new eigenvalues and eigenvectors *)
dmin := d[N];
IF (dmin < small) THEN dmin := small; END;
(* The ratio of the smallest to largest eigenvalue *)
(* determines whether the system is ill conditioned. *)
IllC := (m2*d[1]) > dmin;
IF (Prin > 2) AND (Scbd > 1.0) THEN
VecPrint('Scale Factors ...',z,N);
END;
IF (Prin > 2) THEN
VecPrint('Eigenvalues of A ...',d,N);
END;
IF (Prin > 2) THEN
MatPrint('Eigenvectors of A ...',v,N,N);
END;
END; (* LOOP *)
IF (Prin > 0) THEN
VecPrint('Final solution is ...', X,N);
TIO.WrLn;
TIO.WrStr('ChiSq reduced to '); TIO.WrLngReal(fx,20,-10);
TIO.WrStr(' after '); TIO.WrCard(nf,1); TIO.WrStr(' function calls.');
TIO.WrLn;
END;
END Praxis;
(*=================== SECTION Trudge =====================================*)
PROCEDURE Search(VAR X,dX : ARRAY OF LONGREAL;
NVar : CARDINAL;
VAR Alpha : LONGREAL;
Funkt : MinProc; (* Zu minimierende Funktion. *)
VAR Func : LONGREAL; (* Funktionswert. *)
Noise : LONGREAL;
exact : BOOLEAN;
VAR Curve : LONGREAL;
VAR IFehl : CARDINAL);
(*----------------------------------------------------------------*)
(* Bestimmt ein lokales Minimum der N-dimensionalen Funktion *)
(* Funkt entlang einer vorgegebenen eindimensionalen Suchrichtung *)
(* *)
(* F(LowF) : Func(X+LowF*DX) Suche den optimalen LowF-Wert. *)
(* LookMx : Maximalzahl der erlaubten Funktionsauswertung. *)
(* *)
(* Bei Eintritt in die Routine m"ussen die folgenden Variablen *)
(* gesetzt sein : *)
(* *)
(* X : Anfangsposition im N-dimensionalen Parameterraum. *)
(* dX : Suchrichtung im N-dimensionalen Parameterraum. *)
(* NVar : Anzahl der zu varierenden Parameter (N). *)
(* Alpha : Abgesch"atzte initiale Schrittweite. *)
(* Func : F(0), Funktionswert mit dem initialen X. *)
(* Noise : Genauigkeit der Funktionsauswertung. Variationen, *)
(* die kleiner als Noise sind, werden als nicht *)
(* signifikant gewertet. *)
(* exact : Wenn FALSE, wird der Test des quadratischen Fits *)
(* unterdr"uckt und angenommen, es sei das Minimum *)
(* erreicht. *)
(* *)
(* Bei verlassen der Routine sind die folgenden Werte gesetzt : *)
(* *)
(* X : Berechnete Minimum (X=X(0) + LowF*dX). *)
(* Alpha : Berechne Schrittweite am aufgefundenen Minumum. *)
(* Func : F(Alpha), Funktionswert am Ort des berechneten *)
(* Minimums. *)
(* Curve : Berechnete Kuvatur entlang der Suchrichtung. *)
(* IFehl : 0 bei normalem verlassen der Routine. *)
(* 1 bei "Uberschreiten der Maximalzahl der erlaubten *)
(* Funktionsauswertungen. *)
(* *)
(* F[1], F[3], F[5], F[7] und F[9] sind Funktionswerte an gleich- *)
(* verteilten Punkten mit Abstand Delta. *)
(* Delta kann um Faktoren von 2 vergr"o3ert oder verkleinert *)
(* werden, die F-Werte werden jeweils nach links oder rechts ver- *)
(* schoben, damit F[5] immer der tiefste Wert ist. *)
(* Jedem F_i ist ein Set_i zugeordnet, das als boolscher *)
(* Indikator anzeigt, ob F_i berechnet wurde. *)
(* *)
(* Autor der Routinen : ??? , Version Feb. 1980. *)
(* Angepasst an Modula-2 : M.Riedl, im Juni 1994. *)
(* *)
(* Die Variabelennamen sind zum gro3en Teil unver"andert, aber *)
(* die Kontrollstruktur der Routine ist deutlich gegen"uber dem *)
(* Original ver"andert, um die Routine in M2 implementieren zu *)
(* k"onnen. *)
(*----------------------------------------------------------------*)
CONST LookMx = 50;
VAR i,k,iret,Looks : CARDINAL;
F : ARRAY [1..9] OF LONGREAL;
Set : ARRAY [1..9] OF BOOLEAN;
back,hold,flag : BOOLEAN;
Delta,tol,Dum : LONGREAL;
Falt,FMin,Fquad,FNext : LONGREAL;
LowF,ANext,Spread : LONGREAL;
Xalt : ARRAY [0..MaxVar-1] OF LONGREAL;
PROCEDURE Half();
BEGIN
Delta:=0.5*Delta; hold:=TRUE;
F[1]:=F[3]; Set[1]:=Set[3];
F[3]:=F[4]; Set[3]:=Set[4]; Set[4]:=FALSE;
F[9]:=F[7]; Set[9]:=Set[7];
F[7]:=F[6]; Set[7]:=Set[6]; Set[6]:=FALSE;
Set[2]:=FALSE; Set[8]:=FALSE;
END Half;
PROCEDURE Double();
BEGIN
Delta:=2.0*Delta;
F[4]:=F[3]; Set[4]:=Set[3];
F[3]:=F[1]; Set[3]:=Set[1]; Set[1]:=FALSE;
F[6]:=F[7]; Set[6]:=Set[7];
F[7]:=F[9]; Set[7]:=Set[9]; Set[9]:=FALSE;
Set[2]:=FALSE; Set[8]:=FALSE;
END Double;
BEGIN
Fquad := 0.0; (* 01.02.2015 *)
Delta := 1.0;
hold := FALSE;
back := FALSE;
FOR i:=0 TO NVar-1 DO Xalt[i]:=X[i]; END;
LowF:=Alpha; Alpha:=0.0;
Falt:=Func; FMin:=Func;
(* FNext:=VAL(LONGREAL,MAX(REAL)); Retundant ? *)
iret := 1; IFehl:=0;
Looks:=0;
LOOP (* Funktionsauswertungen *)
IF debug THEN
Write2("/,S,E20.12,/",'LowF : ',LowF);
END;
FOR i:=0 TO NVar-1 DO X[i]:=Xalt[i] + LowF*dX[i]; END;
IF (LowF = 0.0) THEN
Func:=Falt;
ELSE
Funkt(X,NVar,Func);
END;
INC(Looks);
IF (Looks > LookMx) THEN
Fehlerflag:="Zuviele Funktionsauswertungen (Trudge.Search) !";
Fehler:=TRUE; IFehl:=1;
EXIT; (* !!! *)
END;
IF Fehler THEN (* FUSCH ??? - Wie soll das erreicht werden *)
Half()
ELSE
hold:=FALSE;
IF (Func < FMin) THEN
FNext:=FMin; FMin:=Func;
ANext:=Alpha; (* Arithmetische Ausnahme f��r Helix Funktion *)
Alpha:=LowF;
ELSE
FNext:=Func;
ANext:=LowF;
END;
flag:=FALSE;
IF debug THEN
Write2("/,S,C1,/",' iret : ',iret);
END;
CASE iret OF
|2: F[7]:=Func; Set[7]:=TRUE;
IF (LowF = Alpha) THEN (* Verschiebe nach links *)
FOR i:=1 TO 7 DO F[i]:=F[i+2]; Set[i]:=Set[i+2]; END;
Set[8]:=FALSE; Set[9]:=FALSE;
END;
|3: F[3]:=Func; Set[3]:=TRUE;
IF (LowF = Alpha) THEN (* Verschiebe nach rechts *)
FOR i:=1 TO 7 DO F[10-i]:=F[8-i]; Set[10-i]:=Set[8-i]; END;
Set[1]:=FALSE; Set[2]:=FALSE;
END;
|4: IF NOT exact THEN EXIT END; (* !!! *)
tol := 0.25*(Falt - FMin);
IF (tol < Noise) THEN tol:=Noise; END;
IF (ABS(Func - Fquad) < tol) THEN EXIT END; (* !!! *)
IF (LowF = Alpha) OR (LowF = ANext) THEN
flag:=TRUE;
ELSE
Half();
END;
ELSE
(* nichts *)
END; (* CASE *)
IF (iret = 1) OR (iret > 4) OR flag THEN
Delta := ABS(Alpha - ANext);
back := FALSE;
FOR i:=1 TO 9 DO Set[i]:=FALSE; END;
F[5]:=FMin; Set[5]:=TRUE;
IF (Alpha < ANext) THEN k:=7; ELSE k:=3; END;
F[k]:=FNext; Set[k]:=TRUE;
END;
END; (* IF Fehler *)
LOOP (* iret *)
IF NOT Set[7] THEN
IF NOT Set[3] THEN
IF back THEN
LowF:=Alpha - Delta;
iret:=3; EXIT; (* !!! *)
ELSE
LowF:=Alpha + Delta;
iret:=2; EXIT; (* !!! *)
END;
ELSIF NOT Set[1] THEN
LowF:=Alpha + Delta;
iret:=5; EXIT; (* !!! *)
END;
(* Extrpoliere unter Ausnutzung der Funktionswerte F[1], *)
(* F[3] und F[5]. Wenn F[9] voraussichtich gr"o3er als F[7] *)
(* sein wird, verg"o3ere die Schrittweite, ist F[4] voraus- *)
(* sichtlich kleiner als F[5], verkleinere die Schrittweite. *)
back := 4.0*F[3] < (F[1] + 3.0*F[5]);
flag := Set[9] OR hold; (* Logik unter dem Galgen. *)
IF (NOT flag AND (5.0*F[3] <= (2.0*F[1] + 3.0*F[5])))
OR flag
THEN
IF Set[4] OR (6.0*F[3] > (F[1] + 0.5*F[5])) THEN
LowF:=Alpha + Delta;
iret:=5; EXIT; (* !!! *)
END;
Delta:=0.5*Delta;
hold:=TRUE;
F[1]:=F[3]; Set[1]:=Set[3];
F[3]:=F[4]; Set[3]:=Set[4]; Set[4]:=FALSE;
F[9]:=F[7]; Set[9]:=Set[7];
F[7]:=F[6]; Set[7]:=Set[6]; Set[6]:=FALSE;
Set[2]:=FALSE;
Set[8]:=FALSE;
ELSE
Half(); (* ??? *)
END;
ELSIF Set[3] THEN
(* Interpoliere unter Ausnutzung der Funktionswerte F[3], *)
(* F[5] und F[7]. Entweder versuche sofort einen quadra- *)
(* tischen Ausgleich oder reduziere zuerst die Schrittweite. *)
(* Ist Amin = 0, kann die Schrittweite um die Faktoren 2,4,8 *)
(* oder 16 reduziert werden, anderenfalls um einen Faktor 2. *)
IF (F[7] > F[3]) THEN Dum:=F[7]; ELSE Dum:=F[3]; END;
Spread:= Dum - F[5];
tol := 0.5*(Falt - FMin);
IF (tol < 1.0E-03*Noise) THEN tol:=1.0E-03*Noise; END;
IF NOT exact THEN tol:=2.0*tol; END;
Dum := F[3] - 2.0*F[5] + F[7];
Curve := Dum / (2.0*Delta*Delta);
IF (ABS(Dum) < Small) THEN Dum:=Small; END; (* 13.06.16 *)
Dum := (F[3] - F[7]) / Dum;
back := Dum < 0.0;
IF (Spread < tol) THEN
Fquad:= F[5] - 0.125*Dum*(F[3] - F[7]);
LowF:=Alpha + 0.5*Dum*Delta;
iret:=4; EXIT; (* !!! *)
END;
IF (ABS(Dum) < 0.75) OR (Alpha = 0.0) THEN
F[1]:=F[4]; Set[1]:=Set[4]; Set[4]:=FALSE;
F[9]:=F[6]; Set[9]:=Set[6]; Set[6]:=FALSE;
Set[2]:=FALSE; Set[3]:=FALSE;
Set[7]:=FALSE; Set[8]:=FALSE;
Delta:=0.25*Delta;
IF (ABS(Dum) < 0.3750) THEN
Set[1]:=FALSE; Set[9]:=FALSE;
Delta:=0.5*Delta;
END;
IF (ABS(Dum) < 0.1875) THEN Delta:=0.5*Delta; END;
ELSE
Half();
END;
ELSIF Set[9] THEN
back := 4.0*F[7] > (F[9] + 3.0*F[5]);
IF (NOT (Set[1] OR hold)) AND
(5.0*F[7] > (2.0*F[9] + 3.0*F[5]))
THEN
Double();
ELSIF (NOT Set[6]) AND (6.0*F[7] <= (F[9] + 5.0*F[5])) THEN
Half();
ELSE
LowF:=Alpha - Delta; iret:=3; EXIT; (* !!! *)
END;
ELSE
LowF:=Alpha - Delta; iret:=3; EXIT; (* !!! *)
END; (* IF Set[3] *)
END; (* LOOP iret *)
END; (* LOOP Funktionsauswertungen *)
FOR i:=0 TO NVar-1 DO X[i]:=Xalt[i] + Alpha*dX[i]; END;
Func:=FMin;
END Search;
PROCEDURE Trudge( Funkt : MinProc;
VAR X : ARRAY OF LONGREAL;
NVar : CARDINAL;
VAR Func : LONGREAL;
TolF,TolR : LONGREAL;
VAR Alpha : LONGREAL;
VAR V : SuchMatrix;
Noise : LONGREAL;
VAR GeoOptParam : GeoOptParameter;
VAR konvergiert : BOOLEAN;
TimLim : LONGCARD); (* Zeitlimit in Sekunden *)
(*
* COMMON/EXPOPT/ F,FNOT,TOLF,R(10),RNOT(10),TOLR
* COMMON/TRUDGE/ DR(10),ALPH(10),CURV(10),V(10,10),ALPHA,CURVE,
* FNOISE,NDIR,NDIRMX
*)
CONST genau = 1.0E-10; (* F"ur Diagonalisierungsroutine HhQL *)
ord = 1; (* s.o. *)
VAR exact : BOOLEAN;
Time,T0 : LONGCARD;
i,j,l,jj,k1,DirK,DirJ : CARDINAL;
IFehl,MaxL : CARDINAL;
AlfSet,MaxC,dum : LONGREAL;
Dot,DotMax,sqrDR : LONGREAL;
DeltaR,DeltaF,Falt : LONGREAL;
Curve : LONGREAL;
Alph,Curv,Xalt,U,P : ARRAY [1..MaxVar] OF LONGREAL;
dX : Deklera.VEKTOR;
G,Q : POINTER TO Deklera.MATRIX;
BEGIN
Fehler:=FALSE;
NEW(G); NEW(Q); (* Q : Eigenvektormatrix von G *)
IF (G = NIL) OR (Q = NIL) THEN
Fehlerflag:="Kein Freispeicher vorhanden (Trudge) !";
Fehler:=TRUE; ErrOut(Fehlerflag);
END;
IF (NVar > MaxVar) THEN
Fehlerflag:=
"Maximalzahl der variable Parameter ueberschritten (Trudge) !";
Fehler:=TRUE;
END;
IF Fehler THEN
IF (G # NIL) THEN DISPOSE(G); END;
IF (Q # NIL) THEN DISPOSE(Q); END;
RETURN;
END;
Falt:=Func; (* 19.12.15 *)
IF (TolR = 0.0) THEN TolR :=5.0E-02; END;
IF (TolF = 0.0) THEN TolF :=1.0E-03; END;
IF (Noise = 0.0) THEN Noise :=5.0E-04; END;
IF (GeoOptParam.LDirK) THEN (* >= *)
DirK:=GeoOptParam.DirK;
IF GeoOptParam.ReStart THEN Funkt(X,NVar,Func); END;
FOR i:=1 TO NVar DO
IF (Alph[i] <= 0.0) THEN Alph[i]:=0.1; END;
END;
DirJ:=GeoOptParam.DirJ;
IF (DirJ <= 0) OR (DirJ > NVar) THEN DirJ := NVar; END;
GeoOptParam.ReStart := GeoOptParam.ReStart AND (DirK # NVar);
ELSE (* Initialisiere V als Einheitsmatrix, setze Standartwerte *)
DirK := 0;
(* DirJ := 0; 31.01.15 - retundant *)
GeoOptParam.ReStart:=FALSE;
Funkt(X,NVar,Func);
FOR i:=1 TO NVar DO
Alph[i] := 0.1;
FOR j:=1 TO NVar DO V[i,j]:=0.0; END;
Curv[i] := 0.0;
V[i,i] := 1.0
END;
DirJ := NVar;
END;
GetUserTime(T0); konvergiert:=FALSE;
LOOP
REPEAT (* Suche entlag der Suchrichtung j nach einem lokalen Minimum. *)
FOR i:=1 TO NVar DO dX[i] := V[i,DirJ]; END;
exact := DirJ <= DirK;
Alpha := Alph[DirJ];
AlfSet := Alpha;
Search(X,dX,NVar,Alpha,Funkt,Func,Noise,exact,Curve,IFehl);
IF (IFehl # 0) THEN ErrOut(Fehlerflag); END;
Alpha := ABS(Alpha);
IF (Alpha < 0.50*AlfSet) THEN AlfSet := 0.5*AlfSet; END;
IF (Alpha < 0.25*AlfSet) THEN AlfSet := 0.5*AlfSet; END;
IF (Alpha > 2.00*AlfSet) THEN AlfSet := 2.0*AlfSet; END;
IF (Alpha > 4.00*AlfSet) THEN AlfSet := 2.0*AlfSet; END;
Curv[DirJ] := Curve;
Alph[DirJ] := AlfSet;
GetUserTime(Time); IF ((Time - T0) > TimLim) THEN EXIT; END;
DEC(DirJ);
UNTIL (DirJ = 0);
IF NOT GeoOptParam.ReStart THEN
FOR i:=1 TO NVar DO Xalt[i] := X[i-1]; END;
GeoOptParam.ReStart:=TRUE ;
Falt := Func;
DirK := 0;
DirJ := NVar;
END;
(* Erzeuge neue konjugierte Suchrichtungen, wenn keine Konvergenz. *)
sqrDR := 0.0;
FOR i:=1 TO NVar DO
dum := Xalt[i];
dX[i] := dum;
sqrDR:=sqrDR + dum*dum;
END;
AlfSet := sqrt(sqrDR);
Alpha := AlfSet;
IF (Alpha # 0.0) THEN
FOR i:=1 TO NVar DO
U[i] := dX[i] / AlfSet;
dX[i] := U[i];
END;
Search(X,dX,NVar,Alpha,Funkt,Func,Noise,TRUE,Curve,IFehl);
IF (IFehl # 0) THEN ErrOut(Fehlerflag); END;
Alpha := ABS(Alpha);
END; (* IF *)
DeltaF := Falt - Func;
Falt := Func;
sqrDR:=0.0;
FOR i:=1 TO NVar DO
dum := X[i-1] - Xalt[i];
sqrDR:=sqrDR + dum*dum;
Xalt[i] := X[i-1];
END;
DeltaR := sqrt(sqrDR / VAL(LONGREAL,NVar));
Funkt(X,NVar,Func);
IF (DeltaF < TolF) OR (DeltaR < TolR) THEN (* Konvergiert ! *)
konvergiert:=TRUE;
EXIT;
END;
IF (GeoOptParam.LDirK) THEN
DirK := 1;
k1 := 0;
FOR i:=1 TO NVar DO X[i-1]:=dX[i]; END;
dX[NVar] := Curve;
P[1] := Alpha;
ELSE (* Berechne k+1 Q-Vektoren, Q_i \cdot V_j = \delta_{i,j}, j=1,k *)
FOR j:=1 TO DirK DO
Dot := 0.0;
FOR i:=1 TO NVar DO Dot:=Dot + dX[i]*V[i,j]; END;
P[j] := Dot;
FOR i:=1 TO NVar DO dX[i]:=dX[i] - Dot*V[i,j]; END;
END;
Dot:=0.0;
FOR i:=1 TO NVar DO Dot:=Dot + dX[i]*U[i]; END;
k1 := DirK;
DirK := VAL(CARDINAL,k1+1);
GeoOptParam.ReStart := DirK # NVar;
FOR i:=1 TO NVar DO
dum := dX[i] / Dot;
Q^[i,DirK] := dum;
FOR j:=1 TO k1 DO Q^[i,j] := P[j]*dum; END;
END;
FOR i:=1 TO NVar DO (* G = \sum_l^{k+1} Q_l Curv_l (Q_l)^\dag *)
dum := Curve*Q^[i,DirK];
FOR j:=1 TO NVar DO G^[i,j] := dum*Q^[j,DirK]; END;
END;
FOR l:=1 TO k1 DO
FOR i:=1 TO NVar DO
dum := Curv[l]*Q^[i,l];
FOR j:=1 TO NVar DO G^[i,j]:=G^[i,j]*dum; END;
END;
END;
(*
Balance(G^,NVar,1,NVar,Scale);
*)
HhQL(G^,Q^,dX,NVar,genau,ord); (* Diagonalisiere G *)
(*
BalBak(NVar,1,NVar,Scale,NVar,Q^);
*)
Stuerz(Q^,NVar);
FOR j:=1 TO DirK DO
jj := NVar+1-j;
P[j] := Alph[j];
Dot := 0.0;
FOR i:=1 TO NVar DO
X[i-1] := Q^[i,jj];
Dot:=Dot + X[i-1]*V[i,j];
END;
Dot := ABS(Dot);
LOOP
IF (Dot > 0.7) THEN EXIT END;
MaxL := j;
DotMax := Dot;
FOR l:=1 TO k1 DO
IF (l # j) THEN
Dot := 0.0;
FOR i:=1 TO NVar DO Dot:=Dot + X[i-1]*V[i,l]; END;
Dot := ABS(Dot);
IF (Dot > 0.7) THEN P[j] := Alph[l]; EXIT END; (* !!! *)
IF (Dot > DotMax) THEN DotMax := Dot; MaxL := l; END;
END;
END;
P[j] := Alph[MaxL];
Dot := 0.0;
FOR i:=1 TO NVar DO Dot:=Dot + X[i-1]*U[i]; END;
IF (ABS(Dot) > DotMax) THEN P[j] := Alpha; END;
EXIT;
END; (* LOOP *)
END; (* FOR j *)
(* Setze neue konjugierte Vektoren in die V-Matrix *)
(* ein au3er dem k-te Vektor, der in R gesichert wird. *)
FOR j:=1 TO k1 DO
jj := NVar-j+1;
FOR i:=1 TO NVar DO V[i,j] := Q^[i,jj]; END;
Alph[j] := P[j];
Curv[j] := dX[jj]
END;
END; (* IF DirK *)
IF (DirK # NVar) THEN
(* Orthogonalisiere die verbleibenden Such- *)
(* richtungen zu den konjungierten Richt. *)
FOR j:=DirK TO NVar DO
Dot := 0.0;
FOR i:=1 TO NVar DO Dot :=Dot + V[i,j]*X[i-1]; END;
FOR i:=1 TO NVar DO V[i,j]:=V[i,j] - Dot *X[i-1]; END;
IF (DirK # 1) THEN
FOR l:=1 TO k1 DO
Dot := 0.0;
FOR i:=1 TO NVar DO Dot:=Dot + V[i,j]*V[i,l]; END;
FOR i:=1 TO NVar DO V[i,j]:=V[i,j] - Dot*V[i,l]; END
END;
END;
END; (* FOR j *)
MaxC:=0.0; (* Setze die G-Matrix f"ur die *)
FOR l:=DirK TO NVar DO (* verbleibenden Suchrichtungen auf. *)
IF (Curv[l] > MaxC) THEN MaxC := Curv[l]; END;
END;
MaxC := 8.0*MaxC;
FOR i:=1 TO NVar DO
FOR j:=1 TO i DO
dum := -MaxC*X[i-1]*X[j-1];
FOR l:=DirK TO NVar DO dum:=dum + Curv[l]*V[i,l]*V[j,l]; END;
FOR l:=1 TO k1 DO dum:=dum - MaxC *V[i,l]*V[j,l]; END;
G^[i,j]:=dum; G^[j,i]:=dum;
END
END; (* FOR i *)
END; (* IF DirK # NVar *)
FOR i:=1 TO NVar DO (* Setze die k'te konjugierte Such- *)
V[i,DirK] := X[i-1]; (* richtung wieder in die V-Matix ein. *)
X[i-1] := Xalt[i];
END;
Alph[DirK] := P[DirK];
Curv[DirK] := dX[NVar+1-DirK];
IF (DirK # NVar) THEN
IF debug THEN
TIO.WrStr("DirK # NVar"); TIO.WrLn;
FOR i:=1 TO NVar DO
FOR j:=1 TO NVar DO
TIO.WrLngReal(G^[i,j],14,-5);
END;
TIO.WrLn;
END;
TIO.WrLn;
END;
HhQL(G^,Q^,dX,NVar,genau,ord); (* Diagonalisiere G *)
Stuerz(Q^,NVar);
k1 := DirK+1; (* Setze die verbeleibenden Such- *)
FOR j:=k1 TO NVar DO (* richtungen in die V-Matrix ein. *)
Curv[j] := dX[j];
FOR i:=1 TO NVar DO V[i,j] := Q^[i,j]; END;
END;
END; (* IF *)
DirJ := NVar;
GetUserTime(Time); IF ((Time - T0) > TimLim) THEN EXIT END;
END; (* LOOP *)
(* Ausschrieben der f"ur einen Restart notwendigen Daten hier *)
(* oder in der rufenden Routine einf"ugen. *)
(* (Unter anderem DirK und DirJ) aber vor allem noch die lokalen *)
(* Felder *)
(* Alph,Curv,Xalt,U,P : ARRAY [1..MaxVar] OF LONGREAL; *)
(* dX : Deklera.VEKTOR; *)
(* G,Q : POINTER TO Deklera.MATRIX; *)
(* und eventuell Falt und andere Parameter *)
GeoOptParam.DirK:=DirK;
GeoOptParam.DirJ:=DirJ;
DISPOSE(G); DISPOSE(Q);
END Trudge;
(*=================== SECTION FletcherPowell =============================*)
PROCEDURE Sqrt(x : LONGREAL) : LONGREAL;
BEGIN (* kleiner Fusch ... *)
IF (ABS(x) <= Small*Small) THEN
RETURN Small;
ELSE
RETURN LongMath.sqrt(ABS(x));
END;
END Sqrt;
(*===========================================================================
*
* \centerline{Beschreibung der Prozedur {\tt NumAbleit}.}
* % ******************************************
* %
* \bigskip
* Die 1. Ableitung einer Funktion F(X,n) mit $n$ Parametern $X_i$
* l\"a\3t sich numerisch nach
*
* $$
* {\partial F \over \partial X}
* \;\approx\;
* \half\; {{F(X^+) \minus F(X^0)} \over \Delta X}
* $$
*
* berechnen, wobei $X^+ \eq X^0 + \Delta X$ ist. Analog dazu l\"a\3t
* sich die * 2. Ableitung nach
*
* $$
* {\partial^2 F \over \partial X^2}
* \;\approx\;
* {E(X^+) \minus E(X^-) \over {\Delta X}^2 }
* $$
*
* berechnen. Wenn aus einem vorhergehenden Lauf die 2. Ableitung
* $\Delta_2$ * zumindest n\"aherungsweise bekannt ist, kann man die
* 1. Ableitung noch um * den Faktor $- \half \Delta_2 \Delta X$
* korrigieren.
*
* $$
* {\partial F \over \partial X}
* \;\approx\;
* \half \; {{F(X^+) \minus F(X^0)} \over \Delta X} -
* \half \; \Delta_2 \Delta X
* $$
* \smallskip
*
* Diese Option kann in der Prozedur {\tt NumAbleit} mit dem Parameter
* {\tt Korrektur} angew\"ahlt werden.
*
*=========================================================================*)
PROCEDURE NumAbleit( NVar : CARDINAL;
E0 : LONGREAL;
VAR X0 : ARRAY OF LONGREAL; (* ==> *)
VAR DelX0 : ARRAY OF LONGREAL; (* ==> *)
FunkNX : MinProc;
VAR Ableit1 : ARRAY OF LONGREAL; (* <== *)
VAR Ableit2 : ARRAY OF LONGREAL; (* <==> *)
ZweiteAbl : BOOLEAN;
Korrektur : BOOLEAN);
(*------------------------------------------------------------------*)
(* Berechnet auf numerischem Wege die 1. und, wenn mit ZweiteAbl *)
(* gefordert, die Diagonalelemente der 2. Ableitung einer Funktion *)
(* mit NVar Parametern. *)
(* F"ur die erste Ableitung werden NVar zus"atzliche Funktionsaus- *)
(* wertungen ben"otigt, f"ur die 2. Ableitung 2*NVar Auswertungen. *)
(* *)
(* NVar : Anzahl der Parameter der Funktion FunkNX. *)
(* E0 : Energier"uckgabewert der Funktion FunkNX mit dem *)
(* Parameter X0. *)
(* X0 : Parametersatz der Funktion FunkNX. *)
(* DelX0 : Schrittweiten im Parametervektor X0 bei der Be- *)
(* rechnung der Ableitungen. *)
(* Ableit1 : Berechnete 1. Ableitung. *)
(* Ableit2 : Berechnete 2. Ableitung oder Eingabeparameter, *)
(* wenn ZweiteAbl=FALSE und Korrektur=TRUE. *)
(* FunkNX : Funktion, deren Ableitungen berechnet werden *)
(* sollen. FuktNX(NVar,X0,E0) berechnet z.B. den *)
(* Funktionswert E0 mit NVar Parametern X0. *)
(* ZweiteAbl : Wenn auf TRUE gesetzt, werden auch die Diagonal- *)
(* elemente der 2. Ableitung berechnet. *)
(* Korrektur : Wenn auf TRUE gesetzt, wird die 1. Ableitung *)
(* mithilfe der 2. Ableitung korrigiert. Dazu mu3 die *)
(* 2. Ableitung (z.B. aus einem vorhergehenden *)
(* Aufruf) bekannt sein. Diese Option kann NICHT *)
(* bei gesetztem Parameter ZweiteAbl verwandt werden. *)
(*------------------------------------------------------------------*)
CONST eps = Small;
VAR index,i : CARDINAL;
Ep,Em : LONGREAL;
Xakt : ARRAY [0..MaxVar-1] OF LONGREAL;
Nenner : LONGREAL;
BEGIN
Em := 0.0; (* Wg. Compilerwarnung *)
FOR i:=0 TO NVar-1 DO Xakt[i] := X0[i]; END;
FOR index:=0 TO NVar-1 DO
Xakt[index]:=X0[index] + DelX0[index];
FunkNX(Xakt,NVar,Ep);
IF Fehler THEN RETURN END;
IF ZweiteAbl THEN (* Berechne E^- f"ur die 2. Ableitung. *)
Xakt[index]:=X0[index] - DelX0[index];
FunkNX(Xakt,NVar,Em);
IF Fehler THEN RETURN END;
END;
Xakt[index]:=X0[index]; (* Setze Xakt zur"uck. *)
Nenner := DelX0[index];
IF (Nenner < eps) THEN Nenner:=eps; END;
IF ZweiteAbl THEN (* Berechne 1. und 2. Ableitung. *)
Ableit1[index] := 0.5*(Ep - Em) / Nenner;
Ableit2[index] := (Ep + Em - 2.0*E0) / (Nenner*Nenner);
ELSE (* Berechne nur die 1. Ableitung. *)
Ableit1[index] := (Ep - E0) / (2.0*Nenner);
IF Korrektur THEN (* Koorigiere mit vorhandener 2. Ableitung. *)
Ableit1[index]:=Ableit1[index] - 0.5*Ableit2[index]*DelX0[index];
END;
END;
END; (* FOR index *)
END NumAbleit;
PROCEDURE Suche( NVar : CARDINAL;
VAR X0 : ARRAY OF LONGREAL; (* <==> *)
VAR X0alt : ARRAY OF LONGREAL; (* <== *)
FunkNX : MinProc;
VAR E0 : LONGREAL; (* ==> *)
VAR Ableit1 : ARRAY OF LONGREAL; (* ==> *)
VAR Ableit2 : ARRAY OF LONGREAL; (* ==> *)
VAR H : SuchMatrix); (* ==> *)
(*----------------------------------------------------------*)
(* Berechnet ein neues X0 derart, da3 der Funktionswert der *)
(* Funktion FunkNX kleiner ist als zuvor. Ben"otigt dazu *)
(* zwei Funktionsauswertungen. *)
(*----------------------------------------------------------*)
VAR i,j : CARDINAL;
sqrtA2,Xii : LONGREAL;
Xakt,Xi : ARRAY [0..MaxVar-1] OF LONGREAL;
XiAbleit1 : ARRAY [0..MaxVar-1] OF LONGREAL;
E1,E2 : LONGREAL;
Alpha : LONGREAL;
Nenner : LONGREAL;
sqrtAbleit2i : LONGREAL;
sqrtAbleit2 : ARRAY [0..MaxVar-1] OF LONGREAL;
BEGIN
FOR i:=0 TO NVar-1 DO
(* Transformiere X0 und Ableit1 in den \Xi-Raum. *)
sqrtA2 := Sqrt(Ableit2[i]);
IF (sqrtA2 < Small) THEN sqrtA2:=Small; END;
XiAbleit1[i]:=Ableit1[i] / sqrtA2;
X0alt[i] := X0[i]; (* Sicher X0 in X0alt. *)
X0[i] := X0[i] * sqrtA2;
END;
FOR i:=1 TO NVar DO (* Berechne Xi. *)
Xii:=0.0;
FOR j:=1 TO NVar DO Xii:=Xii - H[i,j]*XiAbleit1[j-1]; END;
Xi[i-1]:=Xii;
END;
FOR i:=0 TO NVar-1 DO (* Vermeide Divisionen durch Null *)
sqrtAbleit2i := Sqrt(Ableit2[i]);
IF (sqrtAbleit2i < Small) THEN sqrtAbleit2i:=Small;
END;
sqrtAbleit2[i] := sqrtAbleit2i;
END;
FOR i:=0 TO NVar-1 DO (* 1. Auswertung. *)
Xakt[i] := (X0[i] + 1.0*Xi[i]) / sqrtAbleit2[i];
END;
FunkNX(Xakt,NVar,E1);
FOR i:=0 TO NVar-1 DO (* 2. Auswertung. *)
Xakt[i] := (X0[i] + 2.0*Xi[i]) / sqrtAbleit2[i];
END;
FunkNX(Xakt,NVar,E2);
Nenner:=(2.0*(E0 - 2.0*E1 + E2));
IF (ABS(Nenner) < Small) THEN Nenner:=Small; END;
Alpha := 1.0 - (E2 - E0) / Nenner;
(* Berechne neues X0 im \Xi-Raum. *)
FOR i:=0 TO NVar-1 DO X0[i]:=X0[i] + Alpha*Xi[i]; END;
(* Transfomiere X0 aus dem \Xi-Raum zur"uck. *)
FOR i:=0 TO NVar-1 DO X0[i]:=X0[i] / sqrtAbleit2[i]; END;
END Suche;
PROCEDURE BerHess( NVar : CARDINAL;
VAR X0 : ARRAY OF LONGREAL; (* ==> *)
VAR X0alt : ARRAY OF LONGREAL; (* ==> *)
VAR Ableit1 : ARRAY OF LONGREAL; (* ==> *)
VAR Ableit2 : ARRAY OF LONGREAL; (* ==> *)
VAR Abl1Alt : ARRAY OF LONGREAL; (* ==> *)
VAR H : SuchMatrix); (* <==> *)
(*--------------------------------------------------------------*)
(* Berechnet eine neue, approximative Hesssche Matrix. *)
(* Parameter : Siehe die anderen Prozeduren in diesem Modul. *)
(* Abl1Alt : 1. Ableitung der Funktion aus dem vorhergehenden *)
(* Zyklus. *)
(*--------------------------------------------------------------*)
VAR i,j : CARDINAL;
dA1xH,HxdA1 : ARRAY [0..MaxVar-1] OF LONGREAL;
DeltaX0,DeltaA1 : ARRAY [0..MaxVar-1] OF LONGREAL;
HBar,dX0xA1,sqrtA2 : LONGREAL;
tmp : LONGREAL;
BEGIN
(* Transformiere X0, X0alt und Ableit1 in den \Xi-Raum *)
(* und berechne DeltaX0 und DeltaA1. *)
FOR i:=0 TO NVar-1 DO
(*
sqrtA2 := Sqrt(Ableit2[i]);
*)
tmp := Ableit2[i];
IF (ABS(tmp) < Small) THEN
sqrtA2 := Small;
ELSE
sqrtA2 := sqrt(tmp);
END;
X0[i] := X0[i] * sqrtA2;
X0alt[i] := X0alt[i] * sqrtA2;
Ableit1[i] := Ableit1[i] / sqrtA2;
DeltaX0[i] := X0[i] - X0alt[i]; (* Sigma *)
DeltaA1[i] := Ableit1[i] - Abl1Alt[i]; (* Y *)
END;
FOR i:=0 TO NVar-1 DO
HxdA1[i]:=0.0; (* H^i | \y^i> *)
dA1xH[i]:=0.0; (* <y^j | H^i *)
FOR j:=0 TO NVar-1 DO
HxdA1[i]:=HxdA1[i] + H[i+1,j+1]*DeltaA1[j];
dA1xH[i]:=dA1xH[i] + H[i+1,j+1]*DeltaA1[j];
END;
END;
dX0xA1:=0.0; (* <\sigma^i | y^i> *)
HBar :=0.0; (* <y^i|H^i|y^i> *)
FOR i:=0 TO NVar-1 DO
dX0xA1:=dX0xA1 + DeltaX0[i]*DeltaA1[i];
HBar :=HBar + DeltaA1[i]*HxdA1[i];
END;
IF (ABS(HBar) < Small) THEN HBar:=Small; END;
IF (ABS(dX0xA1) < Small) THEN dX0xA1:=Small; END;
(* Bereche aus den oben erhaltenen Informationen und der *)
(* alten Hessschen Matrix die Hesssche Matrix des aktuellen *)
(* Optimierungszyklus (Gleichung (7) in FP-Artikel) *)
FOR i:=1 TO NVar DO
FOR j:=1 TO NVar DO
H[i,j]:=H[i,j] + ((DeltaX0[i-1]*DeltaX0[j-1]) / dX0xA1)
- ((HxdA1 [i-1]*dA1xH [j-1]) / HBar );
END;
END;
FOR i:=0 TO NVar-1 DO
(* Transformiere X0, X0alt und Ableit1 aus dem \Xi-Raum zur"uck. *)
sqrtA2 := Sqrt(Ableit2[i]);
X0[i] := X0[i] / sqrtA2;
X0alt[i] := X0alt[i] / sqrtA2;
Ableit1[i] := Ableit1[i] * sqrtA2;
END;
END BerHess;
PROCEDURE Konvergenz( NVar : CARDINAL;
VAR Ableit1 : ARRAY OF LONGREAL;
genau : LONGREAL) : BOOLEAN;
VAR i : CARDINAL;
Sum : LONGREAL;
BEGIN
Sum:=0.0; FOR i:=0 TO NVar-1 DO Sum:=Sum + Ableit1[i]*Ableit1[i]; END;
Sum := Sqrt(Sum / VAL(LONGREAL,NVar));
IF debug THEN
Write2("//,S,E20.6,//",'RMS-Kraft : ',Sum);
END;
RETURN (Sum < genau);
END Konvergenz;
PROCEDURE FletcherPowell( NVar : CARDINAL;
VAR X0 : ARRAY OF LONGREAL;
VAR DelX0 : ARRAY OF LONGREAL;
E0 : LONGREAL;
FunkNX : MinProc;
VAR H : SuchMatrix;
genau : LONGREAL;
Restart : BOOLEAN;
MaxZyklus : CARDINAL);
VAR Zyklus : CARDINAL;
Ableit1,Ableit2 : ARRAY [0..MaxVar-1] OF LONGREAL;
X0alt,Abl1Alt : ARRAY [0..MaxVar-1] OF LONGREAL;
PROCEDURE WriteFPDaten(); (* Ausgabe von Iterationsdaten. *)
VAR i : CARDINAL;
dX,dA1 : LONGREAL;
BEGIN
Write2("//,S,C2,//",'Optimierungsdaten in Zyklus : ',Zyklus);
Write5("//,13X,2(A,10X),2(A,5X),4X,S,//",
"X","dX","1. Ableitung","dA1","2. Ableitung");
FOR i:=0 TO NVar-1 DO
dX := X0[i] - X0alt[i];
dA1 := Ableit1[i] - Abl1Alt[i];
IF (Zyklus = 0) THEN dX:=0.0; dA1:=0.0; END;
Write1("C3,2X",i+1);
Write5("2(F14.8,2X,F11.8,2X),F14.8,/",
X0[i],dX,Ableit1[i],dA1,Ableit2[i]);
END;
END WriteFPDaten;
CONST MinAbw = 0.0010;
MaxAbw = 0.0160;
VAR i,j : CARDINAL;
X0A2 : ARRAY [0..MaxVar-1] OF LONGREAL;
DeltaXMax,dX : LONGREAL;
DelXA2Max : LONGREAL;
BEGIN
IF (MaxZyklus = 1) THEN MaxZyklus:=5; END;
IF NOT Restart THEN
FOR i:=1 TO NVar DO (* Initialisiere H mit E-Matrix. *)
FOR j:=1 TO NVar DO H[i,j]:=0.0; END;
H[i,i]:=1.0;
END;
END;
FOR i:=0 TO NVar-1 DO X0A2[i]:=X0[i]; END;
(* Berechne numerische erste und zeite Ableitung *)
NumAbleit(NVar,E0,X0,DelX0,FunkNX,Ableit1,Ableit2,TRUE,FALSE);
Zyklus:=0;
IF debug THEN WriteFPDaten(); END;
LOOP
INC(Zyklus);
(* Berechne verbessertes X0. *)
Suche(NVar,X0,X0alt,FunkNX,E0,Ableit1,Ableit2,H);
IF Konvergenz(NVar,Ableit1,genau) THEN EXIT END; (* !!! *)
FunkNX(X0,NVar,E0); (* Berechne neues E0(X0). *)
FOR i:=0 TO NVar-1 DO Abl1Alt[i]:=Ableit1[i]; END;
DeltaXMax:=0.0; DelXA2Max:=0.0;
FOR i:=0 TO NVar-1 DO
dX := X0[i] - X0A2[i];
IF (ABS(dX) > DelXA2Max) THEN DelXA2Max:=ABS(dX); END;
dX := X0[i] - X0alt[i];
IF (ABS(dX) > DeltaXMax) THEN DeltaXMax:=ABS(dX); END;
END;
IF debug THEN
Write2("//,S,F14.8,//",'Maximale Abweichung in X : ',DeltaXMax);
END;
(*
* Hier ist was falsch - Das kann nicht richtig sein
* DelXA2Max:=2.0*MaxAbw;
*)
DelXA2Max:=2.0*DeltaXMax;
IF (DeltaXMax < MinAbw) OR (DelXA2Max > MaxAbw) THEN
(* DeltaXMax = \sqrt{{\sum_i { X0_i - X0_i^{alt} }^2 \over NVar }} *)
(* Berechne neue 1. und 2. Ableitung. Wenn DelXA2Max > als MaxAbw *)
(* ist es nicht mehr gerechtfertigt, die 2. Ableitung als konstant *)
(* anzusehen. Ist DeltaXMax < MinAbw ist KEIN wesentlicher Fort- *)
(* schritt erricht. *)
FOR i:=0 TO NVar-1 DO X0A2[i]:=X0[i]; END;
(* Berechne numerische erste und zeite Ableitung *)
NumAbleit(NVar,E0,X0,DelX0,FunkNX,Ableit1,Ableit2,TRUE,FALSE);
IF debug THEN
Write1("//,S,//",'Neue 2. Ableitung berechnet !');
END;
ELSE (* Berechne nur die neue 1. Ableitung mit Korrektur *)
NumAbleit(NVar,E0,X0,DelX0,FunkNX,Ableit1,Ableit2,FALSE,TRUE);
END;
IF debug THEN WriteFPDaten(); END;
BerHess(NVar,X0,X0alt,Ableit1,Ableit2,Abl1Alt,H);
IF debug THEN
Write2("//,S,C2,/",'Hessche Matrix in Zyklus : ',Zyklus);
MATaus(Output,NVar,NVar,H,80,12,1);
END;
IF (Zyklus = MaxZyklus) THEN EXIT; END;
END; (* LOOP *)
IF (Zyklus = MaxZyklus) THEN
Write5("//,S,S,/,S,C2,S,//",'Maximalzahl der Zyklen "uberschritten, ',
'die Fletcher-Powell ','Optimierung wurde erfolglos nach ',Zyklus,
' Zyklen abgebrochen !');
ELSE
Write3("//,S,C2,S,//",
'Die Fletcher-Powell Optimierung konvergierte nach ',Zyklus,
' Zyklen !');
END;
END FletcherPowell;
BEGIN
debug:=FALSE;
END OptimLib1.