/* Bivariate nonlinear branch-and-bound optimization with interval estimators. */
/* Last edited on 2003-09-21 16:49:55 by stolfi */

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

#define D 2

/* INTERNAL PROTOS */

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

void bb2_initialize_plot
  ( char *outName,
    bool epsfmt, 
    Interval xr, 
    Interval yr,
    FILE **psfP, 
    double *xsc, 
    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 factors {xsc} and {ysc} to be applied
    to the X and Y coordinates before plotting. The client plot window
    will contain the rectangle {{xsc*xr} × (ysc*yr)}. */

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

void bb2_print_box(Interval *xr, Interval *fr, bool final);
  /* A reporting function that prints the domain box {xr[0..D-1]} and 
    its image {*fr} to {stderr}. */

/* IMPLEMENTATIONS */

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

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

    void draw_box(Interval *xr, Interval fr, bool final)
      { bb2_print_box(xr, &fr, final);
        Interval x = ia_scale(xr[0], xscale, 1.0);
        Interval y = ia_scale(xr[1], 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[0..D-1]} and get its width {xw[0..D-1]}: */
    Interval xr[D];
    Float xw[D], tol[D];
    int i;
    for(i = 0; i < D; i++)
      { xr[i] = (Interval){0.0, 1.0};
        xw[i] = xr[i].hi - xr[i].lo;
        tol[i] = 0.01 * xw[i];
      }

    /* Get the approximate function value range {fr}: */
    /* Interval fr = F_ia(xr); */
    
    /* Get the true global minimum point {sr} and value {mr}: */
    Interval sr[D];
    F_sol(xr, sr);
    Interval mr = F_ia(sr);
    
    /* Generate EPS and PS output: */
    int epsfmt;
    for (epsfmt = 0; epsfmt < 2; epsfmt++)
      { 
        bb2_initialize_plot(outName, epsfmt, xr[0], xr[1], &psfile, &xscale, &yscale);
        
        /* Draw axes and plot of {F}: */
        /* (TO BE WRITTEN) */

        /* Plot the true minimum: */
        Float xsol = xscale * (sr[0].lo + sr[0].hi)/2;
        Float ysol = yscale * (sr[1].lo + sr[1].hi)/2;
        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(2, F_ia, &(xr[0]), &(tol[0]), &gr, draw_box); 
        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");
        bb2_print_box(sr, &mr, FALSE);

         /* Terminate plot: */
        bb2_finalize_plot(&psfile, epsfmt);
      }
    return 0;
  }

void bb2_initialize_plot
  ( char *outName,
    bool epsfmt, 
    Interval xr, 
    Interval yr,
    FILE **psfP, 
    double *xsc, 
    double *ysc
  )
  {
    FILE *psfile;
  
    /* 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.01*xw; xr.hi += 0.01*xw; xw = xr.hi - xr.lo;
    yr.lo -= 0.01*yw; yr.hi += 0.01*yw; yw = yr.hi - yr.lo;

    /* Actual figure dimensions: */
    double hw = 420.0, vw = 420.0;   /* In pt. */
    
    /* Compute {xscale,yscale} so that {(xscale*xw)/(yscale*yw) = hw/vw}: */
    double xscale, yscale;
    { double hyq = hw*yw;
      double vxq = vw*xw;
      if (hyq < vxq)
        { xscale = hyq/vxq; yscale = 1.0; }
      else
        { xscale = 1.0; yscale = vxq/hyq; }
    }

    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,
            xscale*xr.lo, xscale*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,
            xscale*xr.lo, xscale*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;
    *xsc = xscale;
    *ysc = yscale;
  }

void bb2_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 bb2_print_box(Interval *xr, Interval *fr, bool final)
  {
    ia_print(stderr, xr[0]); 
    fprintf(stderr, " × ");
    ia_print(stderr, xr[1]); 
    fprintf(stderr, " = ");
    ia_print(stderr, *fr); 
    if (final) { fprintf(stderr, " *"); }
    fprintf(stderr, "\n");
  }