67 class RemoveFilterEntryPred {
70 RemoveFilterEntryPred(
75 :f_with_boundary_(f_with_boundary),
76 theta_with_boundary_(theta_with_boundary),
83 entry.
f >= f_with_boundary_
85 entry.
theta >= theta_with_boundary_
90 <<
"\nRemoving from the filter the redundant point:"
91 <<
"\n f_with_boundary = " << entry.
f
92 <<
"\n theta_with_boundary = " << entry.
theta
93 <<
"\n iteration added = " << entry.
iter
113 namespace MoochoPack {
117 {
return (x < y) ? x : y; }
120 {
return (x > y) ? x : y; }
124 ,
const std::string obj_iq_name
125 ,
const std::string grad_obj_iq_name
139 gamma_theta_(gamma_theta),
142 gamma_alpha_(gamma_alpha),
146 theta_small_fact_(theta_small_fact),
147 theta_max_(theta_max),
149 back_track_frac_(back_track_frac),
152 grad_obj_f_(grad_obj_iq_name),
158 "Null nlp passed to LineSearchFilter_Step constructor"
161 #if defined(FILTER_DEBUG_OUT)
163 fout <<
"<FilterDebugDocument>" << std::endl;
170 #if defined(FILTER_DEBUG_OUT)
172 fout <<
"</FilterDebugDocument>" << std::endl;
214 IterQuantityAccess<value_type>
216 &alpha_iq = s.alpha();
218 IterQuantityAccess<VectorMutable>
220 *c_iq = m > 0 ? &s.c() : NULL,
225 if (!s.d().updated_k(0) || !x_iq.updated_k(0))
232 if (!alpha_iq.updated_k(0) || alpha_iq.get_k(0) > 1 || alpha_iq.get_k(0) <= 0)
244 ? Gf_iq.get_k(0).inner_product( s.d().get_k(0) )
245 : dyn_cast<NLPInterfacePack::NLPObjGrad>(*nlp_).calc_Gf_prod(
246 x_iq.get_k(0),s.d().get_k(0)
251 out <<
"\ntheta_k = ||c_k||1/c_k.dim() = " << theta_k << std::endl;
253 const value_type theta_small = theta_small_fact_ *
MAX(1.0,theta_k);
263 out <<
"\nBeginning Filter line search method.\n"
264 <<
"\nCurrent Filter\n"
265 <<
"-----------------------------------------------------\n"
266 <<
"|" << setw(25) <<
"f_with_boundary "
267 <<
"|" << setw(25) <<
"theta_with_boundary "
269 <<
"-----------------------------------------------------\n";
271 IterQuantityAccess<Filter_T>& filter_iq =
filter_(s);
273 if (filter_iq.updated_k(-1)) {
274 Filter_T& filter = filter_iq.get_k(-1);
275 if (!filter.empty()) {
277 Filter_T::iterator entry = filter.begin();
278 entry != filter.end();
282 out <<
"|" << setw(25) << entry->
f
283 <<
" " << setw(25) << entry->
theta
288 out <<
"Filter is empty.\n";
292 out <<
"Filter is empty.\n";
296 out <<
"\n Iteration Status\n";
297 out <<
"----------------------------------------------------------------------------------------------------------\n";
299 out <<
"|" << setw(w) <<
"alpha_k "
300 <<
"|" << setw(w) <<
"f_kp1 "
301 <<
"|" << setw(w) <<
"theta_kp1 "
302 <<
"|" << setw(w) <<
"pt. status "
303 <<
"|" << setw(40) <<
"comment "
306 out <<
"----------------------------------------------------------------------------------------------------------"
311 bool augment_filter =
false;
312 bool accepted =
false;
313 while (alpha_k > alpha_min && !accepted)
327 out <<
"|" << setw(w) <<
" --- "
328 <<
" " << setw(w) <<
" --- "
329 <<
" " << setw(w) <<
" --- "
330 <<
" " << setw(w) <<
" failed "
331 <<
" " << setw(40) <<
" nan_or_inf in calc"
349 out <<
"|" << setw(w) << alpha_k
350 <<
" " << setw(w) << f_iq.get_k(+1)
351 <<
" " << setw(w) << theta_kp1;
359 out <<
" " << setw(w) <<
"failed"
360 <<
" " << setw(40) <<
"Unacceptable to filter"
367 if (theta_kp1 > theta_max_) {
371 out <<
" " << setw(w) <<
"failed"
372 <<
" " << setw(40) <<
"theta_kp1 > theta_max"
391 { out <<
" " << setw(w) <<
"accepted"; }
393 { out <<
" " << setw(w) <<
"failed"; }
395 out <<
" " << setw(40) <<
"Switch Cond. Holds (Armijo)" <<
"|" << std::endl;
402 { augment_filter =
true; }
408 { out <<
" " << setw(w) <<
"accepted"; }
410 { out <<
" " << setw(w) <<
"failed"; }
412 out <<
" " << setw(40) <<
"Fraction Reduction (! Switch Cond )" <<
"|" << std::endl;
422 alpha_k = alpha_k*back_track_frac_;
434 out <<
"\nPoint was accepted - augmenting the filter ...\n";
436 out <<
"Point was accepted - did NOT augment filter.\n";
439 if (augment_filter) {
458 ,
"FilterLineSearchFailure : Should go to restoration"
462 if( static_cast<int>(olevel) >= static_cast<int>(
PRINT_VECTORS) )
464 out <<
"\nx_kp1 =\n" << x_iq.get_k(+1);
466 { out <<
"\nc_kp1 =\n" << c_iq->get_k(+1); }
468 { out <<
"\nh_kp1 =\n" << h_iq->get_k(+1); }
471 #if defined(FILTER_DEBUG_OUT)
473 fout <<
" <FilterIteration iter=\"" << s.
k() <<
"\">" << std::endl;
474 fout <<
" <SelectedPoint alpha=\"" << alpha_k
475 <<
"\" f=\"" << f_iq.get_k(+1)
476 <<
"\" theta=\"" << theta_kp1
477 <<
"\" />" << std::endl;
480 fout <<
" <Filter>" << std::endl;
482 IterQuantityAccess<Filter_T>& filter_iq =
filter_(s);
483 if (filter_iq.updated_k(0))
485 Filter_T& current_filter = filter_iq.get_k(0);
487 Filter_T::iterator entry = current_filter.begin();
488 entry != current_filter.end();
492 fout <<
" <FilterPoint iter=\"" << entry->
iter
493 <<
"\" f=\"" << entry->
f
494 <<
"\" theta=\"" << entry->
theta <<
">" << std::endl;
499 fout <<
" <FilterNotUpdated/>" << std::endl;
502 fout <<
" </Filter>" << std::endl;
506 fout <<
" <AlphaCurve>" << std::endl;
508 for (
int i = 0; i < 10 || alpha_tmp > alpha_k; ++i)
514 fout <<
" <AlphaPoint "
515 <<
"alpha=\"" << alpha_tmp <<
"\" "
516 <<
"f=\"" << f_iq.get_k(+1) <<
"\" "
517 <<
"theta=\"" << theta <<
">" << std::endl;
520 alpha_tmp=alpha_tmp*back_track_frac_;
526 fout <<
" </AlphaCurve>" << std::endl;
528 fout <<
" </FilterIteration" << std::endl;
539 ,
poss_type assoc_step_poss, std::ostream&
out,
const std::string& L
543 << L <<
"*** Filter line search method\n"
544 << L <<
"# Assumes initial d_k & alpha_k (0-1) is known and\n"
545 << L <<
"# x_kp1, f_kp1, c_kp1 and h_kp1 are calculated for that alpha_k\n"
546 << L <<
"Gf_t_dk = <Gf,dk>\n"
547 << L <<
"theta_k = norm_1(c_k)/c_k.dim()\n"
548 << L <<
"theta_small = theta_small_fact*max(1.0,theta_k)\n"
549 << L <<
"if f_min != F_MIN_UNBOUNDED then\n"
550 << L <<
" gamma_f_used = gamma_f * (f_k-f_min)/theta_k\n"
553 << L <<
" gamma_f_used = gamma_f\n"
555 << L <<
"if Gf_t_dk < 0 then\n"
556 << L <<
" alpha_min = min(gamma_theta, gamma_f_used*theta_k/(-Gf_t_dk))\n"
557 << L <<
" if theta_k <= theta_small then\n"
558 << L <<
" alpha_min = min(alpha_min, delta_*(theta_k^s_theta)/((-Gf_t_dk)^s_f))\n"
561 << L <<
" alpha_min = gamma_theta\n"
563 << L <<
"alpha_min = alpha_min*gamma_alpha\n"
564 << L <<
"# Start the line search\n"
565 << L <<
"accepted = false\n"
566 << L <<
"augment = false\n"
567 << L <<
"while alpha > alpha_min and accepted = false then\n"
568 << L <<
" accepted = true"
569 << L <<
" if any values in x_kp1, f_kp1, c_kp1, h_kp1 are nan or inf then\n"
570 << L <<
" accepted = false\n"
572 << L <<
" # Check filter\n"
573 << L <<
" if accepted = true then\n"
574 << L <<
" theta_kp1 = norm_1(c_kp1)/c_kp1.dim()\n"
575 << L <<
" for each pt in the filter do\n"
576 << L <<
" if theta_kp1 >= pt.theta and f_kp1 >= pt.f then\n"
577 << L <<
" accepted = false\n"
578 << L <<
" break for\n"
582 << L <<
" #Check Sufficient Decrease\n"
583 << L <<
" if accepted = true then"
584 << L <<
" # if switching condition is satisfied, use Armijo on f\n"
585 << L <<
" if theta_k < theta_small and Gf_t_dk < 0 and\n"
586 << L <<
" ((-Gf_t_dk)^s_f)*alpha_k < delta*theta_k^s_theta then\n"
587 << L <<
" if f_kp1 <= f_k + eta_f*alpha_k*Gf_t_dk then\n"
588 << L <<
" accepted = true\n"
591 << L <<
" # Verify factional reduction\n"
592 << L <<
" if theta_kp1 < (1-gamma_theta)*theta_k or f_kp1 < f_k - gamma_f_used*theta_k then\n"
593 << L <<
" accepted = true\n"
594 << L <<
" augment = true\n"
598 << L <<
" if accepted = false then\n"
599 << L <<
" alpha_k = alpha_k*back_track_frac\n"
600 << L <<
" x_kp1 = x_k + alpha_k * d_k\n"
601 << L <<
" f_kp1 = f(x_kp1)\n"
602 << L <<
" c_kp1 = c(x_kp1)\n"
603 << L <<
" h_kp1 = h(x_kp1)\n"
605 << L <<
"end while\n"
606 << L <<
"if accepted = true then\n"
607 << L <<
" if augment = true then\n"
608 << L <<
" Augment the filter (use f_with_boundary and theta_with_boundary)\n"
611 << L <<
" goto the restoration phase\n"
616 const IterQuantityAccess<VectorMutable>& x,
617 const IterQuantityAccess<value_type>& f,
618 const IterQuantityAccess<VectorMutable>* c,
619 const IterQuantityAccess<VectorMutable>* h,
620 const bool throw_excpt
639 const AlgorithmState& s
642 bool accepted =
true;
644 const IterQuantityAccess<Filter_T>& filter_iq =
filter_(s);
646 if (filter_iq.updated_k(-1))
648 const Filter_T ¤t_filter = filter_iq.get_k(-1);
651 Filter_T::const_iterator entry = current_filter.begin();
652 entry != current_filter.end();
656 if (f >= entry->
f && theta >= entry->
theta)
670 const IterQuantityAccess<value_type>& f_iq
673 bool accepted =
false;
676 double f_kp1 = f_iq.get_k(+1);
677 double f_k = f_iq.get_k(0);
678 double lhs = f_k - f_kp1;
679 double rhs = -eta_f_*alpha_k*Gf_t_dk;
690 const IterQuantityAccess<value_type>& f_iq,
696 bool accepted =
false;
697 if (theta_kp1 <= (1-gamma_theta_)*theta_k
698 || f_iq.get_k(+1) <= f_iq.get_k(0)-gamma_f_used*theta_k )
707 const VectorMutable& d,
709 IterQuantityAccess<VectorMutable> &x,
710 IterQuantityAccess<value_type>& f,
711 IterQuantityAccess<VectorMutable>* c,
712 IterQuantityAccess<VectorMutable>* h,
717 VectorMutable& x_kp1 = x.set_k(+1);
719 Vp_StV( &x_kp1, alpha, d);
724 nlp.unset_quantities();
725 nlp.set_f( &f.set_k(+1) );
726 if (c) nlp.set_c( &c->set_k(+1) );
728 if (c) nlp.calc_c( x_kp1,
false );
729 nlp.unset_quantities();
744 alpha_min =
MIN(gamma_theta_, gamma_f_used*theta_k/(-Gf_t_dk));
745 if (theta_k <= theta_small)
748 alpha_min =
MIN(alpha_min, switch_bound);
753 alpha_min = gamma_theta_;
755 return alpha_min * gamma_alpha_;
759 const IterQuantityAccess<VectorMutable>* c,
760 const IterQuantityAccess<VectorMutable>* h,
769 const Vector &c_k = c->get_k(k);
770 theta += ( c_k.norm_1() / c_k.dim() );
776 const IterQuantityAccess<value_type> &f,
784 out <<
"\nf_min==F_MIN_UNBOUNDED: Setting gamma_f_used = gamma_f = "<<gamma_f()<<
".\n";
788 if( f_k < f_min() ) {
790 out <<
"\nWarning!!! f_k = "<<f_k<<
" < f_min = "<<f_min()<<
"! Setting gamma_f_used = gamma_f = "<<gamma_f()<<
".\n";
793 const value_type gamma_f_used = gamma_f() * ( f_k - f_min() ) / theta_k;
795 out <<
"\nf_min = "<<f_min()<<
"!=F_MIN_UNBOUNDED: Setting gamma_f_used = gamma_f * (f_k - f_min)/theta_k = "
796 <<gamma_f()<<
" * ("<<f_k<<
" - "<<f_min()<<
")/"<<theta_k<<
" = "<<gamma_f_used <<
".\n";
807 if (theta_k < theta_small && Gf_t_dk < 0) {
808 if (
pow(-Gf_t_dk, s_f_)*alpha_k - delta_*
pow(theta_k, s_theta_) > 0) {
819 IterQuantityAccess<Filter_T>& filter_iq =
filter_(s);
821 if (!filter_iq.updated_k(0))
823 if (filter_iq.updated_k(-1))
826 filter_iq.set_k(0,-1);
850 f_with_boundary = f-gamma_f_used*theta,
851 theta_with_boundary = (1.0-gamma_theta())*theta;
854 out <<
"\nAugmenting the filter with the point:"
855 <<
"\n f_with_boundary = f_kp1 - gamma_f_used*theta_kp1 = "<<f<<
" - "<<gamma_f_used<<
"*"<<theta<<
" = " <<f_with_boundary
856 <<
"\n theta_with_boundary = (1-gamma_theta)*theta_kp1 = (1-"<<gamma_theta()<<
")*"<<theta<<
" = "<<theta_with_boundary
864 current_filter.remove_if(
865 RemoveFilterEntryPred(
866 f_with_boundary, theta_with_boundary,
873 current_filter.push_front(
FilterEntry(f_with_boundary, theta_with_boundary, s.
k()));
void UpdateFilter(IterationPack::AlgorithmState &s) const
AbstractLinAlgPack::size_type size_type
bool is_null(const boost::shared_ptr< T > &p)
const std::string FILTER_IQ_STRING
bool CheckFractionalReduction(const IterQuantityAccess< value_type > &f_iq, const value_type gamma_f_used, const value_type theta_kp1, const value_type theta_k) const
void pow(DVectorSlice *vs_lhs, const DVectorSlice &vs_rhs1, const DVectorSlice &vs_rhs2)
vs_lhs = pow(vs_rhs1,vs_rhs2)
value_type CalculateAlphaMin(const value_type gamma_f_used, const value_type Gf_t_dk, const value_type theta_k, const value_type theta_small) const
bool CheckArmijo(const value_type Gf_t_dk, const value_type alpha_k, const IterQuantityAccess< value_type > &f_iq) const
~LineSearchFilter_Step()
Destructor.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
void Vp_StV(VectorMutable *v_lhs, const value_type &alpha, const Vector &v_rhs)
v_lhs = alpha * v_rhs + v_lhs
CastIQMember< Filter_T > filter_
rSQP Algorithm control class.
value_type CalculateGammaFUsed(const IterQuantityAccess< value_type > &f, const value_type theta_k, const EJournalOutputLevel olevel, std::ostream &out) const
LineSearchFilter_Step(Teuchos::RCP< NLPInterfacePack::NLP > nlp, const std::string obj_iq_name="f", const std::string grad_obj_iq_name="Gf", const value_type &gamma_theta=1e-5, const value_type &gamma_f=1e-5, const value_type &f_min=F_MIN_UNBOUNDED, const value_type &gamma_alpha=5e-2, const value_type &delta=1e-4, const value_type &s_theta=1.1, const value_type &s_f=2.3, const value_type &theta_small_fact=1e-4, const value_type &theta_max=1e10, const value_type &eta_f=1e-4, const value_type &back_track_frac=0.5)
Constructor.
value_type MIN(value_type x, value_type y)
bool ValidatePoint(const IterQuantityAccess< VectorMutable > &x, const IterQuantityAccess< value_type > &f, const IterQuantityAccess< VectorMutable > *c, const IterQuantityAccess< VectorMutable > *h, const bool throw_excpt) const
CastIQMember< VectorMutable > grad_obj_f_
ITeration quantity access for objective gradient.
bool do_step(Algorithm &algo, poss_type step_poss, IterationPack::EDoStepType type, poss_type assoc_step_poss)
std::list< FilterEntry > Filter_T
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
const f_int f_dbl_prec const f_int f_int const f_int f_int const f_dbl_prec f_int f_int f_dbl_prec w[]
virtual std::ostream & journal_out() const
Return a reference to a std::ostream to be used to output debug information and the like...
value_type CalculateTheta_k(const IterQuantityAccess< VectorMutable > *c, const IterQuantityAccess< VectorMutable > *h, int k) const
Thrown if a line search failure occurs.
EJournalOutputLevel
enum for journal output.
T_To & dyn_cast(T_From &from)
Reduced space SQP state encapsulation interface.
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
void print_algorithm_step(const Algorithm &algo, Algorithm::poss_type step_poss, EDoStepType type, Algorithm::poss_type assoc_step_poss, std::ostream &out)
Prints to 'out' the algorithm step.
Teuchos::RCP< NLPInterfacePack::NLP > nlp_
Abstacts a set of iteration quantities for an iterative algorithm.
void AugmentFilter(const value_type gamma_f_used, const value_type f_with_boundary, const value_type theta_with_boundary, IterationPack::AlgorithmState &s, const EJournalOutputLevel olevel, std::ostream &out) const
AlgorithmTracker & track()
TypeTo as(const TypeFrom &t)
bool assert_print_nan_inf(const value_type &val, const char name[], bool throw_excpt, std::ostream *out)
This function asserts if a value_type scalare is a NaN or Inf and optionally prints out these entires...
CastIQMember< value_type > obj_f_
Iteration quantity access for objective value.
const f_dbl_prec const f_int const f_int const f_int f_dbl_prec rhs[]
bool ShouldSwitchToArmijo(const value_type Gf_t_dk, const value_type alpha_k, const value_type theta_k, const value_type theta_small) const
static value_type F_MIN_UNBOUNDED
bool CheckFilterAcceptability(const value_type f, const value_type theta, const AlgorithmState &s) const
AbstractLinAlgPack::value_type value_type
Acts as the central hub for an iterative algorithm.
NLPAlgoState & rsqp_state()
<<std aggr>="">> members for algo_cntr
void UpdatePoint(const VectorMutable &d, value_type alpha, IterQuantityAccess< VectorMutable > &x, IterQuantityAccess< value_type > &f, IterQuantityAccess< VectorMutable > *c, IterQuantityAccess< VectorMutable > *h, NLP &nlp) const
virtual size_type m() const
Return the number of general equality constraints.
NLPAlgo & rsqp_algo(Algorithm &algo)
Convert from a Algorithm to a NLPAlgo.
value_type MAX(value_type x, value_type y)