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

#define TestEnergy_C_COPYRIGHT \
  ""
  
#define PROG_INFO \
  "" \
  " "

  /*
  Tests an energy function.  Takes as input two ".st" files with same
  topology but different coordinate vectors, and evaluates a given
  energy function for all configurations that interpolate between the
  two, at "nSteps" equal intervals.  
  
  The output is a ".plot" file that can be examined with "gnuplot".
  Optionally, the program writes the ".st" file with minimum energy.
*/

#include <ParseEnergyParams.h>
#include <MixedEnergy.h>
#include <Triangulation.h>


#define _GNU_SOURCE
#include <stdio.h>
#include <argparser.h>
#include <r4.h>


#include <Triangulation.h>
//  Coords_t;


TYPE
double bool_vec_t = ARRAY OF bool_t;

TYPE
typedef struct Options_t {
      char *tpFile;              /* Topology file (minus ".tp") */
      char *aFile, bFile;        /* Configuration files (minus ".st") */
      char *outFile;             /* Output file names (minus ".plot"/".st") */
      eFunction: MixedEnergy.T;  /* Energy function */
      uint nSteps;          /* Number of steps to take */
      bool_t showMin;          /* TRUE to write interpolated ".st" files */
    }
    
Options_t *GetOptions(int argc, char **argv);

int main(int argc, char **argv)
  <* FATAL Wr.Failure, Thread.Alerted, OSError.E);
  double eMin = (double.nel-1);
    double sMin;
  { Options_t *o = GetOptions(argc, argv);
    char *topo_cmt = jsprintf("Created by %s on %s", PROG_NAME, Today());
    /* Random_t coins = MakeRandomSource(4615); */
    ??? tpData = Triangulation.ReadToMa(o->tpFile);
    ??? top = tpData.top;
    ??? ac = Triangulation.ReadState(o->aFile)^;
    ??? bc = Triangulation.ReadState(o->bFile)^;
    ??? NV = top->node.nel;
    ??? c = r4_vec_new(top->node.nel)^;
    ??? cDs = r4_vec_new(top->node.nel)^;
    ??? eDc = r4_vec_new(top->node.nel)^;
    ??? vTot = bool_vec_new(top->node.nel)^;
    ??? vSome = bool_vec_new(top->node.nel)^;
    ??? pwr = FileWr.Open(o->outFile & ".plot");
    {
      assert((ac.nel) == NV);
      assert((bc.nel) == NV);
      WriteHeader(pwr);
      o->eFunction.defTop(top);
      Triangulation.GetVariableNodes(top, vTot);
      for (i = 0; i < NV; i++)
      {
        cDs[i] = r4_Sub(bc[i], ac[i]);
        vSome[i] = r4_Norm(cDs[i])!=0.0
      }
      GradTest("initial", o->eFunction, ac, vTot, eDc);
      GradTest("final",   o->eFunction, bc, vTot, eDc);

      for (i = 0; i <= (o->nSteps); i++)
      {
        ??? s = ((double)i)/((double)o->nSteps);
        {
          fprintf(pwr, " "); 
          fprintf(pwr, " " & FLR(s, 6, 4)); 
          /* Define coordinates */
          for (i = 0; i < top->node.nel; i++)
          {
            c[i] = r4_Mix(1.0 - s, ac[i], s, bc[i])
          }
          
          /* Evaluate with all or only some nodes variable: */
          ??? eTot = EvalTest(pwr;
          with (o->eFunction, c, cDs, vTot,  eDc),
          double eSome = EvalTest(pwr, o->eFunction, c, cDs, vSome, eDc)
        ){
            /* Compare: */
            fprintf(pwr, "  "); 
            fprintf(pwr, " " & ELR(eTot-eSome, 15, 8));
            if (eTot < eMin)
            {
              eMin = eTot;
              sMin = s
            }
          }
          fprintf(pwr, "\n");
          fflush(pwr);
        }
      }
      fclose(pwr);
      if (o->showMin)
      {
        ??? name = o->outFile & "-min";
        ??? cmt = 
            "  Topology file == " & o->tpFile & "\n" \
            "  Configuration files: " & o->aFile & ";
        with (\n  " & o->bFile & "\n" \
            "  Energy function to test: " & o->eFunction.name() & "\n" \
            "  Minimum energy: " & FLR(eMin, 16, 8) & "\n" \
            "  Position of minimum: " & Fmt.LongReal(sMin)
        DO
          for (i = 0; i < NV; i++)
          {
            c[i] = r4_Mix(1.0 - sMin, ac[i], sMin, bc[i])
          }
          WriteState(name, top, c, cmt)          
        }
      }
    
    return 0;
  }
  
PROCEDURE GradTest(
    char *cfgName;                 /* Configuration name. */
    eFunction: MixedEnergy.T;      /* Energy function. */
    Coords_t *c;                 /* Node coordinates. */
    bool_vec_t *v;             /* Mutable nodes. */
    Coords_t *eDc;               /* Work area for gradient. */
 )
  
      
    double e, ep, cOld;
    eDcNum: r4_t;
  CONST DiffStep == 0.0001;
  { fprintf(stderr, "Gradient test for " & cfgName & " configuration\n"); 
    fprintf(stderr, "\n"); 
    int NV = (c.nel);
    ??? eDcXXX = NEW(REF Coords_t;
    with ( NV)^  /* Garbage gradient. */
  ){
      assert((v.nel) == NV);
      eFunction.defVar(v);
      eFunction.eval(c, e, TRUE, eDc);
      fprintf(stderr, " energy == " & FLR(e, 9, 4) & "\n");
      
      fprintf(stderr, "  vNum V              Gradient (algebrical)                         Gradient (numerical)                 Difference\n"); 
      fprintf(stderr, " ----- -   -------------------------------------------   -------------------------------------------   ----------------\n");
      /* Compute numerical gradient: */
      for (i = 0; i < NV; i++)
      {
        fprintf(stderr, " " & Fmt.Pad(Fmt.Int(i), 5));
        if (v[i])
        {
          fprintf(stderr, " T");
          fprintf(stderr, "   " & FLR4(eDc[i], 10,7));
          for (k = 0; k < 4; k++)
          {
            cOld = c[i][k];
            c[i][k] = cOld + DiffStep;
            eFunction.eval(c, ep, FALSE, eDcXXX);
            eDcNum[k] = (ep - e)/DiffStep;
            c[i][k] = cOld
          }
          fprintf(stderr, "   " & FLR4(eDcNum, 10,7));
          fprintf(stderr, "   " & ELR(r4_Dist(eDcNum, eDc[i]), 16,8));
        }
        else
        {
          fprintf(stderr, " F");
          fprintf(stderr, "   " & FLR4(eDc[i], 10,7));
        }
        fprintf(stderr, "\n")
      }
    }
  } /* END GradTest */

double EvalTest(
    pwr: Wr.T;                     /* File ".plot" */ 
    eFunction: MixedEnergy.T;      /* Energy function */
    Coords_t *c;                 /* Coordinates */
    Coords_t *cDs;          /* Direction vector */
    bool_vec_t *v;             /* Which nodes are variable */
    Coords_t *eDc;               /* Work area for gradient */
 )
  
  double *e, ep;
      double eDs;   /* Directional derivative */
  CONST DiffStep == 0.0001;
  { ??? NV = (c.nel);
    {
      assert((v.nel) == NV);
      eFunction.defVar(v);
      eFunction.eval(c, e, TRUE, eDc);
      fprintf(pwr, "  "); 
      fprintf(pwr, " " & FLR(e, 12, 7)); 
      
      /* Compute directional derivative from gradient: */
      eDs = 0.0;
      for (i = 0; i < NV; i++){ eDs = eDs + r4_dot(eDc[i], cDs[i]); }
      fprintf(pwr, " " & FLR(eDs, 12, 7));
      
      /* Check if setting "grad==FALSE" has any effect on "e": */
      eFunction.eval(c, ep, FALSE, eDc);
      assert(ep == e);
      
      /* Compute directional derivative numerically: */
      ??? ds = DiffStep;
    {
        for (i = 0; i < NV; i++){ c[i] = r4_Mix(1.0, c[i], ds, cDs[i]); }
        eFunction.eval(c, ep, FALSE, eDc);
        eDs = (ep - e)/ds
      }
      fprintf(pwr, " " & FLR(eDs, 12, 7));
    }
    return e;
  } /* END EvalTest */

void WriteHeader(pwr: Wr.T)
  
  { fprintf(pwr, "#"); 
    fprintf(pwr, " " & Fmt.Pad("s", 6));
    fprintf(pwr, "  "); 
    fprintf(pwr, " " & Fmt.Pad("etot", 12));
    fprintf(pwr, " " & Fmt.Pad("etotDc*cDs", 12));
    fprintf(pwr, " " & Fmt.Pad("etotDs", 12));
    fprintf(pwr, "  "); 
    fprintf(pwr, " " & Fmt.Pad("evar", 12));
    fprintf(pwr, " " & Fmt.Pad("evarDc*cDs", 12));
    fprintf(pwr, " " & Fmt.Pad("evarDs", 12));
    fprintf(pwr, "  "); 
    fprintf(pwr, " " & Fmt.Pad("etot-evar", 15));
    fprintf(pwr, "\n");

    fprintf(pwr, "#");
    fprintf(pwr, " " & Fmt.Pad("",  6, '-', Fmt.Align.Right));
    fprintf(pwr, "  "); 
    fprintf(pwr, " " & Fmt.Pad("", 12, '-', Fmt.Align.Right));
    fprintf(pwr, " " & Fmt.Pad("", 12, '-', Fmt.Align.Right));
    fprintf(pwr, " " & Fmt.Pad("", 12, '-', Fmt.Align.Right));
    fprintf(pwr, "  "); 
    fprintf(pwr, " " & Fmt.Pad("", 12, '-', Fmt.Align.Right));
    fprintf(pwr, " " & Fmt.Pad("", 12, '-', Fmt.Align.Right));
    fprintf(pwr, " " & Fmt.Pad("", 12, '-', Fmt.Align.Right));
    fprintf(pwr, "  "); 
    fprintf(pwr, " " & Fmt.Pad("", 15, '-', Fmt.Align.Right));
    fprintf(pwr, "\n");
  } /* END WriteHeader */

Options_t *GetOptions(int argc, char **argv)
  { 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, "-tpFile");                               
        o->tpFile = argparser_get_next(pp);                         

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

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

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

        if ((argparser_keyword_present(pp, "-nSteps")))
        {
          o->nSteps = argparser_get_next_int(pp, 1, 1000);
        }
        else
        {
          o->nSteps = 100;
        }

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

        o->showMin = argparser_keyword_present(pp, "-showMin");

        argparser_finish(pp);                                       
      
  -----------------------------------
#define _HELP \
  fprintf(stderr, "Usage: TestEnergy \\\n" \
   "  -tpFile NAME -aFile NAME -bFile NAME  \\\n" \
   "  -outFile <name> \\\n" \
   "  [ -nSteps <num> ] \\\n" \
   "  [ -showMin ] \\\n" \
   ParseEnergyParams.Help \
   "\n");
          
END¦                                              
      }
    }
    return o;
  } /* END GetOptions */

PROCEDURE FLR(x: double; wd, pr: uint): char *==
  { return Fmt.Pad(Fmt.LongReal(x, Fmt.Style.Fix, pr), wd);
  } /* END FLR */

PROCEDURE ELR(x: double; wd, pr: uint): char *==
  { return Fmt.Pad(Fmt.LongReal(x, Fmt.Style.Sci, pr), wd);
  } /* END ELR */

PROCEDURE FLR4(x: r4_t; wd, pr: uint): char *==
  { return 
      Fmt.Pad(Fmt.LongReal(x[0], Fmt.Style.Fix, pr), wd) & " " \
      Fmt.Pad(Fmt.LongReal(x[1], Fmt.Style.Fix, pr), wd) & " " \
      Fmt.Pad(Fmt.LongReal(x[2], Fmt.Style.Fix, pr), wd) & " " \
      Fmt.Pad(Fmt.LongReal(x[3], Fmt.Style.Fix, pr), wd)
  } /* END FLR4 */

/* end TestEnergy */