hmmer-3.3/easel/esl_minimizer.c
/* bracket():
 * SRE, Wed Jul 27 11:43:32 2005 [St. Louis]
 *
 * Purpose:   Bracket a minimum. 
 *
 *            The minimization is quasi-one-dimensional, 
 *            starting from an initial <n>-dimension vector <ori>
 *            in the <n>-dimensional direction <d>.
 *            
 *            Caller passes a ptr to the objective function <*func()>,
 *            and a void pointer to any necessary conditional 
 *            parameters <prm>. The objective function will
 *            be evaluated at a point <x> by calling
 *            <(*func)(x, n, prm)>. The caller's function
 *            is responsible to casting <prm> to whatever it's
 *            supposed to be, which might be a ptr to a structure,
 *            for example; typically, for a parameter optimization
 *            problem, this holds the observed data.
 *            
 *            The routine works in scalar multipliers relative
 *            to origin <ori> and direction <d>; that is, a new <n>-dimensional
 *            point <b> is defined as <ori> + <bx><d>, for a scalar <bx>.
 *            
 *            The routine identifies a triplet <ax>, <bx>, <cx> such
 *            that $a < b < c$ and such that a minimum is known to
 *            exist in the $(a,b)$ interval because $f(b) < f(a),
 *            f(c)$. Also, the <a..b> and <b...c> intervals are in
 *            a golden ratio; the <b..c> interval is 1.618 times larger
 *            than <a..b>.
 *
 *            Since <d> is usually in the direction of the gradient,
 *            the points <ax>,<bx>,<cx> might be expected to be $\geq 0$;
 *            however, when <ori> is already close to the minimum, 
 *            it is often faster to bracket the minimum using
 *            a negative <ax>. The caller might then try to be "clever"
 *            and assume that the minimum is in the <bx..cx> interval
 *            when <ax> is negative, rather than the full <ax..cx>
 *            interval. That cleverness can fail, though, if <ori>
 *            is already in fact the minimum, because the line minimizer
 *            in brent() assumes a non-inclusive interval. Use
 *            <ax..cx> as the bracket.
 *            
 * Args:      cfg       - OPTIONAL configuration of tunable params, or NULL
 *            ori       - n-dimensional starting vector
 *            d         - n-dimensional direction to minimize along
 *            n         - # of dimensions
 *            firststep - bx is initialized to this scalar multiplier
 *            *func()   - objective function to minimize
 *            prm       - void * to any constant data that *func() needs
 *            wrk       - workspace: 1 allocated n-dimensional vector
 *            ret_ax    - RETURN:  ax < bx < cx scalar bracketing triplet
 *            ret_bx    - RETURN:    ...ax may be negative
 *            ret_cx    - RETURN:    
 *            ret_fa    - RETURN:  function evaluated at a,b,c
 *            ret_fb    - RETURN:    ... f(b) < f(a),f(c)
 *            ret_fc    - RETURN:
 *            opt_dat   - optRETURN: stats collection, if caller provided for it
 *
 * Returns:   <eslOK> on success.
 *
 * Throws:    <eslENORESULT> if it fails to converge.
 *
 * Xref:      STL9/130.
 */


Bracket a minimum.

The minimization is quasi-one-dimensional,
starting from an initial -dimension vector 
in the -dimensional direction .

Caller passes a ptr to the objective function <*func()>,
and a void pointer to any necessary conditional
parameters . The objective function will
be evaluated at a point  by calling
<(*func)(x, n, prm)>. The caller's function
is responsible to casting  to whatever it's
supposed to be, which might be a ptr to a structure,
for example; typically, for a parameter optimization
problem, this holds the observed data.

The routine works in scalar multipliers relative
to origin  and direction ; that is, a new -dimensional
point  is defined as  + , for a scalar .



Since  is usually in the direction of the gradient,
the points ,, might be expected to be $\geq 0$;
however, when  is already close to the minimum,
it is often faster to bracket the minimum using
a negative . The caller might then try to be "clever"
and assume that the minimum is in the  interval
when  is negative, rather than the full 
interval. That cleverness can fail, though, if 
is already in fact the minimum, because the line minimizer
in brent() assumes a non-inclusive interval. Use
 as the bracket.



 on success.

 if it fails to converge.


Bracket a minimum . The minimization is quasi - one - dimensional starting from an initial - dimension vector in the - dimensional direction . Caller passes a ptr to the objective function < * func () > and a void pointer to any necessary conditional parameters . The objective function will be evaluated at a point by calling < ( * func ) ( x n prm ) > . The caller ' s function is responsible to casting to whatever it ' s supposed to be which might be a ptr to a structure for example ; typically for a parameter optimization problem this holds the observed data . The routine works in scalar multipliers relative to origin and direction ; that is a new - dimensional point is defined as + for a scalar . Since is usually in the direction of the gradient the points might be expected to be $ \ geq 0$ ; however when is already close to the minimum it is often faster to bracket the minimum using a negative . The caller might then try to be " clever " and assume that the minimum is in the interval when is negative rather than the full interval . That cleverness can fail though if is already in fact the minimum because the line minimizer in brent () assumes a non - inclusive interval . Use as the bracket . on success . if it fails to converge .
static int
bracket(ESL_MIN_CFG *cfg, double *ori, double *d, int n, double firststep,
	double (*func)(double *, int, void *), void *prm, 
	double *wrk, 
	double *ret_ax, double *ret_bx, double *ret_cx,
	double *ret_fa, double *ret_fb, double *ret_fc, ESL_MIN_DAT *dat)
{
  int    max_iterations = cfg ? cfg->brack_maxiter : eslMIN_BRACK_MAXITER;
  double ax,bx,cx;		
  double fa,fb,fc;		
  double swapper;
  int    niter;
  ax = 0.;  
  fa = (*func)(ori, n, prm);
  bx = firststep;
  esl_vec_DCopy(ori, n, wrk);
  esl_vec_DAddScaled(wrk, d, bx, n);
  fb = (*func)(wrk, n, prm);
  if (fb > fa)
    {
      swapper = ax; ax = bx; bx = swapper;
      swapper = fa; fa = fb; fb = swapper;
    }
  cx = bx + (bx-ax)*1.618;
  esl_vec_DCopy(ori, n, wrk);
  esl_vec_DAddScaled(wrk, d, cx, n);
  fc = (*func)(wrk, n, prm);
  niter = 0;
  while (fc <= fb)
    {
      ax = bx; bx = cx;
      fa = fb; fb = fc;
      cx = bx+(bx-ax)*1.618;
      esl_vec_DCopy(ori, n, wrk);
      esl_vec_DAddScaled(wrk, d, cx, n);
      fc = (*func)(wrk, n, prm);
      if (ax != bx && bx != cx && fa == fb && fb == fc) break;
      niter++;
      if (niter > max_iterations) ESL_EXCEPTION(eslENORESULT, "Failed to bracket a minimum.");
    }
  if (ax > cx)
    {
      swapper = ax; ax = cx; cx = swapper;
      swapper = fa; fa = fc; fc = swapper;
    }
  if (dat)
    {
      dat->brack_n[dat->niter]  = niter;
      dat->nfunc[dat->niter]   += niter + 3;
      dat->brack_ax[dat->niter] = ax;
      dat->brack_bx[dat->niter] = bx;
      dat->brack_cx[dat->niter] = cx;
      dat->brack_fa[dat->niter] = fa;
      dat->brack_fb[dat->niter] = fb;
      dat->brack_fc[dat->niter] = fc;
    }
  ESL_DPRINTF2(("bracket(), in %2d iterations, brackets a minimum: a,b,c = %10.4g %10.4g %10.4g\n", niter, ax, bx, cx));
  ESL_DPRINTF2(("bracket():                                     fa,fb,fc = %10.4g %10.4g %10.4g\n",        fa, fb, fc));
  *ret_ax = ax;  *ret_bx = bx;  *ret_cx = cx;
  *ret_fa = fa;  *ret_fb = fb;  *ret_fc = fc;
  return eslOK;
}
static int bracket ( ESL_MIN_CFG * cfg , double * ori , double * d , int n , double firststep , double ( * func ) ( double * , int , void * ) , void * prm , double * wrk , double * ret_ax , double * ret_bx , double * ret_cx , double * ret_fa , double * ret_fb , double * ret_fc , ESL_MIN_DAT * dat ) { int max_iterations = cfg ? cfg -> brack_maxiter : eslMIN_BRACK_MAXITER ; double ax , bx , cx ; double fa , fb , fc ; double swapper ; int niter ; ax = 0. ; fa = ( * func ) ( ori , n , prm ) ; bx = firststep ; esl_vec_DCopy ( ori , n , wrk ) ; esl_vec_DAddScaled ( wrk , d , bx , n ) ; fb = ( * func ) ( wrk , n , prm ) ; if ( fb > fa ) { swapper = ax ; ax = bx ; bx = swapper ; swapper = fa ; fa = fb ; fb = swapper ; } cx = bx + ( bx - ax ) * 1.618 ; esl_vec_DCopy ( ori , n , wrk ) ; esl_vec_DAddScaled ( wrk , d , cx , n ) ; fc = ( * func ) ( wrk , n , prm ) ; niter = 0 ; while ( fc <= fb ) { ax = bx ; bx = cx ; fa = fb ; fb = fc ; cx = bx + ( bx - ax ) * 1.618 ; esl_vec_DCopy ( ori , n , wrk ) ; esl_vec_DAddScaled ( wrk , d , cx , n ) ; fc = ( * func ) ( wrk , n , prm ) ; if ( ax != bx && bx != cx && fa == fb && fb == fc ) break ; niter ++ ; if ( niter > max_iterations ) ESL_EXCEPTION ( eslENORESULT , " " ) ; } if ( ax > cx ) { swapper = ax ; ax = cx ; cx = swapper ; swapper = fa ; fa = fc ; fc = swapper ; } if ( dat ) { dat -> brack_n [ dat -> niter ] = niter ; dat -> nfunc [ dat -> niter ] += niter + 3 ; dat -> brack_ax [ dat -> niter ] = ax ; dat -> brack_bx [ dat -> niter ] = bx ; dat -> brack_cx [ dat -> niter ] = cx ; dat -> brack_fa [ dat -> niter ] = fa ; dat -> brack_fb [ dat -> niter ] = fb ; dat -> brack_fc [ dat -> niter ] = fc ; } ESL_DPRINTF2 ( ( " \n " , niter , ax , bx , cx ) ) ; ESL_DPRINTF2 ( ( " \n " , fa , fb , fc ) ) ; * ret_ax = ax ; * ret_bx = bx ; * ret_cx = cx ; * ret_fa = fa ; * ret_fb = fb ; * ret_fc = fc ; return eslOK ; }
bracket():
 SRE, Wed Jul 27 11:43:32 2005 [St. Louis]
bracket () : SRE Wed Jul 27 11 : 43 : 32 2005 [ St . Louis ]
/* scalar multipliers */ /* f() evaluations at those points */ /* Set and evaluate our first two points f(a) and f(b), which
   * are initially at 0.0 and <firststep>.
   */ /* always start w/ ax at the origin, ax=0 */ /* In principle, we usually know that the minimum m lies to the
   * right of a, m>=a, because d is likely to be a gradient.  You
   * might think we want 0 = a < b < c.  In practice, there's problems
   * with that. It's far easier to identify bad points (f(x) > f(a))
   * than to identify good points (f(x) < f(a)), because letting f(x)
   * blow up to infinity is fine as far as bracketing is concerned.
   * It can be almost as hard to identify a point b that f(b) < f(a)
   * as it is to find the minimum in the first place!
   * Counterintuitively, in cases where f(b)>f(a), it's better
   * to just swap the a,b labels and look for c on the wrong side
   * of a! This often works immediately, if f(a) was reasonably
   * close to the minimum and f(b) and f(c) are both terrible.
   */ /* Make our first guess at c.
   * Remember, we don't know that b>a any more, and c might go negative.
   * We'll either have:      a..b...c with a=0;
   *                or:  c...b..a     with b=0.
   * In many cases, we'll immediately be done.
   */ /* We're not satisfied until fb < fa, fc; 
   * throughout the routine, we guarantee that fb < fa;
   * so we just check fc.
   */ /* Slide over, discarding the a point; choose 
       * new c point even further away.
       */ /* This is a rare instance. We've reach the minimum
       * by trying to bracket it. Also check that not all
       * three points are the same.
       */ /* We're about to return. Assure the caller that the points
   * are in order a < b < c, not the other way.
   */ /* Optional verbosity (depending on compile-time eslDEBUGLEVEL) */ /* Return. */