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