ROL
Loading...
Searching...
No Matches
algorithm/ROL_Problem_Def.hpp
Go to the documentation of this file.
1// @HEADER
2// *****************************************************************************
3// Rapid Optimization Library (ROL) Package
4//
5// Copyright 2014 NTESS and the ROL contributors.
6// SPDX-License-Identifier: BSD-3-Clause
7// *****************************************************************************
8// @HEADER
9
10#ifndef ROL_PROBLEM_DEF_HPP
11#define ROL_PROBLEM_DEF_HPP
12
13#include <iostream>
14
15#include "ROL_AugmentedLagrangianObjective2.hpp"
16#include "ROL_AugmentedLagrangianPenalty.hpp"
17#include "ROL_IdentityOperator.hpp"
18#include "ROL_PartitionedVector.hpp"
19#include "ROL_PolyhedralProjection.hpp"
20#include "ROL_Projection_Partitioned.hpp"
21#include "ROL_ZeroProjection.hpp"
22
23namespace ROL {
24
25template<typename Real>
26Problem<Real>::Problem( const Ptr<Objective<Real>> &obj,
27 const Ptr<Vector<Real>> &x,
28 const Ptr<Vector<Real>> &g)
29 : isFinalized_(false), hasBounds_(false),
30 hasEquality_(false), hasInequality_(false),
31 hasLinearEquality_(false), hasLinearInequality_(false),
32 hasProximableObjective_(false),
33 cnt_econ_(0), cnt_icon_(0), cnt_linear_econ_(0), cnt_linear_icon_(0),
34 obj_(nullPtr), nobj_(nullPtr), xprim_(nullPtr), xdual_(nullPtr), bnd_(nullPtr),
35 con_(nullPtr), mul_(nullPtr), res_(nullPtr), proj_(nullPtr),
36 problemType_(TYPE_U) {
37 INPUT_obj_ = obj;
38 INPUT_nobj_ = nullPtr;
39 INPUT_xprim_ = x;
40 INPUT_bnd_ = nullPtr;
41 INPUT_con_.clear();
42 INPUT_linear_con_.clear();
43 if (g==nullPtr) INPUT_xdual_ = x->dual().clone();
44 else INPUT_xdual_ = g;
45}
46
47template<typename Real>
48void Problem<Real>::addBoundConstraint(const Ptr<BoundConstraint<Real>> &bnd) {
49 ROL_TEST_FOR_EXCEPTION(isFinalized_,std::invalid_argument,
50 ">>> ROL::Problem: Cannot add bounds after problem is finalized!");
51
52 INPUT_bnd_ = bnd;
53 hasBounds_ = true;
54}
55
56template<typename Real>
57void Problem<Real>::removeBoundConstraint() {
58 ROL_TEST_FOR_EXCEPTION(isFinalized_,std::invalid_argument,
59 ">>> ROL::Problem: Cannot remove bounds after problem is finalized!");
60
61 INPUT_bnd_ = nullPtr;
62 hasBounds_ = false;
63}
64
65template<typename Real>
66void Problem<Real>::addConstraint( std::string name,
67 const Ptr<Constraint<Real>> &econ,
68 const Ptr<Vector<Real>> &emul,
69 const Ptr<Vector<Real>> &eres,
70 bool reset) {
71 ROL_TEST_FOR_EXCEPTION(isFinalized_,std::invalid_argument,
72 ">>> ROL::Problem: Cannot add constraint after problem is finalized!");
73
74 if (reset) INPUT_con_.clear();
75
76 bool isNameUsed = false;
77 for (const auto& map : {INPUT_con_,INPUT_linear_con_,INPUT_proj_}) {
78 if (map.count(name)) {
79 isNameUsed = true;
80 break;
81 }
82 }
83 ROL_TEST_FOR_EXCEPTION(isNameUsed,std::invalid_argument,
84 ">>> ROL::Problem: Constraint names must be distinct!");
85
86 INPUT_con_.insert({name,ConstraintData<Real>(econ,emul,eres)});
87 hasEquality_ = true;
88 cnt_econ_++;
89}
90
91template<typename Real>
92void Problem<Real>::addConstraint( std::string name,
93 const Ptr<Constraint<Real>> &icon,
94 const Ptr<Vector<Real>> &imul,
95 const Ptr<BoundConstraint<Real>> &ibnd,
96 const Ptr<Vector<Real>> &ires,
97 bool reset) {
98 ROL_TEST_FOR_EXCEPTION(isFinalized_,std::invalid_argument,
99 ">>> ROL::Problem: Cannot add constraint after problem is finalized!");
100
101 if (reset) INPUT_con_.clear();
102
103 bool isNameUsed = false;
104 for (const auto& map : {INPUT_con_,INPUT_linear_con_,INPUT_proj_}) {
105 if (map.count(name)) {
106 isNameUsed = true;
107 break;
108 }
109 }
110 ROL_TEST_FOR_EXCEPTION(isNameUsed,std::invalid_argument,
111 ">>> ROL::Problem: Constraint names must be distinct!");
112
113 INPUT_con_.insert({name,ConstraintData<Real>(icon,imul,ires,ibnd)});
114 hasInequality_ = true;
115 cnt_icon_++;
116}
117
118template<typename Real>
119void Problem<Real>::addConstraint( std::string name,
120 const Ptr<Constraint<Real>> &pcon,
121 const Ptr<Vector<Real>> &pmul,
122 const Ptr<Projection<Real>> &proj,
123 const Ptr<Vector<Real>> &pres,
124 bool reset) {
125 ROL_TEST_FOR_EXCEPTION(isFinalized_,std::invalid_argument,
126 ">>> ROL::Problem: Cannot add constraint after problem is finalized!");
127
128 if (reset) INPUT_proj_.clear();
129
130 bool isNameUsed = false;
131 for (const auto& map : {INPUT_con_,INPUT_linear_con_,INPUT_proj_}) {
132 if (map.count(name)) {
133 isNameUsed = true;
134 break;
135 }
136 }
137 ROL_TEST_FOR_EXCEPTION(isNameUsed,std::invalid_argument,
138 ">>> ROL::Problem: Constraint names must be distinct!");
139
140 if (proj == nullPtr) {
141 INPUT_con_.insert({name,ConstraintData<Real>(pcon,pmul,pres)});
142 hasEquality_ = true;
143 cnt_econ_++;
144 }
145 else
146 INPUT_proj_.insert({name,ConstraintData<Real>(pcon,pmul,pres,nullPtr,proj)});
147}
148
149template<typename Real>
150void Problem<Real>::removeConstraint(std::string name) {
151 ROL_TEST_FOR_EXCEPTION(isFinalized_,std::invalid_argument,
152 ">>> ROL::Problem: Cannot remove constraint after problem is finalized!");
153
154 auto it = INPUT_con_.find(name);
155 if (it != INPUT_con_.end()) {
156 if (it->second.bounds==nullPtr) cnt_econ_--;
157 else cnt_icon_--;
158 INPUT_con_.erase(it);
159 }
160 if (cnt_econ_==0) hasEquality_ = false;
161 if (cnt_icon_==0) hasInequality_ = false;
162 it = INPUT_proj_.find(name);
163 if (it != INPUT_proj_.end()) INPUT_proj_.erase(it);
164}
165
166template<typename Real>
167void Problem<Real>::addLinearConstraint( std::string name,
168 const Ptr<Constraint<Real>> &linear_econ,
169 const Ptr<Vector<Real>> &linear_emul,
170 const Ptr<Vector<Real>> &linear_eres,
171 bool reset) {
172 ROL_TEST_FOR_EXCEPTION(isFinalized_,std::invalid_argument,
173 ">>> ROL::Problem: Cannot add linear constraint after problem is finalized!");
174
175 if (reset) INPUT_linear_con_.clear();
176
177 bool isNameUsed = false;
178 for (const auto& map : {INPUT_con_,INPUT_linear_con_,INPUT_proj_}) {
179 if (map.count(name)) {
180 isNameUsed = true;
181 break;
182 }
183 }
184 ROL_TEST_FOR_EXCEPTION(isNameUsed,std::invalid_argument,
185 ">>> ROL::Problem: Constraint names must be distinct!");
186
187 INPUT_linear_con_.insert({name,ConstraintData<Real>(linear_econ,linear_emul,linear_eres)});
188 hasLinearEquality_ = true;
189 cnt_linear_econ_++;
190}
191
192template<typename Real>
193void Problem<Real>::addLinearConstraint( std::string name,
194 const Ptr<Constraint<Real>> &linear_icon,
195 const Ptr<Vector<Real>> &linear_imul,
196 const Ptr<BoundConstraint<Real>> &linear_ibnd,
197 const Ptr<Vector<Real>> &linear_ires,
198 bool reset) {
199 ROL_TEST_FOR_EXCEPTION(isFinalized_,std::invalid_argument,
200 ">>> ROL::Problem: Cannot add linear constraint after problem is finalized!");
201
202 if (reset) INPUT_linear_con_.clear();
203
204 bool isNameUsed = false;
205 for (const auto& map : {INPUT_con_,INPUT_linear_con_,INPUT_proj_}) {
206 if (map.count(name)) {
207 isNameUsed = true;
208 break;
209 }
210 }
211 ROL_TEST_FOR_EXCEPTION(isNameUsed,std::invalid_argument,
212 ">>> ROL::Problem: Constraint names must be distinct!");
213
214 INPUT_linear_con_.insert({name,ConstraintData<Real>(linear_icon,linear_imul,linear_ires,linear_ibnd)});
215 hasLinearInequality_ = true;
216 cnt_linear_icon_++;
217}
218
219template<typename Real>
220void Problem<Real>::removeLinearConstraint(std::string name) {
221 ROL_TEST_FOR_EXCEPTION(isFinalized_,std::invalid_argument,
222 ">>> ROL::Problem: Cannot remove linear inequality after problem is finalized!");
223
224 auto it = INPUT_linear_con_.find(name);
225 if (it!=INPUT_linear_con_.end()) {
226 if (it->second.bounds==nullPtr) cnt_linear_econ_--;
227 else cnt_linear_icon_--;
228 INPUT_linear_con_.erase(it);
229 }
230 if (cnt_linear_econ_==0) hasLinearEquality_ = false;
231 if (cnt_linear_icon_==0) hasLinearInequality_ = false;
232}
233
234template<typename Real>
235void Problem<Real>::setProjectionAlgorithm(ParameterList &list) {
236 ROL_TEST_FOR_EXCEPTION(isFinalized_,std::invalid_argument,
237 ">>> ROL::Problem: Cannot set polyhedral projection algorithm after problem is finalized!");
238
239 ppa_list_ = list;
240}
241
242template<typename Real>
243void Problem<Real>::addProximableObjective(const Ptr<Objective<Real>> &nobj) {
244 ROL_TEST_FOR_EXCEPTION(isFinalized_,std::invalid_argument,
245 ">>> ROL::Problem: Cannot add regularizer after problem is finalized!");
246
247 INPUT_nobj_ = nobj;
248 hasProximableObjective_ = true;
249}
250
251template<typename Real>
252void Problem<Real>::removeProximableObjective() {
253 ROL_TEST_FOR_EXCEPTION(isFinalized_,std::invalid_argument,
254 ">>> ROL::Problem: Cannot remove regularizer after problem is finalized!");
255
256 INPUT_nobj_ = nullPtr;
257 hasProximableObjective_ = false;
258}
259
260template<typename Real>
261void Problem<Real>::addConstraintGroup(const std::string &group_name,
262 const std::vector<std::string> &constraint_names) {
263 ROL_TEST_FOR_EXCEPTION(isFinalized_,std::invalid_argument,
264 ">>> ROL::Problem: Cannot add constraint group after problem is finalized!");
265
266 bool isNameUsed;
267 for (const auto& name : constraint_names) {
268 ROL_TEST_FOR_EXCEPTION(grouped_constraint_names_.count(name),std::invalid_argument,
269 ">>> ROL::Problem: Cannot include the same constraint in two augmented Lagrangian groups!");
270 isNameUsed = false;
271 for (const auto& map : {INPUT_con_,INPUT_linear_con_,INPUT_proj_}) {
272 if (map.count(name)) {
273 isNameUsed = true;
274 break;
275 }
276 }
277 ROL_TEST_FOR_EXCEPTION(!isNameUsed,std::invalid_argument,
278 ">>> ROL::Problem: Constraint names must be distinct!");
279 }
280 for (const auto& name : constraint_names) {
281 grouped_constraint_names_.insert(name);
282 groups_[group_name].push_back(name);
283 }
284}
285
286template<typename Real>
287void Problem<Real>::removeConstraintGroup(const std::string &group_name) {
288 ROL_TEST_FOR_EXCEPTION(isFinalized_,std::invalid_argument,
289 ">>> ROL::Problem: Cannot remove augmented Lagrangian group after problem is finalized!");
290 auto it = groups_.find(group_name);
291 ROL_TEST_FOR_EXCEPTION(it==groups_.end(),std::invalid_argument,
292 ">>> ROL::Problem::removeConstraintGroup: Group name does not exist!");
293 const std::vector<std::string>& constraint_names = it->second;
294 for (const auto& name : constraint_names)
295 grouped_constraint_names_.erase(name);
296 groups_.erase(it);
297}
298
299
300template<typename Real>
301void Problem<Real>::finalize(bool lumpConstraints, bool printToStream, std::ostream &outStream, bool useSlackVariables) {
302 if (!isFinalized_) {
303 // ========================================================================
304 // useSlackVariables = true -> the usual finalize
305 // ========================================================================
306 if (useSlackVariables) {
307 std::unordered_map<std::string,ConstraintData<Real>> con, lcon, icon;
308 bool hasEquality = hasEquality_;
309 bool hasLinearEquality = hasLinearEquality_;
310 bool hasInequality = hasInequality_;
311 bool hasLinearInequality = hasLinearInequality_;
312 bool hasProximableObjective = hasProximableObjective_;
313 con.insert(INPUT_con_.begin(),INPUT_con_.end());
314 if (lumpConstraints) {
315 con.insert(INPUT_linear_con_.begin(),INPUT_linear_con_.end());
316 hasEquality = (hasEquality || hasLinearEquality);
317 hasInequality = (hasInequality || hasLinearInequality);
318 hasLinearEquality = false;
319 hasLinearInequality = false;
320 }
321 else {
322 lcon.insert(INPUT_linear_con_.begin(),INPUT_linear_con_.end());
323 }
324 // Transform optimization problem
325 //std::cout << hasBounds_ << " " << hasEquality << " " << hasInequality << " " << hasLinearEquality << " " << hasLinearInequality << std::endl;
326 nobj_ = nullPtr;
327 if (hasProximableObjective){
328 if (!hasEquality && !hasInequality && !hasBounds_ && !hasLinearEquality && !hasLinearInequality){
329 problemType_ = TYPE_P;
330 obj_ = INPUT_obj_;
331 nobj_ = INPUT_nobj_;
332 xprim_ = INPUT_xprim_;
333 xdual_ = INPUT_xdual_;
334 bnd_ = nullPtr;
335 con_ = nullPtr;
336 mul_ = nullPtr;
337 res_ = nullPtr;
338 }
339 else {
340 throw Exception::NotImplemented(">>> ROL::TypeP - with constraints is not supported");
341 }
342 }
343 else {
344 if (!hasLinearEquality && !hasLinearInequality) {
345 proj_ = nullPtr;
346 if (!hasEquality && !hasInequality && !hasBounds_ ) {
347 problemType_ = TYPE_U;
348 obj_ = INPUT_obj_;
349 xprim_ = INPUT_xprim_;
350 xdual_ = INPUT_xdual_;
351 bnd_ = nullPtr;
352 con_ = nullPtr;
353 mul_ = nullPtr;
354 res_ = nullPtr;
355 }
356 else if (!hasEquality && !hasInequality && hasBounds_) {
357 problemType_ = TYPE_B;
358 obj_ = INPUT_obj_;
359 xprim_ = INPUT_xprim_;
360 xdual_ = INPUT_xdual_;
361 bnd_ = INPUT_bnd_;
362 con_ = nullPtr;
363 mul_ = nullPtr;
364 res_ = nullPtr;
365 }
366 else if (hasEquality && !hasInequality && !hasBounds_) {
367 ConstraintAssembler<Real> cm(con,INPUT_xprim_,INPUT_xdual_);
368 problemType_ = TYPE_E;
369 obj_ = INPUT_obj_;
370 xprim_ = INPUT_xprim_;
371 xdual_ = INPUT_xdual_;
372 bnd_ = nullPtr;
373 con_ = cm.getConstraint();
374 mul_ = cm.getMultiplier();
375 res_ = cm.getResidual();
376 }
377 else {
378 ConstraintAssembler<Real> cm(con,INPUT_xprim_,INPUT_xdual_,INPUT_bnd_);
379 problemType_ = TYPE_EB;
380 obj_ = INPUT_obj_;
381 if (cm.hasInequality()) {
382 obj_ = makePtr<SlacklessObjective<Real>>(INPUT_obj_);
383 }
384 xprim_ = cm.getOptVector();
385 xdual_ = cm.getDualOptVector();
386 bnd_ = cm.getBoundConstraint();
387 con_ = cm.getConstraint();
388 mul_ = cm.getMultiplier();
389 res_ = cm.getResidual();
390 }
391 }
392 else {
393 if (!hasBounds_ && !hasLinearInequality) {
394 ConstraintAssembler<Real> cm(lcon,INPUT_xprim_,INPUT_xdual_);
395 xfeas_ = cm.getOptVector()->clone(); xfeas_->set(*cm.getOptVector());
396 rlc_ = makePtr<ReduceLinearConstraint<Real>>(cm.getConstraint(),xfeas_,cm.getResidual());
397 proj_ = nullPtr;
398 if (!hasEquality && !hasInequality) {
399 problemType_ = TYPE_U;
400 obj_ = rlc_->transform(INPUT_obj_);
401 xprim_ = xfeas_->clone(); xprim_->zero();
402 xdual_ = cm.getDualOptVector();
403 bnd_ = nullPtr;
404 con_ = nullPtr;
405 mul_ = nullPtr;
406 res_ = nullPtr;
407 }
408 else {
409 for (auto it = con.begin(); it != con.end(); ++it) {
410 icon.insert(std::pair<std::string,ConstraintData<Real>>(it->first,
411 ConstraintData<Real>(rlc_->transform(it->second.constraint),
412 it->second.multiplier,it->second.residual,it->second.bounds)));
413 }
414 Ptr<Vector<Real>> xtmp = xfeas_->clone(); xtmp->zero();
415 ConstraintAssembler<Real> cm1(icon,xtmp,cm.getDualOptVector());
416 xprim_ = cm1.getOptVector();
417 xdual_ = cm1.getDualOptVector();
418 con_ = cm1.getConstraint();
419 mul_ = cm1.getMultiplier();
420 res_ = cm1.getResidual();
421 if (!hasInequality) {
422 problemType_ = TYPE_E;
423 obj_ = rlc_->transform(INPUT_obj_);
424 bnd_ = nullPtr;
425 }
426 else {
427 problemType_ = TYPE_EB;
428 obj_ = makePtr<SlacklessObjective<Real>>(rlc_->transform(INPUT_obj_));
429 bnd_ = cm1.getBoundConstraint();
430 }
431 }
432 }
433 else if ((hasBounds_ || hasLinearInequality) && !hasEquality && !hasInequality) {
434 ConstraintAssembler<Real> cm(lcon,INPUT_xprim_,INPUT_xdual_,INPUT_bnd_);
435 problemType_ = TYPE_B;
436 obj_ = INPUT_obj_;
437 if (cm.hasInequality()) {
438 obj_ = makePtr<SlacklessObjective<Real>>(INPUT_obj_);
439 }
440 xprim_ = cm.getOptVector();
441 xdual_ = cm.getDualOptVector();
442 bnd_ = cm.getBoundConstraint();
443 con_ = nullPtr;
444 mul_ = nullPtr;
445 res_ = nullPtr;
446 proj_ = PolyhedralProjectionFactory<Real>(*xprim_,*xdual_,bnd_,
447 cm.getConstraint(),*cm.getMultiplier(),*cm.getResidual(),ppa_list_);
448 }
449 else {
450 ConstraintAssembler<Real> cm(con,lcon,INPUT_xprim_,INPUT_xdual_,INPUT_bnd_);
451 problemType_ = TYPE_EB;
452 obj_ = INPUT_obj_;
453 if (cm.hasInequality()) {
454 obj_ = makePtr<SlacklessObjective<Real>>(INPUT_obj_);
455 }
456 xprim_ = cm.getOptVector();
457 xdual_ = cm.getDualOptVector();
458 con_ = cm.getConstraint();
459 mul_ = cm.getMultiplier();
460 res_ = cm.getResidual();
461 bnd_ = cm.getBoundConstraint();
462 proj_ = PolyhedralProjectionFactory<Real>(*xprim_,*xdual_,bnd_,
463 cm.getLinearConstraint(),*cm.getLinearMultiplier(),
464 *cm.getLinearResidual(),ppa_list_);
465 }
466 }
467 }
468 isFinalized_ = true;
469 if (printToStream) {
470 outStream << std::endl;
471 outStream << " ROL::Problem::finalize" << std::endl;
472 outStream << " Problem Summary:" << std::endl;
473 outStream << " Has Proximable Objective? .......... " << (hasProximableObjective ? "yes" : "no") << std::endl;
474 outStream << " Has Bound Constraint? .............. " << (hasBounds_ ? "yes" : "no") << std::endl;
475 outStream << " Has Equality Constraint? ........... " << (hasEquality ? "yes" : "no") << std::endl;
476 if (hasEquality) {
477 int cnt = 0;
478 for (auto it = con.begin(); it != con.end(); ++it) {
479 if (it->second.bounds==nullPtr) {
480 if (cnt==0) {
481 outStream << " Names: ........................... ";
482 cnt++;
483 }
484 else {
485 outStream << " ";
486 }
487 outStream << it->first << std::endl;
488 }
489 }
490 outStream << " Total: ........................... " << cnt_econ_+(lumpConstraints ? cnt_linear_econ_ : 0) << std::endl;
491 }
492 outStream << " Has Inequality Constraint? ......... " << (hasInequality ? "yes" : "no") << std::endl;
493 if (hasInequality) {
494 int cnt = 0;
495 for (auto it = con.begin(); it != con.end(); ++it) {
496 if (it->second.bounds!=nullPtr) {
497 if (cnt==0) {
498 outStream << " Names: ........................... ";
499 cnt++;
500 }
501 else {
502 outStream << " ";
503 }
504 outStream << it->first << std::endl;
505 }
506 }
507 outStream << " Total: ........................... " << cnt_icon_+(lumpConstraints ? cnt_linear_icon_ : 0) << std::endl;
508 }
509 if (!lumpConstraints) {
510 outStream << " Has Linear Equality Constraint? .... " << (hasLinearEquality ? "yes" : "no") << std::endl;
511 if (hasLinearEquality) {
512 int cnt = 0;
513 for (auto it = lcon.begin(); it != lcon.end(); ++it) {
514 if (it->second.bounds==nullPtr) {
515 if (cnt==0) {
516 outStream << " Names: ........................... ";
517 cnt++;
518 }
519 else {
520 outStream << " ";
521 }
522 outStream << it->first << std::endl;
523 }
524 }
525 outStream << " Total: ........................... " << cnt_linear_econ_ << std::endl;
526 }
527 outStream << " Has Linear Inequality Constraint? .. " << (hasLinearInequality ? "yes" : "no") << std::endl;
528 if (hasLinearInequality) {
529 int cnt = 0;
530 for (auto it = lcon.begin(); it != lcon.end(); ++it) {
531 if (it->second.bounds!=nullPtr) {
532 if (cnt==0) {
533 outStream << " Names: ........................... ";
534 cnt++;
535 }
536 else {
537 outStream << " ";
538 }
539 outStream << it->first << std::endl;
540 }
541 }
542 outStream << " Total: ........................... " << cnt_linear_icon_ << std::endl;
543 }
544 }
545 outStream << std::endl;
546 }
547 }
548 // ========================================================================
549 // useSlackVariables = false -> process constraints into groups
550 // ========================================================================
551 else{
552 // Ungrouped constraints
553 std::unordered_set<std::string> ungrouped_constraints;
554 for (const auto& kv : INPUT_con_) ungrouped_constraints.insert(kv.first);
555 for (const auto& kv : INPUT_linear_con_) ungrouped_constraints.insert(kv.first);
556 for (const auto& kv : INPUT_proj_) ungrouped_constraints.insert(kv.first);
557 for (const auto& kv : groups_) {
558 for (const auto& constraint_name : kv.second)
559 ungrouped_constraints.erase(constraint_name);
560 }
561 // Convert ungrouped inequalities to their own groups.
562 Ptr<ConstraintData<Real>> constraint_data;
563 Ptr<Projection<Real>> projection;
564 bool is_equality = false;
565 for (const std::string& constraint_name : ungrouped_constraints) {
566 constraint_data = getConstraintData(constraint_name);
567 is_equality = (dynamicPtrCast<ZeroProjection<Real>>(constraint_data->projection) != nullPtr);
568 if (!is_equality) {
569 std::vector<std::string> constraint_name_vector(1,constraint_name);
570 addConstraintGroup(constraint_name, constraint_name_vector);
571 }
572 else if (INPUT_linear_con_.count(constraint_name))
573 ungrouped_linear_equality_constraint_names_.push_back(constraint_name);
574 else
575 ungrouped_equality_constraint_names_.push_back(constraint_name);
576 }
577 problemType_ = TYPE_EB;
578 obj_ = INPUT_obj_;
579 nobj_ = INPUT_nobj_;
580 xprim_ = INPUT_xprim_;
581 xdual_ = INPUT_xdual_;
582 bnd_ = INPUT_bnd_;
583 con_ = nullPtr;
584 mul_ = nullPtr;
585 res_ = nullPtr;
586 isFinalized_ = true;
587 }
588 }
589 else {
590 if (printToStream) {
591 outStream << std::endl;
592 outStream << " ROL::Problem::finalize" << std::endl;
593 outStream << " Problem already finalized!" << std::endl;
594 outStream << std::endl;
595 }
596 }
597}
598
599template<typename Real>
600const Ptr<Objective<Real>>& Problem<Real>::getObjective() {
601 finalize();
602 return obj_;
603}
604
605template<typename Real>
606const Ptr<Objective<Real>>& Problem<Real>::getProximableObjective(){
607 finalize();
608 return nobj_;
609}
610
611template<typename Real>
612const Ptr<Vector<Real>>& Problem<Real>::getPrimalOptimizationVector() {
613 finalize();
614 return xprim_;
615}
616
617template<typename Real>
618const Ptr<Vector<Real>>& Problem<Real>::getDualOptimizationVector() {
619 finalize();
620 return xdual_;
621}
622
623template<typename Real>
624const Ptr<BoundConstraint<Real>>& Problem<Real>::getBoundConstraint() {
625 finalize();
626 return bnd_;
627}
628
629template<typename Real>
630const Ptr<Constraint<Real>>& Problem<Real>::getConstraint() {
631 finalize();
632 return con_;
633}
634
635template<typename Real>
636const Ptr<Vector<Real>>& Problem<Real>::getMultiplierVector() {
637 finalize();
638 return mul_;
639}
640
641template<typename Real>
642const Ptr<Vector<Real>>& Problem<Real>::getResidualVector() {
643 finalize();
644 return res_;
645}
646
647template<typename Real>
648const Ptr<PolyhedralProjection<Real>>& Problem<Real>::getPolyhedralProjection() {
649 finalize();
650 return proj_;
651}
652
653template<typename Real>
654EProblem Problem<Real>::getProblemType() {
655 finalize();
656 return problemType_;
657}
658
659template<typename Real>
660Ptr<ConstraintData<Real>> Problem<Real>::getConstraintData(const std::string &constraint_name) {
661 Ptr<ConstraintData<Real>> constraint_data;
662 bool isConstraintMissing = true;
663 for (const auto& map : {INPUT_con_,INPUT_linear_con_,INPUT_proj_}) {
664 if (const auto& it = map.find(constraint_name); it != map.end()) {
665 constraint_data = makePtr<ConstraintData<Real>>(it->second);
666 if (constraint_data->projection == nullPtr) {
667 if (constraint_data->bounds != nullPtr)
668 constraint_data->projection = makePtr<PolyhedralProjection<Real>>(constraint_data->bounds);
669 else
670 constraint_data->projection = makePtr<ZeroProjection<Real>>();
671 }
672 isConstraintMissing = false;
673 break;
674 }
675 }
676 ROL_TEST_FOR_EXCEPTION(isConstraintMissing,std::invalid_argument,
677 ">>> ROL::Problem::getConstraintData: Invalid constraint name!");
678 return constraint_data;
679}
680
681template<typename Real>
682const std::vector<std::string>& Problem<Real>::getUngroupedEqualityConstraintNames() {
683 return ungrouped_equality_constraint_names_;
684}
685
686template<typename Real>
687const std::vector<std::string>& Problem<Real>::getUngroupedLinearEqualityConstraintNames() {
688 return ungrouped_linear_equality_constraint_names_;
689}
690
691template<typename Real>
692const std::vector<std::string> Problem<Real>::getGroupNames() {
693 std::vector<std::string> group_names;
694 for (const auto& kv : groups_) group_names.push_back(kv.first);
695 return group_names; // return a vector to preserve group order
696}
697
698template<typename Real>
699const Ptr<ConstraintData<Real>> Problem<Real>::getGroupConstraintData(const std::string &group_name) {
700 ROL_TEST_FOR_EXCEPTION(groups_.count(group_name)==0,std::invalid_argument,
701 ">>> ROL::Problem::getGroupConstraintData A constraint group with the provided name is missing!");
702 if (groups_[group_name].size() == 1) {
703 std::string& constraint_name = groups_[group_name][0];
704 return getConstraintData(constraint_name);
705 }
706 Ptr<ConstraintData<Real>> constraint_data;
707 std::vector<Ptr<Constraint<Real>>> constraints;
708 std::vector<Ptr<Projection<Real>>> projections;
709 std::vector<Ptr<Vector<Real>>> con_vectors;
710 std::vector<Ptr<Vector<Real>>> mul_vectors;
711 for (const std::string& constraint_name : groups_[group_name]) {
712 constraint_data = getConstraintData(constraint_name);
713 constraints.push_back(constraint_data->constraint);
714 mul_vectors.push_back(constraint_data->multiplier);
715 con_vectors.push_back(constraint_data->residual);
716 projections.push_back(constraint_data->projection);
717 }
718 Ptr<Constraint<Real>> partitioned_constraint = makePtr<Constraint_Partitioned<Real>>(constraints);
719 Ptr<Vector<Real>> partitioned_mul_vector = makePtr<PartitionedVector<Real>>(mul_vectors);
720 Ptr<Vector<Real>> partitioned_con_vector = makePtr<PartitionedVector<Real>>(con_vectors);
721 Ptr<Projection<Real>> partitioned_projection = makePtr<Projection_Partitioned<Real>>(projections);
722 Ptr<ConstraintData<Real>> new_constraint_data = makePtr<ConstraintData<Real>>(partitioned_constraint,
723 partitioned_mul_vector,
724 partitioned_con_vector,
725 nullPtr,
726 partitioned_projection);
727 return new_constraint_data;
728}
729
730
731template<typename Real>
732Real Problem<Real>::checkLinearity(bool printToStream, std::ostream &outStream) const {
733 std::ios_base::fmtflags state(outStream.flags());
734 if (printToStream) {
735 outStream << std::setprecision(8) << std::scientific;
736 outStream << std::endl;
737 outStream << " ROL::Problem::checkLinearity" << std::endl;
738 }
739 const Real one(1), two(2), eps(1e-2*std::sqrt(ROL_EPSILON<Real>()));
740 Real tol(std::sqrt(ROL_EPSILON<Real>())), cnorm(0), err(0), maxerr(0);
741 Ptr<Vector<Real>> x = INPUT_xprim_->clone(); x->randomize(-one,one);
742 Ptr<Vector<Real>> y = INPUT_xprim_->clone(); y->randomize(-one,one);
743 Ptr<Vector<Real>> z = INPUT_xprim_->clone(); z->zero();
744 Ptr<Vector<Real>> xy = INPUT_xprim_->clone();
745 Real alpha = two*static_cast<Real>(rand())/static_cast<Real>(RAND_MAX)-one;
746 xy->set(*x); xy->axpy(alpha,*y);
747 Ptr<Vector<Real>> c1, c2;
748 for (auto it = INPUT_linear_con_.begin(); it != INPUT_linear_con_.end(); ++it) {
749 c1 = it->second.residual->clone();
750 c2 = it->second.residual->clone();
751 it->second.constraint->update(*xy,UpdateType::Temp);
752 it->second.constraint->value(*c1,*xy,tol);
753 cnorm = c1->norm();
754 it->second.constraint->update(*x,UpdateType::Temp);
755 it->second.constraint->value(*c2,*x,tol);
756 c1->axpy(-one,*c2);
757 it->second.constraint->update(*y,UpdateType::Temp);
758 it->second.constraint->value(*c2,*y,tol);
759 c1->axpy(-alpha,*c2);
760 it->second.constraint->update(*z,UpdateType::Temp);
761 it->second.constraint->value(*c2,*z,tol);
762 c1->axpy(alpha,*c2);
763 err = c1->norm();
764 maxerr = std::max(err,maxerr);
765 if (printToStream) {
766 outStream << " Constraint " << it->first;
767 outStream << ": ||c(x+alpha*y) - (c(x)+alpha*(c(y)-c(0)))|| = " << err << std::endl;
768 if (err > eps*cnorm) {
769 outStream << " Constraint " << it->first << " may not be linear!" << std::endl;
770 }
771 }
772 }
773 if (printToStream) {
774 outStream << std::endl;
775 }
776 outStream.flags(state);
777 return maxerr;
778}
779
780template<typename Real>
781void Problem<Real>::checkVectors(bool printToStream, std::ostream &outStream) const {
782 const Real one(1);
783 Ptr<Vector<Real>> x, y;
784 // Primal optimization space vector
785 x = INPUT_xprim_->clone(); x->randomize(-one,one);
786 y = INPUT_xprim_->clone(); y->randomize(-one,one);
787 if (printToStream) {
788 outStream << std::endl << " Check primal optimization space vector" << std::endl;
789 }
790 INPUT_xprim_->checkVector(*x,*y,printToStream,outStream);
791
792 // Dual optimization space vector
793 x = INPUT_xdual_->clone(); x->randomize(-one,one);
794 y = INPUT_xdual_->clone(); y->randomize(-one,one);
795 if (printToStream) {
796 outStream << std::endl << " Check dual optimization space vector" << std::endl;
797 }
798 INPUT_xdual_->checkVector(*x,*y,printToStream,outStream);
799
800 // Check constraint space vectors
801 for (auto it = INPUT_con_.begin(); it != INPUT_con_.end(); ++it) {
802 // Primal constraint space vector
803 x = it->second.residual->clone(); x->randomize(-one,one);
804 y = it->second.residual->clone(); y->randomize(-one,one);
805 if (printToStream) {
806 outStream << std::endl << " " << it->first << ": Check primal constraint space vector" << std::endl;
807 }
808 it->second.residual->checkVector(*x,*y,printToStream,outStream);
809
810 // Dual optimization space vector
811 x = it->second.multiplier->clone(); x->randomize(-one,one);
812 y = it->second.multiplier->clone(); y->randomize(-one,one);
813 if (printToStream) {
814 outStream << std::endl << " " << it->first << ": Check dual constraint space vector" << std::endl;
815 }
816 it->second.multiplier->checkVector(*x,*y,printToStream,outStream);
817 }
818
819 // Check constraint space vectors
820 for (auto it = INPUT_linear_con_.begin(); it != INPUT_linear_con_.end(); ++it) {
821 // Primal constraint space vector
822 x = it->second.residual->clone(); x->randomize(-one,one);
823 y = it->second.residual->clone(); y->randomize(-one,one);
824 if (printToStream) {
825 outStream << std::endl << " " << it->first << ": Check primal linear constraint space vector" << std::endl;
826 }
827 it->second.residual->checkVector(*x,*y,printToStream,outStream);
828
829 // Dual optimization space vector
830 x = it->second.multiplier->clone(); x->randomize(-one,one);
831 y = it->second.multiplier->clone(); y->randomize(-one,one);
832 if (printToStream) {
833 outStream << std::endl << " " << it->first << ": Check dual linear constraint space vector" << std::endl;
834 }
835 it->second.multiplier->checkVector(*x,*y,printToStream,outStream);
836 }
837}
838
839template<typename Real>
840void Problem<Real>::checkDerivatives(bool printToStream, std::ostream &outStream, const Ptr<Vector<Real>> &x0, Real scale) const {
841 const Real one(1);
842 Ptr<Vector<Real>> x, d, v, g, c, w;
843 // Objective check
844 x = x0;
845 if (x == nullPtr) { x = INPUT_xprim_->clone(); x->randomize(-one,one); }
846 d = INPUT_xprim_->clone(); d->randomize(-scale,scale);
847 v = INPUT_xprim_->clone(); v->randomize(-scale,scale);
848 g = INPUT_xdual_->clone(); g->randomize(-scale,scale);
849 if (printToStream)
850 outStream << std::endl << " Check objective function" << std::endl << std::endl;
851 INPUT_obj_->checkGradient(*x,*g,*d,printToStream,outStream);
852 INPUT_obj_->checkHessVec(*x,*g,*d,printToStream,outStream);
853 INPUT_obj_->checkHessSym(*x,*g,*d,*v,printToStream,outStream);
854 //TODO: Proximable Objective Check
855 // Constraint check
856 for (auto it = INPUT_con_.begin(); it != INPUT_con_.end(); ++it) {
857 c = it->second.residual->clone(); c->randomize(-scale,scale);
858 w = it->second.multiplier->clone(); w->randomize(-scale,scale);
859 if (printToStream)
860 outStream << std::endl << " " << it->first << ": Check constraint function" << std::endl << std::endl;
861 it->second.constraint->checkApplyJacobian(*x,*v,*c,printToStream,outStream);
862 it->second.constraint->checkAdjointConsistencyJacobian(*w,*v,*x,printToStream,outStream);
863 it->second.constraint->checkApplyAdjointHessian(*x,*w,*v,*g,printToStream,outStream);
864 }
865
866 // Linear constraint check
867 for (auto it = INPUT_linear_con_.begin(); it != INPUT_linear_con_.end(); ++it) {
868 c = it->second.residual->clone(); c->randomize(-scale,scale);
869 w = it->second.multiplier->clone(); w->randomize(-scale,scale);
870 if (printToStream)
871 outStream << std::endl << " " << it->first << ": Check constraint function" << std::endl << std::endl;
872 it->second.constraint->checkApplyJacobian(*x,*v,*c,printToStream,outStream);
873 it->second.constraint->checkAdjointConsistencyJacobian(*w,*v,*x,printToStream,outStream);
874 it->second.constraint->checkApplyAdjointHessian(*x,*w,*v,*g,printToStream,outStream);
875 }
876}
877
878template<typename Real>
879void Problem<Real>::check(bool printToStream, std::ostream &outStream, const Ptr<Vector<Real>> &x0, Real scale) const {
880 checkVectors(printToStream,outStream);
881 if (hasLinearEquality_ || hasLinearInequality_)
882 checkLinearity(printToStream,outStream);
883 checkDerivatives(printToStream,outStream,x0,scale);
884// if (hasProximableObjective)
885// checkProximableObjective(printToStream, outStream);
886}
887
888template<typename Real>
889bool Problem<Real>::isFinalized() const {
890 return isFinalized_;
891}
892
893template<typename Real>
894void Problem<Real>::edit() {
895 isFinalized_ = false;
896 rlc_ = nullPtr;
897 proj_ = nullPtr;
898}
899
900template<typename Real>
901void Problem<Real>::finalizeIteration() {
902 if (rlc_ != nullPtr) {
903 if (!hasInequality_) {
904 rlc_->project(*INPUT_xprim_,*xprim_);
905 INPUT_xprim_->plus(*rlc_->getFeasibleVector());
906 }
907 else {
908 Ptr<Vector<Real>> xprim = dynamic_cast<PartitionedVector<Real>&>(*xprim_).get(0)->clone();
909 xprim->set(*dynamic_cast<PartitionedVector<Real>&>(*xprim_).get(0));
910 rlc_->project(*INPUT_xprim_,*xprim);
911 INPUT_xprim_->plus(*rlc_->getFeasibleVector());
912 }
913 }
914}
915
916} // namespace ROL
917
918#endif // ROL_NEWOPTIMIZATIONPROBLEM_DEF_HPP
ROL_AugmentedLagrangianObjective2.hpp
ROL_AugmentedLagrangianPenalty.hpp
ROL_IdentityOperator.hpp
obj_
const Ptr< Obj > obj_
Definition ROL_Objective_SerialSimOpt.hpp:60
ROL_PartitionedVector.hpp
ROL_PolyhedralProjection.hpp
ROL_Projection_Partitioned.hpp
ROL_ZeroProjection.hpp
ROL::ConstraintAssembler
Provides a wrapper for multiple constraints.
Definition ROL_ConstraintAssembler.hpp:45
ROL::ConstraintAssembler::getLinearConstraint
const Ptr< Constraint< Real > > & getLinearConstraint() const
Definition ROL_ConstraintAssembler_Def.hpp:283
ROL::ConstraintAssembler::getMultiplier
const Ptr< Vector< Real > > & getMultiplier() const
Definition ROL_ConstraintAssembler_Def.hpp:273
ROL::ConstraintAssembler::getLinearResidual
const Ptr< Vector< Real > > & getLinearResidual() const
Definition ROL_ConstraintAssembler_Def.hpp:293
ROL::ConstraintAssembler::getDualOptVector
const Ptr< Vector< Real > > & getDualOptVector() const
Definition ROL_ConstraintAssembler_Def.hpp:303
ROL::ConstraintAssembler::getBoundConstraint
const Ptr< BoundConstraint< Real > > & getBoundConstraint() const
Definition ROL_ConstraintAssembler_Def.hpp:308
ROL::ConstraintAssembler::getConstraint
const Ptr< Constraint< Real > > & getConstraint() const
Definition ROL_ConstraintAssembler_Def.hpp:268
ROL::ConstraintAssembler::getLinearMultiplier
const Ptr< Vector< Real > > & getLinearMultiplier() const
Definition ROL_ConstraintAssembler_Def.hpp:288
ROL::ConstraintAssembler::getOptVector
const Ptr< Vector< Real > > & getOptVector() const
Definition ROL_ConstraintAssembler_Def.hpp:298
ROL::ConstraintAssembler::getResidual
const Ptr< Vector< Real > > & getResidual() const
Definition ROL_ConstraintAssembler_Def.hpp:278
ROL::Constraint
Defines the general constraint operator interface.
Definition ROL_Constraint.hpp:52
ROL::Objective
Provides the interface to evaluate objective functions.
Definition ROL_Objective.hpp:44
ROL::Problem::Problem
Problem(const Ptr< Objective< Real > > &obj, const Ptr< Vector< Real > > &x, const Ptr< Vector< Real > > &g=nullPtr)
Default constructor for OptimizationProblem.
Definition function/prox/ROL_Problem_Def.hpp:18
ROL::Vector
Defines the linear algebra or vector space interface.
Definition ROL_Vector.hpp:51
ROL::TYPE_U
@ TYPE_U
Definition ROL_Types.hpp:224
ROL::TYPE_E
@ TYPE_E
Definition ROL_Types.hpp:227
ROL::TYPE_EB
@ TYPE_EB
Definition ROL_Types.hpp:228
ROL::TYPE_P
@ TYPE_P
Definition ROL_Types.hpp:225
ROL::TYPE_B
@ TYPE_B
Definition ROL_Types.hpp:226
Generated by doxygen 1.9.8