GlobiPack Package Browser (Single Doxygen Collection) Version of the Day
GlobiPack_ArmijoPolyInterpLineSearch_def.hpp
Go to the documentation of this file.
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 // Redistribution and use in source and binary forms, with or without
00012 // modification, are permitted provided that the following conditions are
00013 // met:
00014 //
00015 // 1. Redistributions of source code must retain the above copyright
00016 // notice, this list of conditions and the following disclaimer.
00017 //
00018 // 2. Redistributions in binary form must reproduce the above copyright
00019 // notice, this list of conditions and the following disclaimer in the
00020 // documentation and/or other materials provided with the distribution.
00021 //
00022 // 3. Neither the name of the Corporation nor the names of the
00023 // contributors may be used to endorse or promote products derived from
00024 // this software without specific prior written permission.
00025 //
00026 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00027 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00028 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00029 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00030 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00031 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00032 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00033 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00034 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00035 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00036 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00037 //
00038 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov) 
00039 // 
00040 // ***********************************************************************
00041 // @HEADER
00042 */
00043 
00044 #ifndef GLOBIPACK_POLY_INTERP_LINE_SEARCH_DEF_HPP
00045 #define GLOBIPACK_POLY_INTERP_LINE_SEARCH_DEF_HPP
00046 
00047 
00048 #include "GlobiPack_ArmijoPolyInterpLineSearch_decl.hpp"
00049 #include "Teuchos_TabularOutputter.hpp"
00050 
00051 
00052 namespace GlobiPack {
00053 
00054 
00055 // Constructor/Initializers/Accessors
00056 
00057 
00058 template<typename Scalar>
00059 ArmijoPolyInterpLineSearch<Scalar>::ArmijoPolyInterpLineSearch()
00060   : eta_(ArmijoPolyInterpLineSearchUtils::eta_default),
00061     minFrac_(ArmijoPolyInterpLineSearchUtils::minFrac_default),
00062     maxFrac_(ArmijoPolyInterpLineSearchUtils::maxFrac_default),
00063     minIters_(ArmijoPolyInterpLineSearchUtils::minIters_default),
00064     maxIters_(ArmijoPolyInterpLineSearchUtils::maxIters_default),
00065     doMaxIters_(ArmijoPolyInterpLineSearchUtils::doMaxIters_default)
00066 {}
00067 
00068 
00069 template<typename Scalar>
00070 Scalar ArmijoPolyInterpLineSearch<Scalar>::eta() const
00071 {
00072   return eta_;
00073 }
00074 
00075 
00076 template<typename Scalar>
00077 Scalar ArmijoPolyInterpLineSearch<Scalar>::minFrac() const
00078 {
00079   return minFrac_;
00080 }
00081 
00082 
00083 template<typename Scalar>
00084 Scalar ArmijoPolyInterpLineSearch<Scalar>::maxFrac() const
00085 {
00086   return maxFrac_;
00087 }
00088 
00089 
00090 template<typename Scalar>
00091 int ArmijoPolyInterpLineSearch<Scalar>::minIters() const
00092 {
00093   return minIters_;
00094 }
00095 
00096 
00097 template<typename Scalar>
00098 int ArmijoPolyInterpLineSearch<Scalar>::maxIters() const
00099 {
00100   return maxIters_;
00101 }
00102 
00103 
00104 template<typename Scalar>
00105 bool ArmijoPolyInterpLineSearch<Scalar>::doMaxIters() const
00106 {
00107   return doMaxIters_;
00108 }
00109 
00110 
00111 // Overridden from ParameterListAcceptor
00112 
00113 
00114 template<class Scalar>
00115 void ArmijoPolyInterpLineSearch<Scalar>::setParameterList(
00116   RCP<ParameterList> const& paramList
00117   )
00118 {
00119   typedef ScalarTraits<Scalar> ST;
00120   namespace AQLSU = ArmijoPolyInterpLineSearchUtils;
00121   using Teuchos::getParameter;
00122   paramList->validateParametersAndSetDefaults(*this->getValidParameters());
00123   eta_ = getParameter<double>(*paramList, AQLSU::eta_name);
00124   minFrac_ = getParameter<double>(*paramList, AQLSU::minFrac_name);
00125   maxFrac_ = getParameter<double>(*paramList, AQLSU::maxFrac_name);
00126   minIters_ = getParameter<int>(*paramList, AQLSU::minIters_name);
00127   maxIters_ = getParameter<int>(*paramList, AQLSU::maxIters_name);
00128   doMaxIters_ = getParameter<bool>(*paramList, AQLSU::doMaxIters_name);
00129   TEUCHOS_ASSERT_INEQUALITY( eta_, >=, ST::zero() );
00130   TEUCHOS_ASSERT_INEQUALITY( eta_, <, ST::one() );
00131   TEUCHOS_ASSERT_INEQUALITY( minFrac_, >=, ST::zero() );
00132   TEUCHOS_ASSERT_INEQUALITY( minFrac_, <, maxFrac_ );
00133   TEUCHOS_ASSERT_INEQUALITY( minIters_, >=, 0 );
00134   TEUCHOS_ASSERT_INEQUALITY( minIters_, <=, maxIters_ );
00135   setMyParamList(paramList);
00136 }
00137 
00138 
00139 template<class Scalar>
00140 RCP<const ParameterList>
00141 ArmijoPolyInterpLineSearch<Scalar>::getValidParameters() const
00142 {
00143   namespace AQLSU = ArmijoPolyInterpLineSearchUtils;
00144   static RCP<const ParameterList> validPL;
00145   if (is_null(validPL)) {
00146     RCP<Teuchos::ParameterList>
00147       pl = Teuchos::rcp(new Teuchos::ParameterList());
00148     pl->set( AQLSU::eta_name, AQLSU::eta_default );
00149     pl->set( AQLSU::minFrac_name, AQLSU::minFrac_default );
00150     pl->set( AQLSU::maxFrac_name, AQLSU::maxFrac_default );
00151     pl->set( AQLSU::minIters_name, AQLSU::minIters_default );
00152     pl->set( AQLSU::maxIters_name, AQLSU::maxIters_default );
00153     pl->set( AQLSU::doMaxIters_name, AQLSU::doMaxIters_default );
00154     validPL = pl;
00155   }
00156   return validPL;
00157 }
00158 
00159 
00160 // Overrridden from LineSearchBase
00161 
00162 
00163 template<typename Scalar>
00164 bool ArmijoPolyInterpLineSearch<Scalar>::requiresBaseDeriv() const
00165 {
00166   return true;
00167 }
00168 
00169 
00170 template<typename Scalar>
00171 bool ArmijoPolyInterpLineSearch<Scalar>::requiresDerivEvals() const
00172 {
00173   return false;
00174 }
00175 
00176 
00177 template<typename Scalar>
00178 bool ArmijoPolyInterpLineSearch<Scalar>::doLineSearch(
00179   const MeritFunc1DBase<Scalar> &phi,
00180   const PointEval1D<Scalar> &point_k,
00181   const Ptr<PointEval1D<Scalar> > &point_kp1,
00182   const Ptr<int> &numIters_out
00183   ) const
00184 {
00185 
00186   using Teuchos::as;
00187   using Teuchos::ptrFromRef;
00188   using Teuchos::TabularOutputter;
00189   typedef Teuchos::TabularOutputter TO;
00190   typedef ScalarTraits<Scalar> ST;
00191   typedef PointEval1D<Scalar> PE1D;
00192   using std::min;
00193   using std::max;
00194 
00195   const RCP<Teuchos::FancyOStream> out = this->getOStream();
00196 
00197   *out << "\nStarting Armijo Quadratic interpolation linesearch ...\n";
00198 
00199 #ifdef TEUCHOS_DEBUG
00200   TEUCHOS_ASSERT_EQUALITY(point_k.alpha, ST::zero());
00201   TEUCHOS_ASSERT_INEQUALITY(point_k.phi, !=, PE1D::valNotGiven());
00202   TEUCHOS_ASSERT_INEQUALITY(point_k.Dphi, !=, PE1D::valNotGiven());
00203   TEUCHOS_ASSERT(!is_null(point_kp1));
00204   TEUCHOS_ASSERT_INEQUALITY(point_kp1->alpha, >, ST::zero());
00205   TEUCHOS_ASSERT_INEQUALITY(point_kp1->phi, !=, PE1D::valNotGiven());
00206   TEUCHOS_ASSERT_EQUALITY(point_kp1->Dphi, PE1D::valNotGiven());
00207 #endif
00208 
00209   const Scalar phi_k = point_k.phi;
00210   Scalar &alpha_k = point_kp1->alpha;
00211   Scalar &phi_kp1 = point_kp1->phi;
00212 
00213   // Loop initialization (technically the first iteration)
00214 
00215   const Scalar Dphi_k = point_k.Dphi;
00216 
00217   // output header
00218 
00219   *out
00220     << "\nDphi_k = " << Dphi_k
00221     << "\nphi_k = " << phi_k
00222     << "\n";
00223   if (minIters_ > 0)
00224     *out << "\nminIters == " << minIters_ << "\n";
00225   if (doMaxIters_)
00226     *out << "\ndoMaxIters == true, maxing out maxIters = "
00227          <<maxIters_<<" iterations!\n"; 
00228   *out << "\n";
00229   
00230   TabularOutputter tblout(out);
00231   
00232   tblout.pushFieldSpec("itr", TO::INT);
00233   tblout.pushFieldSpec("alpha_k", TO::DOUBLE);
00234   tblout.pushFieldSpec("phi_kp1", TO::DOUBLE);
00235   tblout.pushFieldSpec("phi_kp1-frac_phi", TO::DOUBLE);
00236   tblout.pushFieldSpec("alpha_interp", TO::DOUBLE);
00237   tblout.pushFieldSpec("alpha_min", TO::DOUBLE);
00238   tblout.pushFieldSpec("alpha_max", TO::DOUBLE);
00239 
00240   tblout.outputHeader();
00241 
00242   // Check that this is a decent direction
00243   TEUCHOS_TEST_FOR_EXCEPTION( !(Dphi_k < ST::zero()), Exceptions::NotDescentDirection,
00244     "ArmijoPolyInterpLineSearch::doLineSearch(): "
00245     "The given descent direction for the given "
00246     "phi Dphi_k="<<Dphi_k<<" >= 0!" );
00247   
00248   // keep memory of the best value
00249   Scalar best_alpha = alpha_k;
00250   Scalar best_phi = phi_kp1;
00251 
00252   // Perform linesearch.
00253   bool success = false;
00254   int numIters = 0;
00255   for ( ; numIters < maxIters_; ++numIters ) {
00256 
00257     // Print out this iteration.
00258 
00259     Scalar frac_phi = phi_k + eta_ * alpha_k * Dphi_k;
00260     tblout.outputField(numIters);
00261     tblout.outputField(alpha_k);
00262     tblout.outputField(phi_kp1);
00263     tblout.outputField(((phi_kp1)-frac_phi));
00264     
00265     if (ST::isnaninf(phi_kp1)) {
00266 
00267       // Cut back the step to minFrac * alpha_k
00268       alpha_k = minFrac_ * alpha_k;
00269       best_alpha = ST::zero();
00270       best_phi = phi_k;
00271     }
00272     else {    
00273 
00274       // Armijo condition
00275       if (phi_kp1 < frac_phi) {
00276         // We have found an acceptable point
00277         success = true;
00278         if (numIters < minIters_) {
00279           // Keep going!
00280         }
00281         else if ( !doMaxIters_ || ( doMaxIters_ && numIters == maxIters_ - 1 ) ) {
00282           tblout.nextRow(true);
00283           break;  // get out of the loop, we are done!
00284         }
00285       }
00286       else {
00287         success = false;
00288       }
00289 
00290       // Select a new alpha to try:
00291       //   alpha_k = ( minFracalpha_k <= quadratic interpolation <= maxFracalpha_k )
00292 
00293       // Quadratic interpolation of alpha_k that minimizes phi.
00294       // We know the values of phi at the initail point and alpha_k and
00295       // the derivate of phi w.r.t alpha at the initial point and
00296       // that's enough information for a quandratic interpolation.
00297       
00298       Scalar alpha_quad =
00299         ( -as<Scalar>(0.5) * Dphi_k * alpha_k * alpha_k )
00300         /
00301         ( (phi_kp1) - phi_k - alpha_k * Dphi_k );
00302 
00303       tblout.outputField(alpha_quad);
00304 
00305       // 2009/01/29: rabartl: ToDo: Call the interpolation and then add
00306       // support for other types of interpolations based on various points.
00307 
00308       const Scalar alpha_min = minFrac_ * alpha_k;
00309       const Scalar alpha_max = maxFrac_ * alpha_k;
00310 
00311       tblout.outputField(alpha_min);
00312       tblout.outputField(alpha_max);
00313 
00314       alpha_k =
00315         min(
00316           max(alpha_min, alpha_quad),
00317           alpha_max
00318           );
00319 
00320     }
00321 
00322     tblout.nextRow(true);
00323 
00324     
00325     // Evaluate the point
00326 
00327     phi_kp1 = computeValue<Scalar>(phi, alpha_k);
00328 
00329     // Save the best point found
00330     if (phi_kp1 < best_phi) {
00331       best_phi = phi_kp1;
00332       best_alpha = alpha_k;
00333     }
00334 
00335   }
00336 
00337   if (!is_null(numIters_out))
00338     *numIters_out = numIters;
00339 
00340   if( success ) {
00341     *out << "\nLine search success!\n";
00342     return true;
00343   }
00344 
00345   // Line search failure.  Return the best point found and let the client
00346   // decide what to do.
00347   alpha_k = best_alpha;
00348   phi_kp1 = computeValue<Scalar>(phi, best_alpha);
00349   *out << "\nLine search failure!\n";
00350   return false; 
00351   
00352 }
00353 
00354 
00355 } // namespace GlobiPack
00356 
00357 
00358 #endif // GLOBIPACK_POLY_INTERP_LINE_SEARCH_DEF_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Friends