/* Last edited on 2005-06-20 20:12:48 by stolfi */

----------------------------------------------------------------------
void fit_global_zmid_zrad
  ( int np, 
    int *ix,        
    double *Z,
    int Sdim, 
    double *Zmid, 
    double *Zrad
  )
  {
    /* fprintf(stderr, "warning: %s not yet implemented\n", __FUNCTION__);  */
    if (DEBUG) { fprintf(stderr, "%s: Sdim = %4d", __FUNCTION__, Sdim); }
    double Zmin = Z[ix[0]];
    double Zmax = Z[ix[np-1]];
    (*Zmid) = (Zmin + Zmax)/2;
    (*Zrad) = 1.05*(Zmax - Zmin)/2;
    if (DEBUG) { fprintf(stderr, " Zmid = %24.16e Zrad = %24.16e\n", (*Zmid), (*Zrad)); }
  }
----------------------------------------------------------------------
void fit_asymp_zref_coef
  ( int ns, 
    double *Z, 
    int np,
    double Expt,
    double *Zref,
    double *Coef
  )
  {
    /* We replace {T[k]} by {R[k] = T[k]**(1/Expt)}, and fit a linear
      model {R[k] ~ KZ*Z[k] + K1}.  Then we take {Coef=KZ**Expt} and 
      {Zref=-K1/KZ}. */

    /* Set up the normal system for the basis {Z,1} and the 
      scalar product { <f|g> = SUM{ W[i]*f(Z[i])g(Z[i]) : i = 0..nk-1 }}. */
    double SZZ = 0; double SZ1 = 0; double S11 = 0;
    double SRZ = 0; double SR1 = 0;
    int k;
    for (k = 0; k < ns; k++)
      { double Zk = Z[k]; 
        double r = ((double)k + 0.5)/((double) np);
        double Rk = exp(log(r)/Expt);
        double xk = ((double)k + 0.5)/((double)ns);
        double Wk = xk*(1.0 - xk);
        SZZ += Wk*Zk*Zk; SZ1 += Wk*Zk; S11 += Wk;
        SRZ += Wk*Rk*Zk; SR1 += Wk*Rk;
      }
    /* Solve normal system {((SZZ, SZ1),(SZ1,S11)) * (CZ,C1) = (SRZ,SR1)}: */
    double D = SZZ*S11 - SZ1*SZ1;
    double DZ = SRZ*S11 - SR1*SZ1;
    double D1 = SZZ*SR1 - SZ1*SRZ;
    double KZ = DZ/D;
    double K1 = D1/D;
    if (DEBUG) { fprintf(stderr, "KZ = %24.16e D1 = %24.16e\n", KZ, K1 ); }
    if (KZ < 0.0) { KZ = 0.001; }
    (*Coef) = exp(log(KZ)*Expt);
    (*Zref) = -K1/KZ;
  }
----------------------------------------------------------------------

    read_params(stdin, &nv, &d, &distr, &tour, &ms, &seed);
    /* Read tour costs : */
    double Z[ms];
    int ns = 0; /* Actual number of sample tours. */
    while (TRUE)
      {
        double w;
        bool_t ok = read_weight_skip_perm(stdin, &w, nv);
        if (! ok) { break; }
        demand(ns < ms, "too many input values"); 
        Z[ns] = w;
        if (ns > 0) { demand(Z[ns] >= Z[ns-1], "tours out of order"); }
        ns++;
      }


#! /usr/bin/gawk -f
# Last edited on 2004-12-27 23:20:49 by stolfi

BEGIN {
  abort = -1;
  
  
  if (nv == "") { arg_error("must define nv"); }
  nread = 0;
  
  split("", Z);
  
  # Global constants for minimizer:
  Eps = (1.0e-4);                   # Square root of relative machine precision.
  Phi = (1.61803398874989484821);   # Golden ratio.
  Ihp = (0.61803398874989484821);   # 1/Phi.
  IhpSqr = (0.3819660112501051518); # 1/Phi^2.
  
  # 
}

(abort >= 0) { exit abort; }

// { Z[FNR] = $1; nread = FNR; next; }

END {
  if (abort >= 0) { exit abort; }
  ns = nread;
}

void 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) { fprintf(stderr, "tol1 = %18.10f tol2 = %18.10f\n", tol1, tol2 ); }

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

    # check stopping criterion: 
    if (debug) { fprintf(stderr, "x-xa =  %18.10f xb-x = %18.10f\n", x-xa,xb-x ); }
      
    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) { fprintf(stderr, "parabolic step\n" ); }
        d = -(Dfx/DDfx)*(DDq/Dq);
      }
    else
      { # A golden-section interpolation: 
        if (debug) { fprintf(stderr, "golden section\n" ); }
        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) { fprintf(stderr, "u =  %18.10f\n", u ); }
    
    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) { fprintf(stderr, "u =  %18.10f fu = %18.10f\n", u, fu ); }
    
    # 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; }
  }
}

void 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;
}

void evalf(x)
{
  # Returns the total squared discrepancy assuming eponent {x}.
  # Warning: changes the global values of {Coef,Zref}.
  
  # Compute {Zref,Coef} for {Expt = x}: 
  fit_asymp_zref_coef(x);
  # Compute the goal function:
  return compute_discrepancy(Zref,Coef,x);
}

void print_status(xa,xb,v,fv,x,fx,w,fw)
{
  fprintf(stderr, "\n" );
  fprintf(stderr, "xa = %18.10f\n", xa );
  fprintf(stderr, "v =  %18.10f fv = %18.10f\n", v, fv );
  fprintf(stderr, "x =  %18.10f fx = %18.10f\n", x, fx );
  fprintf(stderr, "w =  %18.10f fw = %18.10f\n", w, fw );
  fprintf(stderr, "xb = %18.10f\n", xb );
  fprintf(stderr, "\n" );
}