# Univariate minimization routine (Brent's method)
# Last edited on 2004-12-28 03:06:28 by stolfi

function evalf(x,  S)
{
  ???
}

function brent_min( \
  xini,tol,dist,xa,xb, \
  u,fu,v,fv,w,fw, \
  d,e,r,Dfx,Dq,DDfx,DDq,fxOld,tol1,tol2,tol3 \
)
{
  # Looks for the minimum point of the globaly defined function 
  # {evalf}.  Requires a guess {xini}, a domain interval 
  # {xa,xb}, the desired accuracy of the answer {tol},
  # and an upper bound for the distamce between {xini} and the 
  # true minimum point. Sets the global variables {x,fx,dfx} 
  # to the solution.
  
  # v, fv: Third-smallest point of current triple
  # w, fw: Second-smallest point of current triple
  # u, fu: Splitting point

  # d, e, r:      Step sizes 
  # Dfx,Dq:       Estimated derivative at {x} is {Dfx/Dq}. 
  # DDfx,DDq:     Estimated second derivative at {x} is {DDfx/DDq}. 
  # fxOld:        Previous value of {fx} 
  
  tol3 = tol/3.0; 

  if ((xini < xa) || (xini > xb)) { prog_error("bad xini"); }
  
  x = xini; fx = evalf(x);
  v = x; fv = fx;
  w = x; fw = fx;

  e = 0.0;
  fxOld = fw;
  while (1)
  {
    if ((fw < fx) || (fw > fv)) { prog_error("bad fw"); }
    tol1 = Eps * abs(x) + tol3;
    tol2 = 2.000001 * tol1;
    if (debug) { printf "tol1 = %18.10f tol2 = %18.10f\n", tol1, tol2 > "/dev/stderr"; }

    if (debug) 
      { print_status(xa,xb,v,fv,x,fx,w,fw); }

    # check stopping criterion: 
    if (debug) { printf "x-xa =  %18.10f xb-x = %18.10f\n", x-xa,xb-x > "/dev/stderr"; }
      
    if ((x - xa <= tol2) && (xb - x <= tol2)) { return; }

    if (abs(e) > tol1)
      { # Fit parabola through {v}, {x}, {w}, sets {Dfx,Dq, DDfx,DDq}. 
        compute_derivatives(v, fv, w, fw, x, fx);
        r = e; e = d;
      }
    else
      { Dfx = 0.0; Dq = 0; DDfx = 0.0; DDq = 0.0;
        r = 0.0;
      }

    # Compute the displacement {d} from {x} to the next probe: 
    if (DDfx < 0.0) { Dfx = -Dfx; DDfx = -DDfx; }
    if ( \
       (abs(Dfx)*DDq < 0.5*abs(DDfx*r)*Dq) && \
       (Dfx*DDq < DDfx*(x - xa)*Dq) && \
       (Dfx*DDq > DDfx*(x - xb)*Dq) \
    ) 
      { # A parabolic-interpolation step: 
        if (debug) { printf "parabolic step\n" > "/dev/stderr"; }
        d = -(Dfx/DDfx)*(DDq/Dq);
      }
    else
      { # A golden-section interpolation: 
        if (debug) { printf "golden section\n" > "/dev/stderr"; }
        e = (xb - x > x - xa ? xb - x : xa - x);
        d = IhpSqr * e;
      }

    # The function must not be evaluated too close to {xa} or {xb}: 
    u = x + d;
    
    if (debug) { printf "u =  %18.10f\n", u > "/dev/stderr"; }
    
    if (u - xa < tol1)
      { u = xa + tol1; prte("A"); }
    else if (xb - u < tol1)
      { u = xb - tol1; prte("B"); }

    # The function must not be evaluated too close to {x}: 
    if (abs(u - x) >= tol1)
      { prte("C"); }
    else if (u > x)
      { u = x + tol1; prte("D"); }
    else
      { u = x - tol1; prte("E"); }  

    fu = evalf(u);

    if (debug) { printf "u =  %18.10f fu = %18.10f\n", u, fu > "/dev/stderr"; }
    
    # update  {xa,xb}: 
    if (fu >= fx)
      { if (u < x) { xa = u; } else { xb = u; } }
    if (fu <= fx)
      { if (u < x) { xb = x; } else{ xa = x; } }

    # Update {x,v,w}: 
    if (u < xa) { prog_error("u < xa"); }
    if (u > xb) { prog_error("u > xb"); }
    if (fu <= fx)
      { if ((v != w) && (w != x)) { v = w; fv = fw; }
        w = x; fw = fx;
        x = u; fx = fu;
      }
    else if (fu <= fw)
      { if (w == x) { prog_error("w == x"); } 
        v = w; fv = fw;
        w = u; fw = fu;
      }
    else if (fu <= fv)
      { if (v == x) { prog_error("v == x"); } 
        if (w == x)
          { w = u; fw = fu; }
        else
          { v = u; fv = fu; }
      }
    else if (v == w)
      { v = u; fv= fu; }
  }
}

function compute_derivatives ( \
  u,fu,v,fv,x,fx,
  uv, ux,xv,du,dv,q \
)
{
  # Computes estimates {Dfx/Dq} and {DDfx/DDq} for the first and second
  # derivatives of the goal function {F} at {x}, by fitting 
  # a quadratic through the three points {(u, fu)}, {(v, fv)},
  # and {(x, fx)}.  

  # The scaling factor {q} will be positive if all three arguments
  # {u}, {v}, {w} are sufficiently distinct, and zero otherwise.
  
  uv = v - u;
  ux = x - u;
  xv = v - x;
  du = xv * (fx - fu);
  dv = ux * (fv - fx);
  q = uv * ux * xv;
  
  Dfx = du*xv + dv*ux;
  DDfx = 2.0 * (dv - du);
  if (q < 0.0) { Dfx = -Dfx; DDfx = -DDfx; q = -q; }
  Dq = q;
  DDq = q;
}

function print_status(xa,xb,v,fv,x,fx,w,fw)
{
  printf "\n" > "/dev/stderr";
  printf "xa = %18.10f\n", xa > "/dev/stderr";
  printf "v =  %18.10f fv = %18.10f\n", v, fv > "/dev/stderr";
  printf "x =  %18.10f fx = %18.10f\n", x, fx > "/dev/stderr";
  printf "w =  %18.10f fw = %18.10f\n", w, fw > "/dev/stderr";
  printf "xb = %18.10f\n", xb > "/dev/stderr";
  printf "\n" > "/dev/stderr";
}

function prte(tag)
{
  if (debug) { printf "[%s]\n", tag > "/dev/stderr"; } 
}

function abs(x)
{
  return (x < 0 ? -x : x);
}

function prog_error(msg)
{
  printf ** %s\n", msg > "/dev/stderr";
  abort = 1;
  exit abort;
}