44#ifndef ROL_PROBLEM_DEF_HPP
45#define ROL_PROBLEM_DEF_HPP
51template<
typename Real>
53 const Ptr<Vector<Real>> &x,
54 const Ptr<Vector<Real>> &g)
55 : isFinalized_(false), hasBounds_(false),
56 hasEquality_(false), hasInequality_(false),
57 hasLinearEquality_(false), hasLinearInequality_(false),
58 cnt_econ_(0), cnt_icon_(0), cnt_linear_econ_(0), cnt_linear_icon_(0),
59 obj_(nullPtr), xprim_(nullPtr), xdual_(nullPtr), bnd_(nullPtr),
60 con_(nullPtr), mul_(nullPtr), res_(nullPtr), proj_(nullPtr),
66 INPUT_linear_con_.clear();
67 if (g==nullPtr) INPUT_xdual_ = x->dual().clone();
68 else INPUT_xdual_ = g;
71template<
typename Real>
72void Problem<Real>::addBoundConstraint(
const Ptr<BoundConstraint<Real>> &bnd) {
73 ROL_TEST_FOR_EXCEPTION(isFinalized_,std::invalid_argument,
74 ">>> ROL::Problem: Cannot add bounds after problem is finalized!");
80template<
typename Real>
81void Problem<Real>::removeBoundConstraint() {
82 ROL_TEST_FOR_EXCEPTION(isFinalized_,std::invalid_argument,
83 ">>> ROL::Problem: Cannot remove bounds after problem is finalized!");
89template<
typename Real>
90void Problem<Real>::addConstraint( std::string name,
91 const Ptr<Constraint<Real>> &econ,
92 const Ptr<Vector<Real>> &emul,
93 const Ptr<Vector<Real>> &eres,
95 ROL_TEST_FOR_EXCEPTION(isFinalized_,std::invalid_argument,
96 ">>> ROL::Problem: Cannot add constraint after problem is finalized!");
98 if (reset) INPUT_con_.clear();
100 auto it = INPUT_con_.find(name);
101 ROL_TEST_FOR_EXCEPTION(it != INPUT_con_.end(),std::invalid_argument,
102 ">>> ROL::Problem: Constraint names must be distinct!");
104 INPUT_con_.insert({name,ConstraintData<Real>(econ,emul,eres)});
109template<
typename Real>
110void Problem<Real>::addConstraint( std::string name,
111 const Ptr<Constraint<Real>> &icon,
112 const Ptr<Vector<Real>> &imul,
113 const Ptr<BoundConstraint<Real>> &ibnd,
114 const Ptr<Vector<Real>> &ires,
116 ROL_TEST_FOR_EXCEPTION(isFinalized_,std::invalid_argument,
117 ">>> ROL::Problem: Cannot add constraint after problem is finalized!");
119 if (reset) INPUT_con_.clear();
121 auto it = INPUT_con_.find(name);
122 ROL_TEST_FOR_EXCEPTION(it != INPUT_con_.end(),std::invalid_argument,
123 ">>> ROL::Problem: Constraint names must be distinct!");
125 INPUT_con_.insert({name,ConstraintData<Real>(icon,imul,ires,ibnd)});
126 hasInequality_ =
true;
130template<
typename Real>
131void Problem<Real>::removeConstraint(std::string name) {
132 ROL_TEST_FOR_EXCEPTION(isFinalized_,std::invalid_argument,
133 ">>> ROL::Problem: Cannot remove constraint after problem is finalized!");
135 auto it = INPUT_con_.find(name);
136 if (it!=INPUT_con_.end()) {
137 if (it->second.bounds==nullPtr) cnt_econ_--;
139 INPUT_con_.erase(it);
141 if (cnt_econ_==0) hasEquality_ =
false;
142 if (cnt_icon_==0) hasInequality_ =
false;
145template<
typename Real>
146void Problem<Real>::addLinearConstraint( std::string name,
147 const Ptr<Constraint<Real>> &linear_econ,
148 const Ptr<Vector<Real>> &linear_emul,
149 const Ptr<Vector<Real>> &linear_eres,
151 ROL_TEST_FOR_EXCEPTION(isFinalized_,std::invalid_argument,
152 ">>> ROL::Problem: Cannot add linear constraint after problem is finalized!");
154 if (reset) INPUT_linear_con_.clear();
156 auto it = INPUT_linear_con_.find(name);
157 ROL_TEST_FOR_EXCEPTION(it != INPUT_linear_con_.end(),std::invalid_argument,
158 ">>> ROL::Problem: Linear constraint names must be distinct!");
160 INPUT_linear_con_.insert({name,ConstraintData<Real>(linear_econ,linear_emul,linear_eres)});
161 hasLinearEquality_ =
true;
165template<
typename Real>
166void Problem<Real>::addLinearConstraint( std::string name,
167 const Ptr<Constraint<Real>> &linear_icon,
168 const Ptr<Vector<Real>> &linear_imul,
169 const Ptr<BoundConstraint<Real>> &linear_ibnd,
170 const Ptr<Vector<Real>> &linear_ires,
172 ROL_TEST_FOR_EXCEPTION(isFinalized_,std::invalid_argument,
173 ">>> ROL::Problem: Cannot add linear constraint after problem is finalized!");
175 if (reset) INPUT_linear_con_.clear();
177 auto it = INPUT_linear_con_.find(name);
178 ROL_TEST_FOR_EXCEPTION(it != INPUT_linear_con_.end(),std::invalid_argument,
179 ">>> ROL::Problem: Linear constraint names must be distinct!");
181 INPUT_linear_con_.insert({name,ConstraintData<Real>(linear_icon,linear_imul,linear_ires,linear_ibnd)});
182 hasLinearInequality_ =
true;
186template<
typename Real>
187void Problem<Real>::removeLinearConstraint(std::string name) {
188 ROL_TEST_FOR_EXCEPTION(isFinalized_,std::invalid_argument,
189 ">>> ROL::Problem: Cannot remove linear inequality after problem is finalized!");
191 auto it = INPUT_linear_con_.find(name);
192 if (it!=INPUT_linear_con_.end()) {
193 if (it->second.bounds==nullPtr) cnt_linear_econ_--;
194 else cnt_linear_icon_--;
195 INPUT_linear_con_.erase(it);
197 if (cnt_linear_econ_==0) hasLinearEquality_ =
false;
198 if (cnt_linear_icon_==0) hasLinearInequality_ =
false;
201template<
typename Real>
202void Problem<Real>::setProjectionAlgorithm(ParameterList &list) {
203 ROL_TEST_FOR_EXCEPTION(isFinalized_,std::invalid_argument,
204 ">>> ROL::Problem: Cannot set polyhedral projection algorithm after problem is finalized!");
209template<
typename Real>
210void Problem<Real>::finalize(
bool lumpConstraints,
bool printToStream, std::ostream &outStream) {
212 std::unordered_map<std::string,ConstraintData<Real>> con, lcon, icon;
213 bool hasEquality = hasEquality_;
214 bool hasLinearEquality = hasLinearEquality_;
215 bool hasInequality = hasInequality_;
216 bool hasLinearInequality = hasLinearInequality_;
217 con.insert(INPUT_con_.begin(),INPUT_con_.end());
218 if (lumpConstraints) {
219 con.insert(INPUT_linear_con_.begin(),INPUT_linear_con_.end());
220 hasEquality = (hasEquality || hasLinearEquality);
221 hasInequality = (hasInequality || hasLinearInequality);
222 hasLinearEquality =
false;
223 hasLinearInequality =
false;
226 lcon.insert(INPUT_linear_con_.begin(),INPUT_linear_con_.end());
230 if (!hasLinearEquality && !hasLinearInequality) {
232 if (!hasEquality && !hasInequality && !hasBounds_) {
235 xprim_ = INPUT_xprim_;
236 xdual_ = INPUT_xdual_;
242 else if (!hasEquality && !hasInequality && hasBounds_) {
245 xprim_ = INPUT_xprim_;
246 xdual_ = INPUT_xdual_;
252 else if (hasEquality && !hasInequality && !hasBounds_) {
253 ConstraintAssembler<Real> cm(con,INPUT_xprim_,INPUT_xdual_);
256 xprim_ = INPUT_xprim_;
257 xdual_ = INPUT_xdual_;
259 con_ = cm.getConstraint();
260 mul_ = cm.getMultiplier();
261 res_ = cm.getResidual();
264 ConstraintAssembler<Real> cm(con,INPUT_xprim_,INPUT_xdual_,INPUT_bnd_);
267 if (cm.hasInequality()) {
268 obj_ = makePtr<SlacklessObjective<Real>>(INPUT_obj_);
270 xprim_ = cm.getOptVector();
271 xdual_ = cm.getDualOptVector();
272 bnd_ = cm.getBoundConstraint();
273 con_ = cm.getConstraint();
274 mul_ = cm.getMultiplier();
275 res_ = cm.getResidual();
279 if (!hasBounds_ && !hasLinearInequality) {
280 ConstraintAssembler<Real> cm(lcon,INPUT_xprim_,INPUT_xdual_);
281 xfeas_ = cm.getOptVector()->clone(); xfeas_->set(*cm.getOptVector());
282 rlc_ = makePtr<ReduceLinearConstraint<Real>>(cm.getConstraint(),xfeas_,cm.getResidual());
284 if (!hasEquality && !hasInequality) {
286 obj_ = rlc_->transform(INPUT_obj_);
287 xprim_ = xfeas_->clone(); xprim_->zero();
288 xdual_ = cm.getDualOptVector();
295 for (
auto it = con.begin(); it != con.end(); ++it) {
296 icon.insert(std::pair<std::string,ConstraintData<Real>>(it->first,
297 ConstraintData<Real>(rlc_->transform(it->second.constraint),
298 it->second.multiplier,it->second.residual,it->second.bounds)));
300 Ptr<Vector<Real>> xtmp = xfeas_->clone(); xtmp->zero();
301 ConstraintAssembler<Real> cm1(icon,xtmp,cm.getDualOptVector());
302 xprim_ = cm1.getOptVector();
303 xdual_ = cm1.getDualOptVector();
304 con_ = cm1.getConstraint();
305 mul_ = cm1.getMultiplier();
306 res_ = cm1.getResidual();
307 if (!hasInequality) {
309 obj_ = rlc_->transform(INPUT_obj_);
314 obj_ = makePtr<SlacklessObjective<Real>>(rlc_->transform(INPUT_obj_));
315 bnd_ = cm1.getBoundConstraint();
319 else if ((hasBounds_ || hasLinearInequality) && !hasEquality && !hasInequality) {
320 ConstraintAssembler<Real> cm(lcon,INPUT_xprim_,INPUT_xdual_,INPUT_bnd_);
323 if (cm.hasInequality()) {
324 obj_ = makePtr<SlacklessObjective<Real>>(INPUT_obj_);
326 xprim_ = cm.getOptVector();
327 xdual_ = cm.getDualOptVector();
328 bnd_ = cm.getBoundConstraint();
332 proj_ = PolyhedralProjectionFactory<Real>(*xprim_,*xdual_,bnd_,
333 cm.getConstraint(),*cm.getMultiplier(),*cm.getResidual(),ppa_list_);
336 ConstraintAssembler<Real> cm(con,lcon,INPUT_xprim_,INPUT_xdual_,INPUT_bnd_);
339 if (cm.hasInequality()) {
340 obj_ = makePtr<SlacklessObjective<Real>>(INPUT_obj_);
342 xprim_ = cm.getOptVector();
343 xdual_ = cm.getDualOptVector();
344 con_ = cm.getConstraint();
345 mul_ = cm.getMultiplier();
346 res_ = cm.getResidual();
347 bnd_ = cm.getBoundConstraint();
348 proj_ = PolyhedralProjectionFactory<Real>(*xprim_,*xdual_,bnd_,
349 cm.getLinearConstraint(),*cm.getLinearMultiplier(),
350 *cm.getLinearResidual(),ppa_list_);
355 outStream << std::endl;
356 outStream <<
" ROL::Problem::finalize" << std::endl;
357 outStream <<
" Problem Summary:" << std::endl;
358 outStream <<
" Has Bound Constraint? .............. " << (hasBounds_ ?
"yes" :
"no") << std::endl;
359 outStream <<
" Has Equality Constraint? ........... " << (hasEquality ?
"yes" :
"no") << std::endl;
362 for (
auto it = con.begin(); it != con.end(); ++it) {
363 if (it->second.bounds==nullPtr) {
365 outStream <<
" Names: ........................... ";
371 outStream << it->first << std::endl;
374 outStream <<
" Total: ........................... " << cnt_econ_+(lumpConstraints ? cnt_linear_econ_ : 0) << std::endl;
376 outStream <<
" Has Inequality Constraint? ......... " << (hasInequality ?
"yes" :
"no") << std::endl;
379 for (
auto it = con.begin(); it != con.end(); ++it) {
380 if (it->second.bounds!=nullPtr) {
382 outStream <<
" Names: ........................... ";
388 outStream << it->first << std::endl;
391 outStream <<
" Total: ........................... " << cnt_icon_+(lumpConstraints ? cnt_linear_icon_ : 0) << std::endl;
393 if (!lumpConstraints) {
394 outStream <<
" Has Linear Equality Constraint? .... " << (hasLinearEquality ?
"yes" :
"no") << std::endl;
395 if (hasLinearEquality) {
397 for (
auto it = lcon.begin(); it != lcon.end(); ++it) {
398 if (it->second.bounds==nullPtr) {
400 outStream <<
" Names: ........................... ";
406 outStream << it->first << std::endl;
409 outStream <<
" Total: ........................... " << cnt_linear_econ_ << std::endl;
411 outStream <<
" Has Linear Inequality Constraint? .. " << (hasLinearInequality ?
"yes" :
"no") << std::endl;
412 if (hasLinearInequality) {
414 for (
auto it = lcon.begin(); it != lcon.end(); ++it) {
415 if (it->second.bounds!=nullPtr) {
417 outStream <<
" Names: ........................... ";
423 outStream << it->first << std::endl;
426 outStream <<
" Total: ........................... " << cnt_linear_icon_ << std::endl;
429 outStream << std::endl;
434 outStream << std::endl;
435 outStream <<
" ROL::Problem::finalize" << std::endl;
436 outStream <<
" Problem already finalized!" << std::endl;
437 outStream << std::endl;
442template<
typename Real>
443const Ptr<Objective<Real>>& Problem<Real>::getObjective() {
448template<
typename Real>
449const Ptr<Vector<Real>>& Problem<Real>::getPrimalOptimizationVector() {
454template<
typename Real>
455const Ptr<Vector<Real>>& Problem<Real>::getDualOptimizationVector() {
460template<
typename Real>
461const Ptr<BoundConstraint<Real>>& Problem<Real>::getBoundConstraint() {
466template<
typename Real>
467const Ptr<Constraint<Real>>& Problem<Real>::getConstraint() {
472template<
typename Real>
473const Ptr<Vector<Real>>& Problem<Real>::getMultiplierVector() {
478template<
typename Real>
479const Ptr<Vector<Real>>& Problem<Real>::getResidualVector() {
484template<
typename Real>
485const Ptr<PolyhedralProjection<Real>>& Problem<Real>::getPolyhedralProjection() {
490template<
typename Real>
491EProblem Problem<Real>::getProblemType() {
496template<
typename Real>
497Real Problem<Real>::checkLinearity(
bool printToStream, std::ostream &outStream)
const {
498 std::ios_base::fmtflags state(outStream.flags());
500 outStream << std::setprecision(8) << std::scientific;
501 outStream << std::endl;
502 outStream <<
" ROL::Problem::checkLinearity" << std::endl;
504 const Real one(1), two(2), eps(1e-2*std::sqrt(ROL_EPSILON<Real>()));
505 Real tol(std::sqrt(ROL_EPSILON<Real>())), cnorm(0), err(0), maxerr(0);
506 Ptr<Vector<Real>> x = INPUT_xprim_->clone(); x->randomize(-one,one);
507 Ptr<Vector<Real>> y = INPUT_xprim_->clone(); y->randomize(-one,one);
508 Ptr<Vector<Real>> z = INPUT_xprim_->clone(); z->zero();
509 Ptr<Vector<Real>> xy = INPUT_xprim_->clone();
510 Real alpha = two*
static_cast<Real
>(rand())/
static_cast<Real
>(RAND_MAX)-one;
511 xy->set(*x); xy->axpy(alpha,*y);
512 Ptr<Vector<Real>> c1, c2;
513 for (
auto it = INPUT_linear_con_.begin(); it != INPUT_linear_con_.end(); ++it) {
514 c1 = it->second.residual->clone();
515 c2 = it->second.residual->clone();
516 it->second.constraint->update(*xy,UpdateType::Temp);
517 it->second.constraint->value(*c1,*xy,tol);
519 it->second.constraint->update(*x,UpdateType::Temp);
520 it->second.constraint->value(*c2,*x,tol);
522 it->second.constraint->update(*y,UpdateType::Temp);
523 it->second.constraint->value(*c2,*y,tol);
524 c1->axpy(-alpha,*c2);
525 it->second.constraint->update(*z,UpdateType::Temp);
526 it->second.constraint->value(*c2,*z,tol);
529 maxerr = std::max(err,maxerr);
531 outStream <<
" Constraint " << it->first;
532 outStream <<
": ||c(x+alpha*y) - (c(x)+alpha*(c(y)-c(0)))|| = " << err << std::endl;
533 if (err > eps*cnorm) {
534 outStream <<
" Constraint " << it->first <<
" may not be linear!" << std::endl;
539 outStream << std::endl;
541 outStream.flags(state);
545template<
typename Real>
546void Problem<Real>::checkVectors(
bool printToStream, std::ostream &outStream)
const {
548 Ptr<Vector<Real>> x, y;
550 x = INPUT_xprim_->clone(); x->randomize(-one,one);
551 y = INPUT_xprim_->clone(); y->randomize(-one,one);
553 outStream << std::endl <<
" Check primal optimization space vector" << std::endl;
555 INPUT_xprim_->checkVector(*x,*y,printToStream,outStream);
558 x = INPUT_xdual_->clone(); x->randomize(-one,one);
559 y = INPUT_xdual_->clone(); y->randomize(-one,one);
561 outStream << std::endl <<
" Check dual optimization space vector" << std::endl;
563 INPUT_xdual_->checkVector(*x,*y,printToStream,outStream);
566 for (
auto it = INPUT_con_.begin(); it != INPUT_con_.end(); ++it) {
568 x = it->second.residual->clone(); x->randomize(-one,one);
569 y = it->second.residual->clone(); y->randomize(-one,one);
571 outStream << std::endl <<
" " << it->first <<
": Check primal constraint space vector" << std::endl;
573 it->second.residual->checkVector(*x,*y,printToStream,outStream);
576 x = it->second.multiplier->clone(); x->randomize(-one,one);
577 y = it->second.multiplier->clone(); y->randomize(-one,one);
579 outStream << std::endl <<
" " << it->first <<
": Check dual constraint space vector" << std::endl;
581 it->second.multiplier->checkVector(*x,*y,printToStream,outStream);
585 for (
auto it = INPUT_linear_con_.begin(); it != INPUT_linear_con_.end(); ++it) {
587 x = it->second.residual->clone(); x->randomize(-one,one);
588 y = it->second.residual->clone(); y->randomize(-one,one);
590 outStream << std::endl <<
" " << it->first <<
": Check primal linear constraint space vector" << std::endl;
592 it->second.residual->checkVector(*x,*y,printToStream,outStream);
595 x = it->second.multiplier->clone(); x->randomize(-one,one);
596 y = it->second.multiplier->clone(); y->randomize(-one,one);
598 outStream << std::endl <<
" " << it->first <<
": Check dual linear constraint space vector" << std::endl;
600 it->second.multiplier->checkVector(*x,*y,printToStream,outStream);
604template<
typename Real>
607 Ptr<Vector<Real>> x, d, v, g, c, w;
610 if (x == nullPtr) { x = INPUT_xprim_->clone(); x->randomize(-one,one); }
611 d = INPUT_xprim_->clone(); d->randomize(-scale,scale);
612 v = INPUT_xprim_->clone(); v->randomize(-scale,scale);
613 g = INPUT_xdual_->clone(); g->randomize(-scale,scale);
615 outStream << std::endl <<
" Check objective function" << std::endl << std::endl;
616 INPUT_obj_->checkGradient(*x,*g,*d,printToStream,outStream);
617 INPUT_obj_->checkHessVec(*x,*g,*d,printToStream,outStream);
618 INPUT_obj_->checkHessSym(*x,*g,*d,*v,printToStream,outStream);
621 for (
auto it = INPUT_con_.begin(); it != INPUT_con_.end(); ++it) {
622 c = it->second.residual->clone(); c->randomize(-scale,scale);
623 w = it->second.multiplier->clone(); w->randomize(-scale,scale);
625 outStream << std::endl <<
" " << it->first <<
": Check constraint function" << std::endl << std::endl;
626 it->second.constraint->checkApplyJacobian(*x,*v,*c,printToStream,outStream);
627 it->second.constraint->checkAdjointConsistencyJacobian(*w,*v,*x,printToStream,outStream);
628 it->second.constraint->checkApplyAdjointHessian(*x,*w,*v,*g,printToStream,outStream);
632 for (
auto it = INPUT_linear_con_.begin(); it != INPUT_linear_con_.end(); ++it) {
633 c = it->second.residual->clone(); c->randomize(-scale,scale);
634 w = it->second.multiplier->clone(); w->randomize(-scale,scale);
636 outStream << std::endl <<
" " << it->first <<
": Check constraint function" << std::endl << std::endl;
637 it->second.constraint->checkApplyJacobian(*x,*v,*c,printToStream,outStream);
638 it->second.constraint->checkAdjointConsistencyJacobian(*w,*v,*x,printToStream,outStream);
639 it->second.constraint->checkApplyAdjointHessian(*x,*w,*v,*g,printToStream,outStream);
643template<
typename Real>
645 checkVectors(printToStream,outStream);
646 if (hasLinearEquality_ || hasLinearInequality_)
647 checkLinearity(printToStream,outStream);
648 checkDerivatives(printToStream,outStream,x0,scale);
651template<
typename Real>
656template<
typename Real>
657void Problem<Real>::edit() {
658 isFinalized_ =
false;
663template<
typename Real>
664void Problem<Real>::finalizeIteration() {
665 if (rlc_ != nullPtr) {
666 if (!hasInequality_) {
667 rlc_->project(*INPUT_xprim_,*xprim_);
668 INPUT_xprim_->plus(*rlc_->getFeasibleVector());
671 Ptr<Vector<Real>> xprim =
dynamic_cast<PartitionedVector<Real>&
>(*xprim_).get(0)->clone();
672 xprim->set(*
dynamic_cast<PartitionedVector<Real>&
>(*xprim_).get(0));
673 rlc_->project(*INPUT_xprim_,*xprim);
674 INPUT_xprim_->plus(*rlc_->getFeasibleVector());
Problem(const Ptr< Objective< Real > > &obj, const Ptr< Vector< Real > > &x, const Ptr< Vector< Real > > &g=nullPtr)
Default constructor for OptimizationProblem.
Defines the linear algebra or vector space interface.