MoochoPack : Framework for Large-Scale Optimization Algorithms  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
MoochoPack_PreProcessBarrierLineSearch_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 #include <iostream>
45 #include <math.h>
46 
47 #include "AbstractLinAlgPack_assert_print_nan_inf.hpp"
48 #include "AbstractLinAlgPack_VectorAuxiliaryOps.hpp"
49 #include "AbstractLinAlgPack_MatrixSymDiagStd.hpp"
50 #include "AbstractLinAlgPack_VectorOut.hpp"
51 #include "AbstractLinAlgPack_LinAlgOpPack.hpp"
52 #include "NLPInterfacePack_NLPBarrier.hpp"
53 #include "MoochoPack_PreProcessBarrierLineSearch_Step.hpp"
54 #include "MoochoPack_IpState.hpp"
55 #include "MoochoPack_moocho_algo_conversion.hpp"
56 #include "IterationPack_print_algorithm_step.hpp"
57 #include "Teuchos_dyn_cast.hpp"
58 #include "Teuchos_Assert.hpp"
59 
60 #define min(a,b) ( (a < b) ? a : b )
61 #define max(a,b) ( (a > b) ? a : b )
62 
63 namespace MoochoPack {
64 
67  const value_type tau_boundary_frac
68  )
69  :
70  tau_boundary_frac_(tau_boundary_frac),
71  barrier_nlp_(barrier_nlp),
72  filter_(FILTER_IQ_STRING)
73 {
75  !barrier_nlp_.get(),
76  std::logic_error,
77  "PreProcessBarrierLineSearch_Step given NULL NLPBarrier."
78  );
79 }
80 
81 
83  Algorithm& _algo, poss_type step_poss, IterationPack::EDoStepType type
84  ,poss_type assoc_step_poss
85  )
86 {
87  using Teuchos::dyn_cast;
88  using IterationPack::print_algorithm_step;
89  using AbstractLinAlgPack::assert_print_nan_inf;
93 
94  NLPAlgo &algo = dyn_cast<NLPAlgo>(_algo);
95  IpState &s = dyn_cast<IpState>(_algo.state());
96  NLP &nlp = algo.nlp();
97 
98  EJournalOutputLevel olevel = algo.algo_cntr().journal_output_level();
99  std::ostream& out = algo.track().journal_out();
100 
101  // print step header.
102  if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) )
103  {
104  using IterationPack::print_algorithm_step;
105  print_algorithm_step( _algo, step_poss, type, assoc_step_poss, out );
106  }
107 
108  const value_type& mu_k = s.barrier_parameter().get_k(0);
109 
110  // if using filter and u changed, clear filter
111  if (filter_.exists_in(s))
112  {
113  if ( s.barrier_parameter().updated_k(-1) )
114  {
115  const value_type mu_km1 = s.barrier_parameter().get_k(-1);
116  if (mu_k != mu_km1)
117  {
118  if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) )
119  {
120  out << "\nBarrier Parameter changed - resetting the filter ...\n";
121  }
122  // reset the filter
123  MoochoPack::Filter_T &filter_k = filter_(s).set_k(0);
124  filter_k.clear();
125  }
126  }
127  }
128 
129  // Update the barrier parameter in the NLP
130  barrier_nlp_->mu(s.barrier_parameter().get_k(0));
131 
132  // Calculate the barrier k terms
133  barrier_nlp_->set_Gf( &(s.grad_barrier_obj().set_k(0)) );
134  barrier_nlp_->calc_Gf(s.x().get_k(0), true);
135 
136  barrier_nlp_->set_f( &(s.barrier_obj().set_k(0)) );
137  barrier_nlp_->calc_f(s.x().get_k(0), true);
138 
139 
140  // Calculate the k+1 terms
141  // Get iteration quantities...
142  value_type& alpha_k = s.alpha().set_k(0);
143  value_type& alpha_vl_k = s.alpha_vl().set_k(0);
144  value_type& alpha_vu_k = s.alpha_vu().set_k(0);
145 
146  const Vector& x_k = s.x().get_k(0);
147  VectorMutable& x_kp1 = s.x().set_k(+1);
148 
149  const Vector& d_k = s.d().get_k(0);
150  const Vector& dvl_k = s.dvl().get_k(0);
151  const Vector& dvu_k = s.dvu().get_k(0);
152 
153  const Vector& vl_k = s.Vl().get_k(0).diag();
154  VectorMutable& vl_kp1 = s.Vl().set_k(+1).diag();
155 
156  const Vector& vu_k = s.Vu().get_k(0).diag();
157  VectorMutable& vu_kp1 = s.Vu().set_k(+1).diag();
158 
159  alpha_k = fraction_to_boundary(
160  tau_boundary_frac_,
161  x_k,
162  d_k,
163  nlp.xl(),
164  nlp.xu()
165  );
166 
167  alpha_vl_k = fraction_to_zero_boundary(
168  tau_boundary_frac_,
169  vl_k,
170  dvl_k
171  );
172 
173  alpha_vu_k = fraction_to_zero_boundary(
174  tau_boundary_frac_,
175  vu_k,
176  dvu_k
177  );
178 
179  TEUCHOS_TEST_FOR_EXCEPT( !( alpha_k <= 1.0 && alpha_vl_k <= 1.0 && alpha_vu_k <= 1.0 ) );
180  TEUCHOS_TEST_FOR_EXCEPT( !( alpha_k >= 0.0 && alpha_vl_k >= 0.0 && alpha_vu_k >= 0.0 ) );
181 
182  x_kp1 = x_k;
183  Vp_StV( &x_kp1, alpha_k, d_k);
184 
185  alpha_vl_k = alpha_vu_k = min(alpha_vl_k, alpha_vu_k);
186 
187  vl_kp1 = vl_k;
188  Vp_StV( &vl_kp1, alpha_vl_k, dvl_k);
189 
190  vu_kp1 = vu_k;
191  Vp_StV( &vu_kp1, alpha_vu_k, dvu_k);
192 
193 
194  IterQuantityAccess<VectorMutable>
195  *c_iq = nlp.m() > 0 ? &s.c() : NULL;
196 
197  if (assert_print_nan_inf(x_kp1, "x", true, NULL))
198  {
199  // Calcuate f and c at the new point.
200  barrier_nlp_->unset_quantities();
201  barrier_nlp_->set_f( &s.barrier_obj().set_k(+1) );
202  if (c_iq) {
203  barrier_nlp_->set_c( &c_iq->set_k(+1) );
204  barrier_nlp_->calc_c( x_kp1, true );
205  }
206  barrier_nlp_->calc_f( x_kp1, false );
207  barrier_nlp_->unset_quantities();
208  }
209 
210  if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) )
211  {
212  out << "\nalpha_vl_k = " << alpha_vl_k
213  << "\nalpha_vu_k = " << alpha_vu_k
214  << "\nalpha_k = " << alpha_k
215  << std::endl;
216  }
217 
218  if( static_cast<int>(olevel) >= static_cast<int>(PRINT_VECTORS) )
219  {
220  out << "\nvl_kp1 = \n" << vl_kp1
221  << "\nvu_kp1 = \n" << vu_kp1
222  << "\nx_kp1 = \n" << x_kp1;
223  }
224 
225  return true;
226 }
227 
228 void PreProcessBarrierLineSearch_Step::print_step(
229  const Algorithm& _algo, poss_type step_poss, IterationPack::EDoStepType type
230  ,poss_type assoc_step_poss, std::ostream& out, const std::string& L
231  ) const
232 {
233  //const NLPAlgo &algo = rsqp_algo(_algo);
234  //const NLPAlgoState &s = algo.rsqp_state();
235  out << L << "*** calculate alpha max by the fraction to boundary rule\n"
236  << L << "ToDo: Complete documentation\n";
237 }
238 
239 namespace {
240 
241 const int local_num_options = 1;
242 
243 enum local_EOptions
244 {
245  TAU_BOUNDARY_FRAC
246 };
247 
248 const char* local_SOptions[local_num_options] =
249 {
250  "tau_boundary_frac"
251 };
252 
253 }
254 
255 
256 PreProcessBarrierLineSearch_StepSetOptions::PreProcessBarrierLineSearch_StepSetOptions(
257  PreProcessBarrierLineSearch_Step* target
258  , const char opt_grp_name[] )
259  :
260  OptionsFromStreamPack::SetOptionsFromStreamNode(
261  opt_grp_name, local_num_options, local_SOptions ),
262  OptionsFromStreamPack::SetOptionsToTargetBase< PreProcessBarrierLineSearch_Step >( target )
263 {}
264 
265 void PreProcessBarrierLineSearch_StepSetOptions::setOption(
266  int option_num, const std::string& option_value )
267 {
268  typedef PreProcessBarrierLineSearch_Step target_t;
269  switch( (local_EOptions)option_num )
270  {
271  case TAU_BOUNDARY_FRAC:
272  target().tau_boundary_frac(std::atof(option_value.c_str()));
273  break;
274  default:
275  TEUCHOS_TEST_FOR_EXCEPT(true); // Local error only?
276  }
277 }
278 
279 } // end namespace MoochoPack
void calc_f(const Vector &x, bool newx=true) const
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
void Vp_StV(VectorMutable *v_lhs, const value_type &alpha, const Vector &v_rhs)
T_To & dyn_cast(T_From &from)
rSQP Algorithm control class.
T * get() const
virtual void unset_quantities()
void mu(const value_type mu)
bool do_step(Algorithm &algo, poss_type step_poss, IterationPack::EDoStepType type, poss_type assoc_step_poss)
virtual std::ostream & journal_out() const
void set_c(VectorMutable *c)
AlgorithmTracker & track()
value_type fraction_to_boundary(const value_type tau, const Vector &x, const Vector &d, const Vector &xl, const Vector &xu)
value_type fraction_to_zero_boundary(const value_type tau, const Vector &x, const Vector &d)
void calc_c(const Vector &x, bool newx=true) const
AlgorithmState & state()
void set_Gf(VectorMutable *Gf)
void set_f(value_type *f)
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)
PreProcessBarrierLineSearch_Step(Teuchos::RCP< NLPInterfacePack::NLPBarrier > barrier_nlp, const value_type tau_boundary_frac=0.99)
Constructor.
void calc_Gf(const Vector &x, bool newx=true) const