/* Defines a target function F for bbopt1.c */
/* Last edited on 2003-09-26 08:51:14 by stolfi */
  
#include <ia.h>
#include <js.h>
#include <values.h>
#include <bbgoal.h>

int F_dim = 1;
char *F_name = "F2-IA";
char *F_descr = "(IA) u = x - 1/7; v = x - 6/7; u^2(v^2 + 1/32) + 1/32";

#define ia_One (Interval){1.0, 1.0}

#define F_XUN (6.0)
#define F_XUD (7.0)

#define F_XVN (1.0)
#define F_XVD (7.0)

#define F_EPS (1.0/512.0)

Float F(Float *x) 
  { Float u = (*x) - F_XUN/F_XUD;
    Float v = (*x) - F_XVN/F_XVD;
    return u*u*(v*v + F_EPS) + F_EPS;
  }

Interval F_ia(Interval *xr) 
  { Interval x = *xr;
    Interval mu = ia_scale(ia_One, F_XUN, F_XUD); 
    Interval u = ia_sub(x, mu);
    Interval uu = ia_sqr(u);
    Interval mv = ia_scale(ia_One, F_XVN, F_XVD); 
    Interval v = ia_sub(x, mv);
    Interval vv = ia_sqr(v);
    Interval a = ia_shift(vv, F_EPS);
    Interval b = ia_mul(uu, a);
    Interval f = ia_shift(b, F_EPS);
    return f;
  }
  
/* 
  Function: 
    u^2(v^2 + EPS) + EPS
  Derivative:  
    2u(v^2 + EPS) + u^2(2 v) = 
    2u(v^2 + v u + EPS) = 
    2u(v^2 + v(v + v0 - u0)) + EPS) = 
    2u(2v^2 + v D + EPS)
  where u = x - u0, v = x - v0, D = v0 - u0. 
  Assuming v0 < u0:
  First minimum at    x = v0 + (sqrt(D^2 - 8 EPS) - D)/2
  Internal maximum at x = v0 - (sqrt(D^2 - 8 EPS) + D)/2
  Second minimum at   x = u0 */

void F_sol(Interval *xr, Interval *sr)
  { Interval C = ia_const(F_EPS, Zero); 
    Interval u0 = ia_scale(ia_One, F_XUN, F_XUD);
    Interval v0 = ia_scale(ia_One, F_XVN, F_XVD);
    Interval D = ia_sub(v0, u0);
    Interval sd = ia_sqrt(ia_sub(ia_sqr(D), ia_scale(C, 8.0, 1.0)));
    Interval xmin1 = ia_add(v0, ia_scale(ia_sub(sd, D), 1.0, 2.0)); 
    Interval xmax = ia_sub(v0, ia_scale(ia_add(sd, D), 1.0, 2.0));
    Interval xmin2 = u0;
    
    auto void do_sol(Interval *xr, Interval *sr);
    void do_sol(Interval *xr, Interval *sr)
      { fprintf(stderr, "  xr  = "); ia_print(stderr, *xr); fprintf(stderr, "\n");
        Float slo, shi;
        if (xr->lo > xr->hi)
          { slo = xr->lo; shi = xr->hi; }
        if (xr->lo >= xmin2.hi )
          { slo = shi = xr->lo; }
        if (xr->hi <= xmin1.lo )
          { slo = shi = xr->hi; }
        else if (xr->lo >= xmax.hi)
          { if (xr->hi >= xmin2.hi)
              { slo = xmin2.lo; shi = xmin2.hi; }
            else if (xr->hi >= xmin2.lo)
              { slo = xmin2.lo; shi = xr->hi; }
            else
              { slo = shi = xr->hi; }
          }
        else if (xr->hi <= xmax.lo)
          { if (xr->lo <= xmin1.lo)
              { slo = xmin1.lo; shi = xmin1.hi; }
            else if (xr->lo <= xmin1.hi)
              { slo = xr->lo; shi = xmin1.hi; }
            else
              { slo = shi = xr->lo; }
          }
        else
          { affirm(FALSE, "bad do_sol call"); }
        *sr = (Interval) { slo, shi };
      }
      
    fprintf(stderr, "u0    = "); ia_print(stderr, u0); fprintf(stderr, "\n");
    fprintf(stderr, "v0    = "); ia_print(stderr, v0); fprintf(stderr, "\n");
    fprintf(stderr, "sd    = "); ia_print(stderr, sd); fprintf(stderr, "\n");
    fprintf(stderr, "xmin1 = "); ia_print(stderr, xmin1); fprintf(stderr, "\n");
    fprintf(stderr, "xmax  = "); ia_print(stderr, xmax); fprintf(stderr, "\n");
    fprintf(stderr, "xmin2 = "); ia_print(stderr, xmin2); fprintf(stderr, "\n");
    if ((xr->lo >= xmax.hi) || (xr->hi <= xmax.lo))
      { do_sol(xr, sr); }
    else
      { /* The interval {xr} may span the maximum and hence inresect both basins. */
        /* Trisect {xr} at {xmax} and find an eclosure for the minimum in each part: */
        Interval sra, srm, srb, yra, yrm, yrb;
        Float yminhi = MAXDOUBLE; /* Upper bound for minimum in {xr}. */
        Interval xra = ia_meet((Interval){xr->lo, xmax.lo}, *xr);
        if (xra.lo <= xra.hi) 
          { do_sol(&xra, &sra); yra = F_ia(&sra); 
            if (yra.hi <= yminhi) { yminhi = yra.hi; } 
          }
        else
          { yra = (Interval){ 1.0, -1.0 }; }
        Interval xrm = ia_meet(xmax, *xr);
        if (xrm.lo <= xrm.hi) 
          { srm = xrm; yrm = F_ia(&srm); 
            if (yrm.hi <= yminhi) { yminhi = yrm.hi; }
          }
        else
          { yrm = (Interval){ 1.0, -1.0 };
          }
        Interval xrb = ia_meet((Interval){xmax.hi, xr->hi}, *xr);
        if (xrb.lo <= xrb.hi) 
          { do_sol(&xrb, &srb); yrb = F_ia(&srb); 
            if (yrb.hi <= yminhi) { yminhi = yrb.hi; }
          }
        else
          { yrb = (Interval){ 1.0, -1.0 }; }

        /* Join the intervals {sra, srm, srb} may contain minimum: */
        *sr = (Interval){ 1.0, -1.0 };
        if ((yra.hi >= yra.lo) && (yra.lo <= yminhi)) { *sr = ia_join(*sr, sra); }
        if ((yrm.hi >= yrm.lo) && (yrm.lo <= yminhi)) { *sr = ia_join(*sr, srm); }
        if ((yrb.hi >= yrb.lo) && (yrb.lo <= yminhi)) { *sr = ia_join(*sr, srb); }
      }
  }

Interval F_plot(Interval *xr)
  { Interval w0 = ia_scale(ia_One, F_XUN + F_XVN, F_XUD + F_XVD);
    Interval w = ia_sub(*xr, w0);
    return ia_shift(ia_sqr(ia_sqr(w)), F_EPS);
  }