/* Univariate nonlinear branch-and-bound optimization with interval estimators. */
/* Last edited on 2003-09-25 20:22:36 by stolfi */

#include <fbox.h>
#include <fboxheap.h>
#include <fboxlist.h>
#include <aa.h>
#include <ia.h>
#include <js.h>
#include <ps.h>
#include <fltgraph.h>
#include <stdlib.h>
#include <stdio.h>
#include <bbopt.h>
#include <bbgoal.h>

/* INTERNAL PROTOS */

int main(int argc, char **argv);

void bb1_initialize_plot
  ( char *outName,
    bool epsfmt, 
    Interval xr, 
    Interval fr,
    FILE **psfP, 
    double *ysc
  );
  /* Opens the output plot file and initializes scales, captions, etc.
    If {epsfmt} is true the file is encapsulated Postscript
    ("{outname}.eps") otherwise it is a plain Postscripr
    ("{outname}.ps"). The file is returned in {*psfP}.
    
    The procedure also computes a scale factor {ysc} to be applied to
    the function values before plotting. The client plot window will
    contain the rectangle {xr × (ysc*fr)}. */

void bb1_finalize_plot(FILE **psfP, bool epsfmt);
   /* Finalizes the plot and closes the plot file {*psfP}, whose
     type is indicated by {epsf}. */

void bb1_print_interval(Interval *xr, Interval *fr, bool final);
  /* A reporting function that prints the domain interval {*xr} and 
    its image {*fr} to {stderr}. */

/* IMPLEMENTATIONS */

int main(int argc, char **argv)
  { affirm(F_dim == 1, "Goal function has wrong domain dim");
    char *outName = txtcat("bbopt1-", F_name);
    FILE *psfile = NULL;
        
    double yscale;

    auto void draw_interval(Interval *xr, Interval fr, bool final);

    void draw_interval(Interval *xr, Interval fr, bool final)
      { bb1_print_interval(xr, &fr, final);
        Interval x = *xr, y = ia_scale(fr, yscale, 1.0);
        if (final)
          { ps_fill_and_draw_rectangle
              (psfile, x.lo, x.hi, y.lo, y.hi, 0.5,0.5,0.5);
          }
        else
          { ps_draw_rectangle(psfile, x.lo, x.hi, y.lo, y.hi); }
      }

    ia_init();
    
    /* Choose the initial domain {xr} and get its width {xw}: */
    Interval xr = (Interval){0.0, 1.0};
    Float xw = xr.hi - xr.lo;
    Float tol = 0.01 * xw;
    
    /* Get the approximate function value range {fr}: */
    Interval fr = F_ia(&xr);

    /* Get the true global minimum point {sr} and value {mr}: */
    Interval sr;
    F_sol(&xr, &sr);
    Interval mr = F_ia(&sr);

    /* Generate EPS and PS output: */
    int epsfmt;
    for (epsfmt = 0; epsfmt < 2; epsfmt++)
      { 
        bb1_initialize_plot(outName, epsfmt, xr, fr, &psfile, &yscale);
        
        /* Draw axes and plot of {F}: */
        auto Float F_plot(Float x);
        Float F_plot(Float x) { return yscale * F(&x); }
        Interval yr = ia_scale(fr, yscale, 1.0);
        ps_set_pen(psfile, 0.0,0.2,0.5, 0.25, 0.0,0.0);
        fltgraph_plot(psfile, F_plot, xr, yr, 100);
        
        /* Plot the true minimum: */
        Float xsol = (sr.lo + sr.hi)/2;
        Float ysol = F_plot(xsol);
        ps_set_pen(psfile, 0.5,0.0,0.0, 0.25, 0.0,0.0);
        ps_draw_dot(psfile, xsol, ysol, 2.0);
        
        /* Optimize and plot node boxes: */
        ps_set_pen(psfile, 0.0,0.5,0.0, 0.10, 0.0,0.0);
        Interval gr;
        FBoxList R = bb_optimize(1, F_ia, &xr, &tol, &gr, draw_interval); 
        affirm(R != NULL, "bb_optimize returned no boxes");
        fprintf(stderr, "minimum found = [%8.4f _ %8.4f]\n", gr.lo, gr.hi);
        
        /* Print the true minimum: */
        fprintf(stderr, "true minimum:\n");
        bb1_print_interval(&sr, &mr, FALSE);

        /* Draw frame: */
        ps_set_pen(psfile, 0.0,0.0,0.0, 0.25, 0.0,0.0);
        ps_draw_frame(psfile);
        
        /* Terminate plot: */
        bb1_finalize_plot(&psfile, epsfmt);
      }
    return 0;
  }

void bb1_initialize_plot
  ( char *outName,
    bool epsfmt, 
    Interval xr, 
    Interval fr,
    FILE **psfP, 
    double *ysc
  )
  { FILE *psfile;
  
    Interval yr = F_plot(&xr);

    /* Ensure that the vertical plot range {yr} includes the Y=0 axis: */
    if (yr.lo > 0.0) { yr.lo = 0.0; }
    if (yr.hi < 0.0) { yr.hi = 0.0; }

    /* Compute widths of {xr} and {yr}: */
    double xw = xr.hi - xr.lo;
    double yw = yr.hi - yr.lo;
    
    /* Add a safety margin around {xr,yr}: */
    xr.lo -= 0.03*xw; xr.hi += 0.03*xw; xw = xr.hi - xr.lo;
    yr.lo -= 0.03*yw; yr.hi += 0.03*yw; yw = yr.hi - yr.lo;

    /* Actual figure dimensions: */
    double hw = 420.0, vw = 420.0;   /* In pt. */
    
    /* Compute the scale factor {yscale} for function values: */ 
    double yscale = (xw*vw)/(yw*hw); /* So that {(yscale*yw)/xw = vw/hw}. */

    if (epsfmt)
      { psfile = open_write(txtcat(outName, ".eps"));
        double mrg = 4.0;
        ps_begin_figure(psfile, 0.0, hw + 2*mrg, 0.0, vw + 2*mrg);
        ps_set_window
          ( psfile,
            xr.lo, xr.hi,
            yscale*yr.lo, yscale*yr.hi,
            mrg, hw + mrg,  mrg, hw + mrg,
            1, 1
          );
        ps_set_pen(psfile, 0.0,0.0,0.0, 0.10, 0.0,0.0);
      }
    else
      { double hctr = ( 8.5 * 72.0)/2.0;
        double vctr = (11.0 * 72.0)/2.0;
        psfile = open_write(txtcat(outName, ".ps"));
        ps_begin_document(psfile, "letter");
        ps_begin_page(psfile, "1");
        ps_set_window
          ( psfile,
            xr.lo, xr.hi,
            yscale*yr.lo, yscale*yr.hi,
            hctr - hw/2, hctr + hw/2, vctr - vw/2, vctr + vw/2,
            1, 1
          );
        ps_set_pen(psfile, 0.0,0.0,0.0, 0.10, 0.0,0.0);
        ps_add_caption(psfile, F_descr);
      }
    *psfP = psfile;
    *ysc = yscale;
  }

void bb1_finalize_plot(FILE **psfP, bool epsfmt)  
  {
    FILE *psfile = *psfP;
    if (epsfmt)
      { ps_end_figure(psfile);
        fclose(psfile);
      }
    else
      { ps_end_page(psfile); 
        ps_end_document(psfile);
        fclose(psfile);
      }
  }

void bb1_print_interval(Interval *xr, Interval *fr, bool final)
  { ia_print(stderr, *xr); 
    fprintf(stderr, " × ");
    ia_print(stderr, *fr); 
    if (final) { fprintf(stderr, " *"); }
    fprintf(stderr, "\n");
  }