44 #ifndef ROL_PROBLEM_DEF_HPP
45 #define ROL_PROBLEM_DEF_HPP
51 template<
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),
71 template<
typename Real>
73 ROL_TEST_FOR_EXCEPTION(isFinalized_,std::invalid_argument,
74 ">>> ROL::Problem: Cannot add bounds after problem is finalized!");
80 template<
typename Real>
82 ROL_TEST_FOR_EXCEPTION(isFinalized_,std::invalid_argument,
83 ">>> ROL::Problem: Cannot remove bounds after problem is finalized!");
89 template<
typename Real>
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)});
109 template<
typename Real>
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;
130 template<
typename Real>
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;
145 template<
typename Real>
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;
165 template<
typename Real>
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;
186 template<
typename Real>
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;
201 template<
typename Real>
203 ROL_TEST_FOR_EXCEPTION(isFinalized_,std::invalid_argument,
204 ">>> ROL::Problem: Cannot set polyhedral projection algorithm after problem is finalized!");
209 template<
typename Real>
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;
442 template<
typename Real>
448 template<
typename Real>
454 template<
typename Real>
460 template<
typename Real>
466 template<
typename Real>
472 template<
typename Real>
478 template<
typename Real>
484 template<
typename Real>
490 template<
typename Real>
496 template<
typename Real>
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();
517 it->second.constraint->value(*c1,*xy,tol);
520 it->second.constraint->value(*c2,*x,tol);
523 it->second.constraint->value(*c2,*y,tol);
524 c1->axpy(-alpha,*c2);
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);
545 template<
typename Real>
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);
604 template<
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);
643 template<
typename Real>
645 checkVectors(printToStream,outStream);
646 if (hasLinearEquality_ || hasLinearInequality_)
647 checkLinearity(printToStream,outStream);
648 checkDerivatives(printToStream,outStream,x0,scale);
651 template<
typename Real>
656 template<
typename Real>
658 isFinalized_ =
false;
663 template<
typename Real>
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());
681 #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.