MoochoPack : Framework for Large-Scale Optimization Algorithms  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
MoochoPack_UpdateBarrierParameter_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 "MoochoPack_UpdateBarrierParameter_Step.hpp"
48 #include "MoochoPack_moocho_algo_conversion.hpp"
49 #include "IterationPack_print_algorithm_step.hpp"
50 #include "Teuchos_dyn_cast.hpp"
51 #include "MoochoPack_IpState.hpp"
52 #include "AbstractLinAlgPack_MatrixSymDiagStd.hpp"
53 #include "AbstractLinAlgPack_VectorStdOps.hpp"
54 #include "AbstractLinAlgPack_VectorAuxiliaryOps.hpp"
55 #include "Teuchos_Assert.hpp"
56 
57 #include "OptionsFromStreamPack_StringToBool.hpp"
58 
59 #define min(a,b) ( (a < b) ? a : b )
60 #define max(a,b) ( (a > b) ? a : b )
61 
62 namespace MoochoPack {
63 
65  const value_type init_barrier_parameter,
66  const value_type tau_mu,
67  const value_type theta_mu,
68  const value_type tau_epsilon,
69  const value_type theta_epsilon,
70  const value_type e_tol_max
71  )
72  :
73  init_barrier_parameter_(init_barrier_parameter),
74  tau_mu_(tau_mu),
75  theta_mu_(theta_mu),
76  tau_epsilon_(tau_epsilon),
77  theta_epsilon_(theta_epsilon),
78  e_tol_max_(e_tol_max)
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 
90  NLPAlgo &algo = dyn_cast<NLPAlgo>(_algo);
91  IpState &s = dyn_cast<IpState>(_algo.state());
92  NLP &nlp = algo.nlp();
93 
94  if (!nlp.is_initialized())
95  {
96  nlp.initialize(false);
97  }
98 
99  EJournalOutputLevel olevel = algo.algo_cntr().journal_output_level();
100  std::ostream& out = algo.track().journal_out();
101 
102  // print step header.
103  if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) )
104  {
105  using IterationPack::print_algorithm_step;
106  print_algorithm_step( _algo, step_poss, type, assoc_step_poss, out );
107  }
108 
109 
111  // Get iteration quantities
113  IterQuantityAccess<value_type>& e_tol_iq = s.e_tol();
114  IterQuantityAccess<value_type>& mu_iq = s.barrier_parameter();
115 
117  // Check values and initialize, if necessary
119  /* if (mu_iq.last_updated() == IterQuantity::NONE_UPDATED)
120  {
121  mu_iq.set_k(0) = init_barrier_parameter_;
122  e_tol_iq.set_k(0) = Calculate_e_tol(mu_iq.get_k(0));
123 
124  if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) )
125  {
126  out << "\nInitializing barrier parameter (mu) and sub problem tolerance (e_tol) ...\n";
127  }
128  }
129  else*/
130  {
132  // if tolerance is satisfied, calculate new barrier_parameter
133  // and e_tol, otherwise update mu and e_tol from last iter
135  const value_type opt_err = s.opt_kkt_err().get_k(0);
136  const value_type feas_err = s.feas_kkt_err().get_k(0);
137  const value_type comp_err_mu = s.comp_err_mu().get_k(0);
138 
139  const value_type mu_km1 = mu_iq.get_k(-1);
140  if (e_tol_iq.last_updated() == IterQuantity::NONE_UPDATED)
141  {
142  // First time through, don't let mu change
143  mu_iq.set_k(0,-1);
144  e_tol_iq.set_k(0) = Calculate_e_tol(mu_iq.get_k(0));
145  }
146  else
147  {
148  const value_type e_tol_km1 = e_tol_iq.get_k(-1);
149  bool sub_prob_converged = (opt_err < e_tol_km1 && feas_err < e_tol_km1 && comp_err_mu < e_tol_km1);
150  if (sub_prob_converged)
151  {
152  // Calculate new mu and e_tol
153  value_type& mu_k = mu_iq.set_k(0);
154  mu_k = min(tau_mu_*mu_km1, pow(mu_km1, theta_mu_));
155  e_tol_iq.set_k(0) = Calculate_e_tol(mu_k);
156 
157  if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) )
158  {
159  out << "\nSub-problem converged!\n"
160  << " Updating barrier parameter (mu) and sub problem tolerance (e_tol) ...\n";
161  }
162  }
163  else
164  {
165  if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) )
166  {
167  out << "\nSub-problem not-converged!\n"
168  << " Keeping existing barrier parameter (mu) and sub problem tolerance (e_tol) ...\n";
169  }
170  mu_iq.set_k(0,-1);
171  e_tol_iq.set_k(0,-1);
172  }
173  }
174 
175  if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) )
176  {
177  out << "\nbarrier_parameter (mu) = " << mu_iq.get_k(0)
178  << "\nsub problem tolerance (e_tol) = " << e_tol_iq.get_k(0) << std::endl;
179  }
180  }
181 
182  return true;
183  }
184 
185 
186 void UpdateBarrierParameter_Step::print_step(
187  const Algorithm& _algo, poss_type step_poss, IterationPack::EDoStepType type
188  ,poss_type assoc_step_poss, std::ostream& out, const std::string& L
189  ) const
190  {
191  //const NLPAlgo &algo = rsqp_algo(_algo);
192  //const NLPAlgoState &s = algo.rsqp_state();
193  out << L << "# Update the interior point barrier parameter (mu)\n"
194  << L << "if (KKTerror < e_tol) then\n"
195  << L << " mu_kp1 = min(tau_mu*mu_k, mu_k^theta_mu)\n"
196  << L << " e_tol_kp1 = min(e_tol_max, tau_epsilon*min(mu_k, mu_k^theta_epsilon))\n"
197  << L << "else\n"
198  << L << " mu_kp1 = mu_k\n"
199  << L << " e_tol_kp1 = e_tol_k\n"
200  << L << "end;\n";
201  }
202 
203 value_type UpdateBarrierParameter_Step::Calculate_e_tol(value_type mu)
204  {
205  value_type e_tol = tau_epsilon_*min(mu, pow(mu, theta_epsilon_));
206  e_tol = min(e_tol_max_, e_tol);
207 
208  return e_tol;
209  }
210 
211 
212 namespace {
213 
214 const int local_num_options = 5;
215 
216 enum local_EOptions
217  {
218  TAU_MU,
219  THETA_MU,
220  TAU_EPSILON,
221  THETA_EPSILON,
222  E_TOL_MAX
223  };
224 
225 const char* local_SOptions[local_num_options] =
226  {
227  "tau_mu",
228  "theta_mu",
229  "tau_epsilon",
230  "theta_epsilon",
231  "e_tol_max"
232  };
233 
234 }
235 
236 
237 UpdateBarrierParameter_StepSetOptions::UpdateBarrierParameter_StepSetOptions(
238  UpdateBarrierParameter_Step* target
239  , const char opt_grp_name[] )
240  :
241  OptionsFromStreamPack::SetOptionsFromStreamNode(
242  opt_grp_name, local_num_options, local_SOptions ),
243  OptionsFromStreamPack::SetOptionsToTargetBase< UpdateBarrierParameter_Step >( target )
244  {
245  }
246 
247 void UpdateBarrierParameter_StepSetOptions::setOption(
248  int option_num, const std::string& option_value )
249  {
250  using OptionsFromStreamPack::StringToBool;
251 
252  typedef UpdateBarrierParameter_Step target_t;
253  switch( (local_EOptions)option_num )
254  {
255  case TAU_MU:
256  target().tau_mu(std::atof(option_value.c_str()));
257  break;
258  case THETA_MU:
259  target().theta_mu(std::atof(option_value.c_str()));
260  break;
261  case TAU_EPSILON:
262  target().tau_epsilon(std::atof(option_value.c_str()));
263  break;
264  case THETA_EPSILON:
265  target().theta_epsilon(std::atof(option_value.c_str()));
266  break;
267  case E_TOL_MAX:
268  target().e_tol_max(std::atof(option_value.c_str()));
269  break;
270  default:
271  TEUCHOS_TEST_FOR_EXCEPT(true); // Local error only?
272  }
273  }
274 
275 } // end namespace MoochoPack
UpdateBarrierParameter_Step(const value_type init_barrier_parameter=0.1, const value_type tau_mu=0.2, const value_type theta_mu=1.5, const value_type tau_epsilon=10, const value_type theta_epsilon=1.1, const value_type e_tol_max=1000)
Constructor.
T_To & dyn_cast(T_From &from)
rSQP Algorithm control class.
bool do_step(Algorithm &algo, poss_type step_poss, IterationPack::EDoStepType type, poss_type assoc_step_poss)
virtual std::ostream & journal_out() const
AlgorithmTracker & track()
AlgorithmState & state()
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)