GlobiPack Package Browser (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
GlobiPack_ArmijoPolyInterpLineSearch_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_POLY_INTERP_LINE_SEARCH_DEF_HPP
45 #define GLOBIPACK_POLY_INTERP_LINE_SEARCH_DEF_HPP
46 
47 
50 
51 
52 namespace GlobiPack {
53 
54 
55 // Constructor/Initializers/Accessors
56 
57 
58 template<typename Scalar>
60  : eta_(ArmijoPolyInterpLineSearchUtils::eta_default),
61  minFrac_(ArmijoPolyInterpLineSearchUtils::minFrac_default),
62  maxFrac_(ArmijoPolyInterpLineSearchUtils::maxFrac_default),
63  minIters_(ArmijoPolyInterpLineSearchUtils::minIters_default),
64  maxIters_(ArmijoPolyInterpLineSearchUtils::maxIters_default),
65  doMaxIters_(ArmijoPolyInterpLineSearchUtils::doMaxIters_default)
66 {}
67 
68 
69 template<typename Scalar>
71 {
72  return eta_;
73 }
74 
75 
76 template<typename Scalar>
78 {
79  return minFrac_;
80 }
81 
82 
83 template<typename Scalar>
85 {
86  return maxFrac_;
87 }
88 
89 
90 template<typename Scalar>
92 {
93  return minIters_;
94 }
95 
96 
97 template<typename Scalar>
99 {
100  return maxIters_;
101 }
102 
103 
104 template<typename Scalar>
106 {
107  return doMaxIters_;
108 }
109 
110 
111 // Overridden from ParameterListAcceptor
112 
113 
114 template<class Scalar>
116  RCP<ParameterList> const& paramList
117  )
118 {
119  typedef ScalarTraits<Scalar> ST;
120  namespace AQLSU = ArmijoPolyInterpLineSearchUtils;
121  using Teuchos::getParameter;
122  paramList->validateParametersAndSetDefaults(*this->getValidParameters());
123  eta_ = getParameter<double>(*paramList, AQLSU::eta_name);
124  minFrac_ = getParameter<double>(*paramList, AQLSU::minFrac_name);
125  maxFrac_ = getParameter<double>(*paramList, AQLSU::maxFrac_name);
126  minIters_ = getParameter<int>(*paramList, AQLSU::minIters_name);
127  maxIters_ = getParameter<int>(*paramList, AQLSU::maxIters_name);
128  doMaxIters_ = getParameter<bool>(*paramList, AQLSU::doMaxIters_name);
129  TEUCHOS_ASSERT_INEQUALITY( eta_, >=, ST::zero() );
130  TEUCHOS_ASSERT_INEQUALITY( eta_, <, ST::one() );
131  TEUCHOS_ASSERT_INEQUALITY( minFrac_, >=, ST::zero() );
132  TEUCHOS_ASSERT_INEQUALITY( minFrac_, <, maxFrac_ );
133  TEUCHOS_ASSERT_INEQUALITY( minIters_, >=, 0 );
134  TEUCHOS_ASSERT_INEQUALITY( minIters_, <=, maxIters_ );
135  setMyParamList(paramList);
136 }
137 
138 
139 template<class Scalar>
142 {
143  namespace AQLSU = ArmijoPolyInterpLineSearchUtils;
144  static RCP<const ParameterList> validPL;
145  if (is_null(validPL)) {
154  validPL = pl;
155  }
156  return validPL;
157 }
158 
159 
160 // Overrridden from LineSearchBase
161 
162 
163 template<typename Scalar>
165 {
166  return true;
167 }
168 
169 
170 template<typename Scalar>
172 {
173  return false;
174 }
175 
176 
177 template<typename Scalar>
179  const MeritFunc1DBase<Scalar> &phi,
180  const PointEval1D<Scalar> &point_k,
181  const Ptr<PointEval1D<Scalar> > &point_kp1,
182  const Ptr<int> &numIters_out
183  ) const
184 {
185 
186  using Teuchos::as;
187  using Teuchos::ptrFromRef;
189  typedef Teuchos::TabularOutputter TO;
190  typedef ScalarTraits<Scalar> ST;
191 #ifdef TEUCHOS_DEBUG
192  typedef PointEval1D<Scalar> PE1D;
193 #endif // TEUCHOS_DEBUG
194  using std::min;
195  using std::max;
196 
197  const RCP<Teuchos::FancyOStream> out = this->getOStream();
198 
199  *out << "\nStarting Armijo Quadratic interpolation linesearch ...\n";
200 
201 #ifdef TEUCHOS_DEBUG
202  TEUCHOS_ASSERT_EQUALITY(point_k.alpha, ST::zero());
203  TEUCHOS_ASSERT_INEQUALITY(point_k.phi, !=, PE1D::valNotGiven());
204  TEUCHOS_ASSERT_INEQUALITY(point_k.Dphi, !=, PE1D::valNotGiven());
205  TEUCHOS_ASSERT(!is_null(point_kp1));
206  TEUCHOS_ASSERT_INEQUALITY(point_kp1->alpha, >, ST::zero());
207  TEUCHOS_ASSERT_INEQUALITY(point_kp1->phi, !=, PE1D::valNotGiven());
208  TEUCHOS_ASSERT_EQUALITY(point_kp1->Dphi, PE1D::valNotGiven());
209 #endif // TEUCHOS_DEBUG
210 
211  const Scalar phi_k = point_k.phi;
212  Scalar &alpha_k = point_kp1->alpha;
213  Scalar &phi_kp1 = point_kp1->phi;
214 
215  // Loop initialization (technically the first iteration)
216 
217  const Scalar Dphi_k = point_k.Dphi;
218 
219  // output header
220 
221  *out
222  << "\nDphi_k = " << Dphi_k
223  << "\nphi_k = " << phi_k
224  << "\n";
225  if (minIters_ > 0)
226  *out << "\nminIters == " << minIters_ << "\n";
227  if (doMaxIters_)
228  *out << "\ndoMaxIters == true, maxing out maxIters = "
229  <<maxIters_<<" iterations!\n";
230  *out << "\n";
231 
232  TabularOutputter tblout(out);
233 
234  tblout.pushFieldSpec("itr", TO::INT);
235  tblout.pushFieldSpec("alpha_k", TO::DOUBLE);
236  tblout.pushFieldSpec("phi_kp1", TO::DOUBLE);
237  tblout.pushFieldSpec("phi_kp1-frac_phi", TO::DOUBLE);
238  tblout.pushFieldSpec("alpha_interp", TO::DOUBLE);
239  tblout.pushFieldSpec("alpha_min", TO::DOUBLE);
240  tblout.pushFieldSpec("alpha_max", TO::DOUBLE);
241 
242  tblout.outputHeader();
243 
244  // Check that this is a decent direction
246  "ArmijoPolyInterpLineSearch::doLineSearch(): "
247  "The given descent direction for the given "
248  "phi Dphi_k="<<Dphi_k<<" >= 0!" );
249 
250  // keep memory of the best value
251  Scalar best_alpha = alpha_k;
252  Scalar best_phi = phi_kp1;
253 
254  // Perform linesearch.
255  bool success = false;
256  int numIters = 0;
257  for ( ; numIters < maxIters_; ++numIters ) {
258 
259  // Print out this iteration.
260 
261  Scalar frac_phi = phi_k + eta_ * alpha_k * Dphi_k;
262  tblout.outputField(numIters);
263  tblout.outputField(alpha_k);
264  tblout.outputField(phi_kp1);
265  tblout.outputField(((phi_kp1)-frac_phi));
266 
267  if (ST::isnaninf(phi_kp1)) {
268 
269  // Cut back the step to minFrac * alpha_k
270  alpha_k = minFrac_ * alpha_k;
271  best_alpha = ST::zero();
272  best_phi = phi_k;
273  }
274  else {
275 
276  // Armijo condition
277  if (phi_kp1 < frac_phi) {
278  // We have found an acceptable point
279  success = true;
280  if (numIters < minIters_) {
281  // Keep going!
282  }
283  else if ( !doMaxIters_ || ( doMaxIters_ && numIters == maxIters_ - 1 ) ) {
284  tblout.nextRow(true);
285  break; // get out of the loop, we are done!
286  }
287  }
288  else {
289  success = false;
290  }
291 
292  // Select a new alpha to try:
293  // alpha_k = ( minFracalpha_k <= quadratic interpolation <= maxFracalpha_k )
294 
295  // Quadratic interpolation of alpha_k that minimizes phi.
296  // We know the values of phi at the initail point and alpha_k and
297  // the derivate of phi w.r.t alpha at the initial point and
298  // that's enough information for a quandratic interpolation.
299 
300  Scalar alpha_quad =
301  ( -as<Scalar>(0.5) * Dphi_k * alpha_k * alpha_k )
302  /
303  ( (phi_kp1) - phi_k - alpha_k * Dphi_k );
304 
305  tblout.outputField(alpha_quad);
306 
307  // 2009/01/29: rabartl: ToDo: Call the interpolation and then add
308  // support for other types of interpolations based on various points.
309 
310  const Scalar alpha_min = minFrac_ * alpha_k;
311  const Scalar alpha_max = maxFrac_ * alpha_k;
312 
313  tblout.outputField(alpha_min);
314  tblout.outputField(alpha_max);
315 
316  alpha_k =
317  min(
318  max(alpha_min, alpha_quad),
319  alpha_max
320  );
321 
322  }
323 
324  tblout.nextRow(true);
325 
326 
327  // Evaluate the point
328 
329  phi_kp1 = computeValue<Scalar>(phi, alpha_k);
330 
331  // Save the best point found
332  if (phi_kp1 < best_phi) {
333  best_phi = phi_kp1;
334  best_alpha = alpha_k;
335  }
336 
337  }
338 
339  if (!is_null(numIters_out))
340  *numIters_out = numIters;
341 
342  if( success ) {
343  *out << "\nLine search success!\n";
344  return true;
345  }
346 
347  // Line search failure. Return the best point found and let the client
348  // decide what to do.
349  alpha_k = best_alpha;
350  phi_kp1 = computeValue<Scalar>(phi, best_alpha);
351  *out << "\nLine search failure!\n";
352  return false;
353 
354 }
355 
356 
357 } // namespace GlobiPack
358 
359 
360 #endif // GLOBIPACK_POLY_INTERP_LINE_SEARCH_DEF_HPP
Scalar phi
The value of the merit function phi(alpha).
bool is_null(const boost::shared_ptr< T > &p)
Scalar Dphi
The value of the derivative of the merit function Dphi(alpha).
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)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Thrown if search direction not a descent direction for the merit function.
virtual bool doLineSearch(const MeritFunc1DBase< Scalar > &phi, const PointEval1D< Scalar > &point_k, const Ptr< PointEval1D< Scalar > > &point_kp1, const Ptr< int > &numIters) 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 validateParametersAndSetDefaults(ParameterList const &validParamList, int const depth=1000)
void setParameterList(RCP< ParameterList > const &paramList)
TypeTo as(const TypeFrom &t)
ArmijoPolyInterpLineSearch()
Construct with default parameters.
Base class for 1D merit fucntions used in globalization methods.
#define TEUCHOS_ASSERT(assertion_test)
#define TEUCHOS_ASSERT_EQUALITY(val1, val2)