47 #include "MoochoPack_MeritFunc_PenaltyParamsUpdateWithMult_AddedStep.hpp"
48 #include "MoochoPack_moocho_algo_conversion.hpp"
49 #include "IterationPack_print_algorithm_step.hpp"
50 #include "ConstrainedOptPack_MeritFuncPenaltyParams.hpp"
51 #include "ConstrainedOptPack_MeritFuncNLPDirecDeriv.hpp"
52 #include "ConstrainedOptPack/src/VectorWithNorms.h"
53 #include "DenseLinAlgPack_DVectorOp.hpp"
54 #include "DenseLinAlgPack_DVectorClass.hpp"
55 #include "DenseLinAlgPack_DVectorOut.hpp"
59 typedef MoochoPack::value_type value_type;
60 inline value_type max(value_type v1, value_type v2)
61 {
return (v1 > v2) ? v1 : v2; }
65 namespace MoochoPack {
68 const merit_func_ptr_t& merit_func, value_type small_mu, value_type min_mu_ratio
69 , value_type mult_factor, value_type kkt_near_sol )
70 : merit_func_(merit_func), near_solution_(false)
71 , small_mu_(small_mu), min_mu_ratio_(min_mu_ratio), mult_factor_(mult_factor)
72 , kkt_near_sol_(kkt_near_sol), norm_inf_mu_last_(0.0)
76 , poss_type step_poss, IterationPack::EDoStepType type
77 , poss_type assoc_step_poss)
79 using DenseLinAlgPack::norm_inf;
81 NLPAlgo &algo = rsqp_algo(_algo);
82 NLPAlgoState &s = algo.rsqp_state();
84 EJournalOutputLevel olevel = algo.algo_cntr().journal_output_level();
85 std::ostream& out = algo.track().journal_out();
88 if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
89 using IterationPack::print_algorithm_step;
90 print_algorithm_step( algo, step_poss, type, assoc_step_poss, out );
93 MeritFuncPenaltyParams
94 *params =
dynamic_cast<MeritFuncPenaltyParams*
>(&merit_func());
96 std::ostringstream omsg;
98 <<
"MeritFunc_PenaltyParamsUpdateWithMult_AddedStep::do_step(...), Error "
99 <<
"The class " <<
typeName(&merit_func()) <<
" does not support the "
100 <<
"MeritFuncPenaltyParams iterface\n";
102 throw std::logic_error( omsg.str() );
105 MeritFuncNLPDirecDeriv
106 *direc_deriv =
dynamic_cast<MeritFuncNLPDirecDeriv*
>(&merit_func());
108 std::ostringstream omsg;
110 <<
"MeritFunc_PenaltyParamsUpdateWithMult_AddedStep::do_step(...), Error "
111 <<
"The class " <<
typeName(&merit_func()) <<
" does not support the "
112 <<
"MeritFuncNLPDirecDeriv iterface\n";
114 throw std::logic_error( omsg.str() );
117 bool perform_update =
true;
119 if( s.mu().updated_k(0) ) {
120 if( (
int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
121 out <<
"\nmu_k is already updated by someone else?\n";
123 const value_type mu_k = s.mu().get_k(0);
124 if( mu_k == norm_inf_mu_last_ ) {
125 if( (
int)olevel >= (
int)PRINT_ALGORITHM_STEPS ) {
126 out <<
"\nmu_k " << mu_k <<
" == norm_inf_mu_last = " << norm_inf_mu_last_
127 <<
"\nso we will take this as a signal to skip the update.\n";
129 perform_update =
false;
132 if( (
int)olevel >= (
int)PRINT_ALGORITHM_STEPS ) {
133 out <<
"\nmu_k " << mu_k <<
" != norm_inf_mu_last = " << norm_inf_mu_last_
134 <<
"\nso we will ignore this and perform the update anyway.\n";
140 if ( s.lambda().updated_k(0) ) {
142 if( (
int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
143 out <<
"\nUpdate the penalty parameter...\n";
147 &lambda_k = s.lambda().get_k(0).cv();
149 if( params->mu().size() != lambda_k.size() )
150 params->resize( lambda_k.size() );
155 max_lambda = norm_inf( lambda_k() ),
156 mult_fact = (1.0 + mult_factor_);
159 if( (
int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
160 out <<
"\nNear solution, forcing mu(j) >= mu_old(j)...\n";
162 DVector::const_iterator lb_itr = lambda_k.begin();
163 DVectorSlice::iterator mu_itr = mu.begin();
164 for( ; lb_itr != lambda_k.end(); ++mu_itr, ++ lb_itr )
165 *mu_itr = max( max( *mu_itr, mult_fact * ::fabs(*lb_itr) ), small_mu_ );
168 if( (
int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
169 out <<
"\nNot near solution, allowing reduction in mu(j) ...\n";
171 DVector::const_iterator lb_itr = lambda_k.begin();
172 DVectorSlice::iterator mu_itr = mu.begin();
173 for( ; lb_itr != lambda_k.end(); ++mu_itr, ++ lb_itr ) {
174 const value_type lb_j = ::fabs(*lb_itr);
176 (3.0 * (*mu_itr) + lb_j) / 4.0
177 , max( mult_fact * lb_j, small_mu_ )
180 value_type kkt_error = s.opt_kkt_err().get_k(0) + s.feas_kkt_err().get_k(0);
181 if(kkt_error <= kkt_near_sol_) {
182 if( (
int)olevel >= (
int)PRINT_ALGORITHM_STEPS ) {
183 out <<
"\nkkt_error = " << kkt_error <<
" <= kkt_near_sol = "
184 << kkt_near_sol_ << std::endl
185 <<
"Switching to forcing mu_k >= mu_km1 in the future\n";
187 near_solution_ =
true;
193 max_mu = norm_inf( mu() ),
194 min_mu = min_mu_ratio_ * max_mu;
195 for(DVectorSlice::iterator mu_itr = mu.begin(); mu_itr != mu.end(); ++mu_itr)
196 *mu_itr = max( (*mu_itr), min_mu );
198 s.mu().set_k(0) = norm_inf_mu_last_ = max_mu;
200 if( (
int)olevel >= (
int)PRINT_ALGORITHM_STEPS ) {
202 <<
"\nmin(|mu(j)|) = " << (*std::min_element( mu.begin(), mu.end() ))
206 if( (
int)olevel >= (int)PRINT_VECTORS ) {
207 out <<
"\nmu = \n" << mu;
211 if( (
int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
212 out <<
"\nDon't have the info to update penalty parameter so just use the last updated...\n";
218 direc_deriv->calc_deriv( s.Gf().get_k(0)(), s.c().get_k(0)(), s.d().get_k(0)() );
220 if( (
int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
221 out <<
"\nmu_k = " << s.mu().get_k(0) <<
"\n";
228 , poss_type step_poss, IterationPack::EDoStepType type, poss_type assoc_step_poss
229 , std::ostream& out,
const std::string& L )
const
232 << L <<
"*** Update the penalty parameter for the merit function to ensure\n"
233 << L <<
"*** a descent direction a directional derivatieve.\n"
234 << L <<
"*** phi is a merit function object that uses the penalty parameter mu.\n"
235 << L <<
"default: near_solution = false\n"
236 << L <<
" small_mu = " << small_mu_ << std::endl
237 << L <<
" min_mu_ratio = " << min_mu_ratio_ << std::endl
238 << L <<
" mult_factor = " << mult_factor_ << std::endl
239 << L <<
" kkt_near_sol = " << kkt_near_sol_ << std::endl
240 << L <<
"perform_update = true\n"
241 << L <<
"if mu_k is already updated then\n"
242 << L <<
" if mu_k == norm_inf_mu_last then\n"
243 << L <<
" *** We will use this as a signal to skip the update\n"
244 << L <<
" perform_update = false\n"
246 << L <<
" *** We will perform the update anyway\n"
248 << L <<
"if perform_update == true then\n"
249 << L <<
" if lambda_k is updated then\n"
250 << L <<
" max_lambda = norm(lambda_k,inf)\n"
251 << L <<
" mult_fact = (1+mult_factor)\n"
252 << L <<
" mu = phi.mu()\n"
253 << L <<
" if near_solution == true\n"
254 << L <<
" for j = 1...m\n"
255 << L <<
" mu(j) = max(max(mu(j),mult_fact*abs(lambda_k(j))),small_mu)\n"
258 << L <<
" for j = 1...m\n"
259 << L <<
" mu(j) = max( ( 3.0 * mu(j) + abs(lambda_k(j)) ) / 4.0\n"
260 << L <<
" , max( 1.001 * abs(lambda_k(j)) , small_mu ) )\n"
262 << L <<
" kkt_error = opt_kkt_err_k + feas_kkt_err_k\n"
263 << L <<
" if kkt_error <= kkt_near_sol then\n"
264 << L <<
" near_solution = true\n"
267 << L <<
" min_mu = min_mu_ratio * norm(mu,inf)\n"
268 << L <<
" for j = 1...m\n"
269 << L <<
" mu(j) = max( mu(j), min_mu )\n"
272 << L <<
" *** Don't have the information to perform the update.\n"
275 << L <<
"phi.calc_deriv(Gf_k,c_k,d_k)\n";
297 return min_mu_ratio_;
317 return kkt_near_sol_;
MeritFunc_PenaltyParamsUpdateWithMult_AddedStep(const merit_func_ptr_t &merit_func, value_type small_mu=1e-6, value_type min_mu_ratio=1e-8, value_type mult_factor=1e-4, value_type kkt_near_sol=1e-1)
value_type max_element(const Vector &v)
value_type mult_factor() const
bool do_step(Algorithm &algo, poss_type step_poss, IterationPack::EDoStepType type, poss_type assoc_step_poss)
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
value_type kkt_near_sol() const
value_type min_mu_ratio() const
value_type small_mu() const
std::string typeName(const T &t)