#define PROG_NAME "makefigs"
#define PROG_DESC "generate figures for IA/AA papers"
#define PROG_VERS "1.0"

/* Copyright © 2003 by the State University of Campinas (UNICAMP). */
/* See the copyright, authorship, and warranty notice at end of file. */
/* Last edited on 2007-12-26 09:50:06 by stolfi */

/* Created jul/2003 by Jorge Stolfi, UNICAMP */

#include <aa.h>
#include <ia.h>
#include <flt.h>
#include <figprocs.h>
#include <aaiafuncs.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

/*** INTERNAL PROTOTYPES ***/

/* Command line options passed to the program: */
typedef struct Options
  { double scale;     /* Overall figure size factor */
    bool_t gPlot;       /* TRUE generates the {g} plots. */
    bool_t ggPlot;      /* TRUE generates the {gg} plots. */
    bool_t normPlot;    /* TRUE generates the normal affine approx plot. */
    bool_t chebPlot;    /* TRUE generates the Chebyshev approx plot. */
    bool_t xyrgPlot;    /* TRUE generates the XY range plot. */
    PlotOptions po;   /* Plot format and style options. */
  } Options;
   

int main(int argc, char **argv);
Options *parse_options(int argc, char **argv);
void get_arg_double(double *varp, int *argnp, int argc, char **argv, char *usage);
void get_arg_string(char **varp, int *argnp, int argc, char **argv, char *usage);
void arg_error(char *msg, char *arg, char *pname, char *usage);

void make_ia_ar_aa_plots(char *prefix);
void make_affine_approx_figure(char *prefix, Float alpha, Float gamma, Float delta);
void make_joint_range_figure(char *prefix);

AAP aa_from_4_coeffs
  ( Float x0, 
    Float x1, AAP e1,
    Float x2, AAP e2,
    Float x3, AAP e3,
    Float x4, AAP e4
  );
  /* Returns an affine form {x0 + x1*e1 + x2*e2 + x3*e3 + x4*e4},
    possibly with a few other roundoff error terms. */
   
RangeData make_joint_range
  ( char *tag,
    Float x0, Float x1, Float x2, Float x3, Float x4,
    Float y0, Float y1, Float y2, Float y3, Float y4
  );
  /* Creates two affine forms {x} and {y}, that use four new shared
    noise varianles, with given coefficients and tag. */
    
/* PARAMETERS FOR JOINT RANGE FIGURE */

Float xyrg_x0 = 20.0, xyrg_y0 = 10.0;
Float xyrg_x1 = -4.0, xyrg_y1 =  2.0;
Float xyrg_x2 =  0.0, xyrg_y2 =  1.0;
Float xyrg_x3 =  2.0, xyrg_y3 =  0.0;
Float xyrg_x4 =  3.0, xyrg_y4 = -1.0;

WinData xyrg_w = (WinData)
  { 
    /* xpt: */    {-3.0, 32.0},
    /* ypt: */    {-3.0, 17.0},
    /* scale: */  3.60 /* mm per client unit. */
  };

/* DEFINITIONS FOR BAD-IA FUNCTION GRAPHS */ 

WinData g_w = (WinData)
  { 
    /* xpt: */    {-2.3, 2.3},
    /* ypt: */    {-0.3, 2.3},
    /* scale: */  27.40,
  }; 
  
GraphData g_p = (GraphData)
  { 
    /* tag: */    "g",
    /* descr: */  "g(x) = sqrt(x^2 - x + 1/2)/sqrt(x^2 + 1/2)",
    /* xfp: */    {-2.1, 2.1},
    /* nfp: */    512,
    /* xia: */    {-2.0, 2.0},
    /* nia: */    16,
    /* fp: */     &g_fp,
    /* ia: */     &g_ia,
    /* aa: */     &g_aa
  };

/* DEFINITIONS FOR ITERATED BAD-IA FUNCTION GRAPH */ 

WinData gg_w = (WinData)
  { 
    /* xpt: */    {-2.3, 2.3},
    /* ypt: */    {-0.3, 2.3},
    /* scale: */  27.40,
  }; 
  
GraphData gg_p = (GraphData)
  { 
    /* tag: */    "gg",
    /* descr: */  "gg(x) = g(g(x))",
    /* xfp: */    {-2.1, 2.1},
    /* nfp: */    512,
    /* xia: */    {-2.0, 2.0},
    /* nia: */    16,
    /* fp: */     &gg_fp,
    /* ia: */     &gg_ia,
    /* aa: */     &gg_aa
  };

/* DEFINITIONS FOR NORMAL APPROXIMATION FUNCTION GRAPH */ 

WinData norm_w = (WinData)
  { 
    /* xpt: */    {-0.2, 2.2},
    /* ypt: */    {-0.2, 1.6},
    /* scale: */  52.50,
  };

GraphData norm_p = (GraphData)
  { 
    /* tag: */    "sr",
    /* descr: */  "sr(x) = sqrt(x)",
    /* xfp: */    {-0.1, 2.1},
    /* nfp: */    512,
    /* xia: */    {-0.0, 2.0},
    /* nia: */    8,
    /* fp: */     &sqrt_fp,
    /* ia: */     &sqrt_ia,
    /* aa: */     &sqrt_aa
  };

/* Normal approximation: */
ApproxData norm_ap = (ApproxData)
  { 
    /* tag: */    "nrm",
    /* xap: */    {0.25, 1.75},
    /* alpha: */  0.50,
    /* gamma: */  0.40,
    /* delta: */  0.15
  };

/* PARAMETERS FOR CHEBYSEHV APPROX FIGURE */

WinData cheb_w = (WinData)
  { 
    /* xpt: */    {-0.2, 2.2},
    /* ypt: */    {-0.2, 1.6},
    /* scale: */  52.50,
  };

GraphData cheb_p = (GraphData)
  { 
    /* tag: */    "sr",
    /* descr: */  "sr(x) = sqrt(x)",
    /* xfp: */    {-0.1, 2.1},
    /* nfp: */    512,
    /* xia: */    {-0.0, 2.0},
    /* nia: */    8,
    /* fp: */     &sqrt_fp,
    /* ia: */     &sqrt_ia,
    /* aa: */     &sqrt_aa
  };

/* Chebyshev approximation: */
ApproxData cheb_ap = (ApproxData)
  { 
    /* tag: */    "chb",
    /* xap: */    {0.25, 1.75},
    /* alpha: */  0.54858377035486353,
    /* gamma: */  0.40928648564717897,
    /* delta: */  0.04643242823589485
  };

/* MAIN PROGRAM */

int main(int argc, char **argv)
  { 
    Options *o = parse_options(argc, argv);
    PlotOptions *po = &(o->po);
    aa_init();
    
    if (o->gPlot)
      { /* Overestimation ({g}) plots: */
        g_w.ypt.lo = -0.24/o->scale;
        g_w.scale *= o->scale;
        figprocs_plot_ia_ar_aa_graphs(po, &g_w, &g_p,  "Example of correlated operands");
      }

    if (o->ggPlot)
      { /* Error explosion ({gg}) plots: */
        gg_w.ypt.lo = -0.24/o->scale;
        gg_w.scale *= o->scale;
        figprocs_plot_ia_ar_aa_graphs(po, &gg_w, &gg_p, "Example of error explosion");
      }

    if (o->normPlot)
      { /* Normal affine approximation figures: */
        
        /* Fix ranges to leave enough space for tick labels: */
        double h = (norm_w.xpt.hi - norm_w.xpt.lo)*norm_w.scale;
        norm_w.xpt.lo = -0.20/o->scale;
        norm_w.scale = h/(norm_w.xpt.hi - norm_w.xpt.lo);
        norm_w.ypt.lo = -0.15/o->scale;

        norm_w.scale *= o->scale;
        figprocs_plot_approx(po, &norm_w, &norm_p, &norm_ap, "A normal approximation");
      }

    if (o->chebPlot)
      { /* Chebyshev approximation: */
      
        /* Fix ranges to leave enough space for tick labels: */
        double h = (cheb_w.xpt.hi - cheb_w.xpt.lo)*cheb_w.scale;
        cheb_w.xpt.lo = -0.20/o->scale;
        cheb_w.scale = h/(cheb_w.xpt.hi - cheb_w.xpt.lo);
        cheb_w.scale *= o->scale;

        cheb_w.ypt.lo = -0.15/o->scale;
        figprocs_plot_approx(po, &cheb_w, &cheb_p, &cheb_ap, "Chebishev approximation");
      }

    /* Joint range figures: */
    if (o->xyrgPlot)
      { xyrg_w.scale *= o->scale;
        RangeData xyrg_r =
          make_joint_range
            ( "xyrg",
              xyrg_x0, xyrg_x1, xyrg_x2, xyrg_x3, xyrg_x4,
              xyrg_y0, xyrg_y1, xyrg_y2, xyrg_y3, xyrg_y4
            );
        figprocs_plot_joint_range(po, &xyrg_w, &xyrg_r, "Example of joint range");
      }
    return 0;
  }

/*** IMPLEMENTATIONS ***/

AAP aa_from_4_coeffs
  ( Float x0, 
    Float x1, AAP e1,
    Float x2, AAP e2,
    Float x3, AAP e3,
    Float x4, AAP e4
  )
  { 
    AAP x = aa_const(x0, Zero);
    x = aa_add(x, aa_scale(e1, x1, One));
    x = aa_add(x, aa_scale(e2, x2, One));
    x = aa_add(x, aa_scale(e3, x3, One));
    x = aa_add(x, aa_scale(e4, x4, One));
    return x;
  }


RangeData make_joint_range
  (
   char *tag,
   Float x0, Float x1, Float x2, Float x3, Float x4,
   Float y0, Float y1, Float y2, Float y3, Float y4
  )
  {
    RangeData jr;
    AAP e1 = aa_from_interval((Interval){-One,One});
    AAP e2 = aa_from_interval((Interval){-One,One});
    AAP e3 = aa_from_interval((Interval){-One,One});
    AAP e4 = aa_from_interval((Interval){-One,One});
    jr.tag = tag;
    jr.x = aa_from_4_coeffs(x0, x1, e1, x2, e2, x3, e3, x4, e4);
    jr.y = aa_from_4_coeffs(y0, y1, e1, y2, e2, y3, e3, y4, e4);
    jr.z = aa_const(Zero,Zero);
    return jr;
  }

#define ARG_ERROR(Msg,Arg) arg_error((Msg),(Arg),argv[0],usage)
  
#define GET_STRING(Var) get_arg_string(&(Var), &argn, argc, argv, usage)
#define GET_DOUBLE(Var) get_arg_double(&(Var), &argn, argc, argv, usage)

Options *parse_options(int argc, char **argv)
  {
    Options *o = (Options *)malloc(sizeof(Options));
    char* usage = "\n  [ -help ] [ -outName PREFIX ] [ -eps | -ps ]";
    int argn;

    /* Defaults: */
    o->po.outName = "plt";
    o->po.eps = FALSE;
    o->po.fontSize = 12.0;

    o->scale = 1.0;
    o->gPlot = FALSE;
    o->ggPlot = FALSE;
    o->normPlot = FALSE;
    o->chebPlot = FALSE;
    o->xyrgPlot = FALSE;

    argn = 1;

    /* Scan command line options. */
    while ((argn < argc) && (argv[argn][0] == '-') && (argv[argn][1] != '\0'))
      {
        char *key = argv[argn];
        if ((key[0] == '-') && (key[1] == '-') && (key[2] != '\0')) { key++; }
        if (strcmp(key, "-help") == 0)
          { fprintf(stderr, "usage: %s %s\n", argv[0], usage); exit(0); }
        else if (strcmp(key, "-outName") == 0)
          { GET_STRING(o->po.outName); }
        else if (strcmp(key, "-fontSize") == 0)
          { GET_DOUBLE(o->po.fontSize); }
        else if (strcmp(key, "-eps") == 0)
          { o->po.eps = TRUE; }
        else if (strcmp(key, "-ps") == 0)
          { o->po.eps = FALSE; }
        else if (strcmp(key, "-color") == 0)
          { o->po.color = TRUE; }
        else if (strcmp(key, "-gray") == 0)
          { o->po.color = FALSE; }
        else if (strcmp(key, "-grey") == 0)
          { o->po.color = FALSE; }
        else if (strcmp(key, "-scale") == 0)
          { GET_DOUBLE(o->scale); }
        else if (strcmp(key, "-gPlot") == 0)
          { o->gPlot = TRUE; }
        else if (strcmp(key, "-ggPlot") == 0)
          { o->ggPlot = TRUE; }
        else if (strcmp(key, "-normPlot") == 0)
          { o->normPlot = TRUE; }
        else if (strcmp(key, "-chebPlot") == 0)
          { o->chebPlot = TRUE; }
        else if (strcmp(key, "-xyrgPlot") == 0)
          { o->xyrgPlot = TRUE; }
        else  
          { ARG_ERROR("unknown option", argv[argn]); }
        ++argn;
      }
    if (argn != argc) { ARG_ERROR("extraneous arguments", argv[argn]); }
    return o;
  }

void get_arg_string(char **varp, int *argnp, int argc, char **argv, char *usage)
   /*
     Stores the next command line argument (as a string) into "*varp" */
  {
    int argn = *argnp;
    if (argn+1 >= argc)
      { ARG_ERROR("missing arg value", argv[argn]); }
    (*varp) = argv[argn+1];
    (*argnp) = argn+1;
  }
  
void get_arg_double(double *varp, int *argnp, int argc, char **argv, char *usage)
   /*
     Stores the next command line argument (as a double) into "*varp" */
  {
    int argn = *argnp;
    char *end;
    if (argn+1 >= argc)
      { ARG_ERROR("missing arg value", argv[argn]); }
    (*varp) = strtod(argv[argn+1], &end);
    if ((*end) != '\0') 
      { ARG_ERROR("invalid numeric argument", argv[argn+1]); }
    (*argnp) = argn+1;
  }

void arg_error(char *msg, char *arg, char *pname, char *usage)
   /*
     Prints "msg", "arg", the "usage" message, and exits.
     Handy while parsing the command line arguments. */
  {
    fprintf(stderr, "%s %s\n", msg, arg);
    fprintf(stderr, "usage: %s %s\n", pname, usage);
    exit(1);
  }