GlobiPack Version of the Day
GlobiPack_ArmijoPolyInterpLineSearch_def.hpp
00001 /*
00002 // @HEADER
00003 // ***********************************************************************
00004 // 
00005 //    GlobiPack: Collection of Scalar 1D globalizaton utilities
00006 //                 Copyright (2009) Sandia Corporation
00007 // 
00008 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00009 // license for use of this work by or on behalf of the U.S. Government.
00010 // 
00011 // This library is free software; you can redistribute it and/or modify
00012 // it under the terms of the GNU Lesser General Public License as
00013 // published by the Free Software Foundation; either version 2.1 of the
00014 // License, or (at your option) any later version.
00015 // 
00016 // This library is distributed in the hope that it will be useful, but
00017 // WITHOUT ANY WARRANTY; without even the implied warranty of
00018 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00019 // Lesser General Public License for more details.
00020 // 
00021 // You should have received a copy of the GNU Lesser General Public
00022 // License along with this library; if not, write to the Free Software
00023 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00024 // USA
00025 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov) 
00026 // 
00027 // ***********************************************************************
00028 // @HEADER
00029 */
00030 
00031 #ifndef GLOBIPACK_POLY_INTERP_LINE_SEARCH_DEF_HPP
00032 #define GLOBIPACK_POLY_INTERP_LINE_SEARCH_DEF_HPP
00033 
00034 
00035 #include "GlobiPack_ArmijoPolyInterpLineSearch_decl.hpp"
00036 #include "Teuchos_TabularOutputter.hpp"
00037 
00038 
00039 namespace GlobiPack {
00040 
00041 
00042 // Constructor/Initializers/Accessors
00043 
00044 
00045 template<typename Scalar>
00046 ArmijoPolyInterpLineSearch<Scalar>::ArmijoPolyInterpLineSearch()
00047   : eta_(ArmijoPolyInterpLineSearchUtils::eta_default),
00048     minFrac_(ArmijoPolyInterpLineSearchUtils::minFrac_default),
00049     maxFrac_(ArmijoPolyInterpLineSearchUtils::maxFrac_default),
00050     minIters_(ArmijoPolyInterpLineSearchUtils::minIters_default),
00051     maxIters_(ArmijoPolyInterpLineSearchUtils::maxIters_default),
00052     doMaxIters_(ArmijoPolyInterpLineSearchUtils::doMaxIters_default)
00053 {}
00054 
00055 
00056 template<typename Scalar>
00057 Scalar ArmijoPolyInterpLineSearch<Scalar>::eta() const
00058 {
00059   return eta_;
00060 }
00061 
00062 
00063 template<typename Scalar>
00064 Scalar ArmijoPolyInterpLineSearch<Scalar>::minFrac() const
00065 {
00066   return minFrac_;
00067 }
00068 
00069 
00070 template<typename Scalar>
00071 Scalar ArmijoPolyInterpLineSearch<Scalar>::maxFrac() const
00072 {
00073   return maxFrac_;
00074 }
00075 
00076 
00077 template<typename Scalar>
00078 int ArmijoPolyInterpLineSearch<Scalar>::minIters() const
00079 {
00080   return minIters_;
00081 }
00082 
00083 
00084 template<typename Scalar>
00085 int ArmijoPolyInterpLineSearch<Scalar>::maxIters() const
00086 {
00087   return maxIters_;
00088 }
00089 
00090 
00091 template<typename Scalar>
00092 bool ArmijoPolyInterpLineSearch<Scalar>::doMaxIters() const
00093 {
00094   return doMaxIters_;
00095 }
00096 
00097 
00098 // Overridden from ParameterListAcceptor
00099 
00100 
00101 template<class Scalar>
00102 void ArmijoPolyInterpLineSearch<Scalar>::setParameterList(
00103   RCP<ParameterList> const& paramList
00104   )
00105 {
00106   typedef ScalarTraits<Scalar> ST;
00107   namespace AQLSU = ArmijoPolyInterpLineSearchUtils;
00108   using Teuchos::getParameter;
00109   paramList->validateParametersAndSetDefaults(*this->getValidParameters());
00110   eta_ = getParameter<double>(*paramList, AQLSU::eta_name);
00111   minFrac_ = getParameter<double>(*paramList, AQLSU::minFrac_name);
00112   maxFrac_ = getParameter<double>(*paramList, AQLSU::maxFrac_name);
00113   minIters_ = getParameter<int>(*paramList, AQLSU::minIters_name);
00114   maxIters_ = getParameter<int>(*paramList, AQLSU::maxIters_name);
00115   doMaxIters_ = getParameter<bool>(*paramList, AQLSU::doMaxIters_name);
00116   TEUCHOS_ASSERT_INEQUALITY( eta_, >=, ST::zero() );
00117   TEUCHOS_ASSERT_INEQUALITY( eta_, <, ST::one() );
00118   TEUCHOS_ASSERT_INEQUALITY( minFrac_, >=, ST::zero() );
00119   TEUCHOS_ASSERT_INEQUALITY( minFrac_, <, maxFrac_ );
00120   TEUCHOS_ASSERT_INEQUALITY( minIters_, >=, 0 );
00121   TEUCHOS_ASSERT_INEQUALITY( minIters_, <=, maxIters_ );
00122   setMyParamList(paramList);
00123 }
00124 
00125 
00126 template<class Scalar>
00127 RCP<const ParameterList>
00128 ArmijoPolyInterpLineSearch<Scalar>::getValidParameters() const
00129 {
00130   namespace AQLSU = ArmijoPolyInterpLineSearchUtils;
00131   static RCP<const ParameterList> validPL;
00132   if (is_null(validPL)) {
00133     RCP<Teuchos::ParameterList>
00134       pl = Teuchos::rcp(new Teuchos::ParameterList());
00135     pl->set( AQLSU::eta_name, AQLSU::eta_default );
00136     pl->set( AQLSU::minFrac_name, AQLSU::minFrac_default );
00137     pl->set( AQLSU::maxFrac_name, AQLSU::maxFrac_default );
00138     pl->set( AQLSU::minIters_name, AQLSU::minIters_default );
00139     pl->set( AQLSU::maxIters_name, AQLSU::maxIters_default );
00140     pl->set( AQLSU::doMaxIters_name, AQLSU::doMaxIters_default );
00141     validPL = pl;
00142   }
00143   return validPL;
00144 }
00145 
00146 
00147 // Overrridden from LineSearchBase
00148 
00149 
00150 template<typename Scalar>
00151 bool ArmijoPolyInterpLineSearch<Scalar>::requiresBaseDeriv() const
00152 {
00153   return true;
00154 }
00155 
00156 
00157 template<typename Scalar>
00158 bool ArmijoPolyInterpLineSearch<Scalar>::requiresDerivEvals() const
00159 {
00160   return false;
00161 }
00162 
00163 
00164 template<typename Scalar>
00165 bool ArmijoPolyInterpLineSearch<Scalar>::doLineSearch(
00166   const MeritFunc1DBase<Scalar> &phi,
00167   const PointEval1D<Scalar> &point_k,
00168   const Ptr<PointEval1D<Scalar> > &point_kp1,
00169   const Ptr<int> &numIters_out
00170   ) const
00171 {
00172 
00173   using Teuchos::as;
00174   using Teuchos::ptrFromRef;
00175   using Teuchos::TabularOutputter;
00176   typedef Teuchos::TabularOutputter TO;
00177   typedef ScalarTraits<Scalar> ST;
00178   typedef PointEval1D<Scalar> PE1D;
00179   using std::min;
00180   using std::max;
00181 
00182   const RCP<Teuchos::FancyOStream> out = this->getOStream();
00183 
00184   *out << "\nStarting Armijo Quadratic interpolation linesearch ...\n";
00185 
00186 #ifdef TEUCHOS_DEBUG
00187   TEUCHOS_ASSERT_EQUALITY(point_k.alpha, ST::zero());
00188   TEUCHOS_ASSERT_INEQUALITY(point_k.phi, !=, PE1D::valNotGiven());
00189   TEUCHOS_ASSERT_INEQUALITY(point_k.Dphi, !=, PE1D::valNotGiven());
00190   TEUCHOS_ASSERT(!is_null(point_kp1));
00191   TEUCHOS_ASSERT_INEQUALITY(point_kp1->alpha, >, ST::zero());
00192   TEUCHOS_ASSERT_INEQUALITY(point_kp1->phi, !=, PE1D::valNotGiven());
00193   TEUCHOS_ASSERT_EQUALITY(point_kp1->Dphi, PE1D::valNotGiven());
00194 #endif
00195 
00196   const Scalar phi_k = point_k.phi;
00197   Scalar &alpha_k = point_kp1->alpha;
00198   Scalar &phi_kp1 = point_kp1->phi;
00199 
00200   // Loop initialization (technically the first iteration)
00201 
00202   const Scalar Dphi_k = point_k.Dphi;
00203 
00204   // output header
00205 
00206   *out
00207     << "\nDphi_k = " << Dphi_k
00208     << "\nphi_k = " << phi_k
00209     << "\n";
00210   if (minIters_ > 0)
00211     *out << "\nminIters == " << minIters_ << "\n";
00212   if (doMaxIters_)
00213     *out << "\ndoMaxIters == true, maxing out maxIters = "
00214          <<maxIters_<<" iterations!\n"; 
00215   *out << "\n";
00216   
00217   TabularOutputter tblout(out);
00218   
00219   tblout.pushFieldSpec("itr", TO::INT);
00220   tblout.pushFieldSpec("alpha_k", TO::DOUBLE);
00221   tblout.pushFieldSpec("phi_kp1", TO::DOUBLE);
00222   tblout.pushFieldSpec("phi_kp1-frac_phi", TO::DOUBLE);
00223   tblout.pushFieldSpec("alpha_interp", TO::DOUBLE);
00224   tblout.pushFieldSpec("alpha_min", TO::DOUBLE);
00225   tblout.pushFieldSpec("alpha_max", TO::DOUBLE);
00226 
00227   tblout.outputHeader();
00228 
00229   // Check that this is a decent direction
00230   TEST_FOR_EXCEPTION( !(Dphi_k < ST::zero()), Exceptions::NotDescentDirection,
00231     "ArmijoPolyInterpLineSearch::doLineSearch(): "
00232     "The given descent direction for the given "
00233     "phi Dphi_k="<<Dphi_k<<" >= 0!" );
00234   
00235   // keep memory of the best value
00236   Scalar best_alpha = alpha_k;
00237   Scalar best_phi = phi_kp1;
00238 
00239   // Perform linesearch.
00240   bool success = false;
00241   int numIters = 0;
00242   for ( ; numIters < maxIters_; ++numIters ) {
00243 
00244     // Print out this iteration.
00245 
00246     Scalar frac_phi = phi_k + eta_ * alpha_k * Dphi_k;
00247     tblout.outputField(numIters);
00248     tblout.outputField(alpha_k);
00249     tblout.outputField(phi_kp1);
00250     tblout.outputField(((phi_kp1)-frac_phi));
00251     
00252     if (ST::isnaninf(phi_kp1)) {
00253 
00254       // Cut back the step to minFrac * alpha_k
00255       alpha_k = minFrac_ * alpha_k;
00256       best_alpha = ST::zero();
00257       best_phi = phi_k;
00258     }
00259     else {    
00260 
00261       // Armijo condition
00262       if (phi_kp1 < frac_phi) {
00263         // We have found an acceptable point
00264         success = true;
00265         if (numIters < minIters_) {
00266           // Keep going!
00267         }
00268         else if ( !doMaxIters_ || ( doMaxIters_ && numIters == maxIters_ - 1 ) ) {
00269           tblout.nextRow(true);
00270           break;  // get out of the loop, we are done!
00271         }
00272       }
00273       else {
00274         success = false;
00275       }
00276 
00277       // Select a new alpha to try:
00278       //   alpha_k = ( minFracalpha_k <= quadratic interpolation <= maxFracalpha_k )
00279 
00280       // Quadratic interpolation of alpha_k that minimizes phi.
00281       // We know the values of phi at the initail point and alpha_k and
00282       // the derivate of phi w.r.t alpha at the initial point and
00283       // that's enough information for a quandratic interpolation.
00284       
00285       Scalar alpha_quad =
00286         ( -as<Scalar>(0.5) * Dphi_k * alpha_k * alpha_k )
00287         /
00288         ( (phi_kp1) - phi_k - alpha_k * Dphi_k );
00289 
00290       tblout.outputField(alpha_quad);
00291 
00292       // 2009/01/29: rabartl: ToDo: Call the interpolation and then add
00293       // support for other types of interpolations based on various points.
00294 
00295       const Scalar alpha_min = minFrac_ * alpha_k;
00296       const Scalar alpha_max = maxFrac_ * alpha_k;
00297 
00298       tblout.outputField(alpha_min);
00299       tblout.outputField(alpha_max);
00300 
00301       alpha_k =
00302         min(
00303           max(alpha_min, alpha_quad),
00304           alpha_max
00305           );
00306 
00307     }
00308 
00309     tblout.nextRow(true);
00310 
00311     
00312     // Evaluate the point
00313 
00314     phi_kp1 = computeValue<Scalar>(phi, alpha_k);
00315 
00316     // Save the best point found
00317     if (phi_kp1 < best_phi) {
00318       best_phi = phi_kp1;
00319       best_alpha = alpha_k;
00320     }
00321 
00322   }
00323 
00324   if (!is_null(numIters_out))
00325     *numIters_out = numIters;
00326 
00327   if( success ) {
00328     *out << "\nLine search success!\n";
00329     return true;
00330   }
00331 
00332   // Line search failure.  Return the best point found and let the client
00333   // decide what to do.
00334   alpha_k = best_alpha;
00335   phi_kp1 = computeValue<Scalar>(phi, best_alpha);
00336   *out << "\nLine search failure!\n";
00337   return false; 
00338   
00339 }
00340 
00341 
00342 } // namespace GlobiPack
00343 
00344 
00345 #endif // GLOBIPACK_POLY_INTERP_LINE_SEARCH_DEF_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Friends