GlobiPack Package Browser (Single Doxygen Collection) Version of the Day
GlobiPack_Brents1DMinimization_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_BRENTS_1D_MINIMIZATION_DEF_HPP
00032 #define GLOBIPACK_BRENTS_1D_MINIMIZATION_DEF_HPP
00033 
00034 
00035 #include "GlobiPack_Brents1DMinimization_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 Brents1DMinimization<Scalar>::Brents1DMinimization()
00047   :rel_tol_(Brents1DMinimizationUtils::rel_tol_default),
00048    bracket_tol_(Brents1DMinimizationUtils::bracket_tol_default),
00049    max_iters_(Brents1DMinimizationUtils::max_iters_default)
00050 {}
00051 
00052 
00053 // Overridden from ParameterListAcceptor (simple forwarding functions)
00054 
00055 
00056 template<typename Scalar>
00057 void Brents1DMinimization<Scalar>::setParameterList(RCP<ParameterList> const& paramList)
00058 {
00059   typedef ScalarTraits<Scalar> ST;
00060   namespace BMU = Brents1DMinimizationUtils;
00061   using Teuchos::getParameter;
00062   paramList->validateParametersAndSetDefaults(*this->getValidParameters());
00063   rel_tol_ = getParameter<double>(*paramList, BMU::rel_tol_name);
00064   bracket_tol_ = getParameter<double>(*paramList, BMU::bracket_tol_name);
00065   max_iters_ = getParameter<int>(*paramList, BMU::max_iters_name);
00066   TEUCHOS_ASSERT_INEQUALITY( rel_tol_, >, ST::zero() );
00067   TEUCHOS_ASSERT_INEQUALITY( bracket_tol_, >, ST::zero() );
00068   TEUCHOS_ASSERT_INEQUALITY( max_iters_, >=, 0 );
00069   setMyParamList(paramList);
00070 }
00071 
00072 
00073 template<typename Scalar>
00074 RCP<const ParameterList> Brents1DMinimization<Scalar>::getValidParameters() const
00075 {
00076   namespace BMU = Brents1DMinimizationUtils;
00077   static RCP<const ParameterList> validPL;
00078   if (is_null(validPL)) {
00079     RCP<Teuchos::ParameterList>
00080       pl = Teuchos::rcp(new Teuchos::ParameterList());
00081     pl->set( BMU::rel_tol_name, BMU::rel_tol_default );
00082     pl->set( BMU::bracket_tol_name, BMU::bracket_tol_default );
00083     pl->set( BMU::max_iters_name, BMU::max_iters_default );
00084     validPL = pl;
00085   }
00086   return validPL;
00087 }
00088 
00089 
00090 // Bracket
00091 
00092 
00093 template<typename Scalar>
00094 bool Brents1DMinimization<Scalar>::approxMinimize(
00095   const MeritFunc1DBase<Scalar> &phi,
00096   const PointEval1D<Scalar> &pointLower,
00097   const Ptr<PointEval1D<Scalar> > &pointMiddle,
00098   const 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   using std::min;
00110   using std::max;
00111   
00112 #ifdef TEUCHOS_DEBUG
00113   TEST_FOR_EXCEPT(is_null(pointMiddle));
00114   TEUCHOS_ASSERT_INEQUALITY(pointLower.alpha, <, pointMiddle->alpha);
00115   TEUCHOS_ASSERT_INEQUALITY(pointMiddle->alpha, <, pointUpper.alpha);
00116   TEUCHOS_ASSERT_INEQUALITY(pointLower.phi, !=, PE1D::valNotGiven());
00117   TEUCHOS_ASSERT_INEQUALITY(pointMiddle->phi, !=, PE1D::valNotGiven());
00118   TEUCHOS_ASSERT_INEQUALITY(pointUpper.phi, !=, PE1D::valNotGiven());
00119 #endif
00120 
00121   const RCP<Teuchos::FancyOStream> out = this->getOStream();
00122 
00123   *out << "\nStarting Brent's 1D minimization algorithm ...\n\n";
00124   
00125   TabularOutputter tblout(out);
00126   
00127   tblout.pushFieldSpec("itr", TO::INT);
00128   tblout.pushFieldSpec("alpha_a", TO::DOUBLE);
00129   tblout.pushFieldSpec("alpha_min", TO::DOUBLE);
00130   tblout.pushFieldSpec("alpha_b", TO::DOUBLE);
00131   tblout.pushFieldSpec("phi(alpha_min)", TO::DOUBLE);
00132   tblout.pushFieldSpec("alpha_b - alpha_a", TO::DOUBLE);
00133   tblout.pushFieldSpec("alpha_min - alpha_avg", TO::DOUBLE);
00134   tblout.pushFieldSpec("tol", TO::DOUBLE);
00135 
00136   tblout.outputHeader();
00137   
00138   const Scalar INV_GOLD2=0.3819660112501051518; // (1/golden-ratio)^2
00139   const Scalar TINY = ST::squareroot(ST::eps());
00140   
00141   const Scalar alpha_l = pointLower.alpha, phi_l = pointLower.phi;
00142   Scalar &alpha_m = pointMiddle->alpha, &phi_m = pointMiddle->phi;
00143   const Scalar alpha_u = pointUpper.alpha, phi_u = pointUpper.phi;
00144 
00145   Scalar d = ST::nan();
00146   Scalar e = ST::nan();
00147   Scalar u = ST::nan();
00148   
00149   Scalar phi_w = min(phi_l, phi_u);
00150 
00151   Scalar alpha_v = ST::nan();
00152   Scalar alpha_w = ST::nan();
00153   Scalar phi_v = ST::nan();
00154 
00155   if (phi_w == phi_l){  
00156     alpha_w  = alpha_l;
00157     alpha_v  = alpha_u;
00158     phi_v = phi_u;
00159   }
00160   else {
00161     alpha_w  = alpha_u;
00162     alpha_v  = alpha_l;
00163     phi_v = phi_l;
00164   }
00165 
00166   Scalar alpha_min = alpha_m;
00167   Scalar phi_min = phi_m;
00168   Scalar alpha_a = alpha_l;
00169   Scalar alpha_b = alpha_u;
00170   
00171   bool foundMin = false;
00172 
00173   int iteration = 0;
00174 
00175   for ( ; iteration <= max_iters_; ++iteration) {
00176 
00177     if (iteration < 2)
00178       e = 2.0 * (alpha_b - alpha_a);
00179 
00180     const Scalar alpha_avg = 0.5 *(alpha_a + alpha_b);
00181     const Scalar tol1 = rel_tol_ * ST::magnitude(alpha_min) + TINY;
00182     const Scalar tol2 = 2.0 * tol1;
00183 
00184     const Scalar step_diff = alpha_min - alpha_avg;
00185     const Scalar step_diff_tol = (tol2 + bracket_tol_ * (alpha_b - alpha_a));
00186 
00187     // 2009/02/11: rabartl: Above, I changed from (tol2-0.5*(alpha_b-alpha_a)) which is
00188     // actually in Brent's netlib code!  This gives a negative tolerence when
00189     // the solution alpha_min is near a minimum so you will max out the iters because
00190     // a possitive number can never be smaller than a negative number.  The
00191     // above convergence criteria makes sense to me.
00192 
00193     tblout.outputField(iteration);
00194     tblout.outputField(alpha_a);
00195     tblout.outputField(alpha_min);
00196     tblout.outputField(alpha_b);
00197     tblout.outputField(phi_min);
00198     tblout.outputField(alpha_b - alpha_a);
00199     tblout.outputField(step_diff);
00200     tblout.outputField(step_diff_tol);
00201     tblout.nextRow();
00202 
00203     // If the difference between current point and the middle of the shrinking
00204     // interval [alpha_a, alpha_b] is relatively small, then terminate the
00205     // algorithm.  Also, terminate the algorithm if this difference is small
00206     // relative to the size of alpha.  Does this make sense?  However, don't
00207     // terminate on the very first iteration because we have to take at least
00208     // one step.
00209     if (
00210       ST::magnitude(step_diff) <= step_diff_tol
00211       && iteration > 0
00212       )
00213     {
00214       foundMin = true;
00215       break;
00216     }
00217     // 2009/02/11: rabartl: Above, I added the iteration > 0 condition because
00218     // the original version that I was given would terminate on the first
00219     // first iteration if the initial guess for alpha happened to be too close
00220     // to the midpoint of the bracketing interval!  Is that crazy or what!
00221 
00222     if (ST::magnitude(e) > tol1 || iteration < 2) {
00223 
00224       const Scalar r = (alpha_min - alpha_w) * (phi_min - phi_v);
00225       Scalar q = (alpha_min - alpha_v) * (phi_min - phi_w);
00226       Scalar p = (alpha_min - alpha_v) * q - (alpha_min - alpha_w) * r;
00227       q = 2.0 * (q - r);
00228       if (q > ST::zero())
00229         p = -p;
00230       q = ST::magnitude(q);
00231       const Scalar etemp = e;
00232       e = d;
00233 
00234       if ( ST::magnitude(p) >= ST::magnitude(0.5 * q * etemp)
00235         || p <= q * (alpha_a - alpha_min)
00236         || p >= q * (alpha_b - alpha_min)
00237         )
00238       {
00239         e = (alpha_min >= alpha_avg ? alpha_a - alpha_min : alpha_b - alpha_min);
00240         d = INV_GOLD2 * e;
00241       }
00242       else {
00243         d = p/q;
00244         u = alpha_min + d;
00245         if (u - alpha_a < tol2 || alpha_b - u < tol2) 
00246            // sign(tol1,alpha_avg-alpha_min)
00247           d = ( alpha_avg - alpha_min > ST::zero()
00248             ? ST::magnitude(tol1)
00249             : -ST::magnitude(tol1) );
00250       }
00251 
00252     }
00253     else {
00254 
00255       e = (alpha_min >= alpha_avg ? alpha_a - alpha_min : alpha_b - alpha_min);
00256       d = INV_GOLD2 * e;
00257 
00258     }
00259     
00260     u = ( ST::magnitude(d) >= tol1
00261       ? alpha_min + d
00262       : alpha_min + (d >= 0 ? ST::magnitude(tol1) : -ST::magnitude(tol1))
00263       ); 
00264     
00265     const Scalar phi_eval_u = computeValue<Scalar>(phi, u);
00266 
00267     if (phi_eval_u <= phi_min) {
00268 
00269       if (u >= alpha_min)
00270         alpha_a = alpha_min;
00271       else
00272         alpha_b = alpha_min;
00273 
00274       alpha_v = alpha_w;
00275       phi_v = phi_w;
00276       alpha_w = alpha_min;
00277       phi_w = phi_min;
00278       alpha_min = u;
00279       phi_min = phi_eval_u;
00280 
00281     }
00282     else {
00283 
00284       if (u < alpha_min)
00285         alpha_a = u;
00286       else
00287         alpha_b = u;
00288 
00289       if (phi_eval_u <= phi_w || alpha_w == alpha_min) {
00290         alpha_v = alpha_w;
00291         phi_v = phi_w;
00292         alpha_w = u;
00293         phi_w = phi_eval_u;
00294       }
00295       else if (phi_eval_u <= phi_v || alpha_v == alpha_min || alpha_v == alpha_w) {
00296         alpha_v = u;
00297         phi_v = phi_eval_u;
00298       }
00299 
00300     }
00301   }
00302 
00303   alpha_m = alpha_min;
00304   phi_m = phi_min;
00305   if (!is_null(numIters))
00306     *numIters = iteration;
00307   
00308   if (foundMin) {
00309     *out <<"\nFound the minimum alpha="<<alpha_m<<", phi(alpha)="<<phi_m<<"\n";
00310   }
00311   else {
00312     *out <<"\nExceeded maximum number of iterations!\n";
00313   }
00314 
00315   *out << "\n";
00316 
00317   return foundMin;
00318 
00319 }
00320 
00321 
00322 } // namespace GlobiPack
00323 
00324 
00325 #endif // GLOBIPACK_BRENTS_1D_MINIMIZATION_DEF_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Friends