MODULE GradMinMin EXPORTS Main;

(* Uses SimplexMin to optimize the coefficients Alpha, Beta, Gamma... *)

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

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

CONST
  NC = 3;

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{2.0, 0.5, 1.20};

    (*
      Some good coefficients determined on previous runs:
      Alpha =   2.075516  Beta =   0.491903  Gamma =   0.924999  nits =   7.69
      Quadratic, smallEnough = 2e-6:
      Alpha =   1.788036  Beta =   0.380159  Gamma =   1.144277  nits = 102.96
      Alpha =   1.871431  Beta =   0.457466  Gamma =   1.209523  nits = 102.21
      Alpha =   1.951707  Beta =   0.475276  Gamma =   1.240586  nits = 103.96
      Quadratic, smallEnough = 2e-4:
      Alpha =    1.99363  Beta =   0.493611  Gamma =   1.209342  nits =  82.10
      Alpha =   2.000577  Beta =   0.506291  Gamma =    1.19862  nits =  81.80
    *)

    RanMin.Minimize(
      f := CallGradMin, 
      x := cbest,
      y := avgbest,
      dx := Coefs{0.01, 0.01, 0.20},
      maxRadius := 2.0,
      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 := CallGradMin, 
      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;
  
CONST 
  N = 5;
  RSimp = 10.0;

TYPE
  Point = ARRAY [0..N-1] OF REAL;
  Vector = Point;
  
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 CallGradMin (READONLY coefs: ARRAY OF REAL): REAL =
  VAR x: Point;
  VAR y: REAL;
  VAR dy: Vector;
  VAR res: REAL := 0.0;
  CONST NTrials = 200;
  CONST MaxFuncCalls = 400;
  BEGIN
    Random.Start(pointCoins);
    FOR i := 0 TO NTrials-1 DO
      GenRandomPoint(
        pos := Point{0.0, ..},
        size := RSimp, 
        x := x, 
        coins := pointCoins
      );
      ncalls := 0;
      converged := FALSE;
      smallEnough := 2.0e-4;
      GradMin.Minimize(
        f := Quadratic, 
        x := x,
        y := y,
        dy := dy,
        minStep := 1.0e-6,
        maxStep := 1.0e+30,
        maxCalls := MaxFuncCalls,
        Alpha := coefs[0], 
        Beta := coefs[1],
        Gamma := coefs[2],
        report := ReportFunction
      );
      PrintPoint ("min  ", x, y, dy);
      IF converged THEN 
        res := res + FLOAT(ncalls)
      ELSE
        res := res + FLOAT(2 * MaxFuncCalls) - Log2(y)
      END;
    END;
    RETURN res/FLOAT(NTrials)
  END CallGradMin;

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, "  Beta = ");
    Wr.PutText(stderr, Fmt.Pad(Fmt.Real(coefs[1]), 10));
    Wr.PutText(stderr, "  Gamma = ");
    Wr.PutText(stderr, Fmt.Pad(Fmt.Real(coefs[2]), 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
  mq = Point{0.0, +0.1, +0.2, +0.3, +0.4};
  aq = Point{2.0, -4.0,  0.5, +3.0,  2.5};
  bq = Point{3.0, +2.0,  1.5, -1.0,  7.5};
  cq = Point{3.0, -3.0, -3.5, +1.0,  0.5};

PROCEDURE Quadratic (
    READONLY x: ARRAY OF REAL; 
    VAR y: REAL; 
    VAR dy: ARRAY OF REAL
  ) =
  VAR yr, ya, yb, yc : REAL := 0.0;
      dyr, dya, dyb, dyc: Vector;
  BEGIN
    INC(ncalls);
    yr := 0.0;
    FOR i := 0 TO N-1 DO
      WITH di = x[i] - mq[i] DO
        ya := ya + di * aq[i];  dya[i] := aq[i];
        yb := yb + di * bq[i];  dyb[i] := bq[i];
        yc := yc + di * cq[i];  dyc[i] := cq[i];
        yr := yr + di * di;     dyr[i] := 2.0 * di;
      END
    END;
    y := yr + ya * ya + yb * yb + yc * yc;
    FOR i := 0 TO N-1 DO
      dy[i] := dyr[i]
        + 2.0 * ya * dya[i]
        + 2.0 * yb * dyb[i]
        + 2.0 * yc * dyc[i]
    END;
  END Quadratic;

PROCEDURE QuadValley(
    READONLY x: ARRAY OF REAL; 
    VAR y: REAL; 
    VAR dy: ARRAY OF REAL
  ) =
  VAR ya: REAL := 0.0;
      dya: Vector;
  BEGIN
    INC(ncalls);
    ya := 0.0;
    FOR i := 0 TO N-1 DO
      WITH di = x[i] - mq[i] DO
        ya := ya + di * aq[i];  dya[i] := aq[i];
      END
    END;
    y := ya * ya;
    FOR i := 0 TO N-1 DO
      dy[i] :=  2.0 * ya * dya[i]
    END;
  END QuadValley;

PROCEDURE CubeValley(
    READONLY x: ARRAY OF REAL; 
    VAR y: REAL; 
    VAR dy: ARRAY OF REAL
  ) =
  VAR t, s, m, dsdt: REAL;
      dt: Vector;
  BEGIN
    INC(ncalls);
    m := 0.0;
    FOR i := 0 TO N-1 DO m := m + aq[i] * aq[i] END;
    m := FLOAT(Math.sqrt(FLOAT(m, LONGREAL)));
    t := 0.0;
    FOR i := 0 TO N-1 DO
      WITH vi = x[i] - mq[i] DO
        t := t + vi * aq[i]/m;  dt[i] := aq[i]/m;
      END
    END;
    s := 10.0 + ABS(t);
    IF t > 0.0 THEN dsdt := 1.0 ELSE dsdt := -1.0 END;
    y := t*t*s;
    FOR i := 0 TO N-1 DO
      dy[i] :=  (2.0 * s + t * dsdt ) * t * dt[i]
    END;
  END CubeValley;

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; 
    READONLY z: ARRAY OF REAL;
    READONLY u: REAL;
    READONLY du: ARRAY OF REAL;
  ) =
  <* FATAL Wr.Failure, Thread.Alerted *>
  VAR M := Text.Length(msg);
  BEGIN 
    Wr.PutText(stderr, msg);
    Wr.PutChar(stderr, ' ');
    Wr.PutText(stderr, "x =  ( ");
    FOR j := 0 TO N-1 DO 
      Wr.PutText(stderr, " ");
      Wr.PutText(stderr, Fmt.Pad(Fmt.Real(z[j]), 10));
    END;
    Wr.PutText(stderr, " )\n");

    FOR i := 0 TO M DO Wr.PutChar(stderr, ' ') END;
    Wr.PutText(stderr, "f =    " & Fmt.Pad(Fmt.Real(u), 12));
    Wr.PutText(stderr, "    calls = " & Fmt.Pad(Fmt.Int(ncalls), 6));
    Wr.PutText(stderr, "\n");
    
    FOR i := 0 TO M DO Wr.PutChar(stderr, ' ') END;
    Wr.PutText(stderr, "df = ( ");
    FOR j := 0 TO N-1 DO 
      Wr.PutText(stderr, " ");
      Wr.PutText(stderr, Fmt.Pad(Fmt.Real(du[j]), 10));
    END;
    Wr.PutText(stderr, " )\n");
    Wr.PutText(stderr, "\n");
  END PrintPoint;

PROCEDURE ReportFunction (
    READONLY x: ARRAY OF REAL; 
    y: REAL;
    READONLY dy: ARRAY OF REAL; 
  ): BOOLEAN =
  BEGIN 
    IF y <= smallEnough THEN
      converged := TRUE;
      RETURN TRUE
    ELSE
      RETURN FALSE
    END;
  END ReportFunction;
  
BEGIN
  Main ()
END GradMinMin.