GlobiPack Package Browser (Single Doxygen Collection)  Version of the Day
GlobiPack_Brents1DMinimization_def.hpp
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 #ifndef GLOBIPACK_BRENTS_1D_MINIMIZATION_DEF_HPP
45 #define GLOBIPACK_BRENTS_1D_MINIMIZATION_DEF_HPP
46 
47 
49 #include "Teuchos_TabularOutputter.hpp"
50 
51 
52 namespace GlobiPack {
53 
54 
55 // Constructor/Initializers/Accessors
56 
57 
58 template<typename Scalar>
60  :rel_tol_(Brents1DMinimizationUtils::rel_tol_default),
61  bracket_tol_(Brents1DMinimizationUtils::bracket_tol_default),
62  max_iters_(Brents1DMinimizationUtils::max_iters_default)
63 {}
64 
65 
66 // Overridden from ParameterListAcceptor (simple forwarding functions)
67 
68 
69 template<typename Scalar>
70 void Brents1DMinimization<Scalar>::setParameterList(RCP<ParameterList> const& paramList)
71 {
72  typedef ScalarTraits<Scalar> ST;
73  namespace BMU = Brents1DMinimizationUtils;
74  using Teuchos::getParameter;
75  paramList->validateParametersAndSetDefaults(*this->getValidParameters());
76  rel_tol_ = getParameter<double>(*paramList, BMU::rel_tol_name);
77  bracket_tol_ = getParameter<double>(*paramList, BMU::bracket_tol_name);
78  max_iters_ = getParameter<int>(*paramList, BMU::max_iters_name);
79  TEUCHOS_ASSERT_INEQUALITY( rel_tol_, >, ST::zero() );
80  TEUCHOS_ASSERT_INEQUALITY( bracket_tol_, >, ST::zero() );
81  TEUCHOS_ASSERT_INEQUALITY( max_iters_, >=, 0 );
82  setMyParamList(paramList);
83 }
84 
85 
86 template<typename Scalar>
87 RCP<const ParameterList> Brents1DMinimization<Scalar>::getValidParameters() const
88 {
89  namespace BMU = Brents1DMinimizationUtils;
90  static RCP<const ParameterList> validPL;
91  if (is_null(validPL)) {
92  RCP<Teuchos::ParameterList>
93  pl = Teuchos::rcp(new Teuchos::ParameterList());
97  validPL = pl;
98  }
99  return validPL;
100 }
101 
102 
103 // Bracket
104 
105 
106 template<typename Scalar>
108  const MeritFunc1DBase<Scalar> &phi,
109  const PointEval1D<Scalar> &pointLower,
110  const Ptr<PointEval1D<Scalar> > &pointMiddle,
111  const PointEval1D<Scalar> &pointUpper,
112  const Ptr<int> &numIters
113  ) const
114 {
115 
116  using Teuchos::as;
117  using Teuchos::TabularOutputter;
118  typedef Teuchos::TabularOutputter TO;
119  typedef ScalarTraits<Scalar> ST;
120  using Teuchos::OSTab;
121 #ifdef TEUCHOS_DEBUG
122  typedef PointEval1D<Scalar> PE1D;
123 #endif // TEUCHOS_DEBUG
124  using std::min;
125  using std::max;
126 
127 #ifdef TEUCHOS_DEBUG
128  TEUCHOS_TEST_FOR_EXCEPT(is_null(pointMiddle));
129  TEUCHOS_ASSERT_INEQUALITY(pointLower.alpha, <, pointMiddle->alpha);
130  TEUCHOS_ASSERT_INEQUALITY(pointMiddle->alpha, <, pointUpper.alpha);
131  TEUCHOS_ASSERT_INEQUALITY(pointLower.phi, !=, PE1D::valNotGiven());
132  TEUCHOS_ASSERT_INEQUALITY(pointMiddle->phi, !=, PE1D::valNotGiven());
133  TEUCHOS_ASSERT_INEQUALITY(pointUpper.phi, !=, PE1D::valNotGiven());
134 #endif // TEUCHOS_DEBUG
135 
136  const RCP<Teuchos::FancyOStream> out = this->getOStream();
137 
138  *out << "\nStarting Brent's 1D minimization algorithm ...\n\n";
139 
140  TabularOutputter tblout(out);
141 
142  tblout.pushFieldSpec("itr", TO::INT);
143  tblout.pushFieldSpec("alpha_a", TO::DOUBLE);
144  tblout.pushFieldSpec("alpha_min", TO::DOUBLE);
145  tblout.pushFieldSpec("alpha_b", TO::DOUBLE);
146  tblout.pushFieldSpec("phi(alpha_min)", TO::DOUBLE);
147  tblout.pushFieldSpec("alpha_b - alpha_a", TO::DOUBLE);
148  tblout.pushFieldSpec("alpha_min - alpha_avg", TO::DOUBLE);
149  tblout.pushFieldSpec("tol", TO::DOUBLE);
150 
151  tblout.outputHeader();
152 
153  const Scalar INV_GOLD2=0.3819660112501051518; // (1/golden-ratio)^2
154  const Scalar TINY = ST::squareroot(ST::eps());
155 
156  const Scalar alpha_l = pointLower.alpha, phi_l = pointLower.phi;
157  Scalar &alpha_m = pointMiddle->alpha, &phi_m = pointMiddle->phi;
158  const Scalar alpha_u = pointUpper.alpha, phi_u = pointUpper.phi;
159 
160  Scalar d = ST::nan();
161  Scalar e = ST::nan();
162  Scalar u = ST::nan();
163 
164  Scalar phi_w = min(phi_l, phi_u);
165 
166  Scalar alpha_v = ST::nan();
167  Scalar alpha_w = ST::nan();
168  Scalar phi_v = ST::nan();
169 
170  if (phi_w == phi_l){
171  alpha_w = alpha_l;
172  alpha_v = alpha_u;
173  phi_v = phi_u;
174  }
175  else {
176  alpha_w = alpha_u;
177  alpha_v = alpha_l;
178  phi_v = phi_l;
179  }
180 
181  Scalar alpha_min = alpha_m;
182  Scalar phi_min = phi_m;
183  Scalar alpha_a = alpha_l;
184  Scalar alpha_b = alpha_u;
185 
186  bool foundMin = false;
187 
188  int iteration = 0;
189 
190  for ( ; iteration <= max_iters_; ++iteration) {
191 
192  if (iteration < 2)
193  e = 2.0 * (alpha_b - alpha_a);
194 
195  const Scalar alpha_avg = 0.5 *(alpha_a + alpha_b);
196  const Scalar tol1 = rel_tol_ * ST::magnitude(alpha_min) + TINY;
197  const Scalar tol2 = 2.0 * tol1;
198 
199  const Scalar step_diff = alpha_min - alpha_avg;
200  const Scalar step_diff_tol = (tol2 + bracket_tol_ * (alpha_b - alpha_a));
201 
202  // 2009/02/11: rabartl: Above, I changed from (tol2-0.5*(alpha_b-alpha_a)) which is
203  // actually in Brent's netlib code! This gives a negative tolerence when
204  // the solution alpha_min is near a minimum so you will max out the iters because
205  // a possitive number can never be smaller than a negative number. The
206  // above convergence criteria makes sense to me.
207 
208  tblout.outputField(iteration);
209  tblout.outputField(alpha_a);
210  tblout.outputField(alpha_min);
211  tblout.outputField(alpha_b);
212  tblout.outputField(phi_min);
213  tblout.outputField(alpha_b - alpha_a);
214  tblout.outputField(step_diff);
215  tblout.outputField(step_diff_tol);
216  tblout.nextRow();
217 
218  // If the difference between current point and the middle of the shrinking
219  // interval [alpha_a, alpha_b] is relatively small, then terminate the
220  // algorithm. Also, terminate the algorithm if this difference is small
221  // relative to the size of alpha. Does this make sense? However, don't
222  // terminate on the very first iteration because we have to take at least
223  // one step.
224  if (
225  ST::magnitude(step_diff) <= step_diff_tol
226  && iteration > 0
227  )
228  {
229  foundMin = true;
230  break;
231  }
232  // 2009/02/11: rabartl: Above, I added the iteration > 0 condition because
233  // the original version that I was given would terminate on the first
234  // first iteration if the initial guess for alpha happened to be too close
235  // to the midpoint of the bracketing interval! Is that crazy or what!
236 
237  if (ST::magnitude(e) > tol1 || iteration < 2) {
238 
239  const Scalar r = (alpha_min - alpha_w) * (phi_min - phi_v);
240  Scalar q = (alpha_min - alpha_v) * (phi_min - phi_w);
241  Scalar p = (alpha_min - alpha_v) * q - (alpha_min - alpha_w) * r;
242  q = 2.0 * (q - r);
243  if (q > ST::zero())
244  p = -p;
245  q = ST::magnitude(q);
246  const Scalar etemp = e;
247  e = d;
248 
249  if ( ST::magnitude(p) >= ST::magnitude(0.5 * q * etemp)
250  || p <= q * (alpha_a - alpha_min)
251  || p >= q * (alpha_b - alpha_min)
252  )
253  {
254  e = (alpha_min >= alpha_avg ? alpha_a - alpha_min : alpha_b - alpha_min);
255  d = INV_GOLD2 * e;
256  }
257  else {
258  d = p/q;
259  u = alpha_min + d;
260  if (u - alpha_a < tol2 || alpha_b - u < tol2)
261  // sign(tol1,alpha_avg-alpha_min)
262  d = ( alpha_avg - alpha_min > ST::zero()
263  ? ST::magnitude(tol1)
264  : -ST::magnitude(tol1) );
265  }
266 
267  }
268  else {
269 
270  e = (alpha_min >= alpha_avg ? alpha_a - alpha_min : alpha_b - alpha_min);
271  d = INV_GOLD2 * e;
272 
273  }
274 
275  u = ( ST::magnitude(d) >= tol1
276  ? alpha_min + d
277  : alpha_min + (d >= 0 ? ST::magnitude(tol1) : -ST::magnitude(tol1))
278  );
279 
280  const Scalar phi_eval_u = computeValue<Scalar>(phi, u);
281 
282  if (phi_eval_u <= phi_min) {
283 
284  if (u >= alpha_min)
285  alpha_a = alpha_min;
286  else
287  alpha_b = alpha_min;
288 
289  alpha_v = alpha_w;
290  phi_v = phi_w;
291  alpha_w = alpha_min;
292  phi_w = phi_min;
293  alpha_min = u;
294  phi_min = phi_eval_u;
295 
296  }
297  else {
298 
299  if (u < alpha_min)
300  alpha_a = u;
301  else
302  alpha_b = u;
303 
304  if (phi_eval_u <= phi_w || alpha_w == alpha_min) {
305  alpha_v = alpha_w;
306  phi_v = phi_w;
307  alpha_w = u;
308  phi_w = phi_eval_u;
309  }
310  else if (phi_eval_u <= phi_v || alpha_v == alpha_min || alpha_v == alpha_w) {
311  alpha_v = u;
312  phi_v = phi_eval_u;
313  }
314 
315  }
316  }
317 
318  alpha_m = alpha_min;
319  phi_m = phi_min;
320  if (!is_null(numIters))
321  *numIters = iteration;
322 
323  if (foundMin) {
324  *out <<"\nFound the minimum alpha="<<alpha_m<<", phi(alpha)="<<phi_m<<"\n";
325  }
326  else {
327  *out <<"\nExceeded maximum number of iterations!\n";
328  }
329 
330  *out << "\n";
331 
332  return foundMin;
333 
334 }
335 
336 
337 } // namespace GlobiPack
338 
339 
340 #endif // GLOBIPACK_BRENTS_1D_MINIMIZATION_DEF_HPP
Scalar phi
The value of the merit function phi(alpha).
Brents1DMinimization()
Construct with default parameters.
Scalar alpha
The value of the unknown alpha.
Represents the evaluation point of the merit function phi(alpha) and/or is derivative Dphi(alpha)...
void setParameterList(RCP< ParameterList > const &paramList)
Base class for 1D merit fucntions used in globalization methods.
RCP< const ParameterList > getValidParameters() const
bool approxMinimize(const MeritFunc1DBase< Scalar > &phi, const PointEval1D< Scalar > &pointLower, const Ptr< PointEval1D< Scalar > > &pointMiddle, const PointEval1D< Scalar > &pointUpper, const Ptr< int > &numIters=Teuchos::null) const
Approximatly mimimize a 1D function.