program Minimo (input*, output); const n = 4; beeping = true; type index = 1..n; yesno = (yes, no); vector = record c: array [index] of real end; pvec = ^ vector; arr = record a: array [index] of pvec end; parr = ^ arr; var unitrac: yesno; (* YES if unimin should be traced *) rsol: pvec; (* Correct solution for MultiMin testing *) procedure Beep (c: char); begin if beeping then write (output, c, ' '); end; procedure UniMin ( function F (pvec): real; rv: pvec; d, eps: real; rxx: pvec; var y: real; var nf: integer); (* Minimizes the function F along the direction rv^. (* (* The function F must be eventually concave (i.e. must be plus infinity (* for infinite x). The parameter d is the (* estimated distance from the initial guess (rxx^) to the minimum. (* The output is returned in rxx^ too; if xin and (* xout are the input and output values of rxx^, then: (* (* 1. F(xout) <= F(xin) (* 2. xout = xin + s*rv^ for some real s (* 3. There exist positive reals e1 and e2 between eps/2 and eps (* such that F(xout-e1*rv) >= F(xout) <= F(xout+e2*rv). (* (* The variable nf gives, on output, the number of times F was evaluated. (* *) label 1,2; var rxf, rxb, rt, rxm, rx: pvec; yf, yb, ym, t, sf, sb, sm, dir, tau, mu, sh, sl: real; i, imax: index; procedure TurnAround; begin rt := rxf; rxf := rxb; rxb := rt; t := yf; yf := yb; yb := t; t := sf; sf := -sb; sb := -t; dir := -dir; Beep ('t') end; procedure Trace (rx: pvec; sb, sf, y: real); var i: index; unc, asp: real; begin if unitrac = yes then begin writeln (output); write (output, 'x', nf:3, ': '); write (output, rx^.c [imax] :16:8); unc := (ln (sf - sb) - ln (eps)) / ln (2.0) - 1; asp := (ln (sf) - ln (-sb)) / ln (2.0); write (output, ' (', unc : 6:2, asp : 6:2, ')'); writeln (output, ' F = ', y : 18) end end; begin (* UniMin *) (* Phase 1 - finds an interval containing a local minimum (* (characterized by three consecutive points xb, x and xf (* such that F(xb) >= F(x) <= F(xf)). *) 1: writeln (output); new (rxb); new (rxf); rx := rxx; d := max (d, eps); dir := 1; sb := -d; sf := d; imax := 1; for i := 1 to n do begin rxf^.c [i] := rx^.c [i] + d * rv^.c [i]; rxb^.c [i] := rx^.c [i] - d * rv^.c [i]; if abs(rv^.c [i]) > abs(rv^.c [imax]) then imax := i end; nf := 0; Trace (rx, sb, sf, y); yb := F(rxb); nf := nf + 1; Trace (rxb, sb, sf, yb); yf := F(rxf); nf := nf + 1; Trace (rxf, sb, sf, yf); if (yb >= y) and (yf >= y) then goto 2; if yb < yf then TurnAround; while true do begin tau := 2.0*((yb - y) * sf - (yf - y) * sb); mu := (yb - y) * sf * sf - (yf - y) * sb * sb; if tau > 0.0 then if mu > 5.0*sf*tau then sm := 5.0*sf else sm := max (eps, mu / tau + sf) else sm := 5.0*sf; if sm < sf then begin Beep ('f'); sm := min (max (sm, eps), max (sf/2, sf-eps)); rt := rxb; for i := 1 to n do rt^.c [i] := rx^.c [i] + sm * dir * rv^.c [i]; rxb := rx; yb := y; sb := - sm; rx := rt; y := F(rx); sf := sf - sm; nf := nf + 1; Trace (rx, sb, sf, y) end else begin Beep ('c'); sm := max (sm, sf + eps); rt := rxb; for i := 1 to n do rt^.c [i] := rx^.c [i] + sm * dir * rv^.c [i]; rxb := rx; yb := y; sb := -sf; rx := rxf; y := yf; rxf := rt; yf := F(rxf); sf := sm - sf; nf := nf + 1; Trace (rxf, sb, sf, yf) end; if yf >= y then goto 2 end; (* Phase 2 - reduces the interval until less than eps. *) 2: new (rxm); Beep ('2'); while ((sf > eps) or (sb < -eps)) and (nf < 100) do begin if sf > -4.0*sb then begin Beep ('s'); sm := -2.0 * sb end else if sb < -4.0 * sf then begin Beep ('s'); sm := -2.0 * sf end else begin tau := 2.0 * ((yb - y)* sf - (yf - y)*sb); mu := (yb - y) * sf * sf - (yf - y) * sb * sb; (* if unitrac = yes then *) (* begin *) (* writeln(output); *) (* writeln(output, sb :12, sf :12, tau :20, mu:12) *) (* end; *) if (tau <= 0) then sm := 0 (* with luck, will stop in 2 iter. *) else begin if mu > sf * tau then sm := sf else if mu < sb * tau then sm := sb else sm := mu / tau; end end; if abs (sm) < eps then begin Beep ('j'); if sf > -sb then sm := eps else sm := -eps end; if sm < 0 then begin TurnAround; sm := -sm end; sm := min (sm, max (sf-eps, sf/2)); for i := 1 to n do rxm^.c [i] := rx^.c [i] + sm * dir * rv^.c [i]; ym := F(rxm); nf := nf + 1; if (ym > y) or ((ym = y) and (sf > -sb)) then begin Beep ('>'); rt := rxf; rxf := rxm; rxm := rt; yf := ym; sf := sm; Trace (rxf, sb, sf, yf) end else begin Beep ('<'); rt := rxb; rxb := rx; rx := rxm; rxm := rt; yb := y; y := ym; sb := -sm; sf := sf - sm; Trace (rx, sb, sf, y) end end; (* return the best x in rxx^ *) for i := 1 to n do rxx^.c[i] := rx^.c[i] end; (* UniMin *) procedure UniTest; var rx, rv: pvec; y, eps: real; i, nf: index; function AFunc (rx: pvec): real; var t, f: real; i: index; begin t := rx^.c [1]; f := 1 - 1/(1 + sqr (5 * t)); AFunc := f; end; function QFunc (rx: pvec): real; var t, f: real; i: index; begin t := rx^.c [1]; f := sqr (5 * t) + 1.0; QFunc := f; end; function WFunc (rx: pvec): real; var t, f: real; i: index; begin t := rx^.c [1]; f := sqr(sqr (5 * t)) + 1.0; WFunc := f; end; begin (* UniTest *) new (rx); new (rv); for i := 1 to n do begin rx^.c [i] := 0.418 * i; rv^.c [i] := i end; eps := 0.000001; writeln (output); write (output, 'Epsilon = ', eps : 16:8); y := WFunc (rx); UniMin (WFunc, rv, 1.0, eps, rx, y, nf); writeln (output); write (output, 'Minimum at '); for i := 1 to n do write (output, rx^.c [i] :16:8); write (output, ' value = ', y : 20:12) end; procedure MultiMin ( function F (pvec): real; (* function to be minimized *) rx: pvec; (* input guess and final minimum *) d: real; (* estimated guess-sol. distance *) eps: real; (* minimum step and tolerance *) var y: real; (* output function value at minimum *) var nf, nmin: integer); (* number of calls to F and UinMin *) (* (* Minimizes a black-box function F of n arguments, using O(n) storage. (* (* The function F should be ultimately concave (i.e., F(x) is +infinity (* for infinite x). The procedure returns when the main loop has been (* unable to advance by more than eps for the last 2n iterations. (* It is guaranteed that F(xout) <= F(xin). *) var rxp, dir, theta: pvec; i, k, kmax: index; nit, nok, nft, kk, nwarm: integer; dd, ndd, di, dimax, ddang, lambd: real; obts: real; onf: integer; (* for Trace *) procedure TracePoint (nx: integer; rx: pvec; y: real; var bts: real); var j: index; ej, ee: real; begin ee := 0.0; writeln (output); write(output, '[MULTI] X', nx:3, ': '); for j := 1 to n do begin write (output, rx^.c[j] :16:8); ej := rsol^.c[j] - rx^.c[j]; ee := ee + ej*ej; end; ee := sqrt (ee); if ee > eps then bts := n* (ln (ee) - ln(eps))/ln(2.0) else bts := 0.0; writeln(output); write(output, ' F = ', y:20:12, ' DIST = ', ee: 20:12, ' (', bts:7:2, ')'); end; (* TracePoint *) procedure IniTrace (rx: pvec; y: real); begin TracePoint (0, rx, y, obts); onf := 0; end; (* IniTrace *) procedure TraceDir (rx, dir: pvec; ddang: real); var j: index; ang, ej, ee: real; begin ang := 0.0; ee:= 0.0; for j := 1 to n do begin ej := rsol^.c[j] - rx^.c[j]; ee := ee + ej*ej; ang := ang + dir^.c[j] * ej; end; ee := sqrt (ee); if ee > eps then ang := ang/ee else ang := 0.0; writeln (output); write(output, ' DIR: '); for j := 1 to n do begin write (output, dir^.c[j] :16:8) end; write (output, ' AIM = ', ang:8:4, ' STIFF = ', ddang:8:4) end; (* TraceDir *) procedure TraceNewPoint (nx: integer; rx: pvec; y: real); var bts: real; begin TracePoint (nx, rx, y, bts); write (output, ' DNF = ', nf-onf:5, ' EFFICIENCY = ', (obts-bts)/(nf-onf):7:3); onf := nf; obts := bts end; (* TraceNewPoint *) begin (* MultiMin *) new(rxp); new(dir); new(theta); for i := 1 to n do begin dir^.c[i] := 0.0; theta^.c[i] := 0.0 end; nit := 1; nok := 0; nf := 0; nmin := 0; y := F(rx); nf := nf + 1; IniTrace (rx, y); dir^.c[1] := 1.0; kmax := 1; dd := d; ddang := 0.0; while nok < 2 do begin TraceDir (rx, dir, ddang); (* warming run *) lambd := max (eps, dd/5.0); nwarm := n-1; if nwarm >= n then begin UniMin (F, dir, lambd, eps, rx, y, nft); nf := nf + nft; nmin := nmin + 1; end; kk := max(2, n-nwarm+1); if kmax >= kk then kk := kk -1; for k := kk to n do begin if k = kmax then begin end else begin theta^.c[k] := dir^.c[kmax]; theta^.c[kmax] := -dir^.c[k]; UniMin (F, theta, lambd, eps, rx, y, nft); nf := nf + nft; nmin := nmin + 1; theta^.c[k] := 0.0; theta^.c[kmax] := 0.0; end end; (* save current position *) for i := 1 to n do begin rxp^.c[i] := rx^.c[i] end; (* model run *) UniMin (F, dir, lambd, eps, rx, y, nft); nf := nf + nft; nmin := nmin + 1; for k := 1 to n do begin if k = kmax then begin end else begin theta^.c[k] := dir^.c[kmax]; theta^.c[kmax] := -dir^.c[k]; UniMin (F, theta, lambd, eps, rx, y, nft); nf := nf + nft; nmin := nmin + 1; theta^.c[k] := 0.0; theta^.c[kmax] := 0.0; end end; TraceNewPoint (nit, rx, y); (* recompute extrapolating direction dir *) ndd := 0.0; ddang := 0.0; dimax := 0.0; for i := 1 to n do begin di := 1.0 * (rx^.c[i] - rxp^.c[i]) + 0.0 * dd * dir^.c[i]; if abs (di) > dimax then begin dimax := abs(di); kmax := i end; ddang := ddang + di * dir^.c[i]; ndd := ndd + di * di; dir^.c[i] := di; end; dd := sqrt(ndd); if dd > eps then begin for i := 1 to n do dir^.c[i] := dir^.c[i] / dd; ddang := ddang / dd; nok := 0 end else begin for i := 1 to n do dir^.c[i] := 0.0; dd := eps; dir^.c[kmax] := 1.0; ddang := 0.0; nok := nok + 1 end; nit := nit + 1 end end; (* MultiMin *) procedure MultiTest (var nf, nmin: integer); var rx: pvec; raxis, rdata: pvec; (* internal parameters for test functions *) eps, y: real; i: index; function NFunc (rx: pvec): real; (* a narrow, quadratic valley *) var f, sc, ppi: real; i: index; begin sc := 0.0; for i := 1 to n do sc := sc + (rx^.c[i] - rsol^.c[i])*raxis^.c[i]; f := 0.0; for i := 1 to n do begin ppi := rx^.c[i] - rsol^.c[i] - sc*raxis^.c[i]; f := f + ppi * ppi end; NFunc := (10.0*f + sc*sc + 1.0) end; (* NFunc *) procedure NFuncInit (rx, rsol: pvec); var i: integer; sz: real; begin new(raxis); writeln(output); writeln(output); write(output, '*** AXIS: '); sz := 0.0; for i := 1 to n do begin rsol^.c[i] := 0.0; rx^.c[i] := 4.18; raxis^.c[i] := 1.0 + random (1.0); sz := sz + sqr(raxis^.c[i]); end; sz := sqrt(sz); for i := 1 to n do begin raxis^.c[i] := raxis^.c[i]/sz; write (output, raxis^.c[i]:16:8) end; writeln(output) end; function VFunc (rx: pvec): real; (* a non-linear least squares fit *) var i, k: integer; w, f, t: real; begin f := 0.0; w := rx^.c[n] * 6.2831852 / (n-1); for i := 1 to n do begin t := 0.0; for k := 1 to n-1 do t := t + rx^.c[k] * sin (w * k * i); f := f + sqr (rdata^.c[i] - t) end; VFunc := f end; (* VFunc *) procedure VFuncInit (rx, rsol: pvec); var i, k: integer; w, t: real; begin new(rdata); for i := 1 to n-1 do rsol^.c[i] := n-i;; rsol^.c[n] := 4.615; w := rsol^.c[n] * 6.2831852 / (n-1); for i := 1 to n do begin t := 0.0; for k := 1 to n-1 do t := t + rsol^.c[k] * sin (w * k * i); rx^.c[i] := rsol^.c[i] + 0.15; rdata^.c[i] := t end end; begin (* MultiTest *) new(rx); new(rsol); NFuncInit (rx, rsol); eps := 0.000001; MultiMin (NFunc, rx, 1.0, eps, y, nf, nmin); writeln (output); write (output, '[MULTI] XFIN :'); for i := 1 to n do write (output, rx^.c[i]:16:8); write (output, ' F = ', y:16:8, ' NMIN = ', nmin:5); writeln (output) end; (* MultiTest *) procedure MultiMain; (* exercises MultiMin thru MultiTest and collects *) (* statistics *) var nftot, nmintot: real; m, nfm, nminm: integer; begin nftot := 0.0; nmintot := 0.0; for m := 1 to 5 do begin MultiTest (nfm, nminm); nftot := nftot + nfm; nmintot := nmintot + nminm; end; writeln(output); writeln (output, nftot/5:10:2, ' evaluations and ', nmintot/5:6:2, ' unidimensional minimizations.') end; (* MultiMain *) begin (* MAIN *) writeln (output, 'full trace (YES or NO)?'); break (output); read (input, unitrac); MultiMain; end.