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