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