<<>> <> <> <> DIRECTORY G3dFunction, IO, Real, RealFns; G3dFunctionImpl: CEDAR MONITOR IMPORTS IO, RealFns EXPORTS G3dFunction ~ BEGIN <> FunctionProc: TYPE ~ G3dFunction.FunctionProc; LineProc: TYPE ~ G3dFunction.LineProc; STREAM: TYPE ~ IO.STREAM; <> nScratchIntervals: NAT ~ 100; scratchIntervals: ARRAY [0..nScratchIntervals) OF Interval ¬ ALL[NIL]; ObtainInterval: ENTRY PROC RETURNS [Interval] ~ { FOR i: NAT IN [0..nScratchIntervals) DO ref: Interval ~ scratchIntervals[i]; IF ref = NIL THEN LOOP; scratchIntervals[i] ¬ NIL; ref.parent ¬ ref.prev ¬ ref.next ¬ ref.l ¬ ref.r ¬ NIL; -- shouldn't be necessary RETURN[ref]; ENDLOOP; RETURN[NEW[IntervalRep]]; }; ReleaseInterval: ENTRY PROC [ref: Interval] ~ { ref.parent ¬ ref.prev ¬ ref.next ¬ ref.l ¬ ref.r ¬ NIL; -- very necessary FOR i: NAT IN [0..nScratchIntervals) DO IF scratchIntervals[i] # NIL THEN LOOP; scratchIntervals[i] ¬ ref; EXIT; ENDLOOP; }; Interval: TYPE ~ REF IntervalRep; 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] ~ { <> <> 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] ~ { <> <> <> <> 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] ~ { <> 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] ~ { <> 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] ~ { <> 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] ~ { <> 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; <> }; 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; }; phi: REAL ¬ 0.5*(RealFns.SqRt[5.0]-1.0); phi2: REAL ¬ 1.0-phi; -- = phi*phi MinimizeMonotonic: PUBLIC PROC [ function: FunctionProc, t0: REAL ¬ 0.0, -- minimum abscissa t1: REAL ¬ 360.0, -- maximum abscissa eps: REAL] -- RETURN if dt < eps <> RETURNS [REAL] ~ { <> <> <> <> <> <> <<(not guaranteed to be the global minimum, however).>> <<>> <> <> <> <> <> <<>> MinLeft: PROC [xa, xb, xd, ya, yb, yd: REAL] RETURNS [REAL] ~ { <> xc, yc: REAL; <> IF xd-xa < eps THEN RETURN[xb]; <> IF ya < yb THEN { -- assume yb> xb, yb: REAL; <> IF xd-xa < eps THEN RETURN[xc]; <> IF yc > yd THEN { -- assume ya > yc, DOWN DOWN case: min is 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 ELSE MinRight[xa, xb, xc, ya, yb, yc]]; -- UP UP or DOWN UP UP case: min }; MinParabola: PROC [x1, x2, x3, y1, y2, y3: REAL] RETURNS [REAL] ~ { <> 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]]]; }; <> Gauss: PUBLIC PROC [x: REAL] RETURNS [REAL] ~ { mean: REAL ~ 0.5; standardDeviation: REAL ~ 0.16; factor1: REAL ~ 1.0/(2.0*standardDeviation*standardDeviation); factor2: REAL ~ 1.0/(standardDeviation*RealFns.SqRt[2.0*3.1415926535]); xx: REAL ~ x-mean; RETURN[factor2*RealFns.Exp[-xx*xx*factor1]]; }; Spike: PUBLIC PROC [x: REAL] RETURNS [REAL] ~ { xx: REAL ~ 1.0-ABS[1.0-x-x]; RETURN[xx*xx]; }; Poisson: PUBLIC PROC [x, a: REAL] RETURNS [REAL] ~ {RETURN[a*x*RealFns.Exp[-a*x]]}; Perlin: PUBLIC PROC [x, a: REAL] RETURNS [REAL] ~ { IF x <= 0.0 OR a <= 0.0 THEN RETURN[0.0]; RETURN[RealFns.Power[0.5, RealFns.Log[0.5, x]*RealFns.Log[0.5, a]]]; }; Wyvill: PUBLIC PROC [x: REAL] RETURNS [REAL] ~ { <> <> SELECT x FROM < 0.0 => RETURN[1.0]; > 1.0 => RETURN[0.0]; ENDCASE => { t1: REAL ~ x*x; x4: REAL ~ t1*t1; RETURN[(-4.0/9.0)*t1*x4+(17.0/9.0)*x4+(-22.0/9.0)*t1+1.0]; }; }; END.