MoochoPack : Framework for Large-Scale Optimization Algorithms  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
MoochoPack_LineSearchDirect_Step.cpp
1 // @HEADER
2 // ***********************************************************************
3 //
4 // Moocho: Multi-functional Object-Oriented arCHitecture for Optimization
5 // Copyright (2003) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov)
38 //
39 // ***********************************************************************
40 // @HEADER
41 
42 #include <ostream>
43 #include <typeinfo>
44 
45 #include "MoochoPack_LineSearchDirect_Step.hpp"
46 #include "MoochoPack_Exceptions.hpp"
47 #include "MoochoPack_moocho_algo_conversion.hpp"
48 #include "IterationPack_print_algorithm_step.hpp"
49 #include "ConstrainedOptPack_MeritFuncCalc1DQuadratic.hpp"
50 #include "ConstrainedOptPack_MeritFuncCalcNLP.hpp"
51 #include "AbstractLinAlgPack_VectorMutable.hpp"
52 #include "AbstractLinAlgPack_VectorStdOps.hpp"
53 #include "AbstractLinAlgPack_VectorOut.hpp"
54 #include "AbstractLinAlgPack_assert_print_nan_inf.hpp"
55 #include "AbstractLinAlgPack_LinAlgOpPack.hpp"
56 #include "Teuchos_Assert.hpp"
57 
58 namespace MoochoPack {
59 
61  const direct_line_search_ptr_t& direct_line_search
62  )
63  :direct_line_search_(direct_line_search)
64 {}
65 
67  Algorithm& _algo, poss_type step_poss, IterationPack::EDoStepType type
68  ,poss_type assoc_step_poss
69  )
70 {
72  using LinAlgOpPack::V_VpV;
73 
74  NLPAlgo &algo = rsqp_algo(_algo);
75  NLPAlgoState &s = algo.rsqp_state();
76  NLP &nlp = algo.nlp();
77 
78  EJournalOutputLevel olevel = algo.algo_cntr().journal_output_level();
79  std::ostream& out = algo.track().journal_out();
80 
81  // print step header.
82  if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
83  using IterationPack::print_algorithm_step;
84  print_algorithm_step( algo, step_poss, type, assoc_step_poss, out );
85  }
86 
88  nlpOutputTempState(
89  Teuchos::rcp(&nlp,false), Teuchos::getFancyOStream(Teuchos::rcp(&out,false)),
91  // Above, we don't want to litter the output with any trace from the NLP.
92  // However, if the user forces the verbosity level to be increased, then we
93  // want to set the stream so that it knows where to print to.
94 
95  const size_type
96  m = nlp.m();
97 
98  // /////////////////////////////////////////
99  // Set references to iteration quantities
100  //
101  // Set k+1 first then go back to get k+0 to ensure
102  // we have backward storage!
103 
104  IterQuantityAccess<value_type>
105  &f_iq = s.f(),
106  &alpha_iq = s.alpha(),
107  &phi_iq = s.phi();
108  IterQuantityAccess<VectorMutable>
109  &x_iq = s.x(),
110  &d_iq = s.d(),
111  &c_iq = s.c();
112 
113  VectorMutable &x_kp1 = x_iq.get_k(+1);
114  const Vector &x_k = x_iq.get_k(0);
115  value_type &f_kp1 = f_iq.get_k(+1);
116  const value_type &f_k = f_iq.get_k(0);
117  VectorMutable *c_kp1 = m ? &c_iq.get_k(+1) : NULL;
118  const Vector *c_k = m ? &c_iq.get_k(0) : NULL;
119  const Vector &d_k = d_iq.get_k(0);
120  value_type &alpha_k = alpha_iq.get_k(0);
121 
122  // /////////////////////////////////////
123  // Compute Dphi_k, phi_kp1 and phi_k
124 
125  const MeritFuncNLP
126  &merit_func_nlp_k = s.merit_func_nlp().get_k(0);
127  if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
128  out << "\nBegin definition of NLP merit function phi.value(f(x),c(x)):\n";
129  merit_func_nlp_k.print_merit_func( out, " " );
130  out << "end definition of the NLP merit funciton\n";
131  }
132  // Dphi_k
133  const value_type
134  Dphi_k = merit_func_nlp_k.deriv();
135  if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
136  out << "\nDphi_k = " << Dphi_k << std::endl;
137  }
139  Dphi_k >= 0, LineSearchFailure
140  ,"LineSearch2ndOrderCorrect_Step::do_step(...) : "
141  "Error, d_k is not a descent direction for the merit function "
142  "since Dphi_k = " << Dphi_k << " >= 0" );
143  // ph_kp1
144  value_type
145  &phi_kp1 = phi_iq.set_k(+1) = merit_func_nlp_k.value(
146  f_kp1
147  ,c_kp1
148  ,NULL // h
149  ,NULL // hl
150  ,NULL // hu
151  );
152  // Must compute phi(x) at the base point x_k since the penalty parameter may have changed.
153  const value_type
154  &phi_k = phi_iq.set_k(0) = merit_func_nlp_k.value(
155  f_k
156  ,c_k
157  ,NULL // h
158  ,NULL // hl
159  ,NULL // hu
160  );
161 
162  // //////////////////////////////////////
163  // Setup the calculation merit function
164 
165  // Here f_kp1, c_kp1 and h_kp1 are updated at the same time the
166  // line search is being performed!
167  nlp.unset_quantities();
168  nlp.set_f( &f_kp1 );
169  if(m) nlp.set_c( c_kp1 );
170  MeritFuncCalcNLP
171  phi_calc( &merit_func_nlp_k, &nlp );
172 
173  // //////////////////////
174  // Do the line search
175 
176  const Vector* xd[2] = { &x_k, &d_k };
177  MeritFuncCalc1DQuadratic
178  phi_calc_1d( phi_calc, 1, xd, &x_kp1 );
179 
180  if( !direct_line_search().do_line_search(
181  phi_calc_1d, phi_k, &alpha_k, &phi_kp1
182  ,( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS)
183  ? &out : static_cast<std::ostream*>(0) )
184  )
185  )
186  {
187  // The line search failed!
188  if( static_cast<int>(olevel) >= static_cast<int>(PRINT_BASIC_ALGORITHM_INFO) )
189  out
190  << "\nThe maximum number of linesearch iterations has been exceeded "
191  << "(k = " << algo.state().k() << ")\n"
192  << "(phi_k - phi_kp1)/phi_k = " << ((phi_k - phi_kp1)/phi_k)
193  << "\nso we will reject the step and declare a line search failure.\n";
195  true, LineSearchFailure
196  ,"LineSearchDirect_Step::do_step(): Line search failure" );
197  }
198  nlp.unset_quantities();
199 
200  if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
201  out << "\nalpha_k = " << alpha_k;
202  out << "\n||x_kp1||inf = " << x_kp1.norm_inf();
203  out << "\nf_kp1 = " << f_kp1;
204  if(m)
205  out << "\n||c_kp1||inf = " << c_kp1->norm_inf();
206  out << "\nphi_kp1 = " << phi_kp1;
207  out << std::endl;
208  }
209 
210  if( (int)olevel >= (int)PRINT_VECTORS ) {
211  out << "\nx_kp1 =\n" << x_kp1;
212  if(m)
213  out <<"\nc_kp1 =\n" << *c_kp1;
214  }
215 
216  return true;
217 }
218 
220  const Algorithm& algo, poss_type step_poss, IterationPack::EDoStepType type, poss_type assoc_step_poss
221  ,std::ostream& out, const std::string& L
222  ) const
223 {
224  out
225  << L << "*** Preform a line search along the full space search direction d_k.\n"
226  << L << "Dphi_k = merit_func_nlp_k.deriv()\n"
227  << L << "if Dphi_k >= 0 then\n"
228  << L << " throw line_search_failure\n"
229  << L << "end\n"
230  << L << "phi_kp1 = merit_func_nlp_k.value(f_kp1,c_kp1,h_kp1,hl,hu)\n"
231  << L << "phi_k = merit_func_nlp_k.value(f_k,c_k,h_k,hl,hu)\n"
232  << L << "begin direct line search (where phi = merit_func_nlp_k): \"" << typeName(direct_line_search()) << "\"\n";
233  direct_line_search().print_algorithm( out, L + " " );
234  out
235  << L << "end direct line search\n"
236  << L << "if maximum number of linesearch iterations are exceeded then\n"
237  << L << " throw line_search_failure\n"
238  << L << "end\n";
239 }
240 
241 } // end namespace MoochoPack
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
void Vp_StV(VectorMutable *v_lhs, const value_type &alpha, const Vector &v_rhs)
rSQP Algorithm control class.
void print_step(const Algorithm &algo, poss_type step_poss, IterationPack::EDoStepType type, poss_type assoc_step_poss, std::ostream &out, const std::string &leading_str) const
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
virtual std::ostream & journal_out() const
Thrown if a line search failure occurs.
void V_VpV(VectorMutable *v_lhs, const V1 &V1_rhs1, const V2 &V2_rhs2)
Reduced space SQP state encapsulation interface.
size_t size_type
AlgorithmTracker & track()
bool do_step(Algorithm &algo, poss_type step_poss, IterationPack::EDoStepType type, poss_type assoc_step_poss)
AlgorithmState & state()
NLPAlgoState & rsqp_state()
<<std aggr>="">> members for algo_cntr
LineSearchDirect_Step(const direct_line_search_ptr_t &direct_line_search=Teuchos::null)
std::string typeName(const T &t)