GlobiPack Package Browser (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
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 
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>
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>
88 {
89  namespace BMU = Brents1DMinimizationUtils;
90  static RCP<const ParameterList> validPL;
91  if (is_null(validPL)) {
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;
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).
bool is_null(const boost::shared_ptr< T > &p)
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
#define TEUCHOS_ASSERT_INEQUALITY(val1, comp, val2)
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.
Brents1DMinimization()
Construct with default parameters.
basic_OSTab< char > OSTab
RCP< const ParameterList > getValidParameters() const
Scalar alpha
The value of the unknown alpha.
Represents the evaluation point of the merit function phi(alpha) and/or is derivative Dphi(alpha)...
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
void setParameterList(RCP< ParameterList > const &paramList)
void validateParametersAndSetDefaults(ParameterList const &validParamList, int const depth=1000)
TypeTo as(const TypeFrom &t)
Base class for 1D merit fucntions used in globalization methods.
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)