GlobiPack Package Browser (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
Brents1DMinimization_UnitTests.cpp
Go to the documentation of this file.
1 /*
2 // @HEADER
3 // ***********************************************************************
4 //
5 // GlobiPack: Collection of Scalar 1D globalizaton utilities
6 // Copyright (2009) Sandia Corporation
7 //
8 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
9 // license for use of this work by or on behalf of the U.S. Government.
10 //
11 // Redistribution and use in source and binary forms, with or without
12 // modification, are permitted provided that the following conditions are
13 // met:
14 //
15 // 1. Redistributions of source code must retain the above copyright
16 // notice, this list of conditions and the following disclaimer.
17 //
18 // 2. Redistributions in binary form must reproduce the above copyright
19 // notice, this list of conditions and the following disclaimer in the
20 // documentation and/or other materials provided with the distribution.
21 //
22 // 3. Neither the name of the Corporation nor the names of the
23 // contributors may be used to endorse or promote products derived from
24 // this software without specific prior written permission.
25 //
26 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37 //
38 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov)
39 //
40 // ***********************************************************************
41 // @HEADER
42 */
43 
44 
45 #include "GlobiPack_Brents1DMinimization.hpp"
46 #include "GlobiPack_TestLagrPolyMeritFunc1D.hpp"
47 #include "Teuchos_Tuple.hpp"
49 
50 
51 namespace {
52 
53 
54 //
55 // Helper code and declarations
56 //
57 
58 
62 using GlobiPack::brents1DMinimization;
64 using GlobiPack::computePoint;
65 using Teuchos::as;
66 using Teuchos::inOutArg;
67 using Teuchos::outArg;
68 using Teuchos::null;
69 using Teuchos::RCP;
70 using Teuchos::rcpFromRef;
71 using Teuchos::Array;
72 using Teuchos::ArrayView;
73 using Teuchos::tuple;
75 using Teuchos::parameterList;
76 using Teuchos::OSTab;
77 
78 
79 template<class Scalar>
80 inline Scalar sqr(const Scalar &x) { return x*x; }
81 
82 
83 template<class Scalar>
84 inline Scalar cube(const Scalar &x) { return x*x*x; }
85 
86 
87 // Set up a quadratic merit function with minimizer at alpha=2.0, phi=3.0;
88 template<class Scalar>
89 const RCP<TestLagrPolyMeritFunc1D<Scalar> > quadPhi()
90 {
92  typedef typename ST::magnitudeType ScalarMag;
93  Array<Scalar> alphaPoints = tuple<Scalar>(0.0, 2.0, 4.0);
94  Array<ScalarMag> phiPoints = tuple<ScalarMag>(6.0, 3.0, 6.0);
95  return testLagrPolyMeritFunc1D<Scalar>(alphaPoints, phiPoints);
96 }
97 
98 //
99 // Set up a cubic merit function with minimizer at alpha=2.0, phi=3.0;
100 //
101 // The function being represented approximated is:
102 //
103 // phi(alpha) = (alpha - 2.0)^2 + 1e-3 * (alpha - 2.0)^3 + 3.0
104 //
105 // This function has the first and second derivatives derivatives:
106 //
107 // Dphi(alpha) = 2.0 * (alpha - 2.0) + 3e-3 * (alpha - 2.0)^2
108 //
109 // D2phi(alpha) = 2.0 + 6e-3 * (alpha - 2.0)
110 //
111 // At alpha=2.0, the function has Dphi=0.0 and D2phi = 2.0 and therefore, this
112 // is a local minimum.
113 //
114 
115 const double cubicMut = 1e-3;
116 
117 template<class Scalar>
118 inline Scalar cubicPhiVal(const Scalar &alpha)
119 { return sqr(alpha - 2.0) + cubicMut * cube(alpha - 2.0) + 3.0; }
120 
121 
122 template<class Scalar>
123 const RCP<TestLagrPolyMeritFunc1D<Scalar> > cubicPhi()
124 {
126  typedef typename ST::magnitudeType ScalarMag;
127  Array<Scalar> alphaPoints =
128  tuple<Scalar>(0.0, 1.0, 3.0, 4.0);
129  Array<ScalarMag> phiPoints =
130  tuple<ScalarMag>(
131  cubicPhiVal(alphaPoints[0]),
132  cubicPhiVal(alphaPoints[1]),
133  cubicPhiVal(alphaPoints[2]),
134  cubicPhiVal(alphaPoints[3])
135  );
136  return testLagrPolyMeritFunc1D<Scalar>(alphaPoints, phiPoints);
137 }
138 
139 
140 double g_tol_scale = 100.0;
141 
142 
144 {
146  "tol", &g_tol_scale, "Floating point tolerance scaling of eps." );
147 }
148 
149 
150 //
151 // Unit tests for Brents1DMinimization
152 //
153 
154 
155 //
156 // Check that object can exactly interplate a quadratic merit function. This
157 // takes more than one iteration due to the golden search bracketing algorithm
158 // which takes time to terminate.
159 //
160 
161 TEUCHOS_UNIT_TEST_TEMPLATE_1_DECL( Brents1DMinimization, quadExact, Scalar )
162 {
163 
165 
166  const RCP<TestLagrPolyMeritFunc1D<Scalar> > phi = quadPhi<Scalar>();
167 
168  RCP<Brents1DMinimization<Scalar> > brentsMin = brents1DMinimization<Scalar>();
169 
170  // Set up to do one iteration and that is it. With the quadratic
171  // interplation that should be enough.
172  const RCP<ParameterList> pl = parameterList();
173  //pl->set("Relative Tol", 1.0);
174  //pl->set("Bracket Tol", 1.0);
175  //pl->set("Max Iterations", 3);
176  brentsMin->setParameterList(pl);
177 
178  brentsMin->setOStream(rcpFromRef(out));
179 
180  const Array<Array<double> > brackets =
181  tuple<Array<double> >(
182  tuple(0.0, 2.0, 4.0), // This is the midpoint and the solution!
183  tuple(0.5, 2.5, 4.5), // This is the midpoint but *not* the solution!
184  tuple(0.0, 1.0, 3.0),
185  tuple(1.0, 3.0, 4.0),
186  tuple(1.9, 2.0, 4.0),
187  tuple(1.9, 3.9, 4.0),
188  tuple(0.0, 2.0, 2.1),
189  tuple(0.0, 0.1, 2.1)
190  );
191 
192  for (int i = 0; i < as<int>(brackets.size()); ++i ) {
193 
194  const ArrayView<const double> bracket = brackets[i]();
195 
196  out << "\ni = "<<i<<": bracket = "<<bracket()<<"\n";
197 
198  OSTab tab(out);
199 
200  PointEval1D<Scalar> p_l = computePoint<Scalar>(*phi, bracket[0]);
201  PointEval1D<Scalar> p_m = computePoint<Scalar>(*phi, bracket[1]);
202  PointEval1D<Scalar> p_u = computePoint<Scalar>(*phi, bracket[2]);
203  int numIters = -1;
204 
205  const bool mimimized = brentsMin->approxMinimize(
206  *phi, p_l, inOutArg(p_m), p_u, outArg(numIters) );
207 
208  TEST_ASSERT(mimimized);
209  //TEST_EQUALITY_CONST(numIters, 1);
210  TEST_FLOATING_EQUALITY(p_m.alpha, as<Scalar>(2.0),
211  as<Scalar>(g_tol_scale*ST::squareroot(ST::eps())));
212  TEST_FLOATING_EQUALITY(p_m.phi, as<Scalar>(3.0),
213  as<Scalar>(g_tol_scale)*ST::eps());
214  TEST_COMPARE(p_l.alpha, <=, p_m.alpha);
215  TEST_COMPARE(p_m.alpha, <=, p_u.alpha);
216  TEST_COMPARE(p_m.phi, <=, p_l.phi);
217  TEST_COMPARE(p_m.phi, <=, p_u.phi);
218 
219  }
220 
221 }
222 
223 TEUCHOS_UNIT_TEST_TEMPLATE_1_INSTANT_REAL_SCALAR_TYPES( Brents1DMinimization, quadExact )
224 
225 
226 //
227 // Check that object can approximately mimimize a cubic function.
228 //
229 
230 TEUCHOS_UNIT_TEST_TEMPLATE_1_DECL( Brents1DMinimization, cubicApprox, Scalar )
231 {
232 
234 
235  const RCP<TestLagrPolyMeritFunc1D<Scalar> > phi = cubicPhi<Scalar>();
236 
237  RCP<Brents1DMinimization<Scalar> > brentsMin = brents1DMinimization<Scalar>();
238 
239  // Set up to do one iteration and that is it. With the quadratic
240  // interplation that should be enough.
241  const RCP<ParameterList> pl = parameterList();
242  pl->set("Relative Tol", as<double>(g_tol_scale*ST::eps()));
243  pl->set("Bracket Tol", as<double>(ST::eps()));
244  brentsMin->setParameterList(pl);
245 
246  brentsMin->setOStream(rcpFromRef(out));
247 
248  const Array<Array<double> > brackets =
249  tuple<Array<double> >(
250  tuple(0.0, 2.0, 4.0), // This is the midpoint and the solution!
251  tuple(0.5, 2.5, 4.5), // This is the midpoint but *not* the solution!
252  tuple(0.0, 1.0, 3.0),
253  tuple(1.0, 3.0, 4.0),
254  tuple(1.9, 2.0, 4.0),
255  tuple(1.9, 3.9, 4.0),
256  tuple(0.0, 2.0, 2.1),
257  tuple(0.0, 0.1, 2.1)
258  );
259 
260  for (int i = 0; i < as<int>(brackets.size()); ++i ) {
261 
262  const ArrayView<const double> bracket = brackets[i]();
263 
264  out << "\ni = "<<i<<": bracket = "<<bracket()<<"\n";
265 
266  OSTab tab(out);
267 
268  PointEval1D<Scalar> p_l = computePoint<Scalar>(*phi, bracket[0]);
269  PointEval1D<Scalar> p_m = computePoint<Scalar>(*phi, bracket[1]);
270  PointEval1D<Scalar> p_u = computePoint<Scalar>(*phi, bracket[2]);
271  int numIters = -1;
272 
273  const bool mimimized = brentsMin->approxMinimize(
274  *phi, p_l, inOutArg(p_m), p_u, outArg(numIters) );
275 
276  TEST_ASSERT(mimimized);
277  //TEST_EQUALITY_CONST(numIters, 1);
278  TEST_FLOATING_EQUALITY(p_m.alpha, as<Scalar>(2.0),
279  as<Scalar>(g_tol_scale*ST::squareroot(ST::eps())));
280  TEST_FLOATING_EQUALITY(p_m.phi, as<Scalar>(3.0),
281  as<Scalar>(g_tol_scale*ST::eps()));
282  TEST_COMPARE(p_l.alpha, <=, p_m.alpha);
283  TEST_COMPARE(p_m.alpha, <=, p_u.alpha);
284  TEST_COMPARE(p_m.phi, <=, p_l.phi);
285  TEST_COMPARE(p_m.phi, <=, p_u.phi);
286 
287  }
288 
289 }
290 
291 TEUCHOS_UNIT_TEST_TEMPLATE_1_INSTANT_REAL_SCALAR_TYPES( Brents1DMinimization, cubicApprox )
292 
293 
294 } // namespace
TEUCHOS_UNIT_TEST_TEMPLATE_1_DECL(Teuchos_Conditions, NumberConditionSerialization, T)
const RCP< TestLagrPolyMeritFunc1D< Scalar > > testLagrPolyMeritFunc1D(const ArrayView< const Scalar > &alpha, const ArrayView< const Scalar > &phi)
Simple concrete class that implements a 1D algorithm to mimimize a 1D function.
static CommandLineProcessor & getCLP()
basic_OSTab< char > OSTab
Represents the evaluation point of the merit function phi(alpha) and/or is derivative Dphi(alpha)...
void setOption(const char option_true[], const char option_false[], bool *option_val, const char documentation[]=NULL)
TypeTo as(const TypeFrom &t)
#define TEST_FLOATING_EQUALITY(v1, v2, tol)
TEST_ASSERT(castedDep1->getValuesAndValidators().size()==2)
#define TEUCHOS_UNIT_TEST_TEMPLATE_1_INSTANT_REAL_SCALAR_TYPES(TEST_GROUP, TEST_NAME)
Lagrange Polynomial Merit Function used in testing.
TEUCHOS_STATIC_SETUP()
#define TEST_COMPARE(v1, comp, v2)