
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]
~ {
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<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 B [xa..xc]
ELSE MinLeft[xb, xc, xd, yb, yc, yd]]; -- DOWN DOWN or DOWN DOWN UP: minB[xb..xd]
};
MinRight: PROC [xa, xc, xd, ya, yc, yd: REAL] RETURNS [REAL] ~ {
xb, yb: REAL;
IF xd-xa < eps THEN RETURN[xc];
IF yc > yd THEN {	-- assume ya > yc, DOWN DOWN case: min is B [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 B [xb..xd]
ELSE MinRight[xa, xb, xc, ya, yb, yc]]; -- UP UP or DOWN UP UP case: min B [xa..xc]
};
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.


   	R  
G3dFunctionImpl.mesa
Copyright ╙ 1988, 1992 by Xerox Corporation. All rights reserved.
Bloomenthal, July 14, 1992 3:32 pm PDT
Definitions
A General Univariate Function Minimizer
Find the value of t within the interval [t0..t1] that minimizes function[t].
Written by Paul Heckbert, August, 1988.
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.
Estimate the minimum value of the function over interval n:
Step through the active leaves to find the one with the minimum estimate:
Step through the leaves to find minimum value:
Step through the leaves eliminating those whose min estimates are greater than minValue:
n.prev _ n.next _ n.l _ n.r _ NIL;	-- assist garbage collector?
since our reals are single precision, a suggested eps is about 1e-4 or 1e-5 times x2-x1
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: find min in interval [xa..xd] given left interior sample
invariant: xb = xa+phi2*(xd-xa)
IF xd-xa< eps THEN RETURN[QuadFit[xa, xb, xd, ya, yb, yd]];
MinRight: find min in interval [xa..xd] given right interior sample
invariant: xc = xa+phi*(xd-xa)
IF xd-xa < eps THEN RETURN QuadFit[xa, xc, xd, ya, yc, yd];
Find x coordinate of minimum of parabola passing through (x1,y1), (x2,y2), (x3,y3)
Functions
The paper gives the function ar6/R6+br4/R4+cr2/R2+1, where a=-4/9, b=17/9, and c=-22/9.
We give it in terms of x, where x=r/R.
 Э╩;  Ц"cedarcode" styleХNewlineDelimiter
Щ ЩJЪЬ	╧eЬ5ЩAJЩ&JШ JЪ╧k	Ь Ш)JШ ЧЪ╨blЬЮЬЮШJЪЮЬШJЪЮЬШJШ ЧJЪЬЮШheadЪ╧lЩJЪЬЮЬШ/JЪЬ	ЮЬШ)JЪЮЬЮЬЮЬЮЬШЧЪа'Щ'JЪЬЮЬШЪ	ЬЮЬЮЬЮЬЮЬШFJШ ЧЪ╧nЬЮЬЮЬЮЬШ1ЪЮЬЮЬЮЬЮШ'JЪЬ$Ш$JЪЮЬЮЬЮЬЮЬШJЪЬЮЬШJЪЬ3ЮЬ╧cШQJЪЮЬШJЪЮЬШЧJЪЮЬЮЬШJШJШ ЧЪбЬЮЬЮЬШ/JЪЬ3ЮЬвШIЪЮЬЮЬЮЬЮШ'JЪЮЬЮЬЮЬЮЬШ'JШJЪЮЬШJЪЮЬШЧJШJШ ЧJЪЬ	ЮЬЮЬШ!ЪЬЮЬЮЬв6ШVJЪЬЮЬвШ*JЪЬЮЬвШ-JЪЬЮЬв0ШNJЪЬЮЬв1ШKJЪЬЮ	Ьв'Ш<JЪЬЮЬв,ШCJЪЬЮ	Ьв-ШAJШJШ ЧЪбЬЮЬШJЪЬШJЪЬЮЬШJЪЬ
ЮЬШJЪЬЮЬШJЪЬ
ЮЬШJЪЬЮЬШJЪЬЮЬШJЪЬЮЬШJЪЬЮЬЮЬШJЪЮЬЮЬШJЪЬШJЩLJЩ'ЪбЬЮЬ$ЮЬЮЬШQJШJЪ	Ь
ЮЬ
ЮЬЮЬЮЬШ5JШJШ
JШ
JШJШJЪЮЬШ
JШЧЪбЬЮЬЮЬШ6JЪ
ЮЬЮЬЮЬЮЬЮЬШ#ЪЮЬШJЪЮЬЮЬШЪЮЬШJЪЬ╧bЬШ(JЪЮЬЮЬ	ЮЬЮЬ	ЮЬЮЬШ-JШЧЧJШЧЪбЬЮЬЮЬШ6JЪ
ЮЬЮЬЮЬЮЬЮЬШ#ЪЮЬШJЪЮЬЮЬШЪЮЬШJЪЬгЬШ(JЪЮЬЮЬ	ЮЬЮЬ	ЮЬЮЬШ-JШЧЧJШЧЪ
бЬЮЬЮЬЮЬЮЬШ?JЩIJЩUJЩJЩ*JЪЬЮ╧sЬЮЬдЬдЮЬЮЬ	ЮЬ
ШWЪЮЬЮЬЮЬЮШJЪЮЬШJЪЮЬ;ШFЧJШЧЪбЬЮЬЮЬЮЬШ5JЩ;JШJШJЪЬЮЬШJЪЬЮЬШJШJЪЮЬЮЬ(Ш=JШЧЪбЬЮЬЮЬШ<JЩIJЪЬЮЬШ"ЪЮЬ#ЮЬЮЬЮШ6JЪЮЬЮЬЮЬ*Ш@JШJЪЮЬЮЬ'Ш@JЪЮЬШЧJШЧЪбЬЮЬЮЬЮЬШ-JЩ.JШЪЮЬ#ЮЬЮЬЮШ6JЪЮЬЮЬЮЬЮЬШ<JЪЮЬШЧJШЧЪ
бЬЮЬЮЬЮЬЮЬ	Ш1JЩXЪЮЬ#ЮЬЮЬЮШ6ЪЮЬЮЬШJЪЮЬ
ЮЬЮЬЮЬШDJЪЮЬ
ЮЬЮЬШ*JЪЬЮЬвЩ?JЪЬШЧJЪЮЬШЧJШЧЪбЬЮЬШ(ЪбЬЮЬШJЪЮЬЮЬЮЬгЬШJЪЮЬЮЬЮЬгЬШJШJШЧJЪЬШJШЧJЪЬ	ЮЬ(Ш5JЪЬ	ЮЬ'вШSJЪЬЮЬ&ШFJШJЪЬ
ЮЬШJЪЬЮЬШ%JШJЪЮЬЮЬЮЬ2ШHJЪЮЬЮЬЮЬ2ШHЪЮШJЪЬЮЬШ"JШ JЪЬЮЬШJШЪЮЬЮЬЮЬШ1ЪЮЬШЪЮЬ
ЮЬЮШЪЮЬLШNJЪЮЬЮЬЮЬЮЬЮЬЮЬШXЧЧJШJЪЮЬШ
JШЧЪЮЬШJЪЬЮЬШJШ3JШ3JШJШJЪЮЬЮЬЮЬШ)JШJШJЪЮЬЮЬШ)JЪЮЬЮЬЮЬШ0JШJШJШ JШJШJЪЮЬ
ЮЬЮЬЮЬШ:JЪЮЬ
ЮЬЮЬШ%JШJШJШ JЪЬЮЬШJШЧЧJЪЮЬШЧJШJШ ЧJЪЬЮЬШ(ЪЬЮЬвШ#JШ ЧЪбЬЮЬЮЬШ JЪЬШJЪЬЮЬвШ'JЪЬЮЬвШ(JЪЬЮЬв╨csвШ%JЪЬWЩWJЪЮЬЮЬШЪЬШJЩJЩXJЩTJЩWJЩPJЩEJЩ3JЩ JЩWJЩ[JЩXJЩVJЩ"JЩ ЧЪ
бЬЮЬЮЬЮЬЮЬШ?JЩAJЪЬЮЬШJЩJЪЮЬЮЬЮЬШJЪЮЬЮЬ"Щ;ЪЮЬ	ЮЬв/ШAJШJЪЮЬ,Ш2JШЧJЪЬвШ-JШЪЮЬЮЬ
ШJЪЮЬ"ве
в╨cmв	ШCJЪЮЬ#ве
вевевжвШQЧJШЧЪ
бЬЮЬЮЬЮЬЮЬШ@JЩCJЪЬЮЬШJЩJЪЮЬЮЬЮЬШJЪЮЬЮЬ!Щ;ЪЮЬ	ЮЬв*жв	ШFJШJЪЮЬ-Ш3JЪЬШЧJШ-JШЪЮЬЮЬвШJЪЮЬ#вевжв	ШKJЪ	ЮЬ$веве
вжв	ШSЧJШЧЪ
бЬЮЬЮЬЮЬЮЬШCJЪЬRЩRJЪЬЮЬШJЪЬЮЬШJЪЮЬ	ЮЬЮЬШJЪЮЬ&Ш,JЪЬШЧJЪЬЮЬШJЪЮЬAШGJШЧЧЪа	Щ	Ъ╨bnЬЮЬЮЬЮЬЮЬЮЬШ/JЪЬЮЬШJЪЬЮЬШJЪЬ	ЮЬ1Ш>JЪЬ	ЮЬ:ШGJЪЬЮЬ
ШJЪЮЬ&Ш,JШJШ ЧЪ
зЬЮЬЮЬЮЬЮЬШ/JЪЬЮЬЮЬ
ШJЪЮЬШJШJШ ЧJЪзЬЮЬЮЬЮЬЮЬЮЬЮЬШSJШ ЪзЬЮЬЮЬЮЬЮЬЮЬШ3JЪЮЬ
ЮЬ	ЮЬЮЬШ)JЪЮЬ>ШDJШЧJШ Ъ
бЬЮЬЮЬЮЬЮЬШ0JЪЬ╧uЬиЬиЬиЬиЬиЬ&ЩWJЩ&ЪЮЬЮШJЪЬ	ЮЬШJЪЬ	ЮЬШЪЮЬШJЪЬЮЬШJЪЬЮЬ	ШJЪЮЬ4Ш:JШЧЧJШЧJШ ЧJЪЮЬШJШ JШ Ч ЕЧ    x  4  