ROL
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 namespace ROL {
16 
17 template<typename Real>
18 Problem<Real>::Problem( const Ptr<Objective<Real>> &obj,
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),
27  problemType_(TYPE_U) {
28  INPUT_obj_ = obj;
29  INPUT_xprim_ = x;
30  INPUT_bnd_ = nullPtr;
31  INPUT_con_.clear();
32  INPUT_linear_con_.clear();
33  if (g==nullPtr) INPUT_xdual_ = x->dual().clone();
34  else INPUT_xdual_ = g;
35 }
36 
37 template<typename Real>
38 void Problem<Real>::addBoundConstraint(const Ptr<BoundConstraint<Real>> &bnd) {
39  ROL_TEST_FOR_EXCEPTION(isFinalized_,std::invalid_argument,
40  ">>> ROL::Problem: Cannot add bounds after problem is finalized!");
41 
42  INPUT_bnd_ = bnd;
43  hasBounds_ = true;
44 }
45 
46 template<typename Real>
48  ROL_TEST_FOR_EXCEPTION(isFinalized_,std::invalid_argument,
49  ">>> ROL::Problem: Cannot remove bounds after problem is finalized!");
50 
51  INPUT_bnd_ = nullPtr;
52  hasBounds_ = false;
53 }
54 
55 template<typename Real>
56 void Problem<Real>::addConstraint( std::string name,
57  const Ptr<Constraint<Real>> &econ,
58  const Ptr<Vector<Real>> &emul,
59  const Ptr<Vector<Real>> &eres,
60  bool reset) {
61  ROL_TEST_FOR_EXCEPTION(isFinalized_,std::invalid_argument,
62  ">>> ROL::Problem: Cannot add constraint after problem is finalized!");
63 
64  if (reset) INPUT_con_.clear();
65 
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!");
69 
70  INPUT_con_.insert({name,ConstraintData<Real>(econ,emul,eres)});
71  hasEquality_ = true;
72  cnt_econ_++;
73 }
74 
75 template<typename Real>
76 void Problem<Real>::addConstraint( std::string name,
77  const Ptr<Constraint<Real>> &icon,
78  const Ptr<Vector<Real>> &imul,
79  const Ptr<BoundConstraint<Real>> &ibnd,
80  const Ptr<Vector<Real>> &ires,
81  bool reset) {
82  ROL_TEST_FOR_EXCEPTION(isFinalized_,std::invalid_argument,
83  ">>> ROL::Problem: Cannot add constraint after problem is finalized!");
84 
85  if (reset) INPUT_con_.clear();
86 
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!");
90 
91  INPUT_con_.insert({name,ConstraintData<Real>(icon,imul,ires,ibnd)});
92  hasInequality_ = true;
93  cnt_icon_++;
94 }
95 
96 template<typename Real>
97 void Problem<Real>::removeConstraint(std::string name) {
98  ROL_TEST_FOR_EXCEPTION(isFinalized_,std::invalid_argument,
99  ">>> ROL::Problem: Cannot remove constraint after problem is finalized!");
100 
101  auto it = INPUT_con_.find(name);
102  if (it!=INPUT_con_.end()) {
103  if (it->second.bounds==nullPtr) cnt_econ_--;
104  else cnt_icon_--;
105  INPUT_con_.erase(it);
106  }
107  if (cnt_econ_==0) hasEquality_ = false;
108  if (cnt_icon_==0) hasInequality_ = false;
109 }
110 
111 template<typename Real>
112 void Problem<Real>::addLinearConstraint( std::string name,
113  const Ptr<Constraint<Real>> &linear_econ,
114  const Ptr<Vector<Real>> &linear_emul,
115  const Ptr<Vector<Real>> &linear_eres,
116  bool reset) {
117  ROL_TEST_FOR_EXCEPTION(isFinalized_,std::invalid_argument,
118  ">>> ROL::Problem: Cannot add linear constraint after problem is finalized!");
119 
120  if (reset) INPUT_linear_con_.clear();
121 
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!");
125 
126  INPUT_linear_con_.insert({name,ConstraintData<Real>(linear_econ,linear_emul,linear_eres)});
127  hasLinearEquality_ = true;
128  cnt_linear_econ_++;
129 }
130 
131 template<typename Real>
132 void Problem<Real>::addLinearConstraint( std::string name,
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,
137  bool reset) {
138  ROL_TEST_FOR_EXCEPTION(isFinalized_,std::invalid_argument,
139  ">>> ROL::Problem: Cannot add linear constraint after problem is finalized!");
140 
141  if (reset) INPUT_linear_con_.clear();
142 
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!");
146 
147  INPUT_linear_con_.insert({name,ConstraintData<Real>(linear_icon,linear_imul,linear_ires,linear_ibnd)});
148  hasLinearInequality_ = true;
149  cnt_linear_icon_++;
150 }
151 
152 template<typename Real>
153 void Problem<Real>::removeLinearConstraint(std::string name) {
154  ROL_TEST_FOR_EXCEPTION(isFinalized_,std::invalid_argument,
155  ">>> ROL::Problem: Cannot remove linear inequality after problem is finalized!");
156 
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);
162  }
163  if (cnt_linear_econ_==0) hasLinearEquality_ = false;
164  if (cnt_linear_icon_==0) hasLinearInequality_ = false;
165 }
166 
167 template<typename Real>
168 void Problem<Real>::setProjectionAlgorithm(ParameterList &list) {
169  ROL_TEST_FOR_EXCEPTION(isFinalized_,std::invalid_argument,
170  ">>> ROL::Problem: Cannot set polyhedral projection algorithm after problem is finalized!");
171 
172  ppa_list_ = list;
173 }
174 
175 template<typename Real>
176 void Problem<Real>::finalize(bool lumpConstraints, bool printToStream, std::ostream &outStream) {
177  if (!isFinalized_) {
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;
190  }
191  else {
192  lcon.insert(INPUT_linear_con_.begin(),INPUT_linear_con_.end());
193  }
194  // Transform optimization problem
195  //std::cout << hasBounds_ << " " << hasEquality << " " << hasInequality << " " << hasLinearEquality << " " << hasLinearInequality << std::endl;
196  if (!hasLinearEquality && !hasLinearInequality) {
197  proj_ = nullPtr;
198  if (!hasEquality && !hasInequality && !hasBounds_) {
199  problemType_ = TYPE_U;
200  obj_ = INPUT_obj_;
201  xprim_ = INPUT_xprim_;
202  xdual_ = INPUT_xdual_;
203  bnd_ = nullPtr;
204  con_ = nullPtr;
205  mul_ = nullPtr;
206  res_ = nullPtr;
207  }
208  else if (!hasEquality && !hasInequality && hasBounds_) {
209  problemType_ = TYPE_B;
210  obj_ = INPUT_obj_;
211  xprim_ = INPUT_xprim_;
212  xdual_ = INPUT_xdual_;
213  bnd_ = INPUT_bnd_;
214  con_ = nullPtr;
215  mul_ = nullPtr;
216  res_ = nullPtr;
217  }
218  else if (hasEquality && !hasInequality && !hasBounds_) {
219  ConstraintAssembler<Real> cm(con,INPUT_xprim_,INPUT_xdual_);
220  problemType_ = TYPE_E;
221  obj_ = INPUT_obj_;
222  xprim_ = INPUT_xprim_;
223  xdual_ = INPUT_xdual_;
224  bnd_ = nullPtr;
225  con_ = cm.getConstraint();
226  mul_ = cm.getMultiplier();
227  res_ = cm.getResidual();
228  }
229  else {
230  ConstraintAssembler<Real> cm(con,INPUT_xprim_,INPUT_xdual_,INPUT_bnd_);
231  problemType_ = TYPE_EB;
232  obj_ = INPUT_obj_;
233  if (cm.hasInequality()) {
234  obj_ = makePtr<SlacklessObjective<Real>>(INPUT_obj_);
235  }
236  xprim_ = cm.getOptVector();
237  xdual_ = cm.getDualOptVector();
238  bnd_ = cm.getBoundConstraint();
239  con_ = cm.getConstraint();
240  mul_ = cm.getMultiplier();
241  res_ = cm.getResidual();
242  }
243  }
244  else {
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());
249  proj_ = nullPtr;
250  if (!hasEquality && !hasInequality) {
251  problemType_ = TYPE_U;
252  obj_ = rlc_->transform(INPUT_obj_);
253  xprim_ = xfeas_->clone(); xprim_->zero();
254  xdual_ = cm.getDualOptVector();
255  bnd_ = nullPtr;
256  con_ = nullPtr;
257  mul_ = nullPtr;
258  res_ = nullPtr;
259  }
260  else {
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)));
265  }
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) {
274  problemType_ = TYPE_E;
275  obj_ = rlc_->transform(INPUT_obj_);
276  bnd_ = nullPtr;
277  }
278  else {
279  problemType_ = TYPE_EB;
280  obj_ = makePtr<SlacklessObjective<Real>>(rlc_->transform(INPUT_obj_));
281  bnd_ = cm1.getBoundConstraint();
282  }
283  }
284  }
285  else if ((hasBounds_ || hasLinearInequality) && !hasEquality && !hasInequality) {
286  ConstraintAssembler<Real> cm(lcon,INPUT_xprim_,INPUT_xdual_,INPUT_bnd_);
287  problemType_ = TYPE_B;
288  obj_ = INPUT_obj_;
289  if (cm.hasInequality()) {
290  obj_ = makePtr<SlacklessObjective<Real>>(INPUT_obj_);
291  }
292  xprim_ = cm.getOptVector();
293  xdual_ = cm.getDualOptVector();
294  bnd_ = cm.getBoundConstraint();
295  con_ = nullPtr;
296  mul_ = nullPtr;
297  res_ = nullPtr;
298  proj_ = PolyhedralProjectionFactory<Real>(*xprim_,*xdual_,bnd_,
299  cm.getConstraint(),*cm.getMultiplier(),*cm.getResidual(),ppa_list_);
300  }
301  else {
302  ConstraintAssembler<Real> cm(con,lcon,INPUT_xprim_,INPUT_xdual_,INPUT_bnd_);
303  problemType_ = TYPE_EB;
304  obj_ = INPUT_obj_;
305  if (cm.hasInequality()) {
306  obj_ = makePtr<SlacklessObjective<Real>>(INPUT_obj_);
307  }
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_);
317  }
318  }
319  isFinalized_ = true;
320  if (printToStream) {
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;
326  if (hasEquality) {
327  int cnt = 0;
328  for (auto it = con.begin(); it != con.end(); ++it) {
329  if (it->second.bounds==nullPtr) {
330  if (cnt==0) {
331  outStream << " Names: ........................... ";
332  cnt++;
333  }
334  else {
335  outStream << " ";
336  }
337  outStream << it->first << std::endl;
338  }
339  }
340  outStream << " Total: ........................... " << cnt_econ_+(lumpConstraints ? cnt_linear_econ_ : 0) << std::endl;
341  }
342  outStream << " Has Inequality Constraint? ......... " << (hasInequality ? "yes" : "no") << std::endl;
343  if (hasInequality) {
344  int cnt = 0;
345  for (auto it = con.begin(); it != con.end(); ++it) {
346  if (it->second.bounds!=nullPtr) {
347  if (cnt==0) {
348  outStream << " Names: ........................... ";
349  cnt++;
350  }
351  else {
352  outStream << " ";
353  }
354  outStream << it->first << std::endl;
355  }
356  }
357  outStream << " Total: ........................... " << cnt_icon_+(lumpConstraints ? cnt_linear_icon_ : 0) << std::endl;
358  }
359  if (!lumpConstraints) {
360  outStream << " Has Linear Equality Constraint? .... " << (hasLinearEquality ? "yes" : "no") << std::endl;
361  if (hasLinearEquality) {
362  int cnt = 0;
363  for (auto it = lcon.begin(); it != lcon.end(); ++it) {
364  if (it->second.bounds==nullPtr) {
365  if (cnt==0) {
366  outStream << " Names: ........................... ";
367  cnt++;
368  }
369  else {
370  outStream << " ";
371  }
372  outStream << it->first << std::endl;
373  }
374  }
375  outStream << " Total: ........................... " << cnt_linear_econ_ << std::endl;
376  }
377  outStream << " Has Linear Inequality Constraint? .. " << (hasLinearInequality ? "yes" : "no") << std::endl;
378  if (hasLinearInequality) {
379  int cnt = 0;
380  for (auto it = lcon.begin(); it != lcon.end(); ++it) {
381  if (it->second.bounds!=nullPtr) {
382  if (cnt==0) {
383  outStream << " Names: ........................... ";
384  cnt++;
385  }
386  else {
387  outStream << " ";
388  }
389  outStream << it->first << std::endl;
390  }
391  }
392  outStream << " Total: ........................... " << cnt_linear_icon_ << std::endl;
393  }
394  }
395  outStream << std::endl;
396  }
397  }
398  else {
399  if (printToStream) {
400  outStream << std::endl;
401  outStream << " ROL::Problem::finalize" << std::endl;
402  outStream << " Problem already finalized!" << std::endl;
403  outStream << std::endl;
404  }
405  }
406 }
407 
408 template<typename Real>
409 const Ptr<Objective<Real>>& Problem<Real>::getObjective() {
410  finalize();
411  return obj_;
412 }
413 
414 template<typename Real>
415 const Ptr<Vector<Real>>& Problem<Real>::getPrimalOptimizationVector() {
416  finalize();
417  return xprim_;
418 }
419 
420 template<typename Real>
421 const Ptr<Vector<Real>>& Problem<Real>::getDualOptimizationVector() {
422  finalize();
423  return xdual_;
424 }
425 
426 template<typename Real>
427 const Ptr<BoundConstraint<Real>>& Problem<Real>::getBoundConstraint() {
428  finalize();
429  return bnd_;
430 }
431 
432 template<typename Real>
433 const Ptr<Constraint<Real>>& Problem<Real>::getConstraint() {
434  finalize();
435  return con_;
436 }
437 
438 template<typename Real>
439 const Ptr<Vector<Real>>& Problem<Real>::getMultiplierVector() {
440  finalize();
441  return mul_;
442 }
443 
444 template<typename Real>
445 const Ptr<Vector<Real>>& Problem<Real>::getResidualVector() {
446  finalize();
447  return res_;
448 }
449 
450 template<typename Real>
451 const Ptr<PolyhedralProjection<Real>>& Problem<Real>::getPolyhedralProjection() {
452  finalize();
453  return proj_;
454 }
455 
456 template<typename Real>
458  finalize();
459  return problemType_;
460 }
461 
462 template<typename Real>
463 Real Problem<Real>::checkLinearity(bool printToStream, std::ostream &outStream) const {
464  std::ios_base::fmtflags state(outStream.flags());
465  if (printToStream) {
466  outStream << std::setprecision(8) << std::scientific;
467  outStream << std::endl;
468  outStream << " ROL::Problem::checkLinearity" << std::endl;
469  }
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();
482  it->second.constraint->update(*xy,UpdateType::Temp);
483  it->second.constraint->value(*c1,*xy,tol);
484  cnorm = c1->norm();
485  it->second.constraint->update(*x,UpdateType::Temp);
486  it->second.constraint->value(*c2,*x,tol);
487  c1->axpy(-one,*c2);
488  it->second.constraint->update(*y,UpdateType::Temp);
489  it->second.constraint->value(*c2,*y,tol);
490  c1->axpy(-alpha,*c2);
491  it->second.constraint->update(*z,UpdateType::Temp);
492  it->second.constraint->value(*c2,*z,tol);
493  c1->axpy(alpha,*c2);
494  err = c1->norm();
495  maxerr = std::max(err,maxerr);
496  if (printToStream) {
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;
501  }
502  }
503  }
504  if (printToStream) {
505  outStream << std::endl;
506  }
507  outStream.flags(state);
508  return maxerr;
509 }
510 
511 template<typename Real>
512 void Problem<Real>::checkVectors(bool printToStream, std::ostream &outStream) const {
513  const Real one(1);
514  Ptr<Vector<Real>> x, y;
515  // Primal optimization space vector
516  x = INPUT_xprim_->clone(); x->randomize(-one,one);
517  y = INPUT_xprim_->clone(); y->randomize(-one,one);
518  if (printToStream) {
519  outStream << std::endl << " Check primal optimization space vector" << std::endl;
520  }
521  INPUT_xprim_->checkVector(*x,*y,printToStream,outStream);
522 
523  // Dual optimization space vector
524  x = INPUT_xdual_->clone(); x->randomize(-one,one);
525  y = INPUT_xdual_->clone(); y->randomize(-one,one);
526  if (printToStream) {
527  outStream << std::endl << " Check dual optimization space vector" << std::endl;
528  }
529  INPUT_xdual_->checkVector(*x,*y,printToStream,outStream);
530 
531  // Check constraint space vectors
532  for (auto it = INPUT_con_.begin(); it != INPUT_con_.end(); ++it) {
533  // Primal constraint space vector
534  x = it->second.residual->clone(); x->randomize(-one,one);
535  y = it->second.residual->clone(); y->randomize(-one,one);
536  if (printToStream) {
537  outStream << std::endl << " " << it->first << ": Check primal constraint space vector" << std::endl;
538  }
539  it->second.residual->checkVector(*x,*y,printToStream,outStream);
540 
541  // Dual optimization space vector
542  x = it->second.multiplier->clone(); x->randomize(-one,one);
543  y = it->second.multiplier->clone(); y->randomize(-one,one);
544  if (printToStream) {
545  outStream << std::endl << " " << it->first << ": Check dual constraint space vector" << std::endl;
546  }
547  it->second.multiplier->checkVector(*x,*y,printToStream,outStream);
548  }
549 
550  // Check constraint space vectors
551  for (auto it = INPUT_linear_con_.begin(); it != INPUT_linear_con_.end(); ++it) {
552  // Primal constraint space vector
553  x = it->second.residual->clone(); x->randomize(-one,one);
554  y = it->second.residual->clone(); y->randomize(-one,one);
555  if (printToStream) {
556  outStream << std::endl << " " << it->first << ": Check primal linear constraint space vector" << std::endl;
557  }
558  it->second.residual->checkVector(*x,*y,printToStream,outStream);
559 
560  // Dual optimization space vector
561  x = it->second.multiplier->clone(); x->randomize(-one,one);
562  y = it->second.multiplier->clone(); y->randomize(-one,one);
563  if (printToStream) {
564  outStream << std::endl << " " << it->first << ": Check dual linear constraint space vector" << std::endl;
565  }
566  it->second.multiplier->checkVector(*x,*y,printToStream,outStream);
567  }
568 }
569 
570 template<typename Real>
571 void Problem<Real>::checkDerivatives(bool printToStream, std::ostream &outStream, const Ptr<Vector<Real>> &x0, Real scale) const {
572  const Real one(1);
573  Ptr<Vector<Real>> x, d, v, g, c, w;
574  // Objective check
575  x = x0;
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);
580  if (printToStream)
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);
585 
586  // Constraint check
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);
590  if (printToStream)
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);
595  }
596 
597  // Linear constraint check
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);
601  if (printToStream)
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);
606  }
607 }
608 
609 template<typename Real>
610 void Problem<Real>::check(bool printToStream, std::ostream &outStream, const Ptr<Vector<Real>> &x0, Real scale) const {
611  checkVectors(printToStream,outStream);
612  if (hasLinearEquality_ || hasLinearInequality_)
613  checkLinearity(printToStream,outStream);
614  checkDerivatives(printToStream,outStream,x0,scale);
615 }
616 
617 template<typename Real>
618 bool Problem<Real>::isFinalized() const {
619  return isFinalized_;
620 }
621 
622 template<typename Real>
623 void Problem<Real>::edit() {
624  isFinalized_ = false;
625  rlc_ = nullPtr;
626  proj_ = nullPtr;
627 }
628 
629 template<typename Real>
631  if (rlc_ != nullPtr) {
632  if (!hasInequality_) {
633  rlc_->project(*INPUT_xprim_,*xprim_);
634  INPUT_xprim_->plus(*rlc_->getFeasibleVector());
635  }
636  else {
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());
641  }
642  }
643 }
644 
645 } // namespace ROL
646 
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.
Definition: ROL_Vector.hpp:46
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.
EProblem
Definition: ROL_Types.hpp:223
const Ptr< Obj > obj_
void setProjectionAlgorithm(ParameterList &parlist)
Set polyhedral projection algorithm.