MODULE BrentUniMin;

IMPORT Math;
IMPORT Wr, Thread, Fmt;
FROM Stdio IMPORT stderr;

PROCEDURE Sqrt(x: REAL): REAL =
  BEGIN
    RETURN FLOAT(Math.sqrt(FLOAT(x, LONGREAL)))
  END Sqrt;

PROCEDURE Minimize(
    ax, bx: REAL;
    f: Function;
    tol: REAL;
    good: GoodProc := NIL;     (* Acceptance criterion. *)
    verbose: BOOLEAN := FALSE; (* TRUE to print diagnostics. *)
  ): REAL =

  CONST 
    Eps = 1.0e-3;      (* Square root of relative machine precision. *) 
  VAR 
    IhpSq: REAL := 0.5*(3.0 - Sqrt(5.0)); (* squared inverse of the golden ratio *)

  VAR 
    a, b: REAL;
    u, v, w, fu, fv, fw: REAL;
    x, fx: REAL;
    d, e, xm, p, q, r, tol1, t2, tol3: REAL;
    fxOld: REAL;

  PROCEDURE Message (msg: TEXT) =
    <* FATAL Wr.Failure, Thread.Alerted *>
    BEGIN
      Wr.PutText(stderr, msg);
      Wr.PutText(stderr, "\n");
    END Message;

    PROCEDURE PrintStatus () =
      <* FATAL Wr.Failure, Thread.Alerted *>
      BEGIN 
        Wr.PutText(stderr, "\n");
        PrintAbscissa("a", a);
        PrintPoint("v", v, fv);
        PrintPoint("x", x, fx);
        PrintPoint("w", w, fw);
        PrintAbscissa("b", b);
        Wr.PutText(stderr, "\n");
      END PrintStatus;

    PROCEDURE PrintAbscissa (msg: TEXT; u: REAL) =
      <* FATAL Wr.Failure, Thread.Alerted *>
      BEGIN 
        Wr.PutText(stderr, "    ");
        Wr.PutText(stderr, msg);
        Wr.PutText(stderr, " = ");
        Wr.PutText(stderr, Fmt.Pad(Fmt.Real(u), 12));
        Wr.PutText(stderr, "\n");
      END PrintAbscissa;

    PROCEDURE PrintPoint (msg: TEXT; u, fu: REAL) =
      <* FATAL Wr.Failure, Thread.Alerted *>
      BEGIN 
        Wr.PutText(stderr, "    ");
        Wr.PutText(stderr, msg);
        Wr.PutText(stderr, " = ");
        Wr.PutText(stderr, Fmt.Pad(Fmt.Real(u), 12));
        Wr.PutText(stderr, " f" & msg & " = ");
        Wr.PutText(stderr, Fmt.Pad(Fmt.Real(fu), 12));
        Wr.PutText(stderr, "\n");
      END PrintPoint;

  BEGIN

    a := ax; b := bx;
    v := a + IhpSq * (b-a); fv := f(v); 
    x := v; fx := fv;
    w := v; fw := fx;

    IF good # NIL AND good(x, fx) THEN RETURN x END;
    fxOld := fx;

    tol3 := tol/3.0;

    e := 0.0;
    
    LOOP (* 20 *)
      IF verbose THEN PrintStatus() END;
      
      xm := 0.5 * (a+b);
      tol1 := Eps * ABS(x) + tol3;
      t2 := 2.0 * tol1;
      
      (* check stopping criterion: *)
      IF ABS(x-xm) <= t2 - 0.5*(b-a) THEN EXIT END;
      IF fx < fxOld THEN
        IF good # NIL AND good(x, fx) THEN EXIT END;
        fxOld := fx;
      END;
      
      IF ABS(e) > tol1 THEN 
        (* Fit parabola through "v", "x", "w"  *)
        r := (x-w)*(fx-fv);
        q := (x-v)*(fx-fw);
        p := (x-v)*q - (x-w)*r;
        q := 2.0*(q-r);
        IF q <= 0.0 THEN q := -q ELSE p := -p END;
        r := e;
        e := d;
      ELSE
        p := 0.0;
        q := 0.0;
        r := 0.0;
      END;
      IF ABS(p) < ABS(0.5*q*r) AND p > q*(a-x) AND p < q*(b-x) THEN 
        (* A parabolic-interpolation step: *)
        IF verbose THEN Message("parabolic step") END;
        d := p/q;
        u := x+d;
        
        (* "f" must not be evaluated too close to "a" or "b": *)
        IF u-a < t2 OR b-u < t2 THEN 
          d := tol1;
          IF x >= xm THEN d := -d END;
        END;
      ELSE
        (* A golden-section step: *)
        IF verbose THEN Message("golden section step") END;
        IF x < xm THEN e := b-x ELSE e := a-x END;
        d := IhpSq * e;
      END;
      
      (* f must not be evaluated too close to x: *)
      IF ABS(d) >= tol1 THEN
        u := x+d
      ELSIF d > 0.0 THEN
        u := x + tol1
      ELSE
        u := x - tol1
      END;  
      fu := f(u);

      (* update  a, b, v, w, and x: *)

      IF fx <= fu THEN
        IF u < x THEN a := u ELSE b := u END
      END;
      IF fx >= fu THEN
        IF u < x THEN b := x ELSE a := x END;
        v := w; fv := fw;
        w := x; fw := fx;
        x := u; fx := fu;
      ELSIF fu <= fw OR w = x THEN
        v := w; fv := fw;
        w := u; fw := fu;
      ELSIF fu <= fv OR v = x OR v = w THEN
        v := u; fv := fu;
      END;
    END; 
    RETURN x
  END Minimize;

BEGIN
END BrentUniMin.