10 #ifndef ROL_PROBLEM_DEF_HPP
11 #define ROL_PROBLEM_DEF_HPP
17 template<
typename Real>
19 const Ptr<Vector<Real>> &x,
20 const Ptr<Vector<Real>> &g)
21 : isFinalized_(false), hasBounds_(false),
22 hasEquality_(false), hasInequality_(false),
23 hasLinearEquality_(false), hasLinearInequality_(false),
24 cnt_econ_(0), cnt_icon_(0), cnt_linear_econ_(0), cnt_linear_icon_(0),
25 obj_(nullPtr), xprim_(nullPtr), xdual_(nullPtr), bnd_(nullPtr),
26 con_(nullPtr), mul_(nullPtr), res_(nullPtr), proj_(nullPtr),
37 template<
typename Real>
39 ROL_TEST_FOR_EXCEPTION(isFinalized_,std::invalid_argument,
40 ">>> ROL::Problem: Cannot add bounds after problem is finalized!");
46 template<
typename Real>
48 ROL_TEST_FOR_EXCEPTION(isFinalized_,std::invalid_argument,
49 ">>> ROL::Problem: Cannot remove bounds after problem is finalized!");
55 template<
typename Real>
57 const Ptr<Constraint<Real>> &econ,
58 const Ptr<Vector<Real>> &emul,
59 const Ptr<Vector<Real>> &eres,
61 ROL_TEST_FOR_EXCEPTION(isFinalized_,std::invalid_argument,
62 ">>> ROL::Problem: Cannot add constraint after problem is finalized!");
64 if (reset) INPUT_con_.clear();
66 auto it = INPUT_con_.find(name);
67 ROL_TEST_FOR_EXCEPTION(it != INPUT_con_.end(),std::invalid_argument,
68 ">>> ROL::Problem: Constraint names must be distinct!");
70 INPUT_con_.insert({name,ConstraintData<Real>(econ,emul,eres)});
75 template<
typename Real>
77 const Ptr<Constraint<Real>> &icon,
78 const Ptr<Vector<Real>> &imul,
79 const Ptr<BoundConstraint<Real>> &ibnd,
80 const Ptr<Vector<Real>> &ires,
82 ROL_TEST_FOR_EXCEPTION(isFinalized_,std::invalid_argument,
83 ">>> ROL::Problem: Cannot add constraint after problem is finalized!");
85 if (reset) INPUT_con_.clear();
87 auto it = INPUT_con_.find(name);
88 ROL_TEST_FOR_EXCEPTION(it != INPUT_con_.end(),std::invalid_argument,
89 ">>> ROL::Problem: Constraint names must be distinct!");
91 INPUT_con_.insert({name,ConstraintData<Real>(icon,imul,ires,ibnd)});
92 hasInequality_ =
true;
96 template<
typename Real>
98 ROL_TEST_FOR_EXCEPTION(isFinalized_,std::invalid_argument,
99 ">>> ROL::Problem: Cannot remove constraint after problem is finalized!");
101 auto it = INPUT_con_.find(name);
102 if (it!=INPUT_con_.end()) {
103 if (it->second.bounds==nullPtr) cnt_econ_--;
105 INPUT_con_.erase(it);
107 if (cnt_econ_==0) hasEquality_ =
false;
108 if (cnt_icon_==0) hasInequality_ =
false;
111 template<
typename Real>
113 const Ptr<Constraint<Real>> &linear_econ,
114 const Ptr<Vector<Real>> &linear_emul,
115 const Ptr<Vector<Real>> &linear_eres,
117 ROL_TEST_FOR_EXCEPTION(isFinalized_,std::invalid_argument,
118 ">>> ROL::Problem: Cannot add linear constraint after problem is finalized!");
120 if (reset) INPUT_linear_con_.clear();
122 auto it = INPUT_linear_con_.find(name);
123 ROL_TEST_FOR_EXCEPTION(it != INPUT_linear_con_.end(),std::invalid_argument,
124 ">>> ROL::Problem: Linear constraint names must be distinct!");
126 INPUT_linear_con_.insert({name,ConstraintData<Real>(linear_econ,linear_emul,linear_eres)});
127 hasLinearEquality_ =
true;
131 template<
typename Real>
133 const Ptr<Constraint<Real>> &linear_icon,
134 const Ptr<Vector<Real>> &linear_imul,
135 const Ptr<BoundConstraint<Real>> &linear_ibnd,
136 const Ptr<Vector<Real>> &linear_ires,
138 ROL_TEST_FOR_EXCEPTION(isFinalized_,std::invalid_argument,
139 ">>> ROL::Problem: Cannot add linear constraint after problem is finalized!");
141 if (reset) INPUT_linear_con_.clear();
143 auto it = INPUT_linear_con_.find(name);
144 ROL_TEST_FOR_EXCEPTION(it != INPUT_linear_con_.end(),std::invalid_argument,
145 ">>> ROL::Problem: Linear constraint names must be distinct!");
147 INPUT_linear_con_.insert({name,ConstraintData<Real>(linear_icon,linear_imul,linear_ires,linear_ibnd)});
148 hasLinearInequality_ =
true;
152 template<
typename Real>
154 ROL_TEST_FOR_EXCEPTION(isFinalized_,std::invalid_argument,
155 ">>> ROL::Problem: Cannot remove linear inequality after problem is finalized!");
157 auto it = INPUT_linear_con_.find(name);
158 if (it!=INPUT_linear_con_.end()) {
159 if (it->second.bounds==nullPtr) cnt_linear_econ_--;
160 else cnt_linear_icon_--;
161 INPUT_linear_con_.erase(it);
163 if (cnt_linear_econ_==0) hasLinearEquality_ =
false;
164 if (cnt_linear_icon_==0) hasLinearInequality_ =
false;
167 template<
typename Real>
169 ROL_TEST_FOR_EXCEPTION(isFinalized_,std::invalid_argument,
170 ">>> ROL::Problem: Cannot set polyhedral projection algorithm after problem is finalized!");
175 template<
typename Real>
178 std::unordered_map<std::string,ConstraintData<Real>> con, lcon, icon;
179 bool hasEquality = hasEquality_;
180 bool hasLinearEquality = hasLinearEquality_;
181 bool hasInequality = hasInequality_;
182 bool hasLinearInequality = hasLinearInequality_;
183 con.insert(INPUT_con_.begin(),INPUT_con_.end());
184 if (lumpConstraints) {
185 con.insert(INPUT_linear_con_.begin(),INPUT_linear_con_.end());
186 hasEquality = (hasEquality || hasLinearEquality);
187 hasInequality = (hasInequality || hasLinearInequality);
188 hasLinearEquality =
false;
189 hasLinearInequality =
false;
192 lcon.insert(INPUT_linear_con_.begin(),INPUT_linear_con_.end());
196 if (!hasLinearEquality && !hasLinearInequality) {
198 if (!hasEquality && !hasInequality && !hasBounds_) {
201 xprim_ = INPUT_xprim_;
202 xdual_ = INPUT_xdual_;
208 else if (!hasEquality && !hasInequality && hasBounds_) {
211 xprim_ = INPUT_xprim_;
212 xdual_ = INPUT_xdual_;
218 else if (hasEquality && !hasInequality && !hasBounds_) {
219 ConstraintAssembler<Real> cm(con,INPUT_xprim_,INPUT_xdual_);
222 xprim_ = INPUT_xprim_;
223 xdual_ = INPUT_xdual_;
225 con_ = cm.getConstraint();
226 mul_ = cm.getMultiplier();
227 res_ = cm.getResidual();
230 ConstraintAssembler<Real> cm(con,INPUT_xprim_,INPUT_xdual_,INPUT_bnd_);
233 if (cm.hasInequality()) {
234 obj_ = makePtr<SlacklessObjective<Real>>(INPUT_obj_);
236 xprim_ = cm.getOptVector();
237 xdual_ = cm.getDualOptVector();
238 bnd_ = cm.getBoundConstraint();
239 con_ = cm.getConstraint();
240 mul_ = cm.getMultiplier();
241 res_ = cm.getResidual();
245 if (!hasBounds_ && !hasLinearInequality) {
246 ConstraintAssembler<Real> cm(lcon,INPUT_xprim_,INPUT_xdual_);
247 xfeas_ = cm.getOptVector()->clone(); xfeas_->set(*cm.getOptVector());
248 rlc_ = makePtr<ReduceLinearConstraint<Real>>(cm.getConstraint(),xfeas_,cm.getResidual());
250 if (!hasEquality && !hasInequality) {
252 obj_ = rlc_->transform(INPUT_obj_);
253 xprim_ = xfeas_->clone(); xprim_->zero();
254 xdual_ = cm.getDualOptVector();
261 for (
auto it = con.begin(); it != con.end(); ++it) {
262 icon.insert(std::pair<std::string,ConstraintData<Real>>(it->first,
263 ConstraintData<Real>(rlc_->transform(it->second.constraint),
264 it->second.multiplier,it->second.residual,it->second.bounds)));
266 Ptr<Vector<Real>> xtmp = xfeas_->clone(); xtmp->zero();
267 ConstraintAssembler<Real> cm1(icon,xtmp,cm.getDualOptVector());
268 xprim_ = cm1.getOptVector();
269 xdual_ = cm1.getDualOptVector();
270 con_ = cm1.getConstraint();
271 mul_ = cm1.getMultiplier();
272 res_ = cm1.getResidual();
273 if (!hasInequality) {
275 obj_ = rlc_->transform(INPUT_obj_);
280 obj_ = makePtr<SlacklessObjective<Real>>(rlc_->transform(INPUT_obj_));
281 bnd_ = cm1.getBoundConstraint();
285 else if ((hasBounds_ || hasLinearInequality) && !hasEquality && !hasInequality) {
286 ConstraintAssembler<Real> cm(lcon,INPUT_xprim_,INPUT_xdual_,INPUT_bnd_);
289 if (cm.hasInequality()) {
290 obj_ = makePtr<SlacklessObjective<Real>>(INPUT_obj_);
292 xprim_ = cm.getOptVector();
293 xdual_ = cm.getDualOptVector();
294 bnd_ = cm.getBoundConstraint();
298 proj_ = PolyhedralProjectionFactory<Real>(*xprim_,*xdual_,bnd_,
299 cm.getConstraint(),*cm.getMultiplier(),*cm.getResidual(),ppa_list_);
302 ConstraintAssembler<Real> cm(con,lcon,INPUT_xprim_,INPUT_xdual_,INPUT_bnd_);
305 if (cm.hasInequality()) {
306 obj_ = makePtr<SlacklessObjective<Real>>(INPUT_obj_);
308 xprim_ = cm.getOptVector();
309 xdual_ = cm.getDualOptVector();
310 con_ = cm.getConstraint();
311 mul_ = cm.getMultiplier();
312 res_ = cm.getResidual();
313 bnd_ = cm.getBoundConstraint();
314 proj_ = PolyhedralProjectionFactory<Real>(*xprim_,*xdual_,bnd_,
315 cm.getLinearConstraint(),*cm.getLinearMultiplier(),
316 *cm.getLinearResidual(),ppa_list_);
321 outStream << std::endl;
322 outStream <<
" ROL::Problem::finalize" << std::endl;
323 outStream <<
" Problem Summary:" << std::endl;
324 outStream <<
" Has Bound Constraint? .............. " << (hasBounds_ ?
"yes" :
"no") << std::endl;
325 outStream <<
" Has Equality Constraint? ........... " << (hasEquality ?
"yes" :
"no") << std::endl;
328 for (
auto it = con.begin(); it != con.end(); ++it) {
329 if (it->second.bounds==nullPtr) {
331 outStream <<
" Names: ........................... ";
337 outStream << it->first << std::endl;
340 outStream <<
" Total: ........................... " << cnt_econ_+(lumpConstraints ? cnt_linear_econ_ : 0) << std::endl;
342 outStream <<
" Has Inequality Constraint? ......... " << (hasInequality ?
"yes" :
"no") << std::endl;
345 for (
auto it = con.begin(); it != con.end(); ++it) {
346 if (it->second.bounds!=nullPtr) {
348 outStream <<
" Names: ........................... ";
354 outStream << it->first << std::endl;
357 outStream <<
" Total: ........................... " << cnt_icon_+(lumpConstraints ? cnt_linear_icon_ : 0) << std::endl;
359 if (!lumpConstraints) {
360 outStream <<
" Has Linear Equality Constraint? .... " << (hasLinearEquality ?
"yes" :
"no") << std::endl;
361 if (hasLinearEquality) {
363 for (
auto it = lcon.begin(); it != lcon.end(); ++it) {
364 if (it->second.bounds==nullPtr) {
366 outStream <<
" Names: ........................... ";
372 outStream << it->first << std::endl;
375 outStream <<
" Total: ........................... " << cnt_linear_econ_ << std::endl;
377 outStream <<
" Has Linear Inequality Constraint? .. " << (hasLinearInequality ?
"yes" :
"no") << std::endl;
378 if (hasLinearInequality) {
380 for (
auto it = lcon.begin(); it != lcon.end(); ++it) {
381 if (it->second.bounds!=nullPtr) {
383 outStream <<
" Names: ........................... ";
389 outStream << it->first << std::endl;
392 outStream <<
" Total: ........................... " << cnt_linear_icon_ << std::endl;
395 outStream << std::endl;
400 outStream << std::endl;
401 outStream <<
" ROL::Problem::finalize" << std::endl;
402 outStream <<
" Problem already finalized!" << std::endl;
403 outStream << std::endl;
408 template<
typename Real>
414 template<
typename Real>
420 template<
typename Real>
426 template<
typename Real>
432 template<
typename Real>
438 template<
typename Real>
444 template<
typename Real>
450 template<
typename Real>
456 template<
typename Real>
462 template<
typename Real>
464 std::ios_base::fmtflags state(outStream.flags());
466 outStream << std::setprecision(8) << std::scientific;
467 outStream << std::endl;
468 outStream <<
" ROL::Problem::checkLinearity" << std::endl;
470 const Real one(1), two(2), eps(1e-2*std::sqrt(ROL_EPSILON<Real>()));
471 Real tol(std::sqrt(ROL_EPSILON<Real>())), cnorm(0), err(0), maxerr(0);
472 Ptr<Vector<Real>> x = INPUT_xprim_->clone(); x->randomize(-one,one);
473 Ptr<Vector<Real>> y = INPUT_xprim_->clone(); y->randomize(-one,one);
474 Ptr<Vector<Real>> z = INPUT_xprim_->clone(); z->zero();
475 Ptr<Vector<Real>> xy = INPUT_xprim_->clone();
476 Real alpha = two*
static_cast<Real
>(rand())/static_cast<Real>(RAND_MAX)-one;
477 xy->set(*x); xy->axpy(alpha,*y);
478 Ptr<Vector<Real>> c1, c2;
479 for (
auto it = INPUT_linear_con_.begin(); it != INPUT_linear_con_.end(); ++it) {
480 c1 = it->second.residual->clone();
481 c2 = it->second.residual->clone();
483 it->second.constraint->value(*c1,*xy,tol);
486 it->second.constraint->value(*c2,*x,tol);
489 it->second.constraint->value(*c2,*y,tol);
490 c1->axpy(-alpha,*c2);
492 it->second.constraint->value(*c2,*z,tol);
495 maxerr = std::max(err,maxerr);
497 outStream <<
" Constraint " << it->first;
498 outStream <<
": ||c(x+alpha*y) - (c(x)+alpha*(c(y)-c(0)))|| = " << err << std::endl;
499 if (err > eps*cnorm) {
500 outStream <<
" Constraint " << it->first <<
" may not be linear!" << std::endl;
505 outStream << std::endl;
507 outStream.flags(state);
511 template<
typename Real>
514 Ptr<Vector<Real>> x, y;
516 x = INPUT_xprim_->clone(); x->randomize(-one,one);
517 y = INPUT_xprim_->clone(); y->randomize(-one,one);
519 outStream << std::endl <<
" Check primal optimization space vector" << std::endl;
521 INPUT_xprim_->checkVector(*x,*y,printToStream,outStream);
524 x = INPUT_xdual_->clone(); x->randomize(-one,one);
525 y = INPUT_xdual_->clone(); y->randomize(-one,one);
527 outStream << std::endl <<
" Check dual optimization space vector" << std::endl;
529 INPUT_xdual_->checkVector(*x,*y,printToStream,outStream);
532 for (
auto it = INPUT_con_.begin(); it != INPUT_con_.end(); ++it) {
534 x = it->second.residual->clone(); x->randomize(-one,one);
535 y = it->second.residual->clone(); y->randomize(-one,one);
537 outStream << std::endl <<
" " << it->first <<
": Check primal constraint space vector" << std::endl;
539 it->second.residual->checkVector(*x,*y,printToStream,outStream);
542 x = it->second.multiplier->clone(); x->randomize(-one,one);
543 y = it->second.multiplier->clone(); y->randomize(-one,one);
545 outStream << std::endl <<
" " << it->first <<
": Check dual constraint space vector" << std::endl;
547 it->second.multiplier->checkVector(*x,*y,printToStream,outStream);
551 for (
auto it = INPUT_linear_con_.begin(); it != INPUT_linear_con_.end(); ++it) {
553 x = it->second.residual->clone(); x->randomize(-one,one);
554 y = it->second.residual->clone(); y->randomize(-one,one);
556 outStream << std::endl <<
" " << it->first <<
": Check primal linear constraint space vector" << std::endl;
558 it->second.residual->checkVector(*x,*y,printToStream,outStream);
561 x = it->second.multiplier->clone(); x->randomize(-one,one);
562 y = it->second.multiplier->clone(); y->randomize(-one,one);
564 outStream << std::endl <<
" " << it->first <<
": Check dual linear constraint space vector" << std::endl;
566 it->second.multiplier->checkVector(*x,*y,printToStream,outStream);
570 template<
typename Real>
573 Ptr<Vector<Real>> x, d, v, g, c, w;
576 if (x == nullPtr) { x = INPUT_xprim_->clone(); x->randomize(-one,one); }
577 d = INPUT_xprim_->clone(); d->randomize(-scale,scale);
578 v = INPUT_xprim_->clone(); v->randomize(-scale,scale);
579 g = INPUT_xdual_->clone(); g->randomize(-scale,scale);
581 outStream << std::endl <<
" Check objective function" << std::endl << std::endl;
582 INPUT_obj_->checkGradient(*x,*g,*d,printToStream,outStream);
583 INPUT_obj_->checkHessVec(*x,*g,*d,printToStream,outStream);
584 INPUT_obj_->checkHessSym(*x,*g,*d,*v,printToStream,outStream);
587 for (
auto it = INPUT_con_.begin(); it != INPUT_con_.end(); ++it) {
588 c = it->second.residual->clone(); c->randomize(-scale,scale);
589 w = it->second.multiplier->clone(); w->randomize(-scale,scale);
591 outStream << std::endl <<
" " << it->first <<
": Check constraint function" << std::endl << std::endl;
592 it->second.constraint->checkApplyJacobian(*x,*v,*c,printToStream,outStream);
593 it->second.constraint->checkAdjointConsistencyJacobian(*w,*v,*x,printToStream,outStream);
594 it->second.constraint->checkApplyAdjointHessian(*x,*w,*v,*g,printToStream,outStream);
598 for (
auto it = INPUT_linear_con_.begin(); it != INPUT_linear_con_.end(); ++it) {
599 c = it->second.residual->clone(); c->randomize(-scale,scale);
600 w = it->second.multiplier->clone(); w->randomize(-scale,scale);
602 outStream << std::endl <<
" " << it->first <<
": Check constraint function" << std::endl << std::endl;
603 it->second.constraint->checkApplyJacobian(*x,*v,*c,printToStream,outStream);
604 it->second.constraint->checkAdjointConsistencyJacobian(*w,*v,*x,printToStream,outStream);
605 it->second.constraint->checkApplyAdjointHessian(*x,*w,*v,*g,printToStream,outStream);
609 template<
typename Real>
611 checkVectors(printToStream,outStream);
612 if (hasLinearEquality_ || hasLinearInequality_)
613 checkLinearity(printToStream,outStream);
614 checkDerivatives(printToStream,outStream,x0,scale);
617 template<
typename Real>
622 template<
typename Real>
624 isFinalized_ =
false;
629 template<
typename Real>
631 if (rlc_ != nullPtr) {
632 if (!hasInequality_) {
633 rlc_->project(*INPUT_xprim_,*xprim_);
634 INPUT_xprim_->plus(*rlc_->getFeasibleVector());
637 Ptr<Vector<Real>> xprim =
dynamic_cast<PartitionedVector<Real>&
>(*xprim_).get(0)->clone();
638 xprim->set(*
dynamic_cast<PartitionedVector<Real>&
>(*xprim_).get(0));
639 rlc_->project(*INPUT_xprim_,*xprim);
640 INPUT_xprim_->plus(*rlc_->getFeasibleVector());
647 #endif // ROL_NEWOPTIMIZATIONPROBLEM_DEF_HPP
void check(bool printToStream=false, std::ostream &outStream=std::cout) const
Run vector, linearity and derivative checks for user-supplied vectors, objective function and constra...
const Ptr< Constraint< Real > > & getConstraint()
Get the equality constraint.
void checkDerivatives(bool printToStream=false, std::ostream &outStream=std::cout) const
Run derivative checks for user-supplied objective function and constraints.
const Ptr< Vector< Real > > & getPrimalOptimizationVector()
Get the primal optimization space vector.
void addBoundConstraint(const Ptr< BoundConstraint< Real >> &bnd)
Add a bound constraint.
void addLinearConstraint(std::string name, const Ptr< Constraint< Real >> &linear_econ, const Ptr< Vector< Real >> &linear_emul, const Ptr< Vector< Real >> &linear_eres=nullPtr, bool reset=false)
Add a linear equality constraint.
void removeBoundConstraint()
Remove an existing bound constraint.
void removeLinearConstraint(std::string name)
Remove an existing linear constraint.
const Ptr< BoundConstraint< Real > > & getBoundConstraint()
Get the bound constraint.
Defines the linear algebra or vector space interface.
virtual void finalize(bool lumpConstraints=false, bool printToStream=false, std::ostream &outStream=std::cout)
Tranform user-supplied constraints to consist of only bounds and equalities. Optimization problem can...
Ptr< Objective< Real > > INPUT_obj_
bool isFinalized() const
Indicate whether or no finalize has been called.
Ptr< Vector< Real > > INPUT_xprim_
const Ptr< Objective< Real > > & getObjective()
Get the objective function.
Problem(const Ptr< Objective< Real >> &obj, const Ptr< Vector< Real >> &x, const Ptr< Vector< Real >> &g=nullPtr)
Default constructor for OptimizationProblem.
std::unordered_map< std::string, ConstraintData< Real > > INPUT_linear_con_
std::unordered_map< std::string, ConstraintData< Real > > INPUT_con_
Ptr< Vector< Real > > INPUT_xdual_
const Ptr< Vector< Real > > & getResidualVector()
Get the primal constraint space vector.
const Ptr< Vector< Real > > & getMultiplierVector()
Get the dual constraint space vector.
void finalizeIteration()
Transform the optimization variables to the native parameterization after an optimization algorithm h...
Real checkLinearity(bool printToStream=false, std::ostream &outStream=std::cout) const
Check if user-supplied linear constraints are affine.
void addConstraint(std::string name, const Ptr< Constraint< Real >> &econ, const Ptr< Vector< Real >> &emul, const Ptr< Vector< Real >> &eres=nullPtr, bool reset=false)
Add an equality constraint.
const Ptr< PolyhedralProjection< Real > > & getPolyhedralProjection()
Get the polyhedral projection object. This is a null pointer if no linear constraints and/or bounds a...
EProblem getProblemType()
Get the optimization problem type (U, B, E, or G).
void removeConstraint(std::string name)
Remove an existing constraint.
Ptr< BoundConstraint< Real > > INPUT_bnd_
virtual void edit()
Resume editting optimization problem after finalize has been called.
const Ptr< Vector< Real > > & getDualOptimizationVector()
Get the dual optimization space vector.
void checkVectors(bool printToStream=false, std::ostream &outStream=std::cout) const
Run vector checks for user-supplied vectors.
void setProjectionAlgorithm(ParameterList &parlist)
Set polyhedral projection algorithm.