MODULE SimMinMin EXPORTS Main;

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

IMPORT Random, SimMin, RanMin, Math, Real;

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

TYPE Coefs =  ARRAY [0..3] OF REAL;

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

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

    (*
      Some good coefficients determined on previous runs:
      
      Coefs{0.896189, 0.523152, 3.285137, 0.066973} = 243.00 (NTries = 40)
      Coefs{0.896199, 0.523154, 3.285119, 0.066967} = 245.48 (NTries = 40)
      Coefs{0.934590, 0.505900, 3.343498, 0.0960650} = 276.81 (NTries = 40)
      Coefs{0.977225, 0.495644, 3.364912, 0.0798877} = 243.80 (NTries = 10)
      Coefs{0.977210, 0.495801, 3.364864,  0.079870} = 220.00
      Coefs{0.977372, 0.495623, 3.364900, +0.079987} = 218.60
      Coefs{1.128386, 0.539610, 3.538044, -0.000098} = 292.00
      
      Coefs{1.128466, 0.539616, 3.538078, -0.000054}
      
      Coefs{1.128466, 0.539616, 3.538078, -0.000054}
      Coefs{1.120000, 0.500000, 3.000000, +0.010000}
      Coefs{0.952430, 0.532413, 2.927809, +1.000000}
      Coefs{1.129000, 0.540000, 3.520000, -0.020000}

      Coefs{0.954293, 0.459268, 1.277617, 0.0};
      Coefs{1.128890, 0.539058, 3.539521, 0.0};
      Coefs{0.954295, 0.459284, 1.277607, 0.0};
      Coefs{0.954291, 0.459284, 1.277606, 0.0};
      Coefs{1.128637, 0.539000, 3.537313, -0.000358};
      Coefs{1.130000, 0.550000, 3.550000, -0.010000};
    *)

    RanMin.Minimize(
      f := CallMin, 
      x := cbest,
      y := avgbest,
      dx := Coefs{0.25, ..},
      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 := CallMin, 
      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;
  M = N + 1;
  RSimp = 10.0;

TYPE
  Point = ARRAY [0..N-1] OF REAL;
  
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 CallMin (READONLY coefs: ARRAY OF REAL): REAL =
  VAR x: ARRAY [0..M-1] OF Point;
  VAR y: ARRAY [0..M-1] OF REAL;
  VAR lo, hi: CARDINAL;
  VAR res: REAL := 0.0;
  CONST NTrials = 40;
  CONST MaxFuncCalls = 400;
  BEGIN
    Random.Start(pointCoins);
    FOR i := 0 TO NTrials-1 DO
      GenCenteredSimplex(
        pos := Point{0.0, ..},
        size := RSimp, 
        x := x, 
        coins := pointCoins
      );
      ncalls := 0;
      converged := FALSE;
      smallEnough := 2.0e-6;
      SimMin.Minimize(
        f := Quadratic, 
        x := x,
        y := y,
        lo := lo, hi := hi,
        minSize := 1.0e-6,
        maxCalls := MaxFuncCalls,
        Alpha := coefs[0], 
        Beta := coefs[1], 
        Gamma := coefs[2],
        Delta := coefs[3],
        report := ReportFunction
      );
      PrintPoint (x[lo], y[lo]);
      IF converged THEN 
        res := res + FLOAT(ncalls)
      ELSE
        res := res + FLOAT(2 * MaxFuncCalls) - Log2(y[lo])
      END;
    END;
    RETURN res/FLOAT(NTrials)
  END CallMin;

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, "  Delta = ");
    Wr.PutText(stderr, Fmt.Pad(Fmt.Real(coefs[3]), 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(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]));
      coefs[3] := MAX(-1.0, MIN(5.0, coefs[3]));
    *)
  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): REAL =
  VAR sa, sb, sc, res : REAL := 0.0;
  BEGIN
    INC(ncalls);
    FOR i := 0 TO N-1 DO
      WITH di = x[i] - mq[i] DO
        sa := sa + di * aq[i];
        sb := sb + di * bq[i];
        sc := sc + di * cq[i];
        res := res + di * di
      END
    END;
    res := res + sa * sa + sb * sb + sc * sc;
    RETURN res
  END Quadratic;

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 PrintPoint (READONLY x: ARRAY OF REAL; y: REAL) =
  <* FATAL Wr.Failure, Thread.Alerted *>
  BEGIN 
    Wr.PutText(stderr, "    calls = " & Fmt.Pad(Fmt.Int(ncalls), 6));
    Wr.PutText(stderr, "  fmin = " & Fmt.Pad(Fmt.Real(y), 12));
    Wr.PutText(stderr, "  ( ");
    FOR j := 0 TO N - 1 DO 
      Wr.PutText(stderr, " ");
      Wr.PutText(stderr, Fmt.Pad(Fmt.Real(x[j]), 10));
    END;
    Wr.PutText(stderr, " )\n");
  END PrintPoint;
  
PROCEDURE ReportFunction (READONLY x: ARRAY OF REAL; y: REAL): BOOLEAN =
  BEGIN 
    IF y <= smallEnough THEN
      converged := TRUE;
      RETURN TRUE
    ELSE
      RETURN FALSE
    END;
  END ReportFunction;
  
BEGIN
  Main ()
END SimMinMin.