IntervalRep:
TYPE ~
RECORD [
-- information about interval, used in min search tree
depth: NAT ¬ 0, -- depth within tree
parent: Interval ¬ NIL, -- parent interval
prev, next: Interval ¬ NIL, -- prev and next in linked list of active leaves
l, r: Interval ¬ NIL, -- ptrs to left and right child intervals, IF any
lt, rt: REAL, -- t at left and right ends of interval
lval, rval: REAL, -- values at left and right ends of interval
estval: REAL -- estimated minimum value over this interval
];
MinimizeFunction:
PUBLIC PROC [
function: FunctionProc,
t0, t1: REAL,
logSlope: REAL,
alpha: REAL,
slopeFac: REAL,
certainFac: REAL,
logEps: REAL,
lineProc: LineProc ¬ NIL,
report: STREAM ¬ NIL]
RETURNS [REAL]
~ {
Find the value of t within the interval [t0..t1] that minimizes function[t].
Written by Paul Heckbert, August, 1988.
NewInterval:
PROC [parent: Interval, t0, t1, v0, v1:
REAL]
RETURNS [Interval] ~ {
n: Interval ¬ ObtainInterval[];
n.depth ¬ IF parent # NIL THEN parent.depth+1 ELSE 0;
n.parent ¬ parent;
n.lt ¬ t0;
n.rt ¬ t1;
n.lval ¬ v0;
n.rval ¬ v1;
RETURN[n];
};
PrevSibling:
PROC [n: Interval]
RETURNS [Interval] ~ {
IF n.parent = NIL THEN RETURN[NIL];
IF n.parent.r = n
THEN RETURN[n.parent.l]
ELSE {
uncle: Interval ¬ PrevSibling[n.parent];
RETURN[IF uncle # NIL THEN uncle.r ELSE NIL];
};
};
NextSibling:
PROC [n: Interval]
RETURNS [Interval] ~ {
IF n.parent = NIL THEN RETURN[NIL];
IF n.parent.l = n
THEN RETURN[n.parent.r]
ELSE {
uncle: Interval ¬ NextSibling[n.parent];
RETURN[IF uncle # NIL THEN uncle.l ELSE NIL];
};
};
EndSlope:
PROC [a, b: Interval, sign:
INT16]
RETURNS [
REAL] ~ {
Estimate the slope at the shared point between sibling intervals a and b.
Actually returns the delta value for an interval of width (a.rt-a.lt), not the slope.
a.rt and b.lt should be equal.
Adds or subtracts slack according to sign.
uncertainty: REAL ¬ maxdiff*RealFns.Power[alpha, IF a # NIL THEN a.depth ELSE b.depth];
IF a =
NIL
OR b =
NIL
THEN RETURN[sign*uncertainty]
ELSE RETURN[slopeFac*(b.rval-a.lval)/2.0+sign*certainFac*uncertainty];
};
EstimateMinVal:
PROC [n: Interval]
RETURNS [
REAL] ~ {
Estimate the minimum value of the function over interval n:
l: Interval ¬ PrevSibling[n];
r: Interval ¬ NextSibling[n];
m0: REAL ¬ EndSlope[l, n, -1];
m1: REAL ¬ EndSlope[n, r, 1];
estcount ¬ estcount+1;
RETURN[n.estval ¬ MIN[n.lval, n.rval, n.lval+m0, n.rval-m1]];
};
FindMinEstimate:
PROC []
RETURNS [minInterval: Interval] ~ {
Step through the active leaves to find the one with the minimum estimate:
minEst: REAL ¬ Real.LargestNumber;
FOR n: Interval ¬ firstActive, n.next
WHILE n #
NIL
DO
IF lineProc # NIL THEN lineProc[n.lt, n.estval, n.rt, n.estval];
searchcount ¬ searchcount+1;
IF n.estval < minEst THEN {minEst ¬ n.estval; minInterval ¬ n;};
ENDLOOP;
};
FindMinValue:
PROC []
RETURNS [min:
REAL] ~ {
Step through the leaves to find minimum value:
min ¬ Real.LargestNumber;
FOR n: Interval ¬ firstActive, n.next
WHILE n #
NIL
DO
IF MIN[n.lval, n.rval] < min THEN min ¬ MIN[n.lval, n.rval];
ENDLOOP;
};
Weed:
PROC [minValue:
REAL]
RETURNS [
NAT ¬ 0] ~ {
Step through the leaves eliminating those whose min estimates are greater than minValue:
FOR n: Interval ¬ firstActive, n.next
WHILE n #
NIL
DO
IF n.estval > minValue
THEN {
IF n.prev # NIL THEN n.prev.next ¬ n.next ELSE firstActive ¬ n.next;
IF n.next # NIL THEN n.next.prev ¬ n.prev;
n.prev ← n.next ← n.l ← n.r ← NIL; -- assist garbage collector?
};
ENDLOOP;
};
FreeIntervals:
PROC [root: Interval] ~ {
Inner:
PROC [i: Interval] ~ {
IF i.l # NIL THEN Inner[i.l];
IF i.r # NIL THEN Inner[i.r];
ReleaseInterval[i];
};
Inner[root];
};
maxdiff: REAL ¬ RealFns.Power[10., logSlope]*(t1-t0);
epsilon: REAL ¬ RealFns.Power[10., logEps]*(t1-t0); -- convert to absolute err tol.
root: Interval ¬ NewInterval[NIL, t0, t1, function[t0], function[t1]];
firstActive: Interval ¬ root;
maxDepth: NAT ~ 10;
estcount, searchcount, iter: NAT ¬ 0;
[] ¬ EstimateMinVal[root];
IF lineProc # NIL THEN lineProc[root.lt, root.lval, root.lt, root.lval];
IF lineProc # NIL THEN lineProc[root.rt, root.rval, root.rt, root.rval];
DO
bogus: NAT ¬ Weed[FindMinValue[]];
n: Interval ¬ FindMinEstimate[];
t: REAL ¬ 0.5*(n.lt+n.rt);
iter ¬ iter+1;
IF n.depth > maxDepth
OR
ABS[n.lt-n.rt] < epsilon
THEN {
IF report #
NIL
THEN
IO.PutFL[report, "%g estimates, %g searches, %g iterations, depth %g, t=%g\n",
LIST[IO.int[estcount], IO.int[searchcount], IO.int[iter], IO.int[n.depth], IO.real[t]]];
FreeIntervals[root];
RETURN[t];
}
ELSE {
val: REAL ¬ function[t];
l: Interval ¬ NewInterval[n, n.lt, t, n.lval, val];
r: Interval ¬ NewInterval[n, t, n.rt, val, n.rval];
ll: Interval ¬ PrevSibling[l];
rr: Interval ¬ NextSibling[r];
IF ll # NIL THEN [] ¬ EstimateMinVal[ll];
[] ¬ EstimateMinVal[l];
[] ¬ EstimateMinVal[r];
IF rr # NIL THEN [] ¬ EstimateMinVal[rr];
IF lineProc # NIL THEN lineProc[t, val, t, val];
n.l ¬ l;
n.r ¬ r;
l.prev ¬ n.prev;
r.prev ¬ l;
IF n.prev # NIL THEN n.prev.next ¬ l ELSE firstActive ¬ l;
IF n.next # NIL THEN n.next.prev ¬ r;
l.next ¬ r;
r.next ¬ n.next;
n.prev ¬ n.next ¬ NIL;
};
ENDLOOP;
};
MinimizeMonotonic:
PUBLIC
PROC [
function: FunctionProc,
t0: REAL ¬ 0.0, -- minimum abscissa
t1: REAL ¬ 360.0, -- maximum abscissa
eps: REAL] -- RETURN if dt < eps
since our reals are single precision, a suggested eps is about 1e-4 or 1e-5 times x2-x1
RETURNS [REAL]
~ {
Written by Paul Heckbert, 1/89.
Golden ratio method: split interval into two or three subintervals at each iteration and
recurse on a new subinterval which is typically 0.618 times as long as the previous.
The algorithm expects the function to have a single minimum within the interval, and be
monotonically decreasing to the left and increasing to the right of the minimum.
Even if there are several minima this algorithm will find one of them
(not guaranteed to be the global minimum, however).
Note: since we fit the function with a parabola when the recursion ends, epsilon should
not be set so small that we've lost precision at this point. If epsilon is going to be set
very small then we should not use min←parabola but simply return a point (xb or xc) from
the middle of the interval. If we use min←parabola then the x returned by minimize is
not guaranteed to lie in [x1..x2].
MinLeft:
PROC [xa, xb, xd, ya, yb, yd:
REAL]
RETURNS [
REAL] ~ {
MinLeft: find min in interval [xa..xd] given left interior sample
xc, yc: REAL;
invariant: xb = xa+phi2*(xd-xa)
IF xd-xa < eps THEN RETURN[xb];
IF xd-xa< eps THEN RETURN[QuadFit[xa, xb, xd, ya, yb, yd]];
IF ya < yb
THEN {
-- assume yb<yd, UP UP case: min is in [xa..xb]
xc ¬ xa+phi2*(xb-xa);
RETURN[MinLeft[xa, xc, xb, ya, function[xc], yb]];
};
xc ¬ xb+phi*(xd-xb); -- make a fourth sample
yc ¬ function[xc];
RETURN[
IF (yb < yc)
THEN MinRight[xa, xb, xc, ya, yb, yc] -- DOWN UP UP: min [xa..xc]
ELSE MinLeft[xb, xc, xd, yb, yc, yd]]; -- DOWN DOWN or DOWN DOWN UP: min[xb..xd]
};
MinRight:
PROC [xa, xc, xd, ya, yc, yd:
REAL]
RETURNS [
REAL] ~ {
MinRight: find min in interval [xa..xd] given right interior sample
xb, yb: REAL;
invariant: xc = xa+phi*(xd-xa)
IF xd-xa < eps THEN RETURN[xc];
IF xd-xa < eps THEN RETURN QuadFit[xa, xc, xd, ya, yc, yd];
IF yc > yd
THEN {
-- assume ya > yc, DOWN DOWN case: min is
[xc..xd]
xb ¬ xc+phi*(xd-xc);
RETURN[MinRight[xc, xb, xd, yc, function[xb], yd]];
};
xb ¬ xa+phi*(xc-xa); -- make a fourth sample
yb ¬ function[xb];
RETURN[
IF yb > yc
THEN MinLeft[xb, xc, xd, yb, yc, yd] -- DOWN DOWN UP case: min [xb..xd]
ELSE MinRight[xa, xb, xc, ya, yb, yc]]; -- UP UP or DOWN UP UP case: min [xa..xc]
};
MinParabola:
PROC [x1, x2, x3, y1, y2, y3:
REAL]
RETURNS [
REAL] ~ {
Find x coordinate of minimum of parabola passing through (x1,y1), (x2,y2), (x3,y3)
t1: REAL ¬ (x2-x1)*(y3-y2);
t2: REAL ¬ (x3-x2)*(y2-y1);
IF t1 = t2 THEN RETURN[x2];
RETURN[0.5*(t2*(x3+x2)-t1*(x2+x1))/(t2-t1)];
};
xc: REAL ¬ t0+phi*(t1-t0);
RETURN[MinRight[t0, xc, t1, function[t0], function[xc], function[t1]]];
};