#define PROG_NAME "OptShape"
#define PROG_DESC "???"
#define PROG_VERS "1.0"

#define OptShape_C_COPYRIGHT \
  ""
  
#define PROG_INFO \
  "" \
  " "


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

#include <Octf.h>
#include <r4.h>
#include <r3.h>
#define _GNU_SOURCE
#include <stdio.h>
#include <string.h>
#define _GNU_SOURCE
#include <stdio.h>
#define _GNU_SOURCE
#include <stdio.h>
#include <argparser.h>
// 


#include <CPUTime.h>
#include <math.h>
#include <LRN.h>
#include <ParseMinimizerParams.h>
#include <ParseEnergyParams.h>
#include <Triangulation.h>
#include <MixedEnergy.h>
#include <Energy.h>
#include <Minimizer.h>
#include <ScreenPlot.h>


#include <Tridimensional.h>
//  Coords3D_t;
#include <Triangulation.h>
//  Coords_t, Topology, OrgV;
#include <Energy.h>
//  Gradient;

#define _GNU_SOURCE
#include <stdio.h>
#include <Octf.h>
//  Clock;

#include <Debug.h>


CONST
double Epsilon = 0.0000000001;                /* aditional constants for 4D */
double From4 = (r4_t){0.0,0.0,0.0,-3.0}; /* 4D viewing parameter as    */
double To4 =   (r4_t){0.0,0.0,0.0,0.0};  /* expected by the "Wire4"    */
double Up4 =   (r4_t){0.0,1.0,0.0,0.0};  /* Interactive 4D Wireframe   */
double Over4 = (r4_t){0.0,0.0,1.0,0.0};  /* Display Program.           */

VAR
  Wa,Wb,Wc,Wd: r4_t;    /* the 4 basis vectors for the 4D viewing matrix */
  Data4Radius = 0.0;  /* radius of the 4D data */
  improvements,writestep: uint = 1;
  writetime: double = 0.0;

TYPE
typedef struct Options_t {
     char *inFileTp;           /* Initial guess file name (minus ".tp") */
     char *inFileSt;           /* Initial guess file name (minus ".st") */
     char *outFile;            /* Output file name prefix */
     eFunction: MixedEnergy.T; /* Energy Mixed function to minimize */
     minimizer: Minimizer.T;   /* Algorithm Minimization method to use */
     uint nPasses;            /* Number Of minimization passes */
     uint maxEvals;           /* Max energy evaluations per pass */
     uint nodesPerPass;    /* Max nodes top optimize in each pass */
     bool_t only3D;        /* TRUE==constrains the model to R^{3}     */
     uint writeEvery;         /* When to write ".std" file */
     uint printEvery;         /* Print energies every this many steps. */
     bool_t showAll;            /* TRUE to display after each pass. */
     bool_t showBest;           /* TRUE to display only the best conf so far.*/
     bool_t wait;               /* TRUE to wait for mouse click af each pass */
     double scale;              /* set the scale for the visualization proc. */
     bool_t complete;        /* for compute overall files */
   }

double EvalRec = RECORD        /* A Record of energy data: */
    c: REF Coords_t;        /* Pointer to Node coordinates */
    double e;          /* Total energy */
    termValue: REF double_vec_t; /* Pointer to Unweighted energy terms */
    eDc: REF Gradient     /* Pointer to Gradient of total energy "e" r.to "c"*/
  }

double double = double; double_vec_t == double_vec_t;
double bool_t = bool_t;  bool_vec_t == ARRAY OF bool_t;
double uint = uint; CARDS == uint_vec_t;
    
Options_t *GetOptions(int argc, char **argv);

int main(int argc, char **argv)
  <* FATAL OSError.E, Thread.Alerted, Wr.Failure);
  { /* compute the 4D viewing matrix */
    CalcV4Matrix();
    Options_t *o = GetOptions(argc, argv);
    char *topo_cmt = jsprintf("Created by %s on %s", PROG_NAME, Today());
    /* Random_t coins = MakeRandomSource(4615); */
    ??? tc = Triangulation.ReadToMa(o->inFileTp);
    ??? rc = Triangulation.ReadState(o->inFileSt);
    {
      fprintf(stderr, "Optimizing from: " & o->inFileTp & ".tp\n");
      if (o->complete)
      {
        ??? gnuWr = FileWr.Open(o->outFile & ".gnu");
    {
          WriteGNUPlotCommands(gnuWr, o->outFile, o->eFunction.term^);
          fclose(gnuWr)
        }
      }

      if (o->only3D)
      {
        for(v = 0; v < rc.nel; v++)
        {
          rc[v,3] = 0.0;
        }
      }
      
      ??? plotWr = FileWr.Open(o->outFile & ".plot");
      {
        if (o->complete)
        {
          WritePlotComments(plotWr, o->eFunction, o->minimizer);
          WritePlotHeader(plotWr, o->eFunction.term^);
        }
        if ((o->printEvery < (uint.nel-1)))
        {
          WritePlotHeader(stderr, o->eFunction.term^);
        }
        Minimize(
            tc, rc, o->minimizer, o->eFunction,
            nPasses = o->nPasses,
            maxEvals = o->maxEvals,
            nodesPerPass = o->nodesPerPass,
            only3D     = o->only3D,
            writeEvery = o->writeEvery,
            printEvery = o->printEvery,
            inFile  = o->inFileTp,
            outFile = o->outFile,
            plotWr = plotWr,
            showAll = o->showAll,
            showBest = o->showBest,
            wait = o->wait,
            scale = o->scale,
            complete = o->complete
         );
        fclose(plotWr)
      }
    
    return 0;
  }
  
PROCEDURE Minimize(
    *tc: Triangulation.TopCom_t;
    *rc: REF Coords_t;
    m: Minimizer.T;       /* Minimization method */
    e: MixedEnergy.T;     /* Mixed Energy function */
    uint nPasses;        /* number of passes  */
    uint maxEvals;       /* Evaluation budget per node minimization */
    uint nodesPerPass; /* Min. this at most this many nodes per pass */
    bool_t only3D;    /* TRUE==constrains the model to R^{3} */
    uint writeEvery; /* When to write ".tp" file */
    uint printEvery; /* Print energies every thsi many steps. */
    char *outFile;        /* output file name prefix */
    char *inFile;          /* input file name prefix */
    plotWr: Wr.T;         /* Energy plots */
    bool_t showAll;        /* TRUE to display after each eval. */
    bool_t showBest;       /* TRUE to display only the best config so far. */
    bool_t wait;           /* TRUE to wait for user clicks at each pass */
    double scale;          /* set the scale for the visualization proc. */
    bool_t complete;       /* TRUE to wait for user clicks at each pass */
 )

  /*
    Minimization consists of "nPasses" passes, each starting at the best
    configuration found in the previous pass, and budgeted for 
    "maxEvals" energy evaluations.  At each pass, up to "nodesPerPass"
    nodes are selected (in a round-robin fashion) to be optimized,
    while the rest is held fixed. */
 
  VAR
    totEvals: uint = 0;         /* Count calls to "e.eval". */
    cpuTime: double = 0.0;      /* Accum minimization CPU time, seconds.*/
    cpuWhen: double = 0.0;      /* When minimization was (re)started. */
    uint passCt;                   /* Counts optimization passes. */
    window: ScreenPlot.T;           /* The optimization movie. */
    nextNode: uint = 0;          /* Next node to optimize */
  { ??? top = tc.top;
    int NV = top->node.nel;
    ??? NT = ((e.term^).nel);
    ??? rMin = NewEvalRec(NV;
    with ( NT)^,          /* Best configuration found so far.
                                               Minimun Configuration */
    double variable = VariableNodes(top)^,   /* Oficially variable nodes  */
    double NB = CountTrue(variable),            /* Number of variable nodes  */
    double NS = MIN(NB, nodesPerPass),       /* Num of vert.to opt. per pass. */
    double NC = 4*NS,                           /* Num of variables to optimize. */
    double selected = bool_vec_new(NV)^,      /* Nodes being minimized */

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

      /* Used by "ReportEval: */
    double minComment = 
        "energy function: " & e.name() & "\n" \
        "minimizer: " & m.name() & "\n" \
        "minimization passes == " & Fmt.Int(nPasses) & "\n" \
        "nodes per pass == " & Fmt.Int(nodesPerPass) & "\n" \
        "max evals per pass == " & Fmt.Int(totEvals) & "\n" \
        "---------------------------------------" \
        "---------------------------------------\n" \
        "\nOptimized from: " & inFile & ".tp\n"
         & "\nCreated by OptShape: " & outFile & "\n"
  ){
      /* Create the plot window: */
      
      if ((showAll) || (showBest))
      {
        ??? c3D = ProjectTo3D(top;
          with (rc^,scale)
      ){
          window = MakeWindow(top,c3D^);
        }
      }
      
      /* Start the clock: */
      cpuWhen = CPUTime.Now();

      /* Tell the energy function what topology 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^ = rc^;

      passCt = 0;
      while (1)
      {
     
       /* 
         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.
        */
        double eOld = rMin.e;
        VAR {
          for (v = 0; v < NV; v++){ selected[v] = TRUE; }
          e.defVar(variable);
          e.eval(rMin.c^, rMin.e, FALSE, rMin.eDc^);
          totEvals++;
          rMin.termValue^ = e.termValue^;
          /* Stop the clock: */
          cpuTime = cpuTime + (CPUTime.Now() - cpuWhen);
          IF
            passCt!=0
            AND
            rMin.e - eOld > 1.0d-6 * MAX(fabs(rMin.e), fabs(eOld))
        )){
            Debug.Line("pass == " & Fmt.Int(passCt DIV NV));
            Debug.LongReal("energy discrepancy == ", LRN.T{rMin.e - eOld});
            Debug.Line("")
          }
          ReportConfiguration(outFile,
            top,plotWr, window, showAll, showBest, wait,
            cpuTime, totEvals, passCt, nPasses, writeEvery, printEvery,
            selected, rMin, rMin, 
            e, minComment,scale,complete
         );
          if ((passCt == 0) || (passCt == nPasses))
          {
            /* Write configuration as ".tp": */
            with (
            double name = outFile \
                ARRAY bool_t OF char *{"-initial", "-final"}[passCt >= nPasses]
          ){
              WriteConfiguration(
                name, top, e, rMin, 
                cpuTime = cpuTime, totEvals = totEvals, passCt = passCt, 
                comment = minComment
             )
            }
          }
          /* Restart the clock: */
          cpuWhen = CPUTime.Now();
        }
        
        
        if (passCt >= nPasses){ EXIT; }

        /* Select the nodes to optimize: */
        assert(NS <= NB);
        if (NS == NB){ assert(nextNode == 0); }
        for (v = 0; v < NV; v++){ selected[v] = FALSE; }
        uint k = 0;
        VAR {
          while (k < NS)
          {
            ??? v = nextNode;
    {
              if (variable[v])
              {
                selected[v] = TRUE;
                cIndex[4*k+0] = 4*v+0;
                cIndex[4*k+1] = 4*v+1;
                cIndex[4*k+2] = 4*v+2;
                cIndex[4*k+3] = 4*v+3;
                k++
              }
              nextNode = (v + 1) % NV
            }
          }
          assert(k==NS);
        }
        assert(CountTrue(selected) == NS);    

        void ReportEval(*rWork, rMin: EvalRec)
          {
            /* Stop the clock: */
            cpuTime = cpuTime + (CPUTime.Now() - cpuWhen);
            INC (totEvals);
            ReportConfiguration(outFile,
              top,plotWr, window, showAll, showBest, wait,
              cpuTime, totEvals, passCt, nPasses, writeEvery, printEvery,
              selected, rWork, rMin, 
              e, minComment,scale,complete
           );
            /* Restart the clock: */
            cpuWhen = CPUTime.Now();
          } ReportEval;

        {
          /* Start the clock: */
          cpuWhen = CPUTime.Now();
          e.defVar(selected);
          DoMinimizeStep(
            m, e, rMin,
            cIndex = SUBARRAY(cIndex, 0, NC),
            maxEvals = maxEvals,
            only3D = only3D,
            reportEval = ReportEval,
            rWork = rWork,
            xm = SUBARRAY(xm, 0, NC),
            fm = fm,
            gm = SUBARRAY(gm, 0, NC)
         )
        }
        
        passCt++
      } /*while (1){*/;
      
      /* Stop the clock: */
      cpuTime = cpuTime + (CPUTime.Now() - cpuWhen);
      
      /* Return minimum: */
      rc^ = rMin.c^
    }
  } /* END Minimize */

PROCEDURE ReportConfiguration(
    char *outFile;            /* output file name prefix */
    *ElemTableRec_t *top;
    plotWr: Wr.T;             /* Energy plot file */
    window: ScreenPlot.T;     /* Dynamic display window */
    bool_t showAll;            /* TRUE to show all on "window" */
    bool_t showBest;           /* T. to show only best config so far on screen*/
    bool_t wait;               /* TRUE to wait for user click */
    double cpuTime;        /* Accumulated minimization CPU time so far */
    uint totEvals;           /* Accumulated energy evaluations so far */
    uint passCt;             /* Passes completed so far */
    <*UNUSED);nPasses: uint;  /* Max passes to perform */
    <*UNUSED);writeEvery:uint;/* When to write the ".st" file. */
    uint printEvery;     /* Print energies every this many steps. */
    bool_vec_t *selected; /* Nodes that are or were selected for opt. */
    EvalRec *rWork;  /* Last evaluated configuration */
    EvalRec *rMin;   /* Best configuration found so far */
    e: MixedEnergy.T;         /* The energy function */
    char *minComment;         /* Comment for ".tp" and ".st" files */
    double scale;          /* set the scale for the visualization proc. */
    bool_t complete;
 )
  {    
    if (complete)
    {
      PlotEnergy(plotWr, cpuTime, totEvals, passCt, rWork, e.weight^);
    }
    if (totEvals % printEvery == 0)
    {
      PlotEnergy(stderr, cpuTime, totEvals, passCt, rMin,  e.weight^)
    }  
    
    /* improvements for iterations */
    if (improvements >= 10 * writestep){  writestep = 10 * writestep; }
    if (((improvements % writestep) == 0) && (complete))
    {
      ??? name = outFile & "-" & Fmt.Int(improvements) & "-e";
      {
        WriteConfiguration(
          name, top, e, rMin, 
          cpuTime = cpuTime, totEvals = totEvals, passCt = passCt, 
          comment = ARRAY bool_t OF char *{"", minComment}[passCt == 0]
        );
      }
    }
    improvements++;

    
    IF(writetime-2.0 * FLOAT(ROUND(writetime) DIV 2,double) <= 0.005) AND
      complete))
      {
      /* improvements for cpu time */
      ??? name = outFile & "-" \
         Fmt.Pad(Fmt.LongReal(cpuTime, Fmt.Style.Fix);
      with(prec=2),5,'0') & "-t"
    ){
        WriteConfiguration(
          name, top, e, rMin, 
          cpuTime = cpuTime, totEvals = totEvals, passCt = passCt, 
          comment = ARRAY bool_t OF char *{"", minComment}[passCt == 0]
       );
      }  
    }
    writetime = cpuTime;

    if ((showAll) || ((showBest) && (rWork.e == rMin.e)))
    {
      ??? c3D = ProjectTo3D(top;
        with (rWork.c^,scale)
        /*,
      double name = outFile & "-" & Fmt.Int(improvements)
        */
    ){
        DisplayConfiguration(window, c3D^, selected, wait);
        /*WriteConfigurationToOpenGL(top, c3D^);*/
        /* here display the 3D configuration with OpenGL routines */
        /* testing
        if (improvements >= 10 * writestep){ writestep = 10 * writestep; }
        if (((improvements % writestep) == 0))
        {
          WriteConfigurationToOpenGL(top, c3D^, name);
        }
        */
      }
    }
  } /* END ReportConfiguration */
  
bool_vec_t *VariableNodes(ElemTableRec_t *top)==
  /* Basically, set "vVar.el[v->num]=TRUE" iff node "v" exists and aren't
     fixed. I can set this boolean vector. */
  { ??? r = bool_vec_new(top->node.nel);
    {
      Triangulation.GetVariableNodes(top, r^);
      return r;
    }
  } /* END VariableNodes */
  
EvalRec *NewEvalRec(NV, NT: uint)==
  { ??? r = NEW(REF EvalRec);
    {
      r^ = EvalRec
      {
      c = NEW(REF Coords_t, NV),
      e = 0.0,
      double_vec_t termValue = double_vec_new(NT),
      eDc = NEW(REF Gradient, NV)
      };
      return r;
    }
  } /* END NewEvalRec */
  
TYPE
double ReportEvalProc = PROCEDURE (*rWork, rMin: EvalRec);

PROCEDURE DoMinimizeStep(
    m: Minimizer.T;             /* Minimization method */
    e: MixedEnergy.T;           /* Energy function */
    EvalRec *rMin;          /* In: initial guess, Out: minimum */
    CARDS *cIndex;     /* Maps minimizer args to node coords */
    uint maxEvals;         /* Evaluation budget */
    bool_t only3D;           /* TRUE==constrains the model to R^{3}  */
    ReportEvalProc reportEval; /* Called after every energy evaluation */
    /* Work areas: */
    EvalRec *rWork;         /* 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 nodes defined by the
    "cIndex" vector.
    Assumes that "m.setSize", "e.setTop", "e.setVar"
    have already been called. */
  { ??? NT = ((e.term^).nel);
    ??? NC = (cIndex.nel);
    ??? grad = m.needsGradient();
    with (         /* TRUE to compute gradients */
    double eOffset = NEW(REF double)^,     /* "eval(variable)-eval(selected)" */
    double termOffset = double_vec_new(NT)^, /* Unweigthed terms of "eOffset" */
    double cMin = rMin.c^;
    double eMin = rMin.e;
    double eDcMin = rMin.eDc^;
    double termMin = rMin.termValue^;
    double cWork = rWork.c^;
    double eWork = rWork.e;
    double eDcWork = rWork.eDc^;
    double termWork = rWork.termValue^
  ){
      void Initialize(
         VAR x: Minimizer.Point; 
         VAR f: Minimizer.Value; 
         VAR g: Minimizer.Gradient
       )
        {
          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; t < NT; t++)
          {
            termOffset[t] = termMin[t] - termWork[t];
          }
          termWork = termMin;
          
          reportEval(rWork, rMin);

          /* Set function and gradient: */
          f = eWork + eOffset;
          for (i = 0 TO (x.nel-1))
          {
            ??? vk = cIndex[i], v == vk DIV 4, k == vk % 4;
    {
              x[i] = cWork[v][k];
              if ((k == 3) && (only3D
              )){ g[i] = 0.0;
                     assert(x[i] == 0.0);
                }
                else
                { g[i] = eDcWork[v][k]; }
            }
          }
        } Initialize;

      void Eval(
          *x: Minimizer.Point; 
          VAR f: Minimizer.Value; 
          VAR g: Minimizer.Gradient
       )
        {
          /* Stuff the argument into the "cWork" vector: */
          for (i = 0 TO (x.nel-1))
          {
            ??? vk = cIndex[i], v == vk DIV 4, k == vk % 4;
    {
              cWork[v][k] = x[i];
            }
          }

          /* 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; t < NT; t++)
          {
            termWork[t] = termWork[t] + termOffset[t]
          }

          /* See if we got something good: */
          if (eWork < eMin)
          {
            eMin = eWork;
            cMin = cWork;
            eDcMin = eDcWork;
            termMin = termWork;
          }

          reportEval(rWork, rMin);

          /* Return energy and gradient to optimizer: */
          f = eWork;
          if (grad)
          {
            for (i = 0 TO (g.nel-1))
            {
              ??? vk = cIndex[i], v == vk DIV 4, k == vk % 4;
    {
                if ((k == 3) && (only3D)){ g[i] = 0.0
                                    }
                                    else
                                    { g[i] = eDcWork[v][k]; }
              }
            }
          }
        } Eval;

      /* Now do a partial one to get us started: */
      {
        Initialize(xm, fm, gm);
        m.minimize(
          eval = Eval,
          x = xm,
          f = fm,
          g = gm,
          dist = Math.pow(FLOAT(NC,double),1.0/3.0),
          /* dist==rough guess between the starting point and the optimun set */
          tol = 0.05,
          /* upper bound tolerance for the size of the optimun set */
          flat = 0.0,
          maxEvals = maxEvals
       );
        /* Sanity check: */
        if (fm != eMin)
        {
          Debug.LongReal("** bug: fm, eMin == ", LRN.T{fm, eMin})
        }
      }
    }
  } /* END DoMinimizeStep */

uint CountTrue(*b: bool_vec_t)
      uint n = 0;
  { for (i = 0 TO (b.nel-1)){ if (b[i]){ n++; } }
    return n;
  } /* END CountTrue */

ScreenPlot.T MakeWindow(ElemTableRec_t *top,   
    *c3: Coords3D_t;       
 )
  { /*
      Allocate window for "2*NV+1" points: 
      "[0..NV-1]" are the mesh nodes,
      "[NV..2*NV-1]" are the marks,
      "[2*NV]" is the origin.
    */
    int NV = top->node.nel;
    ??? NE = top->NE;
    ??? NP = 2*NV + 1;
    ??? w = NEW(ScreenPlot.T).init(NP);
    ??? org = NP-1;
    {
      w.pause();
      w.setCoords(0, c3);
      
      /* @{Edge->?}s* */
      for (i = 0; i < NE; i++)
      {
        ??? e = top->edge[i];
          ??? a = e.pa;
        Node_t u = OrgV(a)->num;
        Node_t v = OrgV(Clock(a))->num;
        {
          w.segment(u, v, width = 0.0)
        }
      }
      
      /* Node marks: */
      for (i = 0; i < NV; i++)
      {
        ??? p = NV+i;
    {
          w.setCoord(p, (r3_t){0.0, ..}); /* For the time being */
          w.point(p, size = 5.0)
        }
      }
      
      /* Origin: */
      w.setCoord(org, (r3_t){0.0, ..}); /* For the time being */
      w.point(org, size = 7.0);
      
      w.resume();
      EVAL w.waitKey();
      return w;
    }
  } /* END MakeWindow */

PROCEDURE DisplayConfiguration(
    w: ScreenPlot.T;
    *c3: Coords3D_t;
    bool_vec_t *marked; 
    bool_t wait;
 )
  { int NV = (c3.nel);
    ??? org = 2*NV;
    {
      w.pause();
      w.setCoords(0, c3);

      /* Marked nodes: */
      uint p = NV;
      VAR {
        for (v = 0; v < NV; v++)
        {
          if (marked[v])
          {
            w.setCoord(p, c3[v]);
            p++
          }
        }
        while (p < org)
        {
          w.setCoord(p, (r3_t){0.0, ..}); p++
        }
      }

      w.resume();
      if (wait){ EVAL w.waitKey() }else{ w.waitDone(); }
    }
  } /* END DisplayConfiguration */

void WriteConfiguration(char *name, ElemTableRec_t *top,
    e: MixedEnergy.T;
    EvalRec *r;
    cpuTime: CPUTime.T;
    uint totEvals;
    uint passCt;
    char *comment;
 )
  { WriteState(name,top, r.c^, cmt = 
        comment & IterationComments(cpuTime, totEvals, passCt, r, e) & "\n");
  } /* END WriteConfiguration */
  
PROCEDURE IterationComments(
    double cpuTime;
    uint *totEvals;
    uint *step;
    EvalRec *r;
    *e: MixedEnergy.T;
 ): char *==
  { ??? wr = NEW(TextWr.T).init();
    {
      WritePlotHeader(wr, e.term^);
      PlotEnergy(wr, cpuTime, totEvals, step, r, e.weight^);
      return TextWr.ToText(wr);
    }
  } /* END IterationComments */
  
PROCEDURE PlotEnergy(
    FILE *wr,
    double *cpuTime;
    uint *totEvals;
    uint *step;
    EvalRec *r;
    *weight: ARRAY OF REAL;
 )
  
  { PutText(wr, " ");
    PutText(wr, Fmt.Pad(Fmt.LongReal(cpuTime,Fmt.Style.Fix, prec = 2), 8));
    PutText(wr, " ");
    PutText(wr, Fmt.Pad(Fmt.Int(totEvals), 8));
    PutText(wr, " ");
    PutText(wr, Fmt.Pad(Fmt.LongReal(r.e, Fmt.Style.Fix, 5), 12));
    PutText(wr, " ");
    PutText(wr, Fmt.Pad(Fmt.Int(step), 8));
    for (i = 0 TO LAST(r.termValue^))
    {
      PutText(wr, " ");
      ??? t = r.termValue[i] * ((double)weight[i]);
      {
        PutText(wr, Fmt.Pad(Fmt.LongReal(t, Fmt.Style.Fix, 5), 12))
      }
    }
    PutText(wr, "\n");
    fflush(wr);
  } /* END PlotEnergy */

void WritePlotComments(FILE *wr, e: MixedEnergy.T; m: Minimizer.T)
  
  { PutText(wr, "#");
    PutText(wr, " energy function: " & e.name());
    PutText(wr, "\n");
              
    PutText(wr, "#");
    PutText(wr, " minimizer: " & m.name());
    PutText(wr, "\n");

    fflush(wr);
  } /* END WritePlotComments */

void WritePlotHeader(FILE *wr, *term: ARRAY OF Energy.T)
  
  { PutText(wr, "#");
    PutText(wr, Fmt.Pad("cpuTime", 8));
    PutText(wr, " ");
    PutText(wr, Fmt.Pad("evals", 8));
    PutText(wr, " ");
    PutText(wr, Fmt.Pad("energy", 12));
    PutText(wr, " ");
    PutText(wr, Fmt.Pad("step", 8));
    for (i = 0 TO (term.nel-1))
    {
      PutText(wr, " ");
      PutText(wr, Fmt.Pad(EnergyTag(term[i].name()), 12));
    }
    PutText(wr, "\n");
    fflush(wr);
  } /* END WritePlotHeader */
  
void WriteGNUPlotCommands(FILE *wr, outFile: char *;
                               *term: ARRAY OF Energy.T)
  
  { PutText(wr, "set terminal X11\n");
    /*PutText(wr, "set output \"" & outFile & ".tex\"\n");*/
    PutText(wr, "set xlabel \"Cpu-Time\"\n");
    PutText(wr, "set ylabel \"Energy\"\n");
    PutText(wr, "set title \"" & outFile & "\"\n");
    PutText(wr, "plot \"" & outFile  \
       ".plot\" using 1:3 title \"TotalE\" with lines, \\\n");
    for (i = 0 TO (term.nel-1))
    {
      ??? col = i + 5;
    {
        PutText(wr, "     \"" & outFile & ".plot\" using 1:" \
          Fmt.Int(col) & " title \"" &EnergyTag(term[i].name()) \
          "\" with linespoints");
        if ((i < (term.nel-1))){ PutText(wr, ", \\"); }
        PutText(wr, "\n");
      }
    }
    PutText(wr, "pause 30\n");
    PutText(wr, "quit\n");
    fflush(wr);
  } /* END WriteGNUPlotCommands */
  
PROCEDURE EnergyTag(char *name): char *==
  { ??? n = Text.FindChar(name;
    with ( '(')
  ){
      if (n == -1)
      {
        return Text.Sub(name, 0, 8);
      }
      else
      {
        return Text.Sub(name, 0, MIN(n, 8));
      }
    }
  } /* END EnergyTag */

PROCEDURE CalcV4Matrix()
  /* 
    This procedure computes the four basis vectors for the 4D viewing matrix,
    Wa,Wb,Wc, and Wd.  Note that the Up vector transforms to Wb, the Over
    vector transforms to Wc, and the line of sight transforms to Wd.  The
    Wa vector is then computed from Wb,Wc and Wd. */
  
  { /* Calculate Wd, the 4th coordinate basis vector and line-of-sight. */

    Wd = r4_Sub(To4,From4);
    if ((r4_Norm(Wd) < Epsilon))
    {
      fprintf(stderr, "4D To Point and From Point are the same\n");
      Process.Exit¦(1);
    }      
    Wd = r4_Scale(1.0/r4_Norm(Wd), Wd);

    /* Calculate Wa, the X-axis basis vector. */

    Wa = r4_cross(Up4,Over4,Wd);
    if ((r4_Norm(Wa) < Epsilon))
    {
      fprintf(stderr, 
                 "4D up, over and view vectors are not perpendicular\n");
      Process.Exit¦(1);
    }      
    Wa = r4_Scale(1.0/r4_Norm(Wa), Wa);

    /* Calculate Wb, the perpendicularized Up vector. */

    Wb = r4_cross(Over4,Wd,Wa);
    if ((r4_Norm(Wb) < Epsilon))
    {
      fprintf(stderr, "Invalid 4D over vector\n");
      Process.Exit¦(1);
    }       
    Wb = r4_Scale(1.0/r4_Norm(Wb), Wb);

    /* Calculate Wc, the perpendicularized Over vector.  Note that the
       resulting vector is already normalized, since Wa, Wb and Wd are all
       unit vectors.
    */
    Wc = r4_cross(Wd,Wa,Wb);
  } /* END CalcV4Matrix */  


Coords3D_t *ProjectTo3D(ElemTableRec_t *top, Coords_t *c,
    double scale;
 )==
  /*
    This procedure project all nodes of a configuration from R4 to R3,
    such that permit to create 3D coordinate needs for the "DisplayConfigura-
    tion" procedure. */

  r3_t ProjectNodeTo3D(co: r4_t)
    /*
       Project one 4D node to 3D. This projection is simples in the sense
       that only permit 4D parallel projection and not perspective.
    */
    r3_t c3 ;
    {
      ??? TempV = r4_Sub(co;
        with (From4),
      double rtemp = 1.0/Data4Radius
    ){  
        c3[0] = scale * rtemp * r4_dot(TempV, Wa);
        c3[1] = scale * rtemp * r4_dot(TempV, Wb);
        c3[2] = scale * rtemp * r4_dot(TempV, Wc); 
        return c3;
      }
    } ProjectNodeTo3D;

    void MaxRadius(
      *ElemTableRec_t *top;
      Coords_t *c;
      *ctr: r4_t;
     )
      {
        for (i = 0; i < top->node.nel; i++)
        {
          assert(OrgV(top->node[i])->num == i);
          Node_t v = OrgV(top->node[i]);
          ??? Temp4 = r4_Sub(c[v->num];
          with (ctr),
            dist  == r4_dot(Temp4, Temp4)
        ){
            if (dist > Data4Radius){ Data4Radius = dist;}
          }
        };
      } MaxRadius;

  Coords_t *3D c3D ;
    nc : uint =  (c.nel);
  { MaxRadius(top,c,To4);
    c3D = NEW(REF Coords3D_t, nc);
    for (i = 0; i < nc; i++)
    {
      c3D^[i] = ProjectNodeTo3D(c[i]);
    }
    return c3D;
  } /* END ProjectTo3D */

Options_t GetOptions ()
  { Options_t *o = (Options_t *)malloc(sizeof(Options_t));
    argparser_t *pp = argparser_new(stderr, argc, argv);
    argparser_set_help(pp, PROG_NAME " version " PROG_VERS ", usage:\n" PROG_HELP);
    argparser_set_info(pp, PROG_INFO);
    argparser_process_help_info_options(pp);
    
    argparser_get_keyword(pp, "-inFileTp");                               
        o->inFileTp = argparser_get_next(pp);           

        argparser_get_keyword(pp, "-inFileSt");                               
        o->inFileSt = argparser_get_next(pp);                         

        argparser_get_keyword(pp, "-outFile");                               
        o->outFile = argparser_get_next(pp); 

        if ((argparser_keyword_present(pp, "-nPasses")))
        {
          o->nPasses = argparser_get_next_int(pp, 1, 99999);
        }
        else
        {
          o->nPasses = 1
        }
        
        argparser_get_keyword(pp, "-maxEvals");
        o->maxEvals = argparser_get_next_int(pp, 1, 999999);

        if (o->nPasses > 1)
        {
          if ((argparser_keyword_present(pp, "-nodesPerPass")))
          {
            o->nodesPerPass = argparser_get_next_int(pp, 1, (uint.nel-1));
          }
          else
          {
            o->nodesPerPass = (uint.nel-1)
          }
        }
        else
        {
          o->nodesPerPass = (uint.nel-1)
        }
        
        if ((argparser_keyword_present(pp, "-writeEvery")))
        {
          o->writeEvery = argparser_get_next_int(pp, 0, 1000000)
        }
        else
         if ((argparser_keyword_present(pp, "-writeAll"))){
          o->writeEvery = 1
        }
        else
        {
          o->writeEvery = (uint.nel-1)
        }

        if ((argparser_keyword_present(pp, "-printEvery")))
        {
          o->printEvery = argparser_get_next_int(pp, 0, 1000000)
        }
        else
        {
          o->printEvery = (uint.nel-1)
        }

        o->showAll = argparser_keyword_present(pp, "-showAll");
        o->showBest = (NOT o->showAll)) && (argparser_keyword_present(pp, "-showBest");
        if ((o->showAll) || (o->showBest))
        {
          o->wait = argparser_keyword_present(pp, "-wait");
          if ((argparser_keyword_present(pp, "-scale")))
          {
            o->scale = pp.getNextLongReal(0.10, 100.00);
          }
          else
          {
            o->scale = 1.0;
          }
          o->wait = FALSE;
        }

        o->only3D   = argparser_keyword_present(pp, "-only3D");
        o->complete = argparser_keyword_present(pp, "-complete");

        o->eFunction = ParseEnergyParams.Parse(pp);
        if (o->eFunction == NULL)
        {
          argparser_error(pp, "no energy specified")
        }
        
        o->minimizer = ParseMinimizerParams.Parse(pp);
        if (o->minimizer == NULL)
        {
          argparser_error(pp, "no minimizer specified")
        }

        argparser_finish(pp);                                       
      
  -----------------------------------
#define _HELP \
  fprintf(stderr, "Usage: \\\n");
         fprintf(stderr, "  OptShape -inFileTp <name>  \\\n");
         fprintf(stderr, "   -inFileSt <name> -outFile <name> \\\n");
         fprintf(stderr, "   -maxEvals <num> \\\n");
         fprintf(stderr, "   [-nPasses <num> [ -nodesPerPass <num> ] ] [-only3D]\\\n");
         fprintf(stderr, "   [-writeEvery <num> | -writeAll ] \\\n");
         fprintf(stderr, "   [-printEvery <num> ] \\\n");
         fprintf(stderr, "   [ [ -showAll | -showBest ]  [-wait] [-scale]\\\n");
         fprintf(stderr, "   [-complete ]\\\n");
         fprintf(stderr, ParseMinimizerParams.Help);
         fprintf(stderr, " \\\n");
         fprintf(stderr, ParseEnergyParams.Help);
         fprintf(stderr, "\n");
         
END¦                                              
      }
    }
    return o;
  } /* END GetOptions */

{
  DoIt()
} OptShape.