ROL
algorithm/ROL_Problem_Def.hpp
Go to the documentation of this file.
1// @HEADER
2// ************************************************************************
3//
4// Rapid Optimization Library (ROL) Package
5// Copyright (2014) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact lead developers:
38// Drew Kouri (dpkouri@sandia.gov) and
39// Denis Ridzal (dridzal@sandia.gov)
40//
41// ************************************************************************
42// @HEADER
43
44#ifndef ROL_PROBLEM_DEF_HPP
45#define ROL_PROBLEM_DEF_HPP
46
47#include <iostream>
48
49namespace ROL {
50
51template<typename Real>
52Problem<Real>::Problem( const Ptr<Objective<Real>> &obj,
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),
61 problemType_(TYPE_U) {
62 INPUT_obj_ = obj;
63 INPUT_xprim_ = x;
64 INPUT_bnd_ = nullPtr;
65 INPUT_con_.clear();
66 INPUT_linear_con_.clear();
67 if (g==nullPtr) INPUT_xdual_ = x->dual().clone();
68 else INPUT_xdual_ = g;
69}
70
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!");
75
76 INPUT_bnd_ = bnd;
77 hasBounds_ = true;
78}
79
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!");
84
85 INPUT_bnd_ = nullPtr;
86 hasBounds_ = false;
87}
88
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,
94 bool reset) {
95 ROL_TEST_FOR_EXCEPTION(isFinalized_,std::invalid_argument,
96 ">>> ROL::Problem: Cannot add constraint after problem is finalized!");
97
98 if (reset) INPUT_con_.clear();
99
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!");
103
104 INPUT_con_.insert({name,ConstraintData<Real>(econ,emul,eres)});
105 hasEquality_ = true;
106 cnt_econ_++;
107}
108
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,
115 bool reset) {
116 ROL_TEST_FOR_EXCEPTION(isFinalized_,std::invalid_argument,
117 ">>> ROL::Problem: Cannot add constraint after problem is finalized!");
118
119 if (reset) INPUT_con_.clear();
120
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!");
124
125 INPUT_con_.insert({name,ConstraintData<Real>(icon,imul,ires,ibnd)});
126 hasInequality_ = true;
127 cnt_icon_++;
128}
129
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!");
134
135 auto it = INPUT_con_.find(name);
136 if (it!=INPUT_con_.end()) {
137 if (it->second.bounds==nullPtr) cnt_econ_--;
138 else cnt_icon_--;
139 INPUT_con_.erase(it);
140 }
141 if (cnt_econ_==0) hasEquality_ = false;
142 if (cnt_icon_==0) hasInequality_ = false;
143}
144
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,
150 bool reset) {
151 ROL_TEST_FOR_EXCEPTION(isFinalized_,std::invalid_argument,
152 ">>> ROL::Problem: Cannot add linear constraint after problem is finalized!");
153
154 if (reset) INPUT_linear_con_.clear();
155
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!");
159
160 INPUT_linear_con_.insert({name,ConstraintData<Real>(linear_econ,linear_emul,linear_eres)});
161 hasLinearEquality_ = true;
162 cnt_linear_econ_++;
163}
164
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,
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 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!");
180
181 INPUT_linear_con_.insert({name,ConstraintData<Real>(linear_icon,linear_imul,linear_ires,linear_ibnd)});
182 hasLinearInequality_ = true;
183 cnt_linear_icon_++;
184}
185
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!");
190
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);
196 }
197 if (cnt_linear_econ_==0) hasLinearEquality_ = false;
198 if (cnt_linear_icon_==0) hasLinearInequality_ = false;
199}
200
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!");
205
206 ppa_list_ = list;
207}
208
209template<typename Real>
210void Problem<Real>::finalize(bool lumpConstraints, bool printToStream, std::ostream &outStream) {
211 if (!isFinalized_) {
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;
224 }
225 else {
226 lcon.insert(INPUT_linear_con_.begin(),INPUT_linear_con_.end());
227 }
228 // Transform optimization problem
229 //std::cout << hasBounds_ << " " << hasEquality << " " << hasInequality << " " << hasLinearEquality << " " << hasLinearInequality << std::endl;
230 if (!hasLinearEquality && !hasLinearInequality) {
231 proj_ = nullPtr;
232 if (!hasEquality && !hasInequality && !hasBounds_) {
233 problemType_ = TYPE_U;
234 obj_ = INPUT_obj_;
235 xprim_ = INPUT_xprim_;
236 xdual_ = INPUT_xdual_;
237 bnd_ = nullPtr;
238 con_ = nullPtr;
239 mul_ = nullPtr;
240 res_ = nullPtr;
241 }
242 else if (!hasEquality && !hasInequality && hasBounds_) {
243 problemType_ = TYPE_B;
244 obj_ = INPUT_obj_;
245 xprim_ = INPUT_xprim_;
246 xdual_ = INPUT_xdual_;
247 bnd_ = INPUT_bnd_;
248 con_ = nullPtr;
249 mul_ = nullPtr;
250 res_ = nullPtr;
251 }
252 else if (hasEquality && !hasInequality && !hasBounds_) {
253 ConstraintAssembler<Real> cm(con,INPUT_xprim_,INPUT_xdual_);
254 problemType_ = TYPE_E;
255 obj_ = INPUT_obj_;
256 xprim_ = INPUT_xprim_;
257 xdual_ = INPUT_xdual_;
258 bnd_ = nullPtr;
259 con_ = cm.getConstraint();
260 mul_ = cm.getMultiplier();
261 res_ = cm.getResidual();
262 }
263 else {
264 ConstraintAssembler<Real> cm(con,INPUT_xprim_,INPUT_xdual_,INPUT_bnd_);
265 problemType_ = TYPE_EB;
266 obj_ = INPUT_obj_;
267 if (cm.hasInequality()) {
268 obj_ = makePtr<SlacklessObjective<Real>>(INPUT_obj_);
269 }
270 xprim_ = cm.getOptVector();
271 xdual_ = cm.getDualOptVector();
272 bnd_ = cm.getBoundConstraint();
273 con_ = cm.getConstraint();
274 mul_ = cm.getMultiplier();
275 res_ = cm.getResidual();
276 }
277 }
278 else {
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());
283 proj_ = nullPtr;
284 if (!hasEquality && !hasInequality) {
285 problemType_ = TYPE_U;
286 obj_ = rlc_->transform(INPUT_obj_);
287 xprim_ = xfeas_->clone(); xprim_->zero();
288 xdual_ = cm.getDualOptVector();
289 bnd_ = nullPtr;
290 con_ = nullPtr;
291 mul_ = nullPtr;
292 res_ = nullPtr;
293 }
294 else {
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)));
299 }
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) {
308 problemType_ = TYPE_E;
309 obj_ = rlc_->transform(INPUT_obj_);
310 bnd_ = nullPtr;
311 }
312 else {
313 problemType_ = TYPE_EB;
314 obj_ = makePtr<SlacklessObjective<Real>>(rlc_->transform(INPUT_obj_));
315 bnd_ = cm1.getBoundConstraint();
316 }
317 }
318 }
319 else if ((hasBounds_ || hasLinearInequality) && !hasEquality && !hasInequality) {
320 ConstraintAssembler<Real> cm(lcon,INPUT_xprim_,INPUT_xdual_,INPUT_bnd_);
321 problemType_ = TYPE_B;
322 obj_ = INPUT_obj_;
323 if (cm.hasInequality()) {
324 obj_ = makePtr<SlacklessObjective<Real>>(INPUT_obj_);
325 }
326 xprim_ = cm.getOptVector();
327 xdual_ = cm.getDualOptVector();
328 bnd_ = cm.getBoundConstraint();
329 con_ = nullPtr;
330 mul_ = nullPtr;
331 res_ = nullPtr;
332 proj_ = PolyhedralProjectionFactory<Real>(*xprim_,*xdual_,bnd_,
333 cm.getConstraint(),*cm.getMultiplier(),*cm.getResidual(),ppa_list_);
334 }
335 else {
336 ConstraintAssembler<Real> cm(con,lcon,INPUT_xprim_,INPUT_xdual_,INPUT_bnd_);
337 problemType_ = TYPE_EB;
338 obj_ = INPUT_obj_;
339 if (cm.hasInequality()) {
340 obj_ = makePtr<SlacklessObjective<Real>>(INPUT_obj_);
341 }
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_);
351 }
352 }
353 isFinalized_ = true;
354 if (printToStream) {
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;
360 if (hasEquality) {
361 int cnt = 0;
362 for (auto it = con.begin(); it != con.end(); ++it) {
363 if (it->second.bounds==nullPtr) {
364 if (cnt==0) {
365 outStream << " Names: ........................... ";
366 cnt++;
367 }
368 else {
369 outStream << " ";
370 }
371 outStream << it->first << std::endl;
372 }
373 }
374 outStream << " Total: ........................... " << cnt_econ_+(lumpConstraints ? cnt_linear_econ_ : 0) << std::endl;
375 }
376 outStream << " Has Inequality Constraint? ......... " << (hasInequality ? "yes" : "no") << std::endl;
377 if (hasInequality) {
378 int cnt = 0;
379 for (auto it = con.begin(); it != con.end(); ++it) {
380 if (it->second.bounds!=nullPtr) {
381 if (cnt==0) {
382 outStream << " Names: ........................... ";
383 cnt++;
384 }
385 else {
386 outStream << " ";
387 }
388 outStream << it->first << std::endl;
389 }
390 }
391 outStream << " Total: ........................... " << cnt_icon_+(lumpConstraints ? cnt_linear_icon_ : 0) << std::endl;
392 }
393 if (!lumpConstraints) {
394 outStream << " Has Linear Equality Constraint? .... " << (hasLinearEquality ? "yes" : "no") << std::endl;
395 if (hasLinearEquality) {
396 int cnt = 0;
397 for (auto it = lcon.begin(); it != lcon.end(); ++it) {
398 if (it->second.bounds==nullPtr) {
399 if (cnt==0) {
400 outStream << " Names: ........................... ";
401 cnt++;
402 }
403 else {
404 outStream << " ";
405 }
406 outStream << it->first << std::endl;
407 }
408 }
409 outStream << " Total: ........................... " << cnt_linear_econ_ << std::endl;
410 }
411 outStream << " Has Linear Inequality Constraint? .. " << (hasLinearInequality ? "yes" : "no") << std::endl;
412 if (hasLinearInequality) {
413 int cnt = 0;
414 for (auto it = lcon.begin(); it != lcon.end(); ++it) {
415 if (it->second.bounds!=nullPtr) {
416 if (cnt==0) {
417 outStream << " Names: ........................... ";
418 cnt++;
419 }
420 else {
421 outStream << " ";
422 }
423 outStream << it->first << std::endl;
424 }
425 }
426 outStream << " Total: ........................... " << cnt_linear_icon_ << std::endl;
427 }
428 }
429 outStream << std::endl;
430 }
431 }
432 else {
433 if (printToStream) {
434 outStream << std::endl;
435 outStream << " ROL::Problem::finalize" << std::endl;
436 outStream << " Problem already finalized!" << std::endl;
437 outStream << std::endl;
438 }
439 }
440}
441
442template<typename Real>
443const Ptr<Objective<Real>>& Problem<Real>::getObjective() {
444 finalize();
445 return obj_;
446}
447
448template<typename Real>
449const Ptr<Vector<Real>>& Problem<Real>::getPrimalOptimizationVector() {
450 finalize();
451 return xprim_;
452}
453
454template<typename Real>
455const Ptr<Vector<Real>>& Problem<Real>::getDualOptimizationVector() {
456 finalize();
457 return xdual_;
458}
459
460template<typename Real>
461const Ptr<BoundConstraint<Real>>& Problem<Real>::getBoundConstraint() {
462 finalize();
463 return bnd_;
464}
465
466template<typename Real>
467const Ptr<Constraint<Real>>& Problem<Real>::getConstraint() {
468 finalize();
469 return con_;
470}
471
472template<typename Real>
473const Ptr<Vector<Real>>& Problem<Real>::getMultiplierVector() {
474 finalize();
475 return mul_;
476}
477
478template<typename Real>
479const Ptr<Vector<Real>>& Problem<Real>::getResidualVector() {
480 finalize();
481 return res_;
482}
483
484template<typename Real>
485const Ptr<PolyhedralProjection<Real>>& Problem<Real>::getPolyhedralProjection() {
486 finalize();
487 return proj_;
488}
489
490template<typename Real>
491EProblem Problem<Real>::getProblemType() {
492 finalize();
493 return problemType_;
494}
495
496template<typename Real>
497Real Problem<Real>::checkLinearity(bool printToStream, std::ostream &outStream) const {
498 std::ios_base::fmtflags state(outStream.flags());
499 if (printToStream) {
500 outStream << std::setprecision(8) << std::scientific;
501 outStream << std::endl;
502 outStream << " ROL::Problem::checkLinearity" << std::endl;
503 }
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);
518 cnorm = c1->norm();
519 it->second.constraint->update(*x,UpdateType::Temp);
520 it->second.constraint->value(*c2,*x,tol);
521 c1->axpy(-one,*c2);
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);
527 c1->axpy(alpha,*c2);
528 err = c1->norm();
529 maxerr = std::max(err,maxerr);
530 if (printToStream) {
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;
535 }
536 }
537 }
538 if (printToStream) {
539 outStream << std::endl;
540 }
541 outStream.flags(state);
542 return maxerr;
543}
544
545template<typename Real>
546void Problem<Real>::checkVectors(bool printToStream, std::ostream &outStream) const {
547 const Real one(1);
548 Ptr<Vector<Real>> x, y;
549 // Primal optimization space vector
550 x = INPUT_xprim_->clone(); x->randomize(-one,one);
551 y = INPUT_xprim_->clone(); y->randomize(-one,one);
552 if (printToStream) {
553 outStream << std::endl << " Check primal optimization space vector" << std::endl;
554 }
555 INPUT_xprim_->checkVector(*x,*y,printToStream,outStream);
556
557 // Dual optimization space vector
558 x = INPUT_xdual_->clone(); x->randomize(-one,one);
559 y = INPUT_xdual_->clone(); y->randomize(-one,one);
560 if (printToStream) {
561 outStream << std::endl << " Check dual optimization space vector" << std::endl;
562 }
563 INPUT_xdual_->checkVector(*x,*y,printToStream,outStream);
564
565 // Check constraint space vectors
566 for (auto it = INPUT_con_.begin(); it != INPUT_con_.end(); ++it) {
567 // Primal constraint space vector
568 x = it->second.residual->clone(); x->randomize(-one,one);
569 y = it->second.residual->clone(); y->randomize(-one,one);
570 if (printToStream) {
571 outStream << std::endl << " " << it->first << ": Check primal constraint space vector" << std::endl;
572 }
573 it->second.residual->checkVector(*x,*y,printToStream,outStream);
574
575 // Dual optimization space vector
576 x = it->second.multiplier->clone(); x->randomize(-one,one);
577 y = it->second.multiplier->clone(); y->randomize(-one,one);
578 if (printToStream) {
579 outStream << std::endl << " " << it->first << ": Check dual constraint space vector" << std::endl;
580 }
581 it->second.multiplier->checkVector(*x,*y,printToStream,outStream);
582 }
583
584 // Check constraint space vectors
585 for (auto it = INPUT_linear_con_.begin(); it != INPUT_linear_con_.end(); ++it) {
586 // Primal constraint space vector
587 x = it->second.residual->clone(); x->randomize(-one,one);
588 y = it->second.residual->clone(); y->randomize(-one,one);
589 if (printToStream) {
590 outStream << std::endl << " " << it->first << ": Check primal linear constraint space vector" << std::endl;
591 }
592 it->second.residual->checkVector(*x,*y,printToStream,outStream);
593
594 // Dual optimization space vector
595 x = it->second.multiplier->clone(); x->randomize(-one,one);
596 y = it->second.multiplier->clone(); y->randomize(-one,one);
597 if (printToStream) {
598 outStream << std::endl << " " << it->first << ": Check dual linear constraint space vector" << std::endl;
599 }
600 it->second.multiplier->checkVector(*x,*y,printToStream,outStream);
601 }
602}
603
604template<typename Real>
605void Problem<Real>::checkDerivatives(bool printToStream, std::ostream &outStream, const Ptr<Vector<Real>> &x0, Real scale) const {
606 const Real one(1);
607 Ptr<Vector<Real>> x, d, v, g, c, w;
608 // Objective check
609 x = x0;
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);
614 if (printToStream)
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);
619
620 // Constraint check
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);
624 if (printToStream)
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);
629 }
630
631 // Linear constraint check
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);
635 if (printToStream)
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);
640 }
641}
642
643template<typename Real>
644void Problem<Real>::check(bool printToStream, std::ostream &outStream, const Ptr<Vector<Real>> &x0, Real scale) const {
645 checkVectors(printToStream,outStream);
646 if (hasLinearEquality_ || hasLinearInequality_)
647 checkLinearity(printToStream,outStream);
648 checkDerivatives(printToStream,outStream,x0,scale);
649}
650
651template<typename Real>
652bool Problem<Real>::isFinalized() const {
653 return isFinalized_;
654}
655
656template<typename Real>
657void Problem<Real>::edit() {
658 isFinalized_ = false;
659 rlc_ = nullPtr;
660 proj_ = nullPtr;
661}
662
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());
669 }
670 else {
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());
675 }
676 }
677}
678
679} // namespace ROL
680
681#endif // ROL_NEWOPTIMIZATIONPROBLEM_DEF_HPP
const Ptr< Obj > obj_
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.
Definition: ROL_Vector.hpp:84
EProblem
Definition: ROL_Types.hpp:257
@ TYPE_U
Definition: ROL_Types.hpp:258
@ TYPE_E
Definition: ROL_Types.hpp:260
@ TYPE_EB
Definition: ROL_Types.hpp:261
@ TYPE_B
Definition: ROL_Types.hpp:259