MoochoPack : Framework for Large-Scale Optimization Algorithms  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
MoochoPack_LineSearchNLE_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_LineSearchNLE_Step.hpp"
46 #include "MoochoPack_Exceptions.hpp"
47 #include "MoochoPack_moocho_algo_conversion.hpp"
48 #include "IterationPack_print_algorithm_step.hpp"
49 #include "ConstrainedOptPack_MeritFuncNLESqrResid.hpp"
50 #include "ConstrainedOptPack_MeritFuncCalcNLE.hpp"
51 #include "ConstrainedOptPack_MeritFuncCalc1DQuadratic.hpp"
52 #include "AbstractLinAlgPack_VectorMutable.hpp"
53 #include "AbstractLinAlgPack_VectorStdOps.hpp"
54 #include "AbstractLinAlgPack_VectorOut.hpp"
55 #include "AbstractLinAlgPack_assert_print_nan_inf.hpp"
56 #include "AbstractLinAlgPack_LinAlgOpPack.hpp"
57 #include "Teuchos_Assert.hpp"
58 
59 namespace MoochoPack {
60 
62  const direct_line_search_ptr_t& direct_line_search
63  )
64  :direct_line_search_(direct_line_search)
65 {}
66 
68  Algorithm& _algo, poss_type step_poss, IterationPack::EDoStepType type
69  ,poss_type assoc_step_poss
70  )
71 {
73  using LinAlgOpPack::V_VpV;
74 
75  NLPAlgo &algo = rsqp_algo(_algo);
76  NLPAlgoState &s = algo.rsqp_state();
77  NLP &nlp = algo.nlp();
78 
79  EJournalOutputLevel olevel = algo.algo_cntr().journal_output_level();
80  std::ostream& out = algo.track().journal_out();
81 
82  // print step header.
83  if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
84  using IterationPack::print_algorithm_step;
85  print_algorithm_step( algo, step_poss, type, assoc_step_poss, out );
86  }
87 
88  const size_type
89  //n = nlp.n(),
90  m = nlp.m();
91 
92  TEUCHOS_TEST_FOR_EXCEPTION( m == 0, std::logic_error, "LineSearchNLE_Step::do_step(...) : Error!" );
93 
94  // /////////////////////////////////////////
95  // Set references to iteration quantities
96  //
97  // Set k+1 first then go back to get k+0 to ensure
98  // we have backward storage!
99 
100  IterQuantityAccess<value_type>
101  &alpha_iq = s.alpha();
102  IterQuantityAccess<VectorMutable>
103  &x_iq = s.x(),
104  &d_iq = s.d(),
105  &c_iq = s.c();
106 
107  VectorMutable &x_kp1 = x_iq.get_k(+1);
108  const Vector &x_k = x_iq.get_k(0);
109  VectorMutable &c_kp1 = c_iq.get_k(+1);
110  const Vector &c_k = c_iq.get_k(0);
111  const Vector &d_k = d_iq.get_k(0);
112  value_type &alpha_k = alpha_iq.get_k(0);
113 
114  // //////////////////////////////////////
115  // Build the merit function
116 
118 
119  // /////////////////////////////////////
120  // Compute Dphi_k, phi_kp1 and phi_k
121 
122  // phi_k, phi_kp1
123  const value_type phi_k = phi_c.value(c_k);
124  value_type phi_kp1 = phi_c.value(c_kp1);
125  if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
126  out << "\nBegin definition of NLP merit function phi_c.value(c(x)):\n";
127  phi_c.print_merit_func( out, " " );
128  out << "end definition of the NLP merit funciton\n";
129  }
130  // Dphi_k
131  phi_c.calc_deriv(c_k);
132  const value_type
133  Dphi_k = phi_c.deriv();
134  if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
135  out << "\nDphi_k = " << Dphi_k << std::endl;
136  }
138  Dphi_k >= 0, LineSearchFailure
139  ,"LineSearchNLE_Step::do_step(...) : "
140  "Error, d_k is not a descent direction for the merit function "
141  "since Dphi_k = " << Dphi_k << " >= 0" );
142 
143  // //////////////////////
144  // Do the line search
145 
146  nlp.unset_quantities();
147  nlp.set_c( &c_kp1 );
148  ConstrainedOptPack::MeritFuncCalcNLE phi_c_calc( &phi_c, &nlp );
149  const Vector* xd[2] = { &x_k, &d_k };
150  MeritFuncCalc1DQuadratic phi_calc_1d( phi_c_calc, 1, xd, &x_kp1 );
151 
152  if( !direct_line_search().do_line_search(
153  phi_calc_1d, phi_k, &alpha_k, &phi_kp1
154  ,( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS)
155  ? &out : static_cast<std::ostream*>(0) )
156  )
157  )
158  {
159  // The line search failed!
160  if( static_cast<int>(olevel) >= static_cast<int>(PRINT_BASIC_ALGORITHM_INFO) )
161  out
162  << "\nThe maximum number of linesearch iterations has been exceeded "
163  << "(k = " << algo.state().k() << ")\n"
164  << "(phi_k - phi_kp1)/phi_k = " << ((phi_k - phi_kp1)/phi_k)
165  << "\nso we will reject the step and declare a line search failure.\n";
167  true, LineSearchFailure
168  ,"LineSearchNLE_Step::do_step(): Line search failure" );
169  }
170 
171  nlp.unset_quantities();
172 
173  if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
174  out << "\nalpha_k = " << alpha_k;
175  out << "\n||x_kp1||inf = " << x_kp1.norm_inf();
176  out << "\n||c_kp1||inf = " << c_kp1.norm_inf();
177  out << "\nphi_kp1 = " << phi_kp1;
178  out << std::endl;
179  }
180 
181  if( (int)olevel >= (int)PRINT_VECTORS ) {
182  out << "\nx_kp1 =\n" << x_kp1;
183  out << "\nc_kp1 =\n" << c_kp1;
184  }
185 
186  return true;
187 }
188 
190  const Algorithm& algo, poss_type step_poss, IterationPack::EDoStepType type, poss_type assoc_step_poss
191  ,std::ostream& out, const std::string& L
192  ) const
193 {
194  out
195  << L << "*** Preform a line search for c(x_k + alpha_k*d_k) along the full space search direction d_k.\n"
196  << L << "ToDo: Fill this in!\n";
197 }
198 
199 } // end namespace MoochoPack
value_type calc_deriv(const Vector &c_k)
void print_merit_func(std::ostream &out, const std::string &leading_str) const
bool do_step(Algorithm &algo, poss_type step_poss, IterationPack::EDoStepType type, poss_type assoc_step_poss)
value_type value(const Vector &c) const
LineSearchNLE_Step(const direct_line_search_ptr_t &direct_line_search=Teuchos::null)
#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
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()
AlgorithmState & state()
NLPAlgoState & rsqp_state()
<<std aggr>="">> members for algo_cntr