MODULE JSUniMinMin EXPORTS Main;

(* Uses SimMin to optimize the coefficient Alpha OF JSUniMin... *)

IMPORT Random, SimMin, JSUniMin, RanMin, Math, Real, Text;

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

CONST
  NC = 1;

TYPE 
  Coefs =  ARRAY [0..NC-1] OF REAL;

VAR
  coefCoins: Random.T := Random.New(); 

PROCEDURE Main () =
  <* FATAL Wr.Failure, Thread.Alerted *>
  VAR coefs: ARRAY [0..NC] OF Coefs;
  VAR avgcalls: ARRAY [0..NC] OF REAL;
  VAR lo, hi: CARDINAL;
  VAR cbest: Coefs;
  VAR avgbest: REAL;
  BEGIN
    cbest := Coefs{0.125};

    (*
      Some good coefficients determined on previous runs:
    *)

    RanMin.Minimize(
      f := CallJSUniMin, 
      x := cbest,
      y := avgbest,
      dx := Coefs{0.01},
      maxRadius := 0.25,
      minRadius := 0.0001,
      maxCalls := 200,
      report := ReportCoefs,
      normalize := NormalizeCoefs,
      verbose := TRUE
    );
    Wr.PutText(stderr, "\n");
    Wr.PutText(stderr, "** Best coefficients found by RanMin:\n");
    PrintCoefs(cbest, avgbest);
    Wr.PutText(stderr, "\n");
    
    (* Minimize by SimMin: *)
    GenRootedSimplex(
      pos := cbest, 
      size := 0.02, 
      x := coefs,
      coins := coefCoins
    );
    SimMin.Minimize(
      f := CallJSUniMin, 
      x := coefs,
      y := avgcalls,
      lo := lo, hi := hi,
      minSize := 0.00001,
      maxCalls := 1000,
      report := ReportCoefs,
      normalize := NormalizeCoefs,
      verbose := TRUE
    );
    Wr.PutText(stderr, "\n");
    Wr.PutText(stderr, "** Best coefficients found by SimMin:\n");
    PrintCoefs(coefs[lo], avgcalls[lo]);

  END Main;
  
VAR
  pointCoins: Random.T := Random.New(); 

VAR 
  ncalls: CARDINAL;    (* Num of calls to F *)
  smallEnough: REAL;   (* Stop minimization when F < smallEnough. *)
  converged: BOOLEAN;  (* TRUE if F became smallenough. *)

PROCEDURE CallJSUniMin (READONLY coefs: ARRAY OF REAL): REAL =
  VAR x, fx, dfx, ddfx: REAL;
  VAR res: REAL := 0.0;
  CONST NTrials = 200;
  CONST MaxFuncCalls = 400;
  BEGIN
    Random.Start(pointCoins);
    FOR i := 0 TO NTrials-1 DO
      ncalls := 0;
      converged := FALSE;
      smallEnough := 2.0e-4;
      x := 2.0 * Random.Real(pointCoins) - 1.0;
      fx := Cubic(x);
      JSUniMin.Minimize(
        f := Cubic, 
        x := x,
        fx := fx,
        dfx := dfx,
        ddfx := ddfx,
        step := 1.0,
        minStep := 1.0e-6,
        maxStep := 1.0e+30,
        maxCalls := MaxFuncCalls,
        Alpha := coefs[0], 
        report := ReportFunction
      );
      PrintPoint ("min  ", x, fx, dfx, ddfx);
      IF converged THEN 
        res := res + FLOAT(ncalls)
      ELSE
        res := res + FLOAT(2 * MaxFuncCalls) - Log2(fx)
      END;
    END;
    RETURN res/FLOAT(NTrials)
  END CallJSUniMin;

PROCEDURE PrintCoefs (READONLY coefs: ARRAY OF REAL; nits: REAL) =
  <* FATAL Wr.Failure, Thread.Alerted *>
  BEGIN 
    Wr.PutText(stderr, "\n");
    Wr.PutText(stderr, "**** Alpha = ");
    Wr.PutText(stderr, Fmt.Pad(Fmt.Real(coefs[0]), 10));
    Wr.PutText(stderr, "  nits = " & 
      Fmt.Pad(Fmt.Real(nits, 2, Fmt.Style.Flo), 6)
    );
    Wr.PutText(stderr, "\n");
    Wr.PutText(stderr, "\n");
  END PrintCoefs;
  
PROCEDURE ReportCoefs (READONLY coefs: ARRAY OF REAL; nits: REAL): BOOLEAN =
  BEGIN 
    PrintCoefs (coefs, nits);
    RETURN FALSE
  END ReportCoefs;
  
PROCEDURE Log2(x: REAL): REAL =
  (* Log base 2 of x, clamped at smallpos *)
  BEGIN
    WITH
      ee = FLOAT(Real.MinPos, LONGREAL),
      xx = MAX(ee, FLOAT(x, LONGREAL)),
      log = (Math.log(xx) - Math.log(ee))/Math.log(2.0D0)
    DO
      RETURN FLOAT(log)
    END
  END Log2;
  
PROCEDURE NormalizeCoefs(<*UNUSED*> VAR coefs: ARRAY OF REAL) =
  BEGIN
    (*
      coefs[0] := MAX( 0.5, MIN(2.0, coefs[0]));
      coefs[1] := MAX( 0.1, MIN(0.9, coefs[1]));
      coefs[2] := MAX( 0.1, MIN(4.0, coefs[2]));
    *)
  END NormalizeCoefs;

CONST
  p = 0.0;
  q = 2.0;

PROCEDURE Quadratic (x: REAL): REAL =
  BEGIN
    INC(ncalls);
    WITH xp = x - p DO
      RETURN xp * xp;
    END;
  END Quadratic;

PROCEDURE AlmostQuadratic (x: REAL): REAL =
  BEGIN
    INC(ncalls);
    WITH xp = x - p, xq = x - q DO
      RETURN xp * xp + 0.00001 * xq * xq * xq * xq;
    END;
  END AlmostQuadratic;

PROCEDURE Cubic (x: REAL): REAL =
  BEGIN
    INC(ncalls);
    WITH xp = x - p, xq = x - q DO
      RETURN xp * xp + 0.0001 * xq * xq * xq;
    END;
  END Cubic;

PROCEDURE GenRootedSimplex(
    READONLY pos: ARRAY OF REAL; 
    size: REAL; 
    VAR x: ARRAY OF ARRAY OF REAL;
    coins: Random.T;
  ) =
  (* Generate a simplex with diameter "size" AND corner at "pos": *)
  VAR sig: REAL;
  BEGIN
    WITH 
      NC = NUMBER(x[0]),
      MC = NUMBER(x)
    DO
      <* ASSERT MC = NC + 1 *>
      IF NC = 3 THEN
        sig := 0.25
      ELSE
        WITH
          fn1 = FLOAT(NC + 1, LONGREAL),
          fn3 = FLOAT(NC - 3, LONGREAL)
        DO
          sig := FLOAT((Math.sqrt(fn1) + 2.0d0)/fn3);
        END;
      END;
      FOR i := 0 TO MC-1 DO
        FOR j := 0 TO NC-1 DO
          IF i = MC-1 THEN 
            x[i,j] := 0.0
          ELSIF j = i THEN 
            x[i,j] := 1.0 
          ELSE 
            x[i,j] := sig 
          END
        END
      END;
      
      (* rotate randomly: *)
      RandomRotateSimplex(x, coins);
      
      (* Scale and displace: *)
      FOR i := 0 TO MC-1 DO
        FOR j := 0 TO NC-1 DO
          x[i,j] := pos[j] + size * x[i,j] 
        END
      END;
    END;
  END GenRootedSimplex;

PROCEDURE GenCenteredSimplex (
    READONLY pos: ARRAY OF REAL; 
    size: REAL; 
    VAR x: ARRAY OF ARRAY OF REAL; 
    coins: Random.T;
  ) =
  (* Generate a simplex with radius "size" surrounding "pos": *)
  VAR tau: REAL;
  BEGIN
    WITH 
      NC = NUMBER(x[0]),
      MC = NUMBER(x)
    DO
      <* ASSERT MC = NC + 1 *>
      IF NC = 3 THEN
        tau := 0.5
      ELSE
        WITH
          fn1 = FLOAT(NC+1, LONGREAL),
          fn3 = FLOAT(NC-3, LONGREAL)
        DO
          tau := FLOAT((fn1 - 2.0D0 * Math.sqrt(fn1))/fn3);
        END
      END;
      FOR i := 0 TO MC-1 DO
        FOR j := 0 TO NC-1 DO
          IF i = MC-1 THEN 
            x[i,j] := -1.0
          ELSIF j = i THEN 
            x[i,j] := +1.0 
          ELSE 
            x[i,j] := - tau 
          END
        END
      END;
      
      (* rotate randomly: *)
      RandomRotateSimplex(x, coins);
      
      (* Scale and displace: *)
      FOR i := 0 TO MC-1 DO
        FOR j := 0 TO NC-1 DO
          x[i,j] := pos[j] + size * x[i,j] 
        END
      END;
    END;
  END GenCenteredSimplex;

PROCEDURE RandomRotateSimplex(
    VAR x: ARRAY OF ARRAY OF REAL; 
    coins: Random.T;
  ) =
  BEGIN
    WITH 
      NC = NUMBER(x[0]),
      MC = NUMBER(x)
    DO
      (* Do some random rotations: *)
      FOR j1 := 0 TO NC-1 DO
        WITH 
          j2 = (j1 + 1) MOD NC,
          t = FLOAT(2.0 * Math.Pi * Random.Real(coins), LONGREAL),
          s = FLOAT(Math.sin(t)),
          c = FLOAT(Math.cos(t))
        DO
          FOR i := 0 TO MC-1 DO
            WITH 
              c1 = c * x[i,j1] - s * x[i,j2],
              c2 = s * x[i,j1] + c * x[i,j2]
            DO
              x[i,j1] := c1; x[i,j2] := c2
            END
          END
        END
      END
    END
  END RandomRotateSimplex;

PROCEDURE GenRandomPoint (
    READONLY pos: ARRAY OF REAL; 
    size: REAL; 
    VAR x: ARRAY OF REAL; 
    coins: Random.T;
  ) =
  (* Generate a random point with mean "pos" and 
     approximate spread radius "size".
  *)
  BEGIN
    WITH 
      N = NUMBER(x)
    DO
      FOR i := 0 TO N-1 DO
        x[i] := pos[i] + size * (2.0 * Random.Real(coins) - 1.0)
      END
    END
  END GenRandomPoint;

PROCEDURE PrintPoint (
    msg: TEXT; 
    x, fx, dfx, ddfx: REAL;
  ) =
  <* FATAL Wr.Failure, Thread.Alerted *>
  BEGIN 
    Wr.PutText(stderr, msg);
    Wr.PutText(stderr, "  ");
    Wr.PutText(stderr, "x = ");
    Wr.PutText(stderr, Fmt.Pad(Fmt.Real(x), 10));
    Wr.PutText(stderr, "  ");
    Wr.PutText(stderr, "fx = ");
    Wr.PutText(stderr, Fmt.Pad(Fmt.Real(fx), 10));
    Wr.PutText(stderr, "  ");
    Wr.PutText(stderr, "dfx = ");
    Wr.PutText(stderr, Fmt.Pad(Fmt.Real(dfx), 10));
    Wr.PutText(stderr, "  ");
    Wr.PutText(stderr, "ddfx = ");
    Wr.PutText(stderr, Fmt.Pad(Fmt.Real(ddfx), 10));
    Wr.PutText(stderr, "\n");
  END PrintPoint;

PROCEDURE ReportFunction (x, fx, dfx, ddfx: REAL): BOOLEAN =
  BEGIN 
    IF fx<= smallEnough THEN
      converged := TRUE;
      RETURN TRUE
    ELSE
      RETURN FALSE
    END;
  END ReportFunction;
  
BEGIN
  Main ()
END JSUniMinMin.