DEFINITION 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. *)
(*------------------------------------------------------------------------*)
(* Most routines are based on other author work - please see a reference *)
(* list in the implementation module and the headings of the procedures *)
(* in this definition module *)
(*------------------------------------------------------------------------*)
(* Implementation : Michael Riedl *)
(* Licence : GNU Lesser General Public License (LGPL) *)
(*------------------------------------------------------------------------*)
(* $Id: OptimLib1.def,v 1.2 2018/05/03 07:19:19 mriedl Exp mriedl $ *)
IMPORT LMathLib;
CONST MaxVar = 32; (* Maximal number of parameters to be optimized *)
TYPE Funktion1 = LMathLib.Funktion1;
FunktionN = LMathLib.FunktionN;
TYPE MinProc = PROCEDURE(VAR ARRAY OF LONGREAL, (* Parametervektor *)
CARDINAL, (* Anzahl Parameter *)
VAR LONGREAL); (* Funktionsergebniss *)
SuchMatrix = ARRAY [1..MaxVar] OF ARRAY [1..MaxVar] OF LONGREAL;
PROCEDURE FMin( A,B : LONGREAL;
Func : Funktion1;
Tol : LONGREAL;
VAR iter : CARDINAL;
VAR Min : LONGREAL);
(* ---------------------------------------------------------------*)
(* An approximation x to the point where Func=F(x) attains *)
(* a minimum on the interval (a,b) is determined. *)
(* *)
(* Parameter *)
(* *)
(* A : left endpoint of initial interval *)
(* B : right endpoint of initial interval *)
(* Func : function procedure which evaluates F(x) for *)
(* any x in the interval (A,B) *)
(* Tol : desired length of the interval of uncertainty *)
(* of the final result (Tol .ge. 0.0d0) *)
(* *)
(* Return value *)
(* *)
(* Min : abcissa approximating the point where f(fmin) *)
(* attains a minimum *)
(* iter : number of iterations taken *)
(* *)
(* The method used is a combination of golden section search *)
(* and successive parabolic interpolation. *)
(* Convergence is never much slower than that for a fibonacci *)
(* search. If F(x) has a continuous second derivative which is *)
(* positive at the minimum (which is not at A or B), then *)
(* convergenc is superlinear, and usually of the order of *)
(* about 1.324.... *)
(* The function F(x) is never evaluated at two points closer *)
(* together than eps*abs(fmin) + (Tol/3), where eps is *)
(* approximately the square root of the relative machine *)
(* precision. If F(x) is a unimodal function and the computed *)
(* values of F(x) are always unimodal when separated by at least *)
(* eps*abs(x) + (Tol/3), then fmin approximates the abcissa of *)
(* the global minimum of F(x) on the interval (A,B) with an *)
(* error less than 3*eps*abs(fmin) + Tol. If F(x) is not *)
(* unimodal, then fmin may approximate a local, but perhaps non- *)
(* global, minimum to the same accuracy. *)
(* *)
(* This function is a slightly modified version of the Algol 60 *)
(* procedure localmin given in Richard Brent, "Algorithms for *)
(* minimization without derivatives", Prentice-Hall (1973). *)
(* *)
(* Fortran source is based on a procedures from the book *)
(* Forsythe, George E.; Malcolm, Michael A.; Moler, Cleve B.; *)
(* "Computer Methods for Mathematical Computations", *)
(* Prentice-Hall, 1977. *)
(*----------------------------------------------------------------*)
PROCEDURE FMinBr( a,b : LONGREAL;
f : Funktion1;
tol : LONGREAL;
VAR iter : CARDINAL;
VAR Min : LONGREAL);
(*----------------------------------------------------------------*)
(* Parameters and Algorithm as described for FMin. *)
(* Translated from "C" source as published in the Internet *)
(*----------------------------------------------------------------*)
PROCEDURE GloMin( a,b,c : LONGREAL;
m : LONGREAL;
e : LONGREAL;
t : LONGREAL;
f : Funktion1;
VAR x : LONGREAL) : LONGREAL;
(*----------------------------------------------------------------*)
(* glomin returns the global minuimum of the function f(x) *)
(* defined on [a,b], The procedure assumes that f \in C^2[a,b] *)
(* and f''(x) <= m \forall x \in [a,b] weaker conditions are *)
(* sufficient: see Section 7). *)
(* e and t are positive tolerance: we assume that f(x) is *)
(* computed with an absolute error bounded by e, i.e. then *)
(* | f1(f(x(1 +- macheps))) - f(x)| <= e, where macheps is the *)
(* relative machine precision. Then x and glomin are returned *)
(* so than min(f) <= min(f) + t + 2e and *)
(* min(f) - e <= glomin = f1(f(x)) <= min(f) + t + e. *)
(* c is an initial guess at x (a or b will do). The number of *)
(* function ealuations required is usually close to the least *)
(* possible, provided t is not unreasonable small *)
(* (see Section 3 to 5) *)
(* *)
(* Parameters: *)
(* *)
(* a,b : The endpoints of the interval with a < b *)
(* c : An initial guess for the global minimizer. *)
(* If no good guess is known, c = a or b is acceptable. *)
(* m : The bound on the second derivative. *)
(* e : a positive tolerance, a bound for the absolute error *)
(* in the evaluation of f(x) for any x in [A,B]. *)
(* t : a positive error tolerance. *)
(* f : f(x), a user-supplied function whose global minimum *)
(* is being sought. *)
(* *)
(* x : the estimated value of the abscissa for which f(x) *)
(* attains its global minimum value in [A,B]. *)
(* *)
(* Return value is the value of f(x) *)
(* *)
(* Algol60 version by R.Brent, adopted to Modula-2 by M.Riedl *)
(*----------------------------------------------------------------*)
PROCEDURE LocalMin( A,B : LONGREAL;
eps : LONGREAL;
T : LONGREAL;
F : Funktion1;
VAR x : LONGREAL) : LONGREAL;
(*----------------------------------------------------------------*)
(* Parameters: *)
(* *)
(* A,B : The endpoints of the interval with A < B *)
(* eps : A positive relative error tolerance. *)
(* eps should be no smaller than twice the relative *)
(* machine precision and preferably not much less *)
(* than the square root of the relative machine *)
(* precision. *)
(* T : a positive error tolerance. *)
(* F : F(x), a user-supplied function whose global minimum *)
(* is being sought. *)
(* *)
(* x : the estimated value of the abscissa for which F(x) *)
(* attains its global minimum value in [A,B]. *)
(* *)
(* Return value is the value of F(x) *)
(* *)
(* [1] Brent, Richard; "Algorithms for finding zeros and extrema *)
(* of functions without calculating derivatives", Computer *)
(* Sience Departmet, Stanfort university, Stanford (US) *)
(* (STAN-CS-71-198) (1971) *)
(*----------------------------------------------------------------*)
PROCEDURE Zero(a,b : LONGREAL;
t : LONGREAL;
f : Funktion1) : LONGREAL;
(*----------------------------------------------------------------*)
(* Zero seeks the root of a function f(x) in an interval [A,B]. *)
(* *)
(* Parameters: *)
(* *)
(* a,b : The endpoints of the interval with a < b *)
(* t : a positive error tolerance. *)
(* f : f(x), a user-supplied function whose zero *)
(* is being sought. *)
(* *)
(* x : the estimated value of the abscissa for which f(x) *)
(* attains its global minimum value in [A,B]. *)
(* *)
(* Return value is the value of x such that f(x) = 0 *)
(* *)
(* If no zero is found in the interval [a,b] the global variable *)
(* Errors.Fehler ist set to true and Errors.Fehlerflag is set *)
(* accordingly *)
(*----------------------------------------------------------------*)
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);
(*----------------------------------------------------------------*)
(* Optimierer nach Nelder-Mead, basierend auf F77 Version von *)
(* R. V. O'Neill mit Modifikationen von John Burkardt. *)
(* *)
(* Simplex function minimisation procedure due to *)
(* nelder+mead(1965), as implemented by o'neill *)
(* (1971, appl.statist. 20, 338-45), with subsequent comments *)
(* by chambers+ertel(1974, 23, 250-1), benyon(1976,25, 97) *)
(* and hill(1978, 27, 380-2) *)
(* *)
(* This routine seeks the minimum value of a user-specified *)
(* function *)
(* *)
(* The function to be minimized must be defined as a procedure *)
(* of the form *)
(* *)
(* PROCEDURE fn(VAR X : ARRAY OF LONGREAL; *)
(* n : INTEGER; *)
(* VAR Fx : LONGREAL); *)
(* *)
(* and passed as the argument FNPorc. *)
(* *)
(* This routine does not include a termination test using *)
(* the fitting of a quadratic surface. *)
(* *)
(* Parameters: *)
(* *)
(* FNProc : input,external *)
(* the name of the function which evaluates *)
(* the function to be minimized. *)
(* N : input *)
(* the number of variables. *)
(* Start(N) : input/output *)
(* on input, a starting point for the iteration. *)
(* on output, this data may have been overwritten. *)
(* XMin(N) : output *)
(* the coordinates of the point which is estimated *)
(* to minimize the function. *)
(* YnewLo : output *)
(* the minimum value of the function. *)
(* ReqMin : input *)
(* the terminating limit for the variance of *)
(* function values. *)
(* Step(n) : input *)
(* determines the size and shape of the initial *)
(* simplex. The relative magnitudes of its elements *)
(* should reflect the units of the variables. *)
(* Konvge : input *)
(* the convergence check is carried out every *)
(* Konvge iterations. *)
(* KCount : input *)
(* the maximum number of function evaluations. *)
(* ICount : output *)
(* the number of function evaluations used. *)
(* NumRes : output *)
(* the number of restarts. *)
(* IFault : output *)
(* error indicator. *)
(* 0: no errors detected. *)
(* 1: ReqMin, N, or Konvge has an illegal value. *)
(* 2: iteration terminated because KCount was *)
(* exceeded without convergence. *)
(*----------------------------------------------------------------*)
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);
(* ---------------------------------------------------------------*)
(* Praxis is for the minimization of a function in several *)
(* variables. The algorithm used is a modification of a conjugate *)
(* gradient method developed by Powell. Changes are due to Brent, *)
(* who gives an ALGOL W program. *)
(* Users who are interested in more of the details should read *)
(* *)
(* - Powell, M. J. D., 1962. An efficient method for finding *)
(* the minimum of a function of several variables without *)
(* calculating derivatives, Computer Journal 7, 155-162. *)
(* - Brent, R. P., 1973. Algorithms for minimization without *)
(* derivatives. Prentice Hall, Englewood Cliffs. *)
(* *)
(* Original Pascal version by Karl Gegenfurtner *)
(* ---------------------------------------------------------------*)
(* Parameter *)
(* *)
(* X : on input a vector containing an initial guess of the *)
(* solution. *)
(* on output X holds the final solution of the system *)
(* N : number of unknown parameters *)
(* least calculated value of the function F *)
(* F : F(X,n,fx) is the function to be minimized *)
(* fx : on input the calculaded value of F(X,n,fx) at *)
(* the initial guess *)
(* on output the calculaded value of F(X,n,fx) at *)
(* the proposed minimum *)
(* T : is the tolerance for the precision of the solution *)
(* H : is a steplength parameter and shoul be set to the *)
(* expected distance to the solution. An exceptional *)
(* large or small value of H leads to slower convergence *)
(* on the first few iterations. *)
(* Ktm : parameter to control the termination condition. Its *)
(* default value is 1 and a value of 4 leads to a very *)
(* cautious stopping criterion *)
(* Praxis returns if the criterion *)
(* 2 * ||x(k)-x(k-1)|| <= sqrt(MachEps) * ||x(k)|| + T *)
(* is fulfilled more than KTM times *)
(* Scbd : is a scaling parameter and should be set to about 10. *)
(* The default is 1 and with that value no scaling is *)
(* done at all. It is only necessary when the parameters *)
(* are scaled very different *)
(* IllC : Set to true if problem is known to be ill-conditioned *)
(* Prin : controls the printout from PRAXIS *)
(* 0 : no printout at all *)
(* 1 : only initial and final values *)
(* 2 : detailed map of the minimization process *)
(* 3 : also prints eigenvalues and eigenvectors of *)
(* the direction matices in use *)
(* (for insiders only). *)
(* MaxF : Input variable - limit the numer of function *)
(* evaluations to MaxF *)
(* nf : on output, number of function evaluations performed *)
(* iFehl : Output variable, indicating errors *)
(* 0 : no error *)
(* 1 : *)
(* 2 : *)
(* 3 : *)
(* ---------------------------------------------------------------*)
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);
(*-----------------------------------------------------------------*)
(* Fletcher-Powell Optimierungsroutine. Benutzt numerische *)
(* 1. und 2. Ableitungen. *)
(* Diese Routine sowie deren Unterroutine Suche sind stark an *)
(* die entsprechenden Routinen in GAUSSIAN 80 angelehnt, der *)
(* Beitrag deren Autoren zu diesem Programmteil sei hiermit *)
(* entsprechend gew"urdigt. *)
(* *)
(* NVar : Anzahl zu varierender Parameter. *)
(* X0 : Initialer Parametersatz, der zu optimieren ist. *)
(* DelX0 : Schrittweite bei der numerischen Berechnung der *)
(* 1. und 2. Ableitungen. *)
(* E0 : Funktionswert der Funktion FunkNX am Punkt X0. *)
(* FunkNX : Zu minimierende Funktion in NVar Variablen. *)
(* H : Hessche Matrix. Diese mu3 bei einem Restart *)
(* initialisiert sein. *)
(* genau : Abbruchkriterium. *)
(* Restart : Offensichtlich. *)
(* MaxZyklus : Maximalzahl der zul"assigen Optimierungszyklen. *)
(* Im ersten Zyklus werden 2 NVar + 3 Funktionsaus- *)
(* wertungen ben"otigt, in jedem weiteren Zyklus *)
(* NVar + 3 Funktionsauswertungen, falls die 2. Ab- *)
(* leitung nicht neu berechnet werden mu3. *)
(* *)
(* Literatur : *)
(* *)
(* Fletcher, R.; Powell, M.J.D.; Comput. J. 6, 163 (1963) *)
(* Collins, J.B.; Schleyer, P.; Binkley, J.S.; Pople, J.A.; *)
(* J. of Chem. Phys. 64, 5142 (1976) *)
(*-----------------------------------------------------------------*)
TYPE GeoOptParameter = RECORD (* LDirK set to reflect DirK < 0 in Original *)
ReStart : BOOLEAN;
LDirK : BOOLEAN;
DirK,DirJ : CARDINAL;
END;
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 *)
(*----------------------------------------------------------------*)
(* Aus exp.f von Hondo VII *)
(* *)
(* Minimiert eine Funktion von NVar Variablen mithilfe einer *)
(* modifizierten Powell-methode durch Suchen entlang *)
(* konjungierter Richtungen. *)
(* *)
(* Wird die Routine mit ReStart = TRUE aufgerufen, m"ussen alle *)
(* Werte in den entsprechenden RECORDS gesetzt worden sein. *)
(* Wird die Routine mit GeoOptParam.LDirK = TRUE aufgerufen, *)
(* m"ussen nur TimLim,TolF,TolR,Noise,NVar und X gesetzt sein. *)
(* Wird die Routine mit GeoOptParam.LDirK = FALSE aber mit *)
(* ReStart = FALSE aufgerufen, m"ussen TimLim, TolF, TolR, Noise, *)
(* NVar, X, und V entsprechend gesetzt sein. In diesem Fall wird *)
(* DirK nach dem ersten Durchgang auf Null gesetzt. *)
(* *)
(* Parameter *)
(* *)
(* Funkt : Prozedurparameter Funkt(X,NVar,Func) - *)
(* zu optimierende Funktion *)
(* X : Eingangseitig die Anfangsposition im *)
(* NVar-dimensionalen Parameterraum. *)
(* Ausgansseitung das berechnete Minimum *)
(* NVar : Anzahl der Parameter der zu optimierenden *)
(* Funktion *)
(* Func : Funktionswert. Die Varibale Func muss beim *)
(* ersten Aufruf mit dem Funktionswert an der *)
(* Stelle der Startwerte von X belegt sein. *)
(* Beim Verlassen der Routine enthaelt Func den *)
(* Funktionswert am gefundenen Optimum *)
(* TolF : Wenn die Funkttionswerte einer neuen Iteration *)
(* nur um weniger als TolF vom vorher gefundenen *)
(* Minimum abweicht --- UND --- *)
(* TolR : die Wurzel der Summe der Quadrate der *)
(* Abweichungnen in den X_i zum vorherigen Iter- *)
(* rationsschritt, gewichtet mit der Anzahl der zu *)
(* optimierenden Parameter, kleiner als TolR ist, *)
(* wird die Optimierung als konvergiert angesehen. *)
(* Vorgaben: TolF = 0.001, TolR = 0.05 *)
(* Alpha : Eingangsseitig die initiale Schrittweite fuer *)
(* die Veraenderung der Parameter X *)
(* Ausgansseitig die berechnete Schrittweite am *)
(* aufgefundenen Minumum. *)
(* V : Spalten der V-Matrix sind die jeweiligen *)
(* Suchrichtungen, wobei Spalten 1 bis DirK *)
(* annaehernd konjungiert sind *)
(* Noise : Genauigkeit der Funktionsauswertung. Variationen *)
(* die kleiner als Noise sind werden als nicht *)
(* signifikant gewertet (Vorgabe = 0.0005) *)
(* GeoOptParam : Steuerparameter *)
(* ReStart : Restart eines vorherigen Laufs *)
(* (nicht implementiert) *)
(* Restart = Restart AND DirK # NVar *)
(* LDirK : Siehe oben *)
(* DirK : Anzahl d. konjungierten Spalten in V *)
(* DirJ : Gibt die Suchrichtung an. Die *)
(* Zeile DirJ in V wird als erste Such- *)
(* richtung gewaeht wenn ein Restart *)
(* durchgefuehrt wird *)
(* konvergiert : Gibt an ob die Optimierung konvergiert ist oder *)
(* nicht (siehe auch TolF,TolR) *)
(* TimLim : Zeitlimit in Sekunden. Wird dies ueberschritten *)
(* wird die Routine abgebrochen, konvergiert wird *)
(* dann auf "falsch" gesetzt *)
(*----------------------------------------------------------------*)
(* Minimizes a function of NDir variables by a modified powell *)
(* method of searches along conjugate directions. columns of V *)
(* matrix are current search directions. Columns 1 to KDir are *)
(* (approximately) conjugate. *)
(* When subroutine is entered with Restart=TRUE then values of *)
(* LdirK and of all quantities in common should have been set *)
(* previously by the calling program. *)
(* When Trudge is called with LdirK negative then only required *)
(* values are TimLim,TolF,TolR,Noise,NDir,X and Noise for the *)
(* search. In this case initial values of KDir,JDir,Restart, *)
(* Alpha and V are ignored. *)
(* If Trudge is called with non-negative LDirK and Restart=FALSE *)
(* then values of TimLim,TolF,TolR,Noise,NDir,X and V are *)
(* required. In this case KDir is reset to zero after first pass *)
(* *)
(* GeoOptParam : control parameter set *)
(* DirK : indicates the conjugate gradient *)
(* direction in which the optimization *)
(* will proceed (Def. = FALSE) *)
(* FALSE : indicates that this is a *)
(* non-restart run. *)
(* TRUE : corresponds to a restart run *)
(* Noise : Accuracy of function values. Variation smaller *)
(* than Noise are not considered to be significant *)
(* (Def. 0.0005) *)
(* TolF : accuracy required of the function (Def. 0.001) *)
(* TolR : accuracy required of conjugate directions *)
(* (Def. 0.05) *)
(* *)
(* For geometry optimization, the values which give better *)
(* results (closer to the ones obtained with gradient methods) *)
(* are: TolF=0.0001, TolR=0.001, Noise=0.00001 *)
(*----------------------------------------------------------------*)
END OptimLib1.