{TOptim1D class
author: Nikolai V. Shokhirev http://www.shokhirev.com/nikolai.html
created: 2003.06.06
last modified(: 2004.01.21
©Nikolai V. Shokhirev, 2001-2003 }
unit uOptimization;
interface
uses
uMatTypes, uDynArrays;
type
TOptim1D = class
private
x1, x2, x3, x4, x0: TFloat;
f1, f2, f3, f4, f0: TFloat;
d1f, d2f, d3f: TFloat;
fXmin, fXmax: TFloat;
fLimitLow, fLimitHigh: TFloat;
fReverse, fBounded: boolean;
fOutOfLimits: boolean;
fIterMax: TInt;
function GetXmax: TFloat;
function GetXmin: TFloat;
procedure SetXmax(const Value: TFloat);
procedure SetXmin(const Value: TFloat);
function GetIterMax: TInt;
procedure SetIterMax(const Value: TInt);
function GetLimitHigh: TFloat;
function GetLimitLow: TFloat;
procedure SetLimitHigh(const Value: TFloat);
procedure SetLimitLow(const Value: TFloat);
function GetBounded: boolean;
function GetOutOfLimits: boolean;
function GetOptimalX: TFloat;
public
constructor Create;
procedure Bracketing(f: TFloatFunction1D; x, step: TFloat);
procedure GoldenSection(f: TFloatFunction1D; eps: TFloat);
procedure BernetMin(f: TFloatFunction1D; eps: TFloat);
property Bounded: boolean read GetBounded;
property OutOfLimits: boolean read GetOutOfLimits;
property Xmin: TFloat read GetXmin write SetXmin;
property Xmax: TFloat read GetXmax write SetXmax;
property LimitLow: TFloat read GetLimitLow write SetLimitLow;
property LimitHigh: TFloat read GetLimitHigh write SetLimitHigh;
property IterMax: TInt read GetIterMax write SetIterMax;
property OptimalX: TFloat read GetOptimalX;
end;
var
GoldenRatio,
GoldenFactor: TFloat;
implementation
uses
math;
{ TOptim }
constructor TOptim1D.Create;
begin
fLimitLow := NegMaxFloat;
fLimitHigh := MaxFloat;
fIterMax := 999;
end;
{ initial bracketing: f(Xmin) > f(x0) < f(Xmax), Xmin < x0 < Xmax
x - search start, step - initial step
LimitLow < LimitHigh - search limits }
procedure TOptim1D.Bracketing(f: TFloatFunction1D; x, step: TFloat);
var
h: TFloat;
iter: TInt;
procedure CheckLimits;
begin
if (x3 < LimitLow) then
begin
fOutOfLimits := true;
Xmin := LimitLow;
Xmax := x1;
end else
if (x3 > LimitHigh) then
begin
fOutOfLimits := true;
Xmin := x1;
Xmax := LimitHigh;
end;
end;
begin
f3 := 0.0;
fReverse := false;
fBounded := false;
fOutOfLimits := false;
if (x < LimitLow) or (x > LimitHigh) then
begin { wrong start point }
fOutOfLimits := true;
Xmin := LimitLow;
Xmax := LimitHigh;
exit;
end;
h := step;
x1 := x;
x2 := x;
f2 := f(x2);
x3 := x2 + step;
CheckLimits;
if OutOfLimits then { wrong initial step }
exit;
f3 := f(x3);
if f3 > f2 then
begin { wrong initial direction, changing to the opposite }
fReverse := true;
Swap(x3, x2);
Swap(f3, f2);
h := -h;
end;
iter := 0;
{ x1 x2 x3
+------+---------------+ }
repeat { step double and search }
inc(iter);
h := h+h;
x1 := x2;
f1 := f2;
x2 := x3;
f2 := f3;
x3 := x2+h;
CheckLimits;
if OutOfLimits then
exit
else
f3 := f(x3);
fBounded := (f3 > f2);
until fBounded or (iter > IterMax);
{ x1 x2 x3 x4
+------+------+------+ }
if fBounded then
begin { half-step back }
x4 := x3;
f4 := f3;
h := h/2.0;
x3 := x3 + h;
f3 := f(x3);
if fReverse then
begin
Swap(x1, x4);
Swap(f1, f4);
Swap(x2, x3);
Swap(f2, f3);
h := -h;
end;
{ estimation of Xmin, x0, f0, Xmax and the derivatives }
d3f := ((f4-f1)-3.0*(f3-f2))/(h*sqr(h));
if (abs(f3-f2) < SqrtMachEps*abs(f3+f2)) then
begin
Xmin := x1;
Xmax := x4;
x0 := (x3+x2)/2.0;
f0 := (9.0*(f3+f2)-(f4+f1))/16.0;
d1f := (27.0*(f3-f2)-(f4-f1))/(24.0*h);
d2f := ((f4+f1)-(f3+f2))/(2.0*sqr(h));
end else
if f2 < f3 then
begin
Xmin := x1;
Xmax := x3;
x0 := x2;
f0 := f2;
d1f := (f3-f1)/(2.0*h);
d2f := (f3-2.0*f2+f1)/sqr(h);
end else
begin
Xmin := x2;
Xmax := x4;
x0 := x3;
f0 := f3;
d1f := (f4-f2)/(2.0*h);
d2f := (f4-2.0*f3+f2)/sqr(h);
end;
end;
end;
{ Golden Section Minimization
Search termination conditions:
|dx| < Sqrt(MachEps*2*f0/(d2f/dx2)) or |dx| < eps
x1 x2 x3 x4
+--------+------+--------+ }
procedure TOptim1D.GoldenSection(f: TFloatFunction1D; eps: TFloat);
var
iter: TInt;
dxIsSmall: boolean;
begin
X4 := Xmax;
X1 := Xmin;
f1 := f(x1);
f4 := f(x4);
x2 := x1+GoldenFactor*(x4-x1);
f2 := f(x2);
x3 := x4-GoldenFactor*(x4-x1);
f3 := f(x3);
fBounded := ((f1>f2) and (f2<f3)) or ((f2>f3) and (f3<f4)) or
((abs(f3-f2) < SqrtMachEps*abs(f3+f2)) and (f1>f2) and (f3<f4));
if not fBounded then
exit; // Bracketing is needed
iter := 0;
repeat
inc(iter);
if f2 < f3 then
begin
x4 := x3;
f4 := f3;
x3 := x2;
f3 := f2;
x2 := x1+GoldenFactor*(x4-x1);
f2 := f(x2);
// already have x1, f1
end else
begin
x1 := x2;
f1 := f2;
x2 := x3;
f2 := f3;
x3 := x4-GoldenFactor*(x4-x1);
f3 := f(x3);
// already have x4, f4
end;
// H = x4-x1; h = x3-x2 = (sqrt5-2)*H
// f0 = ((sqrt5+2)*(f3+f2)-(sqrt5-2)*(f4+f1))/8 ~(f2+f3)/2
// d2f := (sqrt5+2)*(f4-f3-f2+f1)/(H*H)=(sqrt5-2)*(f4-f3-f2+f1)/(h*h)
// 2*f0/d2f ~ (f3+f2)*h*h*(sqrt5+2)/(f4-f3-f2+f1)
// (sqrt5+2) ~ 4.2360679774997896964091736687313
// done: h/2 < sqrt(MachEps*(f3+f2)*h*h*(sqrt5+2)/(f4-f3-f2+f1))
x0 := (x3+x2)/2.0;
dxIsSmall := (sqrt(abs(f4-f3-f2+f1)) < 4.0*SqrtMachEps*sqrt(abs(f3+f2)+MachEps))
or (abs(x3-x2)< abs(eps));
until dxIsSmall or (iter > Itermax);
f0 := ((sqrt(5.0)+2.0)*(f3+f2)-(sqrt(5.0)-2.0)*(f4+f1))/8.0;
// for bracketing
Xmin := x2;
Xmax := x3;
end;
procedure TOptim1D.BernetMin(f: TFloatFunction1D; eps: TFloat);
var
dx: TFloat;
dxIsSmall, MinExists, MinFound: boolean;
iter: TInt;
procedure OptimalAbscissa;
begin
x0 := x2-(sqr(x2-x1)*(f2-f3)-sqr(x2-x3)*(f2-f1))/((x2-x1)*(f2-f3)-(x2-x3)*(f2-f1))/2.0;
if (abs(x2-x1) > SqrtMinFloat) and (abs(x3-x2) > SqrtMinFloat) then
d2f := 2.0*(-f1/(x2-x1)/(x1-x3)+ f2/(x2-x3)/(x2-x1)+f3/(x1-x3)/(x2-x3) );
end;
begin
x1 := Xmin;
x3 := Xmax;
x2 := (x3+x1)/2.0;
f1 := f(x1);
f2 := f(x2);
f3 := f(x3);
fBounded := ((f1>f2) and (f2<f3));
if not Bounded then
exit; // Bracketing is needed
iter := 0;
repeat
inc(iter);
OptimalAbscissa;
f0 := f(x0);
dxIsSmall := abs(x0-x2) < eps;
if (x0 < x2) then
begin
MinExists := (f1>f0) and (f0<f2);
if MinExists then
begin // min is in [x1,x2], use x0 as x2
x3 := x2;
f3 := f2;
x2 := x0;
f2 := f0;
end else
begin // min is NOT in [x1,x2], reduce the interval (bracketing)
x1 := x0;
f1 := f0;
x2 := (x3+x1)/2.0;
f2 := f(x2);
end;
end else
begin
MinExists := (f2>f0) and (f0<f3);
if MinExists then
begin // min is in [x2,x3], use x0 as x2
x1 := x2;
f1 := f2;
x2 := x0;
f2 := f0;
end else
begin // min is NOT in [x2,x3], reduce the interval (bracketing)
x3 := x0;
f3 := f0;
x2 := (x3+x1)/2.0;
f2 := f(x2);
end;
end;
dx := 4.0*SqrtMachEps*sqrt(abs(f2/d2f));
MinFound := (f(x0+dx)>f0) and (f(x0-dx)>f0);
until MinFound or dxIsSmall or (iter > Itermax);
f0 := f2;
end;
function TOptim1D.GetXmax: TFloat;
begin
result := fXmax;
end;
function TOptim1D.GetXmin: TFloat;
begin
result := fXmin;
end;
procedure TOptim1D.SetXmax(const Value: TFloat);
begin
fXmax := Value;
end;
procedure TOptim1D.SetXmin(const Value: TFloat);
begin
fXmin := Value;
end;
function TOptim1D.GetBounded: boolean;
begin
result := fBounded;
end;
function TOptim1D.GetIterMax: TInt;
begin
result := fIterMax;
end;
procedure TOptim1D.SetIterMax(const Value: TInt);
begin
fIterMax := Value;
end;
function TOptim1D.GetLimitHigh: TFloat;
begin
result := fLimitHigh;
end;
function TOptim1D.GetLimitLow: TFloat;
begin
result := fLimitLow;
end;
function TOptim1D.GetOutOfLimits: boolean;
begin
result := fOutOfLimits;
end;
procedure TOptim1D.SetLimitHigh(const Value: TFloat);
begin
fLimitHigh := Value;
end;
procedure TOptim1D.SetLimitLow(const Value: TFloat);
begin
fLimitLow := Value;
end;
function TOptim1D.GetOptimalX: TFloat;
begin
result := x0;
end;
initialization
{ GoldenRatio = 1.618033988749895 }
GoldenRatio := (sqrt(5.0)+1.0)/2.0;
{ gold = 0.381966011250105 }
GoldenFactor := 0.5*(3.0-sqrt(5.0));
end.
Generated by Lore's Source to HTML Converter(http://www.newty.de/lsc/index.html)