/*$Id: morethuente.c,v 1.1 2002/08/11 01:15:06 benson Exp $*/

#include "morethuente.h"  /*I "tao_solver.h
" I*/

static int TaoSetOptions_BoundLineSearch(TAO_SOLVER,void*);
static int TaoApply_BoundLineSearch(TAO_SOLVER,TaoVec*,TaoVec*,TaoVec*,TaoVec*,
			       double*,double*,int*,void*);

/*------------------------------------------------------------*/
#undef __FUNCT__  
#define __FUNCT__ "TaoSetOptions_LineSearch"
static int TaoSetOptions_BoundLineSearch(TAO_SOLVER tao, void*linectx)
{
  TAO_LINESEARCH *ctx = (TAO_LINESEARCH *)linectx;
  int            info;

  TaoFunctionBegin;
  info = TaoOptionsHead("Modified More-Thuente line line search options for bound constrained minimization");CHKERRQ(info);
  info = TaoOptionsTail();CHKERRQ(info);
  info = TaoSetOptions_LineSearch(tao,ctx);
  TaoFunctionReturn(0);
}
/*------------------------------------------------------------*/

/* ---------------------------------------------------------- */
#undef __FUNCT__  
#define __FUNCT__ "TaoApply_BoundLineSearch"
/* @ 
   TaoApply_BoundLineSearch - This routine performs a line search algorithm.

   Input Parameters:
+  tao - TAO_SOLVER context
.  X - current iterate (on output X contains new iterate, X + step*S)
.  S - search direction
.  f - objective function evaluated at X
.  G - gradient evaluated at X
.  W - work vector
-  step - initial estimate of step length

   Output parameters:
+  f - objective function evaluated at new iterate, X + step*S
.  G - gradient evaluated at new iterate, X + step*S
.  X - new iterate
.  gnorm - 2-norm of G
-  step - final step length


   Info is set to one of:
+   0 - the line search succeeds; the sufficient decrease
   condition and the directional derivative condition hold

   negative number if an input parameter is invalid
.   -1 -  step < 0 
.   -2 -  ftol < 0 
.   -3 -  rtol < 0 
.   -4 -  gtol < 0 
.   -5 -  stepmin < 0 
.   -6 -  stepmax < stepmin
-   -7 -  maxfev < 0

   positive number > 1 if the line search otherwise terminates
+    2 -  Relative width of the interval of uncertainty is 
         at most rtol.
.    3 -  Maximum number of function evaluations (maxfev) has 
         been reached.
.    4 -  Step is at the lower bound, stepmin.
.    5 -  Step is at the upper bound, stepmax.
.    6 -  Rounding errors may prevent further progress. 
         There may not be a step that satisfies the 
         sufficient decrease and curvature conditions.  
         Tolerances may be too small.
-    7 -  Search direction is not a descent direction.

   Notes:
   This algorithm is is a modification of the algorithm by More' and Thuente.
   The modifications concern bounds.  This algorithm steps in the direction
   passed into this routine.  This point get projected back into the feasible set.
   In the context of bound constrained
   optimization, there may not be a point in the piecewise linear 
   path that satisfies the Wolfe conditions.  When the active set
   is changing, decrease in the objective function may be sufficient
   to terminate this line search.

   Note: Much of this coded is identical to the More' Thuente line search.
   Variations to the code are commented.

   Notes:
   This routine is used within the following TAO boiund constrained minimization 
   solvers: Newton linesearch (tao_tron) and limited memory variable metric
   (tao_blmvm).

   Level: advanced

.keywords: TAO_SOLVER, linesearch
@ */
static int TaoApply_BoundLineSearch(TAO_SOLVER tao,TaoVec* X,TaoVec* G,TaoVec* S,TaoVec* W,double *f,
                        double *step,int *info2,void*ctx)
{
  TAO_LINESEARCH *neP = (TAO_LINESEARCH *) tao->linectx;
  TaoVec *XL,*XU;
  double    zero = 0.0, two = 2.0, p5 = 0.5, p66 = 0.66, xtrapf = 4.0;
  double    finit, width, width1, dginit,fm, fxm, fym, dgm, dgxm, dgym;
  double    dgx, dgy, dg, fx, fy, stx, sty, dgtest, ftest1=0.0,ftest2=0.0;
  double    bstepmin1,bstepmin2,bstepmax;
  double    dg1,dg2;
  int       info, i, stage1;

#if defined(PETSC_USE_COMPLEX)
  PetscScalar    cdginit, cdg, cstep = 0.0;
#endif

  TaoFunctionBegin;
  /* neP->stepmin - lower bound for step */
  /* neP->stepmax - upper bound for step */
  /* neP->rtol 	  - relative tolerance for an acceptable step */
  /* neP->ftol 	  - tolerance for sufficient decrease condition */
  /* neP->gtol 	  - tolerance for curvature condition */
  /* neP->nfev 	  - number of function evaluations */
  /* neP->maxfev  - maximum number of function evaluations */

  /* Check input parameters for errors */
  if (*step < zero) {
    info = PetscLogInfo((tao,"TaoApply_BoundLineSearch:Line search error: step (%g) < 0\n",*step)); CHKERRQ(info);
    *info2 = -1; TaoFunctionReturn(0);
  } else if (neP->ftol < zero) {
    info = PetscLogInfo((tao,"TaoApply_BoundLineSearch:Line search error: ftol (%g) < 0\n",neP->ftol)); CHKERRQ(info);
    *info2 = -2; TaoFunctionReturn(0);
  } else if (neP->rtol < zero) {
    info = PetscLogInfo((tao,"TaoApply_BoundLineSearch:Line search error: rtol (%g) < 0\n",neP->rtol)); CHKERRQ(info);
    *info2 = -3; TaoFunctionReturn(0);
  } else if (neP->gtol < zero) {
    info = PetscLogInfo((tao,"TaoApply_BoundLineSearch:Line search error: gtol (%g) < 0\n",neP->gtol)); CHKERRQ(info);
    *info2 = -4; TaoFunctionReturn(0);
  } else if (neP->stepmin < zero) {
    info = PetscLogInfo((tao,"TaoApply_BoundLineSearch:Line search error: stepmin (%g) < 0\n",neP->stepmin)); CHKERRQ(info);
    *info2 = -5; TaoFunctionReturn(0);
  } else if (neP->stepmax < neP->stepmin) {
    info = PetscLogInfo((tao,"TaoApply_BoundLineSearch:Line search error: stepmax (%g) < stepmin (%g)\n",
			 neP->stepmax,neP->stepmin)); CHKERRQ(info);
    *info2 = -6; TaoFunctionReturn(0);
  } else if (neP->maxfev < zero) {
    info = PetscLogInfo((tao,"TaoApply_BoundLineSearch:Line search error: maxfev (%d) < 0\n",neP->maxfev)); CHKERRQ(info);
    *info2 = -7; TaoFunctionReturn(0);
  }

  /* Next 6 lines added by S. Benson to accomodate bounds on the variables */
  /* Compute step length needed to make all variables equal a bound */ 
  /* Compute the smallest steplength that will make one nonbinding variable equal the bound */ 
  info = TaoGetVariableBounds(tao,&XL,&XU); CHKERRQ(info);
  info = S->Negate(); CHKERRQ(info);
  info = S->BoundGradientProjection(S,XL,X,XU); CHKERRQ(info);
  info = S->Negate(); CHKERRQ(info);
  info = X->StepBoundInfo(XL,XU,S,&bstepmin1,&bstepmin2,&bstepmax); CHKERRQ(info);
  neP->stepmax=TaoMin(bstepmax,1.0e+15);

  /* Check that search direction is a descent direction */

#if defined(PETSC_USE_COMPLEX)
  info = G->Dot(S,&cdginit);CHKERRQ(info); dginit = TaoReal(cdginit);
#else
  info = G->Dot(S,&dginit);CHKERRQ(info);
#endif


  if (dginit >= zero) {
    info = PetscLogInfo((tao,"TaoApply_BoundLineSearch:Search direction not a descent direction\n")); CHKERRQ(info);
    *info2 = 7; TaoFunctionReturn(0);
  }

  /* Initialization */
  neP->bracket = 0;
  *info2	  = 0;
  stage1  = 1;
  finit   = *f;
  dgtest  = neP->ftol * dginit;
  width   = neP->stepmax - neP->stepmin;
  width1  = width * two;
  info = W->CopyFrom(X);CHKERRQ(info);
  /* Variable dictionary:  
     stx, fx, dgx - the step, function, and derivative at the best step
     sty, fy, dgy - the step, function, and derivative at the other endpoint 
                   of the interval of uncertainty
     step, f, dg - the step, function, and derivative at the current step */

  stx = zero;
  fx  = finit;
  dgx = dginit;
  sty = zero;
  fy  = finit;
  dgy = dginit;
 
  neP->nfev = 0;
  for (i=0; i< neP->maxfev; i++) {
    /* Set min and max steps to correspond to the interval of uncertainty */
    if (neP->bracket) {
      neP->stepmin = TaoMin(stx,sty); 
      neP->stepmax = TaoMax(stx,sty); 
    } else {
      neP->stepmin = stx;
      neP->stepmax = *step + xtrapf * (*step - stx);
    }

    /* Force the step to be within the bounds */
    *step = TaoMax(*step,neP->stepmin);
    *step = TaoMin(*step,neP->stepmax);

    /* If an unusual termination is to occur, then let step be the lowest
       point obtained thus far */
    if (((neP->bracket) && (*step <= neP->stepmin || *step >= neP->stepmax)) ||
        ((neP->bracket) && (neP->stepmax - neP->stepmin <= neP->rtol * neP->stepmax)) ||
        (neP->nfev >= neP->maxfev - 1) || (neP->infoc == 0)) 
      *step = stx;
    
    info = PetscLogInfo((tao,"TaoApply_LineSearch: Try step size: %6.4e\n",*step)); CHKERRQ(info);
#if defined(PETSC_USE_COMPLEX)
    cstep = *step;
    info = X->Waxpby(cstep,S,1.0,W);CHKERRQ(info);
#else
    info = X->Waxpby(*step,S,1.0,W);CHKERRQ(info); 	/* X = W + step*S */
#endif
    /* This if statement added by S. Benson.   --  Project the solution, if necessary, to keep feasible */
    if (*step >= bstepmin1){
      info=X->Median(XL,X,XU);CHKERRQ(info);
    }
    info = TaoComputeMeritFunctionGradient(tao,X,f,G);CHKERRQ(info);
    neP->nfev++;
#if defined(PETSC_USE_COMPLEX)
    info = G->Dot(S,&cdg);CHKERRQ(info); dg = TaoReal(cdg);
    return 1;
#else
    //    info = G->Dot(S,&dg);CHKERRQ(info);	        /* dg = G^T S */
     /* Next 3 lines  added by S. Benson.   to compute an estimate of decrease in the actual direction and curvature */
    info = G->Dot(W,&dg1);CHKERRQ(info);	        /* dg = G^T S */
    info = G->Dot(X,&dg2);CHKERRQ(info);	        /* dg = G^T S */
    dg=dg2-dg1;
#endif
    ftest1 = finit + (*step) * dgtest;

    /* This test added by S. Benson.   --  Be satisfied with Armijo condition if the Projection X->Median(XL,X,XU) changed X. */
    ftest2 = finit + (*step) * dgtest * neP->ftol;
    

    /* Convergence testing */
    if (((neP->bracket) && (*step <= neP->stepmin||*step >= neP->stepmax)) || (!neP->infoc)) {
      info = PetscLogInfo((tao,"TaoApply_LineSearch:Rounding errors may prevent further progress.  May not be a step satisfying\n")); CHKERRQ(info);
      info = PetscLogInfo((tao,"TaoApply_BoundLineSearch:sufficient decrease and curvature conditions. Tolerances may be too small.\n")); CHKERRQ(info);
      *info2 = 6; break;
    }
    if ((*step == neP->stepmax) && (*f <= ftest1) && (dg <= dgtest)) {
      info = PetscLogInfo((tao,"TaoApply_BoundLineSearch:Step is at the upper bound, stepmax (%g)\n",neP->stepmax)); CHKERRQ(info);
      *info2 = 5; break;
    }
    if ((*step == neP->stepmin) && (*f >= ftest1) && (dg >= dgtest)) {
      info = PetscLogInfo((tao,"TaoApply_BoundLineSearch:Step is at the lower bound, stepmin (%g)\n",neP->stepmin)); CHKERRQ(info);
      *info2 = 4; break;
    }
    if (neP->nfev >= neP->maxfev) {
      info = PetscLogInfo((tao,"TaoApply_BoundLineSearch:Number of line search function evals (%d) > maximum (%d)\n",neP->nfev,neP->maxfev)); CHKERRQ(info);
      *info2 = 3; break;
    }
    if ((neP->bracket) && (neP->stepmax - neP->stepmin <= neP->rtol*neP->stepmax)){
      info = PetscLogInfo((tao,"TaoApply_BoundLineSearch:Relative width of interval of uncertainty is at most rtol (%g)\n",neP->rtol)); CHKERRQ(info);
      *info2 = 2; break;
    }
    if ((*f <= ftest1) && (TaoAbsDouble(dg) <= neP->gtol*(-dginit))) {
      info = PetscLogInfo((tao,"TaoApply_BoundLineSearch:Line search success: Sufficient decrease and directional deriv conditions hold, %6.4e, %6.4e\n",TaoAbsDouble(dg), -dginit )); CHKERRQ(info);
      *info2 = 0; break;
    }

    /* This test added by S. Benson.   --  Be satisfied with Armijo condition if the Projection X->Median(XL,X,XU) changed X. */
    if ( (*f <= ftest2) && ((*step) < bstepmin2) ) {
      info = PetscLogInfo((tao,"TaoApply_BoundLineSearch:Line search success: Sufficient decrease new bounds apply\n")); CHKERRQ(info);
      *info2 = 0; break;
    }

    /* In the first stage, we seek a step for which the modified function
        has a nonpositive value and nonnegative derivative */
    if ((stage1) && (*f <= ftest1) && (dg >= dginit * TaoMin(neP->ftol, neP->gtol)))
      stage1 = 0;

    /* A modified function is used to predict the step only if we
       have not obtained a step for which the modified function has a 
       nonpositive function value and nonnegative derivative, and if a
       lower function value has been obtained but the decrease is not
       sufficient */

    if ((stage1) && (*f <= fx) && (*f > ftest1)) {
      fm   = *f - *step * dgtest;	/* Define modified function */
      fxm  = fx - stx * dgtest;	        /* and derivatives */
      fym  = fy - sty * dgtest;
      dgm  = dg - dgtest;
      dgxm = dgx - dgtest;
      dgym = dgy - dgtest;

      /* Update the interval of uncertainty and compute the new step */
      info = TaoStep_LineSearch(tao,&stx,&fxm,&dgxm,&sty,&fym,&dgym,step,&fm,&dgm);CHKERRQ(info);

      fx  = fxm + stx * dgtest;	/* Reset the function and */
      fy  = fym + sty * dgtest;	/* gradient values */
      dgx = dgxm + dgtest; 
      dgy = dgym + dgtest; 
    } else {
      /* Update the interval of uncertainty and compute the new step */
      info = TaoStep_LineSearch(tao,&stx,&fx,&dgx,&sty,&fy,&dgy,step,f,&dg);CHKERRQ(info);
    }

   /* Force a sufficient decrease in the interval of uncertainty */
   if (neP->bracket) {
     if (TaoAbsDouble(sty - stx) >= p66 * width1) *step = stx + p5*(sty - stx);
       width1 = width;
       width = TaoAbsDouble(sty - stx);
     }
   }

  /* Finish computations */
  info = PetscLogInfo((tao,"TaoApply_BoundLineSearch:%d function evals in line search, step = %10.4e\n",neP->nfev,*step)); CHKERRQ(info);
  /*
  info = G->Norm2(dginitt);CHKERRQ(info);
  */
  TaoFunctionReturn(0);
}

/* ---------------------------------------------------------- */

#undef __FUNCT__  
#define __FUNCT__ "TaoCreateMoreThuenteBoundLineSearch"
/*@C
   TaoCreateMoreThuenteBoundLineSearch - Create a line search

   Input Parameters:
+  tao - TAO_SOLVER context
.  fftol - the sufficient descent parameter , greater than 0.
-  ggtol - the curvature tolerance, greater than 0, less than 1.


   Note:
   If either fftol or ggtol is 0, default parameters will be used.

   Note:
   This algorithm is is a modification of the algorithm by More' and Thuente.
   The modifications concern bounds.  This algorithm steps in the direction
   passed into this routine.  This point get projected back into the feasible set.
   It tries to satisfy the Wolfe conditions, but in the context of bound constrained
   optimization, there may not be a point in the piecewise linear 
   path that satisfies the Wolfe conditions.  When the active set
   is changing, decrease in the objective function may be sufficient
   to terminate this line search.

   Note:
   This routine is used within the following TAO bound constrained
   minimization solvers: Newton trust region (tao_tron) and limited memory variable 
   metric (tao_blmvm).

   Level: developer

.keywords: TAO_SOLVER, linesearch
@*/
int TaoCreateMoreThuenteBoundLineSearch(TAO_SOLVER tao, double fftol, double ggtol)
{
  int info;
  TAO_LINESEARCH *neP;

  TaoFunctionBegin;

  info = TaoNew(TAO_LINESEARCH,&neP); CHKERRQ(info);
  info = PetscLogObjectMemory(tao,sizeof(TAO_LINESEARCH)); CHKERRQ(info);

  if (fftol<=0) fftol=0.001;
  if (ggtol<=0) ggtol=0.9;

  neP->ftol		  = fftol;
  neP->rtol		  = 1.0e-10;
  neP->gtol		  = ggtol;
  neP->stepmin		  = 1.0e-30;
  neP->stepmax		  = 1.0e+20;
  neP->nfev		  = 0; 
  neP->bracket		  = 0; 
  neP->infoc              = 1;
  neP->maxfev		  = 30;

  info = TaoSetLineSearch(tao,0,
			  TaoSetOptions_BoundLineSearch,
			  TaoApply_BoundLineSearch,
			  TaoView_LineSearch,
			  TaoDestroy_LineSearch,
			  (void *) neP);CHKERRQ(info);

  TaoFunctionReturn(0);
}




/*


*/
