GlobiPack Package Browser (Single Doxygen Collection) Version of the Day
Brents1DMinimization_UnitTests.cpp
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 
00032 #include "GlobiPack_Brents1DMinimization.hpp"
00033 #include "GlobiPack_TestLagrPolyMeritFunc1D.hpp"
00034 #include "Teuchos_Tuple.hpp"
00035 #include "Teuchos_UnitTestHarness.hpp"
00036 
00037 
00038 namespace {
00039 
00040 
00041 //
00042 // Helper code and declarations
00043 //
00044 
00045 
00046 using GlobiPack::TestLagrPolyMeritFunc1D;
00047 using GlobiPack::testLagrPolyMeritFunc1D;
00048 using GlobiPack::Brents1DMinimization;
00049 using GlobiPack::brents1DMinimization;
00050 using GlobiPack::PointEval1D;
00051 using GlobiPack::computePoint;
00052 using Teuchos::as;
00053 using Teuchos::inOutArg;
00054 using Teuchos::outArg;
00055 using Teuchos::null;
00056 using Teuchos::RCP;
00057 using Teuchos::rcpFromRef;
00058 using Teuchos::Array;
00059 using Teuchos::ArrayView;
00060 using Teuchos::tuple;
00061 using Teuchos::ParameterList;
00062 using Teuchos::parameterList;
00063 using Teuchos::OSTab;
00064 
00065 
00066 template<class Scalar>
00067 inline Scalar sqr(const Scalar &x) { return x*x; }
00068 
00069 
00070 template<class Scalar>
00071 inline Scalar cube(const Scalar &x) { return x*x*x; }
00072 
00073 
00074 // Set up a quadratic merit function with minimizer at alpha=2.0, phi=3.0;
00075 template<class Scalar>
00076 const RCP<TestLagrPolyMeritFunc1D<Scalar> > quadPhi()
00077 {
00078   typedef Teuchos::ScalarTraits<Scalar> ST;
00079   typedef typename ST::magnitudeType ScalarMag;
00080   Array<Scalar> alphaPoints = tuple<Scalar>(0.0, 2.0, 4.0);
00081   Array<ScalarMag> phiPoints = tuple<ScalarMag>(6.0, 3.0, 6.0);
00082   return testLagrPolyMeritFunc1D<Scalar>(alphaPoints, phiPoints);
00083 }
00084 
00085 //
00086 // Set up a cubic merit function with minimizer at alpha=2.0, phi=3.0;
00087 //
00088 // The function being represented approximated is:
00089 //
00090 //   phi(alpha) = (alpha - 2.0)^2 + 1e-3 * (alpha - 2.0)^3 + 3.0
00091 //
00092 // This function has the first and second derivatives derivatives:
00093 //
00094 //   Dphi(alpha) = 2.0 * (alpha - 2.0) + 3e-3 * (alpha - 2.0)^2
00095 //
00096 //   D2phi(alpha) = 2.0 + 6e-3 * (alpha - 2.0)
00097 //
00098 // At alpha=2.0, the function has Dphi=0.0 and D2phi = 2.0 and therefore, this
00099 // is a local minimum.
00100 //
00101 
00102 const double cubicMut = 1e-3;
00103 
00104 template<class Scalar>
00105 inline Scalar cubicPhiVal(const Scalar &alpha)
00106 { return sqr(alpha - 2.0) + cubicMut * cube(alpha - 2.0) + 3.0; }
00107 
00108 
00109 template<class Scalar>
00110 const RCP<TestLagrPolyMeritFunc1D<Scalar> > cubicPhi()
00111 {
00112   typedef Teuchos::ScalarTraits<Scalar> ST;
00113   typedef typename ST::magnitudeType ScalarMag;
00114   Array<Scalar> alphaPoints =
00115     tuple<Scalar>(0.0, 1.0, 3.0, 4.0);
00116   Array<ScalarMag> phiPoints =
00117     tuple<ScalarMag>(
00118       cubicPhiVal(alphaPoints[0]),
00119       cubicPhiVal(alphaPoints[1]),
00120       cubicPhiVal(alphaPoints[2]),
00121       cubicPhiVal(alphaPoints[3])
00122       );
00123   return testLagrPolyMeritFunc1D<Scalar>(alphaPoints, phiPoints);
00124 }
00125 
00126 
00127 double g_tol_scale = 100.0;
00128 
00129 
00130 TEUCHOS_STATIC_SETUP()
00131 {
00132   Teuchos::UnitTestRepository::getCLP().setOption(
00133     "tol", &g_tol_scale, "Floating point tolerance scaling of eps." );
00134 }
00135 
00136 
00137 //
00138 // Unit tests for Brents1DMinimization
00139 //
00140 
00141 
00142 //
00143 // Check that object can exactly interplate a quadratic merit function.  This
00144 // takes more than one iteration due to the golden search bracketing algorithm
00145 // which takes time to terminate.
00146 //
00147 
00148 TEUCHOS_UNIT_TEST_TEMPLATE_1_DECL( Brents1DMinimization, quadExact, Scalar )
00149 {
00150 
00151   typedef Teuchos::ScalarTraits<Scalar> ST;
00152 
00153   const RCP<TestLagrPolyMeritFunc1D<Scalar> > phi = quadPhi<Scalar>();
00154   
00155   RCP<Brents1DMinimization<Scalar> > brentsMin = brents1DMinimization<Scalar>();
00156 
00157   // Set up to do one iteration and that is it.  With the quadratic
00158   // interplation that should be enough.
00159   const RCP<ParameterList> pl = parameterList();
00160   //pl->set("Relative Tol", 1.0);
00161   //pl->set("Bracket Tol", 1.0);
00162   //pl->set("Max Iterations", 3);
00163   brentsMin->setParameterList(pl);
00164 
00165   brentsMin->setOStream(rcpFromRef(out));
00166 
00167   const Array<Array<double> > brackets =
00168     tuple<Array<double> >(
00169       tuple(0.0, 2.0, 4.0), // This is the midpoint and the solution!
00170       tuple(0.5, 2.5, 4.5), // This is the midpoint but *not* the solution!
00171       tuple(0.0, 1.0, 3.0),
00172       tuple(1.0, 3.0, 4.0),
00173       tuple(1.9, 2.0, 4.0),
00174       tuple(1.9, 3.9, 4.0),
00175       tuple(0.0, 2.0, 2.1),
00176       tuple(0.0, 0.1, 2.1)
00177       );
00178 
00179   for (int i = 0; i < as<int>(brackets.size()); ++i ) {
00180 
00181     const ArrayView<const double> bracket = brackets[i]();
00182 
00183     out << "\ni = "<<i<<": bracket = "<<bracket()<<"\n";
00184 
00185     OSTab tab(out);
00186     
00187     PointEval1D<Scalar> p_l = computePoint<Scalar>(*phi, bracket[0]);
00188     PointEval1D<Scalar> p_m = computePoint<Scalar>(*phi, bracket[1]);
00189     PointEval1D<Scalar> p_u = computePoint<Scalar>(*phi, bracket[2]);
00190     int numIters = -1;
00191     
00192     const bool mimimized = brentsMin->approxMinimize(
00193       *phi, p_l, inOutArg(p_m), p_u, outArg(numIters) );
00194     
00195     TEST_ASSERT(mimimized);
00196     //TEST_EQUALITY_CONST(numIters, 1);
00197     TEST_FLOATING_EQUALITY(p_m.alpha, as<Scalar>(2.0),
00198       as<Scalar>(g_tol_scale*ST::squareroot(ST::eps())));
00199     TEST_FLOATING_EQUALITY(p_m.phi, as<Scalar>(3.0),
00200       as<Scalar>(g_tol_scale)*ST::eps());
00201     TEST_COMPARE(p_l.alpha, <=, p_m.alpha);
00202     TEST_COMPARE(p_m.alpha, <=, p_u.alpha);
00203     TEST_COMPARE(p_m.phi, <=, p_l.phi);
00204     TEST_COMPARE(p_m.phi, <=, p_u.phi);
00205     
00206   }
00207   
00208 }
00209 
00210 TEUCHOS_UNIT_TEST_TEMPLATE_1_INSTANT_REAL_SCALAR_TYPES( Brents1DMinimization, quadExact )
00211 
00212 
00213 //
00214 // Check that object can approximately mimimize a cubic function.
00215 //
00216 
00217 TEUCHOS_UNIT_TEST_TEMPLATE_1_DECL( Brents1DMinimization, cubicApprox, Scalar )
00218 {
00219 
00220   typedef Teuchos::ScalarTraits<Scalar> ST;
00221 
00222   const RCP<TestLagrPolyMeritFunc1D<Scalar> > phi = cubicPhi<Scalar>();
00223   
00224   RCP<Brents1DMinimization<Scalar> > brentsMin = brents1DMinimization<Scalar>();
00225 
00226   // Set up to do one iteration and that is it.  With the quadratic
00227   // interplation that should be enough.
00228   const RCP<ParameterList> pl = parameterList();
00229   pl->set("Relative Tol", as<double>(g_tol_scale*ST::eps()));
00230   pl->set("Bracket Tol", as<double>(ST::eps()));
00231   brentsMin->setParameterList(pl);
00232 
00233   brentsMin->setOStream(rcpFromRef(out));
00234 
00235   const Array<Array<double> > brackets =
00236     tuple<Array<double> >(
00237       tuple(0.0, 2.0, 4.0), // This is the midpoint and the solution!
00238       tuple(0.5, 2.5, 4.5), // This is the midpoint but *not* the solution!
00239       tuple(0.0, 1.0, 3.0),
00240       tuple(1.0, 3.0, 4.0),
00241       tuple(1.9, 2.0, 4.0),
00242       tuple(1.9, 3.9, 4.0),
00243       tuple(0.0, 2.0, 2.1),
00244       tuple(0.0, 0.1, 2.1)
00245       );
00246 
00247   for (int i = 0; i < as<int>(brackets.size()); ++i ) {
00248 
00249     const ArrayView<const double> bracket = brackets[i]();
00250 
00251     out << "\ni = "<<i<<": bracket = "<<bracket()<<"\n";
00252 
00253     OSTab tab(out);
00254     
00255     PointEval1D<Scalar> p_l = computePoint<Scalar>(*phi, bracket[0]);
00256     PointEval1D<Scalar> p_m = computePoint<Scalar>(*phi, bracket[1]);
00257     PointEval1D<Scalar> p_u = computePoint<Scalar>(*phi, bracket[2]);
00258     int numIters = -1;
00259     
00260     const bool mimimized = brentsMin->approxMinimize(
00261       *phi, p_l, inOutArg(p_m), p_u, outArg(numIters) );
00262     
00263     TEST_ASSERT(mimimized);
00264     //TEST_EQUALITY_CONST(numIters, 1);
00265     TEST_FLOATING_EQUALITY(p_m.alpha, as<Scalar>(2.0),
00266       as<Scalar>(g_tol_scale*ST::squareroot(ST::eps())));
00267     TEST_FLOATING_EQUALITY(p_m.phi, as<Scalar>(3.0),
00268       as<Scalar>(g_tol_scale*ST::eps()));
00269     TEST_COMPARE(p_l.alpha, <=, p_m.alpha);
00270     TEST_COMPARE(p_m.alpha, <=, p_u.alpha);
00271     TEST_COMPARE(p_m.phi, <=, p_l.phi);
00272     TEST_COMPARE(p_m.phi, <=, p_u.phi);
00273     
00274   }
00275   
00276 }
00277 
00278 TEUCHOS_UNIT_TEST_TEMPLATE_1_INSTANT_REAL_SCALAR_TYPES( Brents1DMinimization, cubicApprox )
00279 
00280 
00281 } // namespace
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Friends