MODULE OptShape EXPORTS Main;

(*
  Adjusts the vertex coordinates of a triangulation
  so as to minimize some energy function.
  
  Created 1994 by Rober Marcone Rosi and J.Stolfi.
*)

IMPORT Wr, FileWr, Text, TextWr, OSError, Fmt, Process, Thread, ParseParams;
IMPORT LRN, LR3, Math, CPUTime;
IMPORT ParseEnergyParams, ParseMinimizerParams;
IMPORT Energy, MixedEnergy, Triang; 
IMPORT Minimizer;
IMPORT ScreenPlot;
FROM Triang IMPORT Coords, Topology, OrgV;
FROM Oct IMPORT Sym;
FROM Energy IMPORT Gradient;
FROM Stdio IMPORT stderr;

IMPORT Debug;
 
TYPE
  LONG = LONGREAL; LONGS = ARRAY OF LONGREAL;
  BOOL = BOOLEAN;  BOOLS = ARRAY OF BOOLEAN;
  CARD = CARDINAL; CARDS = ARRAY OF CARDINAL;

TYPE
  Options = RECORD
      inName: TEXT;             (* Initial guess file name (minus ".top") *)
      outName: TEXT;            (* Output file name prefix *)
      eFunction: MixedEnergy.T; (* Energy function to minimize *)
      minimizer: Minimizer.T;   (* Minimization method to use. *)
      nPasses: CARD;            (* Number OF minimization passes. *)
      maxEvals: CARD;           (* Max energy evaluations per pass. *)
      verticesPerPass: CARD;    (* Max vertices to optimze in each pass. *)
      writeEvery: CARD;         (* When to write ".st" file. *)
      showAll: BOOL;            (* TRUE to display after each pass. *)
      showBest: BOOL;           (* TRUE to display only the best config so far. *)
      wait: BOOL;               (* TRUE to wait for mouse click after each pass. *)
    END;

TYPE
  EvalRec = RECORD           (* A configuration and its energy data: *)
      c: REF Coords;         (* Vertex coordinates *)
      e: LONG;               (* Total energy *)
      termValue: REF LONGS;  (* Unweighted energy terms *)
      eDc: REF Gradient      (* Gradient of "e" relto "c" *)
    END;
    
PROCEDURE DoIt() =
  <* FATAL OSError.E, Thread.Alerted, Wr.Failure *>
  BEGIN
    WITH 
      o = GetOptions(), 
      tc = Triang.Read(o.inName)
    DO
      WITH gnuWr = FileWr.Open(o.outName & ".gnuplot") DO
        WriteGNUPlotCommands(gnuWr, o.outName, o.eFunction.term^);
        Wr.Close(gnuWr)
      END;
      (* Write triangles for animation: *)
      IF o.writeEvery < LAST(CARD) THEN 
        WITH
          su = FileWr.Open(o.outName & ".su"),
          ma = FileWr.Open(o.outName & ".ma")
        DO
          Triang.WriteRLuisSuMa(su, ma, tc.top, all := FALSE, comment := tc.comment)
        END
      END;
      WITH 
        plotWr = FileWr.Open(o.outName & ".plot"),
        stWr = FileWr.Open(o.outName & ".st")
      DO
        WritePlotComment(plotWr, o.eFunction, o.minimizer);
        WritePlotHeader(plotWr, o.eFunction.term^);
        Minimize(
          tc, o.minimizer, o.eFunction,
          nPasses := o.nPasses,
          maxEvals := o.maxEvals, 
          verticesPerPass := o.verticesPerPass,
          writeEvery := o.writeEvery,
          outName := o.outName,
          plotWr := plotWr,
          stWr := stWr,
          showAll := o.showAll,
          showBest := o.showBest,
          wait := o.wait
        );
        Wr.Close(plotWr);
        Wr.Close(stWr)
      END
    END
  END DoIt;
  
PROCEDURE Minimize(
    READONLY tc: Triang.Configuration;  (* In: initial guess, Out: minimum *)
    m: Minimizer.T;         (* Minimization method. *)
    e: MixedEnergy.T;       (* Energy function. *)
    nPasses: CARD;          (* Number of passes. *)
    maxEvals: CARD;         (* Total evaluation budget per pass. *)
    verticesPerPass: CARD;  (* Minimize this at most this many vertices per pass. *)
    writeEvery: CARD;       (* When to write the ".st" file. *)
    outName: TEXT;          (* Output file name prefix. *)
    plotWr: Wr.T;           (* Energy plots. *)
    stWr: Wr.T;             (* Configuration file for animation. *)
    showAll: BOOL;          (* TRUE to display after each eval. *)
    showBest: BOOL;         (* TRUE to display only the best config so far. *)
    wait: BOOL;             (* TRUE to wait for user clicks at each pass *)
  ) =
  (*
    Minimization consists of "nPasses" passes, each starting at the best
    configuration found in the the previous pass, and budgeted for 
    "maxEvals" energy evaluations.  At each pass, up to "verticesPerPass"
    vertices are selected (in a round-robin fashion) to be optimized,
    while the rest is held fixed.
  *)
  VAR
    totEvals: CARD := 0;            (* Counts calls to "e.eval". *)
    cpuTime: LONG := 0.0d0;         (* Accum minimization CPU time, seconds. *)
    cpuWhen: LONG;                  (* When minimization was (re)started. *)
    passCt: CARD;                   (* Counts optimization passes. *)
    window: ScreenPlot.T;           (* The optimization movie. *)
    nextVertex: CARD := 0;          (* Next vertex to optimize *)
  BEGIN
    WITH
      top = tc.top,
      NV = top.NV,
      NT = NUMBER(e.term^),
      
      rMin = NewEvalRec(NV, NT)^,             (* Best configuration found so far. *) 
      
      variable = VariableVertices(top)^,      (* Oficially variable vertices. *)
      NB = CountTrue(variable),               (* Num of variable vertices. *)
      NS = MIN(NB, verticesPerPass),          (* Num of vertices to optimize per pass. *)
      NC = 3*NS,                              (* Num of variables to optimize. *)
      selected = NEW(REF BOOLS, NV)^,         (* Vertices being minimized. *)

      (* The following are work areas for "DoMinimizeStep": *) 
      
      rWork = NewEvalRec(NV, NT)^,            (* Probe configuration. *)
      cIndex = NEW(REF CARDS, NC)^,           (* Maps minimizer vars to coords. *)
      xm = NEW(REF Minimizer.Point, 3*NV)^,   (* Argument vector for minimizer. *)
      fm = NEW(REF Minimizer.Value)^,         (* Function value for minimizer. *)
      gm = NEW(REF Minimizer.Gradient, 3*NV)^,(* Gradient for minimizer. *) 

      (* Used by "ReportEval: *)
      minComment = 
        "energy function: " & e.name() & "\n" &
        "minimization algorithm: " & m.name() & "\n" &
        "minimization passes = " & Fmt.Int(nPasses) & "\n" &
        "vertices per pass = " & Fmt.Int(verticesPerPass) & "\n" &
        "max evals per pass = " & Fmt.Int(maxEvals) & "\n" &
        "--------------------------------------------------------\n" &
        tc.comment
    DO
      
      (* Create the plot window: *)
      IF showAll OR showBest THEN
        window := MakeWindow(top, tc.c^)
      END;

      (* Start the clock: *)
      cpuWhen := CPUTime.Now();

      (* Tell the energy function what topolgy we are evaluationg: *)
      e.defTop(tc.top);
      
      (* Tell the minimizer how many variables are we going to minimize over; *)
      EVAL m.setSize(NC);

      (* Initialize the minimum configuration: *)
      rMin.c^ := tc.c^;

      passCt := 0;
      LOOP

        (* 
          Do a full evaluation once in a while, to get the total energy right.
          Ideally, a full evaluation should be necessary only once at the beginning,
          but doing it at every pass reduces the effect of rounding errors, and
          and provides a good check against programming errors.
        *)

        VAR eOld: LONG := rMin.e;
        BEGIN
          FOR v := 0 TO NV-1 DO selected[v] := TRUE END;
          e.defVar(variable);
          e.eval(rMin.c^, rMin.e, FALSE, rMin.eDc^);
          INC(totEvals);
          rMin.termValue^ := e.termValue^;
          (* Stop the clock: *)
          cpuTime := cpuTime + (CPUTime.Now() - cpuWhen);
          IF passCt # 0 AND rMin.e - eOld > 1.0d-6 * MAX(ABS(rMin.e), ABS(eOld)) THEN
            Debug.Line("pass = " & Fmt.Int(passCt DIV NV));
            Debug.LongReal("energy discrepancy = ", LRN.T{rMin.e - eOld});
            Debug.Line("")
          END;
          ReportConfiguration(
            plotWr, stWr, window, showAll, showBest, wait,
            cpuTime, totEvals, passCt, nPasses, writeEvery,
            selected, rMin, rMin, 
            e, minComment
          );
          IF passCt = 0 OR passCt = nPasses THEN
            (* Write configuration as ".top": *)
            WITH
              name = outName & 
                ARRAY BOOL OF TEXT{"-initial", "-final"}[passCt >= nPasses]
            DO
              WriteConfiguration(
                name, top, e, rMin, 
                cpuTime := cpuTime, totEvals := totEvals, passCt := passCt, 
                comment := minComment
              )
            END
          END;
          (* Restart the clock: *)
          cpuWhen := CPUTime.Now();
        END;
        
        
        IF passCt >= nPasses THEN EXIT END;

        (* Select the vertices to optimize: *)
        <* ASSERT NS <= NB *>
        IF NS = NB THEN <* ASSERT nextVertex = 0 *> END;
        FOR v := 0 TO NV-1 DO selected[v] := FALSE END;
        VAR k: CARD := 0;
        BEGIN
          WHILE k < NS DO
            WITH v = nextVertex DO
              IF variable[v] THEN
                selected[v] := TRUE;
                cIndex[3*k+0] := 3*v+0;
                cIndex[3*k+1] := 3*v+1;
                cIndex[3*k+2] := 3*v+2;
                INC(k)
              END;
              nextVertex := (v + 1) MOD NV
            END
          END;
          <* ASSERT k=NS *>
        END;
        <* ASSERT CountTrue(selected) = NS *>    

        (* Now optimize the selected vertices. *)

        PROCEDURE ReportEval(READONLY rWork, rMin: EvalRec) =
          BEGIN
            (* Stop the clock: *)
            cpuTime := cpuTime + (CPUTime.Now() - cpuWhen);
            INC (totEvals);
            ReportConfiguration(
              plotWr, stWr, window, showAll, showBest, wait,
              cpuTime, totEvals, passCt, nPasses, writeEvery, 
              selected, rWork, rMin, 
              e, minComment
            );
            (* Restart the clock: *)
            cpuWhen := CPUTime.Now();
          END ReportEval;

        BEGIN
          (* Start the clock: *)
          cpuWhen := CPUTime.Now();
          e.defVar(selected);
          DoMinimizeStep(
            m, e, rMin,
            cIndex := SUBARRAY(cIndex, 0, NC),
            maxEvals := maxEvals,
            reportEval := ReportEval,
            rWork := rWork,
            xm := SUBARRAY(xm, 0, NC),
            fm := fm,
            gm := SUBARRAY(gm, 0, NC)
          )
        END;
        
        INC(passCt)
      END (*LOOP*);
      
      (* Stop the clock: *)
      cpuTime := cpuTime + (CPUTime.Now() - cpuWhen);
      
      (* Return minimum: *)
      tc.c^ := rMin.c^
    END
  END Minimize;
  
PROCEDURE ReportConfiguration(
    plotWr: Wr.T;              (* Energy plot file *)
    stWr: Wr.T;                (* Animation ".st" file *)
    window: ScreenPlot.T;      (* Dynamic display window *)
    showAll: BOOL;             (* TRUE to show all on "window" *)
    showBest: BOOL;            (* TRUE to show only best config so far on "window" *)
    wait: BOOL;                (* TRUE to wait for user click *)
    cpuTime: LONGREAL;         (* Accumulated minimization CPU time so far *)
    totEvals: CARD;            (* Accumulated energy evaluations so far *)
    passCt: CARD;              (* Passes completed so far *)
    nPasses: CARD;             (* Max passes to perform *)
    writeEvery: CARD;          (* When to write the ".st" file. *)
    READONLY selected: BOOLS;  (* Vertices that are or were selected for optimization *)
    READONLY rWork: EvalRec;   (* Last evaluated configuration *)
    READONLY rMin: EvalRec;    (* Best configuration found so far *)
    e: MixedEnergy.T;          (* The energy function *)
    minComment: TEXT;          (* Comment for ".top" and ".st" files *)
  ) =
  BEGIN
    PlotEnergy(plotWr, cpuTime, totEvals, passCt, rWork, e.weight^);
    IF totEvals MOD writeEvery = 1 OR passCt >= nPasses THEN
      (* Write best configuration as ".st": *)
      WriteCoords(
        stWr, e, rMin, 
        cpuTime := cpuTime, totEvals := totEvals, passCt := passCt, 
        comment := ARRAY BOOL OF TEXT{"", minComment}[passCt = 0]
      );
    END;
    IF showAll OR (showBest AND rWork.e = rMin.e) THEN
      DisplayConfiguration(window, rWork.c^, selected, wait)
    END;
  END ReportConfiguration;

PROCEDURE VariableVertices(READONLY top: Topology): REF BOOLS =
  BEGIN
    WITH 
      r = NEW(REF BOOLS, top.NV)
    DO
      Triang.GetVariableVertices(top, r^);
      RETURN r
    END
  END VariableVertices;
  
PROCEDURE NewEvalRec(NV, NT: CARD): REF EvalRec =
  BEGIN
    WITH 
      r = NEW(REF EvalRec)
    DO
      r^ := EvalRec{
        c := NEW(REF Coords, NV),
        e := 0.0d0,
        termValue := NEW(REF LONGS, NT),
        eDc := NEW(REF Gradient, NV)
      };
      RETURN r
    END
  END NewEvalRec;
  
TYPE
  ReportEvalProc = PROCEDURE (READONLY rWork, rMin: EvalRec);

PROCEDURE DoMinimizeStep(
    m: Minimizer.T;             (* Minimization method *)
    e: MixedEnergy.T;           (* Energy function *)
    VAR rMin: EvalRec;          (* In: initial guess, Out: minimum *)
    READONLY cIndex: CARDS;     (* Maps minimizer args to vertex coords *)
    maxEvals: CARD;             (* Evaluation budget *)
    reportEval: ReportEvalProc; (* Called after every energy evaluation *)
    (* Work areas: *)
    VAR rWork: EvalRec;         (* Probe configuration *)
    VAR xm: Minimizer.Point;    (* Argument vector for minimizer *)
    VAR fm: Minimizer.Value;    (* Function value for minimizer *)
    VAR gm: Minimizer.Gradient; (* Gradient vector for minimizer *)
  ) =
  (*
    Performs an energy minimization for the vertices defined by the
    "cIndex" vector.
    Assumes that "m.setSize", "e.setTop", "e.setVar"
    have already been called.
  *)
  BEGIN
    WITH
      NT = NUMBER(e.term^),
      NC = NUMBER(cIndex),
      grad = m.needsGradient(),         (* TRUE to compute gradients *)
      eOffset = NEW(REF LONG)^,         (* "eval(variable)-eval(selected)" *)
      termOffset = NEW(REF LONGS, NT)^, (* Unweigthed terms of "eOffset" *)
      
      cMin = rMin.c^, 
      eMin = rMin.e,
      eDcMin = rMin.eDc^,
      termMin = rMin.termValue^,
      
      cWork = rWork.c^, 
      eWork = rWork.e,
      eDcWork = rWork.eDc^,
      termWork = rWork.termValue^
      
    DO

      PROCEDURE Initialize(
          VAR x: Minimizer.Point; 
          VAR f: Minimizer.Value; 
          VAR g: Minimizer.Gradient
        ) =
        BEGIN
          cWork := cMin;
          e.eval(cWork, eWork, grad, eDcWork);
          termWork := e.termValue^;

          (* Save offsets between full and partial evaluations: *)
          eOffset := eMin - eWork;
          eWork := eMin;
          FOR t := 0 TO NT-1 DO
            termOffset[t] := termMin[t] - termWork[t];
          END;
          termWork := termMin;
          
          reportEval(rWork, rMin);

          (* Set function and gradient: *)
          f := eWork + eOffset;
          FOR i := 0 TO LAST(x) DO 
            WITH vk = cIndex[i], v = vk DIV 3, k = vk MOD 3 DO
              x[i] := cWork[v][k];
              g[i] := eDcWork[v][k]
            END
          END;
        END Initialize;
        
      VAR nEvals: CARD := 0;

      PROCEDURE Eval(
          READONLY x: Minimizer.Point; 
          VAR f: Minimizer.Value; 
          VAR g: Minimizer.Gradient
        ): BOOLEAN =
        BEGIN
          IF nEvals >= maxEvals THEN RETURN TRUE END;

          (* Stuff the argument into the "cWork" vector: *)
          FOR i := 0 TO LAST(x) DO 
            WITH vk = cIndex[i], v = vk DIV 3, k = vk MOD 3 DO
              cWork[v][k] := x[i];
            END
          END;

          (* Evaluate the partial energy for "cWork": *)
          e.eval(cWork, eWork, grad, eDcWork);
          termWork := e.termValue^;

          (* Correct to full energy: *)
          eWork := eWork + eOffset;
          FOR t := 0 TO NT-1 DO
            termWork[t] := termWork[t] + termOffset[t]
          END;

          (* See if we got something good: *)
          IF eWork < eMin THEN
            eMin := eWork;
            cMin := cWork;
            eDcMin := eDcWork;
            termMin := termWork;
          END;

          reportEval(rWork, rMin);

          (* Return energy and gradient to optimizer: *)
          f := eWork;
          IF grad THEN
            FOR i := 0 TO LAST(g) DO 
              WITH vk = cIndex[i], v = vk DIV 3, k = vk MOD 3 DO
                g[i] := eDcWork[v][k]
              END
            END
          END;
          INC(nEvals);
          RETURN FALSE;
        END Eval;

      (* Now do a partial one to get us started: *)
      BEGIN
        Initialize(xm, fm, gm);
        m.minimize(
          eval := Eval,
          x := xm,
          fx := fm,
          dfx := gm,
          dist := Math.sqrt(FLOAT(NC, LONG)),
          tol := 0.01d0,
          check := NIL
        );
        (* Sanity check: *)
        IF fm # eMin THEN
          Debug.LongReal("** bug: fm, eMin = ", LRN.T{fm, eMin})
        END;
      END
    END
  END DoMinimizeStep;
          
PROCEDURE CountTrue(READONLY b: BOOLS): CARD = 
  VAR n: CARD := 0;
  BEGIN
    FOR i := 0 TO LAST(b) DO IF b[i] THEN INC(n) END END;
    RETURN n
  END CountTrue;

PROCEDURE MakeWindow(
    READONLY top: Topology;   
    READONLY c: Coords;       
  ): ScreenPlot.T =
  BEGIN
    (*
      Allocate window for "2*NV+1" points: 
      "[0..NV-1]" are the mesh vertices,
      "[NV..2*NV-1]" are the marks,
      "[2*NV]" is the origin.
    *)
    WITH
      NV = top.NV,
      NE = top.NE,
      NP = 2*NV + 1,
      w = NEW(ScreenPlot.T).init(NP),
      org = NP-1
    DO
      w.pause();
      w.setCoords(0, c);
      
      (* Edges* *)
      FOR i := 0 TO NE-1 DO 
        WITH 
          e = top.edge[i],
          u = OrgV(e).num, v = OrgV(Sym(e)).num
        DO
          w.segment(u, v, width := 0.0)
        END
      END;
      
      (* Vertex marks: *)
      FOR i := 0 TO NV-1 DO
        WITH p = NV+i DO
          w.setCoord(p, LR3.T{0.0d0, ..}); (* For the time being *)
          w.point(p, size := 5.0)
        END
      END;
      
      (* Origin: *)
      w.setCoord(org, LR3.T{0.0d0, ..}); (* For the time being *)
      w.point(org, size := 7.0);
      
      w.resume();
      EVAL w.waitKey();
      RETURN w
    END
  END MakeWindow;
  
PROCEDURE DisplayConfiguration(
    w: ScreenPlot.T;
    READONLY c: Coords;
    READONLY marked: BOOLS; 
    wait: BOOL;
  ) =
  BEGIN
    WITH 
      NV = NUMBER(c),
      org = 2*NV
    DO
      w.pause();
      w.setCoords(0, c);

      (* Marked vertices: *)
      VAR p: CARD := NV;
      BEGIN
        FOR v := 0 TO NV-1 DO
          IF marked[v] THEN 
            w.setCoord(p, c[v]);
            INC(p)
          END
        END;
        WHILE p < org DO 
          w.setCoord(p, LR3.T{0.0d0, ..}); INC(p)
        END;
      END;

      w.resume();
      IF wait THEN EVAL w.waitKey() ELSE w.waitDone() END;
    END;
  END DisplayConfiguration;
  
PROCEDURE WriteCoords(
    wr: Wr.T;
    e: MixedEnergy.T;
    READONLY r: EvalRec;
    cpuTime: CPUTime.T;
    totEvals: CARD;
    passCt: CARD;
    comment: TEXT;
  ) =
  (* 
    Writes only coords values to file "wr" in ".st" format
  *)    
  BEGIN
    Triang.WriteRLuisSt(
      wr, r.c^, r.eDc^,
      prec := 5,
      comment := 
        IterationComment(cpuTime, totEvals, passCt, r, e) & "\n" & 
        comment
    )
  END WriteCoords;

PROCEDURE WriteConfiguration(
    name: TEXT;
    READONLY top: Topology;
    e: MixedEnergy.T;
    READONLY r: EvalRec;
    cpuTime: CPUTime.T;
    totEvals: CARD;
    passCt: CARD;
    comment: TEXT;
  ) =
  BEGIN
    Triang.Write(
      name, top, r.c^,
      comment := 
        IterationComment(cpuTime, totEvals, passCt, r, e) & "\n" &
        comment
    )
  END WriteConfiguration;
  
PROCEDURE IterationComment(
    cpuTime: LONG;
    READONLY totEvals: CARD;
    READONLY passCt: CARD;
    READONLY r: EvalRec;
    READONLY e: MixedEnergy.T;
  ): TEXT =
  BEGIN
    WITH
      wr = NEW(TextWr.T).init()
    DO
      WritePlotHeader(wr, e.term^);
      PlotEnergy(wr, cpuTime, totEvals, passCt, r, e.weight^);
      RETURN TextWr.ToText(wr)
    END
  END IterationComment;
  
PROCEDURE PlotEnergy(
    wr: Wr.T;
    READONLY cpuTime: LONG;
    READONLY totEvals: CARD;
    READONLY passCt: CARD;
    READONLY r: EvalRec;
    READONLY weight: ARRAY OF REAL;
  ) =
  <* FATAL Thread.Alerted, Wr.Failure *>
  BEGIN
    Wr.PutText(wr, " ");
    Wr.PutText(wr, Fmt.Pad(Fmt.LongReal(cpuTime, Fmt.Style.Fix, prec := 2), 8));
    Wr.PutText(wr, " ");
    Wr.PutText(wr, Fmt.Pad(Fmt.Int(totEvals), 8));
    Wr.PutText(wr, " ");
    Wr.PutText(wr, Fmt.Pad(Fmt.LongReal(r.e, Fmt.Style.Fix, 5), 10));
    Wr.PutText(wr, " ");
    Wr.PutText(wr, Fmt.Pad(Fmt.Int(passCt), 8));
    FOR i := 0 TO LAST(r.termValue^) DO 
      Wr.PutText(wr, " ");
      WITH
        t = r.termValue[i] * FLOAT(weight[i], LONG)
      DO
        Wr.PutText(wr, Fmt.Pad(Fmt.LongReal(t, Fmt.Style.Fix, 5), 10))
      END;
    END;
    Wr.PutText(wr, "\n");
    Wr.Flush(wr);
  END PlotEnergy;

PROCEDURE WritePlotComment(wr: Wr.T; e: MixedEnergy.T; m: Minimizer.T) =
  <* FATAL Thread.Alerted, Wr.Failure *>
  BEGIN
    Wr.PutText(wr, "#");
    Wr.PutText(wr, "energy function: " & e.name());
    Wr.PutText(wr, "\n");
              
    Wr.PutText(wr, "#");
    Wr.PutText(wr, "minimizer: " & m.name());
    Wr.PutText(wr, "\n");

    Wr.Flush(wr);
  END WritePlotComment;

PROCEDURE WritePlotHeader(wr: Wr.T; READONLY term: ARRAY OF Energy.T) =
  <* FATAL Thread.Alerted, Wr.Failure *>
  BEGIN
    Wr.PutText(wr, "#");
    Wr.PutText(wr, Fmt.Pad("cpuTime", 8));
    Wr.PutText(wr, " ");
    Wr.PutText(wr, Fmt.Pad("evals", 8));
    Wr.PutText(wr, " ");
    Wr.PutText(wr, Fmt.Pad("energy", 10));
    Wr.PutText(wr, " ");
    Wr.PutText(wr, Fmt.Pad("steps", 8));
    FOR i := 0 TO LAST(term) DO 
      Wr.PutText(wr, " ");
      Wr.PutText(wr, Fmt.Pad(EnergyTag(term[i].name()), 10));
    END;
    Wr.PutText(wr, "\n");
    Wr.Flush(wr);
  END WritePlotHeader;
  
PROCEDURE WriteGNUPlotCommands(
    wr: Wr.T; 
    outName: TEXT;
    READONLY term: ARRAY OF Energy.T;
  ) =
  <* FATAL Thread.Alerted, Wr.Failure *>
  BEGIN
    Wr.PutText(wr, "set terminal X11\n");
    Wr.PutText(wr, "set title \"" & outName & "\"\n");
    Wr.PutText(wr, 
      "plot \"" & outName & ".plot\" using 1:3 title \"etot\" with lines, \\\n"
    );
    FOR i := 0 TO LAST(term) DO 
      WITH col = i + 5 DO
        Wr.PutText(wr, 
          "     \"" & outName & ".plot\" using 1:" & 
          Fmt.Int(col) & " title \"" &
          EnergyTag(term[i].name()) & 
          "\" with linespoints"
        );
        IF i < LAST(term) THEN Wr.PutText(wr, ", \\") END;
        Wr.PutText(wr, "\n");
      END
    END;
    Wr.PutText(wr, "pause 300\n");
    Wr.PutText(wr, "quit\n");
    Wr.Flush(wr);
  END WriteGNUPlotCommands;
  
PROCEDURE EnergyTag(name: TEXT): TEXT =
  BEGIN
    WITH
      n = Text.FindChar(name, '(')
    DO
      IF n = -1 THEN 
        RETURN Text.Sub(name, 0, 8)
      ELSE
        RETURN Text.Sub(name, 0, MIN(n, 8))
      END
    END
  END EnergyTag;

PROCEDURE GetOptions (): Options =
  <* FATAL Thread.Alerted, Wr.Failure *>
  VAR o: Options;
  BEGIN
    WITH pp = NEW(ParseParams.T).init(stderr) DO
      TRY
        pp.getKeyword("-inName");                               
        o.inName := pp.getNext();                         

        pp.getKeyword("-outName");                               
        o.outName := pp.getNext(); 

        IF pp.keywordPresent("-nPasses") THEN
          o.nPasses := pp.getNextInt(1, 99999);
        ELSE
          o.nPasses := 1
        END;
        
        pp.getKeyword("-maxEvals");
        o.maxEvals := pp.getNextInt(1, 999999);

        IF o.nPasses > 1 THEN 
          IF pp.keywordPresent("-verticesPerPass") THEN
            o.verticesPerPass := pp.getNextInt(1, LAST(CARD));
          ELSE
            o.verticesPerPass := LAST(CARD)
          END
        ELSE
          o.verticesPerPass := LAST(CARD)
        END;
        
        IF pp.keywordPresent("-writeEvery") THEN
          o.writeEvery := pp.getNextInt(0, 1000000)
        ELSIF pp.keywordPresent("-writeAll") THEN
          o.writeEvery := 1
        ELSE
          o.writeEvery := LAST(CARD)
        END;

        o.showAll := pp.keywordPresent("-showAll");
        o.showBest := (NOT o.showAll) AND pp.keywordPresent("-showBest");
        IF o.showAll OR o.showBest THEN
          o.wait := pp.keywordPresent("-wait")
        ELSE
          o.wait := FALSE
        END;

        o.eFunction := ParseEnergyParams.Parse(pp);
        IF o.eFunction = NIL THEN 
          pp.error("no energy specified")
        END;
        
        o.minimizer := ParseMinimizerParams.Parse(pp);
        IF o.minimizer = NIL THEN
          pp.error("no minimizer specified")
        END;

        pp.finish();                                       
      EXCEPT                                                            
      | ParseParams.Error =>                                              
          Wr.PutText(stderr, "Usage: OptShape \\\n");
          Wr.PutText(stderr, "  -inName <name>   -outName <name> \\\n");
          Wr.PutText(stderr, "  -maxEvals <num> \\\n");
          Wr.PutText(stderr, "  [ -nPasses <num> [ -verticesPerPass <num> ] ] \\\n");
          Wr.PutText(stderr, "  [ -writeEvery <num> | -writeAll ] \\\n");
          Wr.PutText(stderr, "  [ { -showAll | -showBest } [ -wait ] ] \\\n");
          Wr.PutText(stderr, ParseMinimizerParams.Help);
          Wr.PutText(stderr, " \\\n");
          Wr.PutText(stderr, ParseEnergyParams.Help);
          Wr.PutText(stderr, "\n");
          Process.Exit (1);                                              
      END;
    END;
    RETURN o
  END GetOptions;

BEGIN
  DoIt()
END OptShape.