GlobiPack Package Browser (Single Doxygen Collection) Version of the Day
GlobiPack_GoldenQuadInterpBracket_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 // 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_GOLDEN_BRACKET_QUAD_INTERP_DEF_HPP
00032 #define GLOBIPACK_GOLDEN_BRACKET_QUAD_INTERP_DEF_HPP
00033 
00034 
00035 #include "GlobiPack_GoldenQuadInterpBracket_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 GoldenQuadInterpBracket<Scalar>::GoldenQuadInterpBracket()
00047 {}
00048 
00049 
00050 // Overridden from ParameterListAcceptor (simple forwarding functions)
00051 
00052 
00053 template<typename Scalar>
00054 void GoldenQuadInterpBracket<Scalar>::setParameterList(RCP<ParameterList> const& paramList)
00055 {
00056   typedef ScalarTraits<Scalar> ST;
00057   paramList->validateParametersAndSetDefaults(*this->getValidParameters());
00058   // ToDo: Add parameters!
00059   setMyParamList(paramList);
00060 }
00061 
00062 
00063 template<typename Scalar>
00064 RCP<const ParameterList> GoldenQuadInterpBracket<Scalar>::getValidParameters() const
00065 {
00066   static RCP<const ParameterList> validPL;
00067   if (is_null(validPL)) {
00068     RCP<Teuchos::ParameterList>
00069       pl = Teuchos::rcp(new Teuchos::ParameterList());
00070     // ToDo: Add parameters!
00071     validPL = pl;
00072   }
00073   return validPL;
00074 }
00075 
00076 
00077 // Bracket
00078 
00079 
00080 template<typename Scalar>
00081 bool GoldenQuadInterpBracket<Scalar>::bracketMinimum(
00082   const MeritFunc1DBase<Scalar> &phi,
00083   const Ptr<PointEval1D<Scalar> > &pointLower,
00084   const Ptr<PointEval1D<Scalar> > &pointMiddle,
00085   const Ptr<PointEval1D<Scalar> > &pointUpper,
00086   const Ptr<int> &numIters
00087   ) const
00088 {
00089 
00090   using Teuchos::as;
00091   using Teuchos::TabularOutputter;
00092   typedef Teuchos::TabularOutputter TO;
00093   typedef ScalarTraits<Scalar> ST;
00094   using Teuchos::OSTab;
00095   typedef PointEval1D<Scalar> PE1D;
00096   
00097 #ifdef TEUCHOS_DEBUG
00098   TEST_FOR_EXCEPT(is_null(pointLower));
00099   TEST_FOR_EXCEPT(is_null(pointUpper));
00100   TEST_FOR_EXCEPT(is_null(pointMiddle));
00101   TEUCHOS_ASSERT_INEQUALITY(pointLower->alpha, <, pointMiddle->alpha);
00102   TEUCHOS_ASSERT_INEQUALITY(pointLower->phi, !=, PE1D::valNotGiven());
00103   TEUCHOS_ASSERT_INEQUALITY(pointMiddle->phi, !=, PE1D::valNotGiven());
00104 #endif
00105 
00106   const RCP<Teuchos::FancyOStream> out = this->getOStream();
00107   
00108   // ToDo: Make these variable!
00109   const Scalar GOLDEN_RATIO = 1.618033988749895;
00110   const Scalar SMALL_DIV = 1e-20;
00111   const Scalar MAX_EXTRAP_FACTOR = 100.0;
00112   const int MAX_TOTAL_ITERS = 30;
00113 
00114   *out << "\nStarting golden quadratic interpolating bracketing of the minimum ...\n\n";
00115   
00116   // Repeatedly evaluate the function along the search direction until
00117   // we know we've bracketed a minimum.
00118  
00119   Scalar &alpha_l = pointLower->alpha, &phi_l = pointLower->phi;
00120   Scalar &alpha_m = pointMiddle->alpha, &phi_m = pointMiddle->phi;
00121   Scalar &alpha_u = pointUpper->alpha = ST::nan(), &phi_u = pointUpper->phi = ST::nan();
00122 
00123   Scalar tmp = ST::nan(), q = ST::nan(), r = ST::nan();
00124 
00125   const Scalar zero = ST::zero();
00126  
00127   // This does a simple backtracking 
00128   alpha_u = zero;
00129   const Scalar goldinv = 1.0/(1.0+GOLDEN_RATIO);
00130   
00131   TabularOutputter tblout(out);
00132   
00133   tblout.pushFieldSpec("itr", TO::INT);
00134   tblout.pushFieldSpec("alpha_l", TO::DOUBLE);
00135   tblout.pushFieldSpec("alpha_m", TO::DOUBLE);
00136   tblout.pushFieldSpec("alpha_u", TO::DOUBLE);
00137   tblout.pushFieldSpec("phi_l", TO::DOUBLE);
00138   tblout.pushFieldSpec("phi_m", TO::DOUBLE);
00139   tblout.pushFieldSpec("phi_u", TO::DOUBLE);
00140   tblout.pushFieldSpec("step type             ", TO::STRING);
00141 
00142   tblout.outputHeader();
00143 
00144   int icount = 0;
00145 
00146   std::string stepType = "";
00147 
00148   //
00149   // A) Find phi_l > phi_m first
00150   //
00151 
00152   tblout.outputField("-");
00153   tblout.outputField(alpha_l);
00154   tblout.outputField(alpha_m);
00155   tblout.outputField("-");
00156   tblout.outputField(phi_l);
00157   tblout.outputField(phi_m);
00158   tblout.outputField("-");
00159   tblout.outputField("start");
00160   tblout.nextRow();
00161 
00162   for (; icount < MAX_TOTAL_ITERS; ++icount) {
00163 
00164     // ToDo: Put in a check for NAN and backtrack if you find it!
00165 
00166     if (phi_l > phi_m) {
00167       break;
00168     }
00169 
00170     stepType = "golden back";
00171     alpha_u = alpha_m;
00172     phi_u = phi_m;
00173     alpha_m = goldinv * (alpha_u + GOLDEN_RATIO*alpha_l);
00174     phi_m = computeValue<Scalar>(phi, alpha_m);
00175 
00176     tblout.outputField(icount);
00177     tblout.outputField(alpha_l);
00178     tblout.outputField(alpha_m);
00179     tblout.outputField(alpha_u);
00180     tblout.outputField(phi_l);
00181     tblout.outputField(phi_m);
00182     tblout.outputField(phi_u);
00183     tblout.outputField(stepType);
00184     tblout.nextRow();
00185 
00186   }
00187   
00188   if (alpha_u == zero) {
00189     // The following factor of gold was reduced to (GOLDEN_RATIO-1) to save
00190     // one function evaluation near convergence.
00191     alpha_u = alpha_m + (GOLDEN_RATIO-1.0) * (alpha_m-alpha_l);
00192     phi_u = computeValue<Scalar>(phi, alpha_u);
00193   }
00194   
00195   //
00196   // B) Quadratic interpolation iterations
00197   //
00198   
00199   bool bracketedMin = false;
00200 
00201   for (; icount < MAX_TOTAL_ITERS; ++icount) {
00202     
00203     if (phi_m < phi_u) {
00204       bracketedMin = true;
00205       break;
00206     }
00207       
00208     // find the extremum alpha_quad of a quadratic model interpolating there
00209     // points
00210     q = (phi_m-phi_l)*(alpha_m-alpha_u);
00211     r = (phi_m-phi_u)*(alpha_m-alpha_l);
00212     
00213     // avoid division by small (q-r) by bounding with signed minimum
00214     tmp = ST::magnitude(q-r);
00215     tmp = (tmp > SMALL_DIV ? tmp : SMALL_DIV);
00216     tmp = (q-r >= 0  ? tmp : -tmp);
00217 
00218     Scalar alpha_quad =
00219       alpha_m - (q*(alpha_m-alpha_u) - r*(alpha_m-alpha_l))/(2.0*tmp);
00220     
00221     // maximum point for which we trust the interpolation
00222     const Scalar alpha_lim = alpha_m + MAX_EXTRAP_FACTOR * (alpha_u-alpha_m);
00223     
00224     // now detect which interval alpha_quad is in and act accordingly
00225     bool skipToNextIter = false;
00226     Scalar phi_quad = ST::nan();
00227     if ( (alpha_m-alpha_quad)*(alpha_quad-alpha_u) > zero ) {  // [alpha_m, alpha_u]
00228       phi_quad = computeValue<Scalar>(phi, alpha_quad);
00229       if (phi_quad < phi_u) {                        // use points [b, alpha_quad, c]
00230         alpha_l = alpha_m;
00231         phi_l = phi_m;
00232         alpha_m = alpha_quad;
00233         phi_m = phi_quad;
00234         skipToNextIter = true;
00235         stepType = "alpha_quad middle";
00236       }
00237       else if (phi_quad > phi_m) {                   // use points [a, b, alpha_quad]
00238         alpha_u = alpha_quad;
00239         phi_u = phi_quad;
00240         skipToNextIter = true;
00241         stepType = "alpha_quad upper";
00242       }
00243       else {
00244         alpha_quad = alpha_u + GOLDEN_RATIO*(alpha_u-alpha_m);
00245         phi_quad = computeValue<Scalar>(phi, alpha_quad);
00246       }
00247     }
00248 
00249     if (!skipToNextIter) {
00250       
00251       if ((alpha_u-alpha_quad)*(alpha_quad-alpha_lim) > zero) {  // [alpha_u, alpha_lim]
00252         phi_quad = computeValue<Scalar>(phi, alpha_quad);
00253         stepType = "[alpha_u, alpha_lim]";
00254         if (phi_quad < phi_u) {
00255           alpha_m = alpha_u;
00256           alpha_u = alpha_quad;
00257           alpha_quad = alpha_u + GOLDEN_RATIO*(alpha_u-alpha_m);
00258           phi_m = phi_u;
00259           phi_u = phi_quad;  
00260           phi_quad = computeValue<Scalar>(phi, alpha_quad);
00261           stepType = "phi_quad < phi_u";
00262         }
00263       }
00264       else if ((alpha_quad-alpha_lim)*(alpha_lim-alpha_u) >= zero ) { // [alpha_lim, inf]
00265         alpha_quad = alpha_lim;
00266         phi_quad = computeValue<Scalar>(phi, alpha_quad);
00267         stepType = "[alpha_lim, inf]";
00268       }
00269       else {                                  // [0,alpha_m]
00270         alpha_quad = alpha_u + GOLDEN_RATIO*(alpha_u-alpha_m);
00271         phi_quad = computeValue<Scalar>(phi, alpha_quad);
00272         stepType = "[0, alpha_m]";
00273       }
00274       
00275       // shift to newest 3 points before loop
00276       alpha_l = alpha_m;
00277       phi_l = phi_m;
00278       alpha_m = alpha_u;
00279       phi_m = phi_u;
00280       alpha_u = alpha_quad;
00281       phi_u = phi_quad;
00282 
00283     }
00284     
00285     tblout.outputField(icount);
00286     tblout.outputField(alpha_l);
00287     tblout.outputField(alpha_m);
00288     tblout.outputField(alpha_u);
00289     tblout.outputField(phi_l);
00290     tblout.outputField(phi_m);
00291     tblout.outputField(phi_u);
00292     tblout.outputField(stepType);
00293     tblout.nextRow();
00294    
00295   }  // end for loop
00296   
00297   if (icount >= MAX_TOTAL_ITERS) {
00298     *out <<"\nExceeded maximum number of iterations!.\n";
00299   }
00300 
00301   if (!is_null(numIters)) {
00302     *numIters = icount;
00303   }
00304 
00305   *out << "\n";
00306  
00307   return bracketedMin;
00308 
00309 }
00310 
00311 
00312 } // namespace GlobiPack
00313 
00314 
00315 #endif // GLOBIPACK_GOLDEN_BRACKET_QUAD_INTERP_DEF_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Friends