src/GSS/OptGSS.C

//------------------------------------------------------------------------
// Copyright (C) 2003:
// R.A.Oliva, Lawrence Berkeley National Laboratory.
// raoliva@lbl.gov
//------------------------------------------------------------------------

#ifdef HAVE_CONFIG_H
#include "OPT++_config.h"
#endif

#ifdef HAVE_STD
#include <cstring>
#include <ctime>
#else
#include <string.h>
#include <time.h>
#endif

#include "OptGSS.h"
#include "NLF.h"
#include "ioformat.h"

using namespace std;
using Teuchos::SerialDenseVector;

namespace OPTPP {

void OptGSS::setParams() {
  // use defaults parameters
  OptGSS_params op;
  setParams(op);
}
void OptGSS::setParams(OptGSS_params op) {
  // use the parameters passed as arguments
  Delta     = op.Delta;
  Delta_tol = op.Delta_tol;
  Phi       = op.Phi;
  Theta     = op.Theta;
  Iter_max  = op.Iter_max;
  printCOPYRIGHT = op.printCOPYRIGHT;
  printXiter  = op.printXiter;
  printGiter  = op.printGiter;

  // default non-user parameters:
  computeGrad = ((nlp1)? true : false);
  mpi_rank = 0;
#ifdef OPTPP_HAVE_MPI
  setpid();
#endif
}

void OptGSS::reset()
{
   int n = nlp->getDim();
   nlp->reset();
   OptimizeClass::defaultReset(n);
   setParams();
   // Still to do - double check parameters in the genSet class.
}

void OptGSS::initOpt()
{
  if (debug_)
    nlp->setDebug();

  if(nlp->hasConstraints()){
    cerr << "Error: OptGSS does not support bound, linear, or nonlinear "
         << "constraints.\n       Please select a different method for "
         << "constrained problems." << endl;
    abort_handler(-1);
  }

  bool debug = nlp->getDebug();

  extras_srched = false; // before optimize()

  nlp->initFcn();  // sets initial point

  if (nlp1)
    nlp1->eval();
  else
    nlp->eval();

  if (debug) {
    *optout << "NLP Initialized in OptGSS::initOpt()\n";
    if (nlp1)
      *optout << "GSS::initOpt() - NLP1 eval()\n";
    else
      *optout << "GSS::initOpt() - NLP0 eval()\n";
    optout->flush();

  } // END isInit==false

  //--
  // Initialize X, fX, gX (if available), and genset.
  //--
  X = nlp->getXc();

  if (nlp1) {
    gX = nlp1->getGrad();
    gset->init(gX);
  }
  else {
    gset->init();
  }

  fX = nlp->getF();
  fprev = fX;

  //--
  // init step-size if not set by user
  //--
  if (Delta==0.0) {

    // set to max(|x_i|) ....
    for (int i=0; i<X.length(); i++) {
      double xi = fabs(X(i));
      if (xi > Delta) Delta = xi;
    }

    // ... or to 1 if all |x_i| == 0
    if (Delta==0.0) {
      Delta = 1.0;
    }

  } // END if (Delta==0)

  //--
  // Print out header with start time and "0th-iteration" info
  //--
  printHeader();
  printIter(0,0);

}

void OptGSS::optimize()
//------------------------------------------------------------------------
// Generating Set Search.
//
// Solves the unconstrained minimization problem
//
//          min F(x),    x = (x_1, x_2, ..., x_N)
//
// using the Generating Set Search Algorithm.
//
// References:
//
// T. Kolda and V. Torczon, Optimization by Direct Search: New Perspectives
// on Some Classical and Modern Methods, SIAM J. of Optimization. 2004.
//
// Description of the algorithm
//
// 1. On each iteration, evaluate f on the search-space.
// 2. If a function decrease is found, advance to that point
//    and (optionally) increase the step-size;
//    Otherwise, decrease the step size.
//
// 3. Stop iterations if step-size becomes too small or
//    either function or gradient convergence tolarances are met.
//    Else, do updates and iterate.
//
// Notes:
//
// (1)
// The search-space consists of points generated by the
// generating set, pruned by the gradient when available,
// plus an optional matrix of extra search points (in columns).
//
// (2)
// Each iteration changes either the step size
// or the current point, but not both, so we split checkConvg()
// into its 3 componenets: step-tol, fcn-tol, gradnorm-tol.
//
// (3)
// The update of the gradient will be skipped when
//   * there is only one iteration,
//               --AND--
//   * the computeGrad flag is false.
// This is to save evaluating the gradient when it might
// not be used by the calling routine (i.e. trustGSS).
// WARNING:
//   In this case, the grad in the NLP1 may not correspond to the
//   current point xc when optimize() terminates; the calling
//   program must do the gradient update.
//
{ // optimize BEGIN

  //--

  SpecOption SpecTmp = nlp->getSpecOption();
  nlp->setSpecOption(NoSpec);

  initOpt();

  bool done = StepCondition(); // should be false initially
  if (done) {
    *optout << "!!! Step tolerance met "
	    << "before iterations begin !!!\n";
    cerr   << "Warning: step tolerance met "
	    << "before iterations begin!\n*******\n";
    ret_code = 1;
    setReturnCode(ret_code);
  }
  else
    ret_code = 0;

  int iter = 0;
  int bestid = -1;

  while (!done) {

    ++iter;

    //--
    // search for a better point on {gen set mesh} union {extra search pts}
    // update current point if we find one.
    bestid = search();

    //--
    // update step size
    //--
    if (bestid==-1) { // search failed

      contractStep();

    }
    else if (bestid <= gset->size()) {   // advanced to a mesh point";

      expandStep();

    }

    // print iteration data, pass best point
    printIter(iter, bestid);

    // if search failed, check reduced step > min step
    if (bestid == -1)

      ret_code = StepCondition();

    else { // search succeeded ==> have a new point

      // check convergence on fcn value if we are iterating
      if (Iter_max > 1)
	ret_code = checkConvg_fcn();

      // gradient update -- see note (2) above.
      if (nlp1)
	if (Iter_max > 1 || computeGrad==true) {
	  nlp1->evalG();
	  gX = nlp1->getGrad();

	  if (ret_code == 0)
	    ret_code = checkConvg_grad();
	}

    } //

    // check iteration condition
    done = (ret_code || iter == Iter_max);

    if (!done) { // update: compute D_k = G_k union H_k
      if (nlp1 && gset->prunes() )
	gset->update(gX);
      else
	gset->update();
    }

  } // END main loop

  iter_taken = iter;
  if (ret_code ==0) { // no conv. criteria was met
    ret_code = -4;
    setReturnCode(ret_code);
    std::strcpy(mesg,"Algorithm terminated - Number of iterations exceeds the specified limit ");
  }

  nlp->setSpecOption(SpecTmp);

} // END optimize() ____________________________________


int OptGSS::checkConvg() {
  // all convergence tests - currently not used.
  int rc;
  if ( rc = StepCondition() ) return rc;
  if ( rc = checkConvg_fcn() ) return rc;
  if ( rc = checkConvg_grad() ) return rc;
  return 0;
}

//--
// Check convergence wrt step size
int OptGSS::StepCondition() {

  if ( Delta > Delta_tol ) return 0;

  std::strcpy(mesg,"Algorithm converged - Norm of last step is less than step tolerance");
  if (mpi_rank==0) {
    *optout << "             \tSteplength = " << e(Delta,12,4)
	    << " Steplength Tolerance: " << e(Delta_tol,12,4) << endl;
  }
  setReturnCode(1);
  return 1;

}

//--
// Checks convergence of function value
int OptGSS::checkConvg_fcn() {

  double ftol = tol.getFTol();
  double fvalue = fX;
  double rftol = ftol*max(1.0,fabs(fprev));
  double deltaf = fprev - fvalue;

  if (deltaf <= rftol) {
    std::strcpy(mesg,"Algorithm converged - Difference of successive fcn values is less than fcn tolerance");
    if (mpi_rank==0) {
      *optout << "checkConvg():\tdeltaf = " << e(deltaf,12,4)
	      << "  ftol = " << e(ftol,12,4) << "\n";
    }
    setReturnCode(2);
    return 2;
  }
  return 0;
}

//--
// Checks convergence wrt norm of the gradient
int OptGSS::checkConvg_grad()
{
  if (!nlp1) return 0; // no gradient

  // Test 3. gradient tolerance

  //  SerialDenseVector<int,double> grad(nlp1->getGrad());
  double gtol = tol.getGTol();
  double rgtol = gtol*max(1.0,fabs(fX));
  //double gnorm = Norm2(gX);
  double gnorm = std::sqrt(gX.dot(gX));
  if (gnorm <= rgtol) {
    std::strcpy(mesg,"Algorithm converged - Norm of gradient is less than gradient tolerance");
    if (mpi_rank==0) {
      *optout << "checkConvg():\tgnorm = " << e(gnorm,12,4)
	      << "  gtol = " << e(rgtol, 12,4) << "\n";
    }
    setReturnCode(3);
    return 3;
  }

  // Test 4. absolute gradient tolerance

  if (gnorm <= gtol) {
    std::strcpy(mesg,"Algorithm converged - Norm of gradient is less than gradient tolerance");
    if (mpi_rank==0) {
      *optout << "checkConvg: gnorm = " << e(gnorm,12,4)
	      << " gtol = " << e(gtol, 12,4) << "\n";
    }
    setReturnCode(4);
    return 4;
  }

  return 0; // Nothing to report
}

// --
// Search Method: Parallel and Serial, For GenSet AND Extras
//--
int OptGSS::search() {

  bool debug = nlp->getDebug();

  // size of search space:
  int searchsize;
  searchsize = gset->nActive() + extras.numCols();

  if (searchsize == 0) {
    if (debug) *optout << "*-*-* empty search set! -- nothing done\n";
    return -1;
  }

  // the size of the active gen-set mesh
  int gssz = gset->nActive();

  // vars to keep best point, its index and its fval
  SerialDenseVector<int,double> bestx(X);
  int          besti = -1;
  double       bestf = fX;

  // vars to keep trial point and its fval
  SerialDenseVector<int,double> xi(X);
  double       fi;

#ifdef OPTPP_HAVE_MPI

  double pll_bestf;
  int    pll_besti;

  // [imin,imax] <==> this processor's set of directions

  int imin = searchsize * mpi_rank / mpi_size + 1;
  int imax = searchsize * (mpi_rank+1) / mpi_size;

#else

  // [imin,imax] <==> all search directions
  int imin = 1;
  int imax = searchsize;

#endif

  if (debug) {

    int gssz = gset->nActive();
    if (imax <= gssz)
      *optout << "Searching genset " << imin << " to " << imax << endl;
    else if (imin > gssz)
      *optout << "Searching extras " << imin-gssz << " to " << imax-gssz << endl;
    else
      *optout << "Searching genset " << imin << " to " << gssz
	      << " and extras 1 to " << imax - gssz << endl;

    optout->flush();
  }

  //--
  // Search loop
  //--

  for (int i = imin; i <= imax;  i++) {

    if (i<=gssz)
      gset->generateActive(i-1, Delta, X, xi); // sets xi = X + Delta * Act_i
    else{
      for(int j=0;j<extras.numRows();j++)
	{xi(j) = extras(j,i-gssz-1);}
      // xi = extras.Column(i-gssz);
    }

    fi = nlp->evalF(xi);

    if (fi < bestf) {  // (fi < fX - forcingfun()) {
      besti = i-1;
      bestf = fi;
      bestx = xi;
    }

  } // END for


#ifdef OPTPP_HAVE_MPI
  // MIN-reduce local best to global best among all processors
  MPI_Allreduce(&bestf, &pll_bestf, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD);
  if (debug) {
    *optout << "First MPI_Allreduce() completed\n";
    *optout << "-- Current fX   = " << fX << endl;
    *optout << "-- local  bestf = " << bestf << endl;
    *optout << "-- global bestf = " << pll_bestf << endl;
    optout->flush();
  }

  // if there is no better point than current pt, all process are done
  if (pll_bestf == fX) { return -1; }

  // if we dont hold the best point unset besti
  if (pll_bestf != bestf)  besti = -1;

  // MAX-reduce best i to catch the global best
  MPI_Allreduce(&besti, &pll_besti, 1, MPI_INT, MPI_MAX, MPI_COMM_WORLD);
  if (debug) {
    *optout << "Second MPI_Allreduce() completed\n";
    *optout << "-- local besti  = " << besti << endl;
    *optout << "-- global besti = " << pll_besti << endl;
    optout->flush();
  }

  // Update X, fX, ...

  besti  = pll_besti;
  bestf  = pll_bestf;
  if (pll_besti > gssz)
    {
      // bestx = extras.Column(pll_besti-gssz);
      for(int j=1;j<=bestx.length();j++)
	{bestx(j) = extras(j,pll_besti-gssz);}
    }
  else {
    gset->generateActive(pll_besti, Delta, X, bestx);
  }

#else

  if (besti == -1) {  // search failed to find better pt

    if (debug)
      *optout << "search() done. No improved pt found.\n";
      //  *optout << "search(" << flag << ") done. No improved pt found.\n";

    return (-1);
  }

#endif

  if (debug) {
    if (besti==-1)
      *optout << "!!! GSS::search() error: besti is zero past return condition.\n";
      //      *optout << "!!! GSS::search(" << flag
      //	   << ") error: besti is zero past return condition.\n";
  }

  if (debug && besti<=gssz) {
    int acid = gset->activeID(besti);
    if ( acid==0 || acid > gset->size() )
      *optout << "!!! GSS:search() error: Invalid gset->ActiveID(" << besti << "):  " << acid << "\n";
  }

  SerialDenseVector<int,double> xc(X);
  xc = nlp->getXc();

  if (debug) {
    if (besti<=gssz) {
      gset->generateActive(besti, Delta, X, xi);
      if (bestx != xi ) {
	*optout << "!!! GSS search() : gset->Active(" << besti << ") != bestx ***\n";
      }
    }
    else {
      // xi = extras.Column(besti-gssz);
      for(int j=0;j<xi.length();j++)
	{xi(j) = extras(j,besti-gssz-1);}
      if (bestx != xi ) {
	*optout << "!!! GSS search() : extras(:," << besti-gssz << ") != bestx ***\n";
      }
    }
  }


  // Search succeeded, so do update:
  fprev = fX;
  X     = bestx;
  fX    = bestf;
  nlp->setX(X);
  nlp->setF(fX);

  // return actual idx of best dir in ordering [genset | extras]
  int idx;
  if (besti<=gssz)
    idx = gset->activeID(besti);
  else
    idx = besti-gssz + gset->size();

  if (debug)
    *optout << "search() done. Best pt idx = " << idx << endl;
    //    *optout << "search(" << flag << ") done. Best pt idx = " << idx << endl;

  return idx;

} // END search()


void OptGSS::printStatus(char *s, bool printSoln) // set Message
{
  *optout << "\n\n=========  " << s << "  ===========\n\n";
  *optout << "Optimization method       = " << method << "\n";
  *optout << "Dimension of the problem  = " << nlp->getDim() << "\n";
  *optout << "Return code               = " << ret_code << " ("
       << mesg << ")\n";
  *optout << "No. iterations taken      = " << iter_taken  << "\n";
  *optout << "No. iterations allowed    = " << Iter_max    << "\n";
  *optout << "No. function evaluations  = " << nlp->getFevals() << "\n";
  *optout << "Last step length          = " << Delta << "\n";
  *optout << "Last function value       = " << nlp->getF() << "\n";
  *optout << "Norm of last point        = " << std::sqrt((nlp->getXc()).dot(nlp->getXc())) << "\n";
  if (nlp1)
    *optout << "Norm of last gradient     = " << std::sqrt((nlp1->getGrad()).dot(nlp1->getGrad())) << "\n";


  if (printSoln) {
    *optout << "\n\n=========  " << "Solution"  << "  ===========\n\n";
    *optout << "   i   \t" << "x" << endl;
    for (int i=0; i<gset->vdim(); i++)
      *optout << d(i,5) << "\t" << e(X(i),12,4) << endl;
    *optout << "\n\n";
  }

  tol.printTol(optout);
}

//--
// Prints iteration information to output file
//
void OptGSS::printIter(int iter, int imp) {

  *optout << d(iter,5) << " " << e(fX,12,4) << "\t"
	  << e(Delta,12,4);
  if (nlp1) {
    *optout << "\t" << e(std::sqrt(gX.dot(gX)),4);
  }

  int nsearched = 0;
  if (iter==1) nsearched = gset->nActive() + extras.numCols();
  *optout << "\t" << d(nsearched,5);

  *optout << "\t" << d(imp,5)
	  << "\t" << d(nlp->getFevals(),8);

  if (printXiter) {
    // print first 3 components of X
    *optout << "\t";
    int k = nlp->getDim();
    if (k>3) k=3;
    for (int i=0; i<k; i++)
      *optout << f(X(i),8,4) << " ";
  }

  if (printGiter && nlp1) {
  // print first 3 components of gX
    *optout << "\t";
    int k = nlp->getDim();
    if (k>3) k=3;
    for (int i=0; i<k; i++)
      *optout << f(gX(i),8,4) << " ";
  }

  // finish output line;
  *optout << endl;
}


//--
// Prints output header before iterations begin
void OptGSS::printHeader() {

  //  time_t t;
  //  char *c;

  //  t = time(NULL);
  //  t = time(0);
  //  c = asctime(localtime(&t));

  if (printCOPYRIGHT) {
    *optout << "************************************************************\n";
    *optout << "OPT++ version " << OPT_GLOBALS::OPT_VERSION << "\n";
    //    *optout << "Job run at " << c << "\n";
    copyright();
    *optout << "************************************************************\n";
  }

  *optout << method << endl
	  << "Iter \t\t F(x)\t    ||step||";
  if (nlp1) {
    *optout << "\t||gX||"
	    << "\t ndir";
      }
  *optout << "\tbesti\t   fevals \t";
  if (printXiter)         *optout << "\t X(1:3)";
  if (nlp1 && printGiter) *optout << "\t gX(1:3)";
  *optout << "\n\n";

}

} // namespace OPTPP