ROL
function/prox/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>
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>
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>
57  ROL_TEST_FOR_EXCEPTION(isFinalized_,std::invalid_argument,
58  ">>> ROL::Problem: Cannot add prox objective after problem is finalized!");
59 
60  INPUT_prox_ = prox;
61  hasProx_ = true;
62 }
63 
64 template<typename Real>
66  ROL_TEST_FOR_EXCEPTION(isFinalized_,std::invalid_argument,
67  ">>> ROL::Problem: Cannot remove prox objective after problem is finalized!");
68 
69  INPUT_prox_ = nullPtr;
70  hasProx_ = false;
71 }
72 
73 template<typename Real>
74 void Problem<Real>::addConstraint( std::string name,
75  const Ptr<Constraint<Real>> &econ,
76  const Ptr<Vector<Real>> &emul,
77  const Ptr<Vector<Real>> &eres,
78  bool reset) {
79  ROL_TEST_FOR_EXCEPTION(isFinalized_,std::invalid_argument,
80  ">>> ROL::Problem: Cannot add constraint after problem is finalized!");
81 
82  if (reset) INPUT_con_.clear();
83 
84  auto it = INPUT_con_.find(name);
85  ROL_TEST_FOR_EXCEPTION(it != INPUT_con_.end(),std::invalid_argument,
86  ">>> ROL::Problem: Constraint names must be distinct!");
87 
88  INPUT_con_.insert({name,ConstraintData<Real>(econ,emul,eres)});
89  hasEquality_ = true;
90  cnt_econ_++;
91 }
92 
93 template<typename Real>
94 void Problem<Real>::addConstraint( std::string name,
95  const Ptr<Constraint<Real>> &icon,
96  const Ptr<Vector<Real>> &imul,
97  const Ptr<BoundConstraint<Real>> &ibnd,
98  const Ptr<Vector<Real>> &ires,
99  bool reset) {
100  ROL_TEST_FOR_EXCEPTION(isFinalized_,std::invalid_argument,
101  ">>> ROL::Problem: Cannot add constraint after problem is finalized!");
102 
103  if (reset) INPUT_con_.clear();
104 
105  auto it = INPUT_con_.find(name);
106  ROL_TEST_FOR_EXCEPTION(it != INPUT_con_.end(),std::invalid_argument,
107  ">>> ROL::Problem: Constraint names must be distinct!");
108 
109  INPUT_con_.insert({name,ConstraintData<Real>(icon,imul,ires,ibnd)});
110  hasInequality_ = true;
111  cnt_icon_++;
112 }
113 
114 template<typename Real>
115 void Problem<Real>::removeConstraint(std::string name) {
116  ROL_TEST_FOR_EXCEPTION(isFinalized_,std::invalid_argument,
117  ">>> ROL::Problem: Cannot remove constraint after problem is finalized!");
118 
119  auto it = INPUT_con_.find(name);
120  if (it!=INPUT_con_.end()) {
121  if (it->second.bounds==nullPtr) cnt_econ_--;
122  else cnt_icon_--;
123  INPUT_con_.erase(it);
124  }
125  if (cnt_econ_==0) hasEquality_ = false;
126  if (cnt_icon_==0) hasInequality_ = false;
127 }
128 
129 template<typename Real>
130 void Problem<Real>::addLinearConstraint( std::string name,
131  const Ptr<Constraint<Real>> &linear_econ,
132  const Ptr<Vector<Real>> &linear_emul,
133  const Ptr<Vector<Real>> &linear_eres,
134  bool reset) {
135  ROL_TEST_FOR_EXCEPTION(isFinalized_,std::invalid_argument,
136  ">>> ROL::Problem: Cannot add linear constraint after problem is finalized!");
137 
138  if (reset) INPUT_linear_con_.clear();
139 
140  auto it = INPUT_linear_con_.find(name);
141  ROL_TEST_FOR_EXCEPTION(it != INPUT_linear_con_.end(),std::invalid_argument,
142  ">>> ROL::Problem: Linear constraint names must be distinct!");
143 
144  INPUT_linear_con_.insert({name,ConstraintData<Real>(linear_econ,linear_emul,linear_eres)});
145  hasLinearEquality_ = true;
146  cnt_linear_econ_++;
147 }
148 
149 template<typename Real>
150 void Problem<Real>::addLinearConstraint( std::string name,
151  const Ptr<Constraint<Real>> &linear_icon,
152  const Ptr<Vector<Real>> &linear_imul,
153  const Ptr<BoundConstraint<Real>> &linear_ibnd,
154  const Ptr<Vector<Real>> &linear_ires,
155  bool reset) {
156  ROL_TEST_FOR_EXCEPTION(isFinalized_,std::invalid_argument,
157  ">>> ROL::Problem: Cannot add linear constraint after problem is finalized!");
158 
159  if (reset) INPUT_linear_con_.clear();
160 
161  auto it = INPUT_linear_con_.find(name);
162  ROL_TEST_FOR_EXCEPTION(it != INPUT_linear_con_.end(),std::invalid_argument,
163  ">>> ROL::Problem: Linear constraint names must be distinct!");
164 
165  INPUT_linear_con_.insert({name,ConstraintData<Real>(linear_icon,linear_imul,linear_ires,linear_ibnd)});
166  hasLinearInequality_ = true;
167  cnt_linear_icon_++;
168 }
169 
170 template<typename Real>
172  ROL_TEST_FOR_EXCEPTION(isFinalized_,std::invalid_argument,
173  ">>> ROL::Problem: Cannot remove linear inequality after problem is finalized!");
174 
175  auto it = INPUT_linear_con_.find(name);
176  if (it!=INPUT_linear_con_.end()) {
177  if (it->second.bounds==nullPtr) cnt_linear_econ_--;
178  else cnt_linear_icon_--;
179  INPUT_linear_con_.erase(it);
180  }
181  if (cnt_linear_econ_==0) hasLinearEquality_ = false;
182  if (cnt_linear_icon_==0) hasLinearInequality_ = false;
183 }
184 
185 template<typename Real>
186 void Problem<Real>::setProjectionAlgorithm(ParameterList &list) {
187  ROL_TEST_FOR_EXCEPTION(isFinalized_,std::invalid_argument,
188  ">>> ROL::Problem: Cannot set polyhedral projection algorithm after problem is finalized!");
189 
190  ppa_list_ = list;
191 }
192 
193 template<typename Real>
194 void Problem<Real>::finalize(bool lumpConstraints, bool printToStream, std::ostream &outStream) {
195  if (!isFinalized_) {
196  std::unordered_map<std::string,ConstraintData<Real>> con, lcon, icon;
197  bool hasEquality = hasEquality_;
198  bool hasLinearEquality = hasLinearEquality_;
199  bool hasInequality = hasInequality_;
200  bool hasLinearInequality = hasLinearInequality_;
201  con.insert(INPUT_con_.begin(),INPUT_con_.end());
202  if (lumpConstraints) {
203  con.insert(INPUT_linear_con_.begin(),INPUT_linear_con_.end());
204  hasEquality = (hasEquality || hasLinearEquality);
205  hasInequality = (hasInequality || hasLinearInequality);
206  hasLinearEquality = false;
207  hasLinearInequality = false;
208  }
209  else {
210  lcon.insert(INPUT_linear_con_.begin(),INPUT_linear_con_.end());
211  }
212  // Transform optimization problem
213  //std::cout << hasBounds_ << " " << hasEquality << " " << hasInequality << " " << hasLinearEquality << " " << hasLinearInequality << std::endl;
214  bool proxCompatible = false;
215  if (!hasLinearEquality && !hasLinearInequality) {
216  proj_ = nullPtr;
217  if (!hasEquality && !hasInequality && !hasBounds_) {
218  problemType_ = TYPE_U;
219  obj_ = INPUT_obj_;
220  xprim_ = INPUT_xprim_;
221  xdual_ = INPUT_xdual_;
222  bnd_ = nullPtr;
223  con_ = nullPtr;
224  mul_ = nullPtr;
225  res_ = nullPtr;
226  if (hasProx_) { prox_ = INPUT_prox_; }
227  else { prox_ = nullPtr; }
228  proxCompatible = true;
229  }
230  else if (!hasEquality && !hasInequality && hasBounds_) {
231  problemType_ = TYPE_B;
232  obj_ = INPUT_obj_;
233  xprim_ = INPUT_xprim_;
234  xdual_ = INPUT_xdual_;
235  bnd_ = INPUT_bnd_;
236  con_ = nullPtr;
237  mul_ = nullPtr;
238  res_ = nullPtr;
239  }
240  else if (hasEquality && !hasInequality && !hasBounds_) {
241  ConstraintAssembler<Real> cm(con,INPUT_xprim_,INPUT_xdual_);
242  problemType_ = TYPE_E;
243  obj_ = INPUT_obj_;
244  xprim_ = INPUT_xprim_;
245  xdual_ = INPUT_xdual_;
246  bnd_ = nullPtr;
247  con_ = cm.getConstraint();
248  mul_ = cm.getMultiplier();
249  res_ = cm.getResidual();
250  }
251  else {
252  ConstraintAssembler<Real> cm(con,INPUT_xprim_,INPUT_xdual_,INPUT_bnd_);
253  problemType_ = TYPE_EB;
254  obj_ = INPUT_obj_;
255  if (cm.hasInequality()) {
256  obj_ = makePtr<SlacklessObjective<Real>>(INPUT_obj_);
257  }
258  xprim_ = cm.getOptVector();
259  xdual_ = cm.getDualOptVector();
260  bnd_ = cm.getBoundConstraint();
261  con_ = cm.getConstraint();
262  mul_ = cm.getMultiplier();
263  res_ = cm.getResidual();
264  }
265  }
266  else {
267  if (!hasBounds_ && !hasLinearInequality) {
268  ConstraintAssembler<Real> cm(lcon,INPUT_xprim_,INPUT_xdual_);
269  xfeas_ = cm.getOptVector()->clone(); xfeas_->set(*cm.getOptVector());
270  rlc_ = makePtr<ReduceLinearConstraint<Real>>(cm.getConstraint(),xfeas_,cm.getResidual());
271  proj_ = nullPtr;
272  if (!hasEquality && !hasInequality) {
273  problemType_ = TYPE_U;
274  obj_ = rlc_->transform(INPUT_obj_);
275  xprim_ = xfeas_->clone(); xprim_->zero();
276  xdual_ = cm.getDualOptVector();
277  bnd_ = nullPtr;
278  con_ = nullPtr;
279  mul_ = nullPtr;
280  res_ = nullPtr;
281  }
282  else {
283  for (auto it = con.begin(); it != con.end(); ++it) {
284  icon.insert(std::pair<std::string,ConstraintData<Real>>(it->first,
285  ConstraintData<Real>(rlc_->transform(it->second.constraint),
286  it->second.multiplier,it->second.residual,it->second.bounds)));
287  }
288  Ptr<Vector<Real>> xtmp = xfeas_->clone(); xtmp->zero();
289  ConstraintAssembler<Real> cm1(icon,xtmp,cm.getDualOptVector());
290  xprim_ = cm1.getOptVector();
291  xdual_ = cm1.getDualOptVector();
292  con_ = cm1.getConstraint();
293  mul_ = cm1.getMultiplier();
294  res_ = cm1.getResidual();
295  if (!hasInequality) {
296  problemType_ = TYPE_E;
297  obj_ = rlc_->transform(INPUT_obj_);
298  bnd_ = nullPtr;
299  }
300  else {
301  problemType_ = TYPE_EB;
302  obj_ = makePtr<SlacklessObjective<Real>>(rlc_->transform(INPUT_obj_));
303  bnd_ = cm1.getBoundConstraint();
304  }
305  }
306  }
307  else if ((hasBounds_ || hasLinearInequality) && !hasEquality && !hasInequality) {
308  ConstraintAssembler<Real> cm(lcon,INPUT_xprim_,INPUT_xdual_,INPUT_bnd_);
309  problemType_ = TYPE_B;
310  obj_ = INPUT_obj_;
311  if (cm.hasInequality()) {
312  obj_ = makePtr<SlacklessObjective<Real>>(INPUT_obj_);
313  }
314  xprim_ = cm.getOptVector();
315  xdual_ = cm.getDualOptVector();
316  bnd_ = cm.getBoundConstraint();
317  con_ = nullPtr;
318  mul_ = nullPtr;
319  res_ = nullPtr;
320  proj_ = PolyhedralProjectionFactory<Real>(*xprim_,*xdual_,bnd_,
321  cm.getConstraint(),*cm.getMultiplier(),*cm.getResidual(),ppa_list_);
322  }
323  else {
324  ConstraintAssembler<Real> cm(con,lcon,INPUT_xprim_,INPUT_xdual_,INPUT_bnd_);
325  problemType_ = TYPE_EB;
326  obj_ = INPUT_obj_;
327  if (cm.hasInequality()) {
328  obj_ = makePtr<SlacklessObjective<Real>>(INPUT_obj_);
329  }
330  xprim_ = cm.getOptVector();
331  xdual_ = cm.getDualOptVector();
332  con_ = cm.getConstraint();
333  mul_ = cm.getMultiplier();
334  res_ = cm.getResidual();
335  bnd_ = cm.getBoundConstraint();
336  proj_ = PolyhedralProjectionFactory<Real>(*xprim_,*xdual_,bnd_,
338  *cm.getLinearResidual(),ppa_list_);
339  }
340  }
341 
342  std::stringstream out;
343  out << ">>> ROL::Problem: Cannot solve ";
344  if (problemType_ == TYPE_B) out << "TypeB";
345  else if (problemType_ == TYPE_E) out << "TypeE";
346  else if (problemType_ == TYPE_EB) out << "TypeG";
347  else out << "TypeU with linear constraints";
348  out << " problems with a prox objective!";
349  ROL_TEST_FOR_EXCEPTION(hasProx_&&!proxCompatible,std::invalid_argument,out.str());
350 
351  isFinalized_ = true;
352  if (printToStream) {
353  outStream << std::endl;
354  outStream << " ROL::Problem::finalize" << std::endl;
355  outStream << " Problem Summary:" << std::endl;
356  outStream << " Has Bound Constraint? .............. " << (hasBounds_ ? "yes" : "no") << std::endl;
357  outStream << " Has Equality Constraint? ........... " << (hasEquality ? "yes" : "no") << std::endl;
358  if (hasEquality) {
359  int cnt = 0;
360  for (auto it = con.begin(); it != con.end(); ++it) {
361  if (it->second.bounds==nullPtr) {
362  if (cnt==0) {
363  outStream << " Names: ........................... ";
364  cnt++;
365  }
366  else {
367  outStream << " ";
368  }
369  outStream << it->first << std::endl;
370  }
371  }
372  outStream << " Total: ........................... " << cnt_econ_+(lumpConstraints ? cnt_linear_econ_ : 0) << std::endl;
373  }
374  outStream << " Has Inequality Constraint? ......... " << (hasInequality ? "yes" : "no") << std::endl;
375  if (hasInequality) {
376  int cnt = 0;
377  for (auto it = con.begin(); it != con.end(); ++it) {
378  if (it->second.bounds!=nullPtr) {
379  if (cnt==0) {
380  outStream << " Names: ........................... ";
381  cnt++;
382  }
383  else {
384  outStream << " ";
385  }
386  outStream << it->first << std::endl;
387  }
388  }
389  outStream << " Total: ........................... " << cnt_icon_+(lumpConstraints ? cnt_linear_icon_ : 0) << std::endl;
390  }
391  if (!lumpConstraints) {
392  outStream << " Has Linear Equality Constraint? .... " << (hasLinearEquality ? "yes" : "no") << std::endl;
393  if (hasLinearEquality) {
394  int cnt = 0;
395  for (auto it = lcon.begin(); it != lcon.end(); ++it) {
396  if (it->second.bounds==nullPtr) {
397  if (cnt==0) {
398  outStream << " Names: ........................... ";
399  cnt++;
400  }
401  else {
402  outStream << " ";
403  }
404  outStream << it->first << std::endl;
405  }
406  }
407  outStream << " Total: ........................... " << cnt_linear_econ_ << std::endl;
408  }
409  outStream << " Has Linear Inequality Constraint? .. " << (hasLinearInequality ? "yes" : "no") << std::endl;
410  if (hasLinearInequality) {
411  int cnt = 0;
412  for (auto it = lcon.begin(); it != lcon.end(); ++it) {
413  if (it->second.bounds!=nullPtr) {
414  if (cnt==0) {
415  outStream << " Names: ........................... ";
416  cnt++;
417  }
418  else {
419  outStream << " ";
420  }
421  outStream << it->first << std::endl;
422  }
423  }
424  outStream << " Total: ........................... " << cnt_linear_icon_ << std::endl;
425  }
426  }
427  outStream << std::endl;
428  }
429  }
430  else {
431  if (printToStream) {
432  outStream << std::endl;
433  outStream << " ROL::Problem::finalize" << std::endl;
434  outStream << " Problem already finalized!" << std::endl;
435  outStream << std::endl;
436  }
437  }
438 }
439 
440 template<typename Real>
441 const Ptr<Objective<Real>>& Problem<Real>::getObjective() {
442  finalize();
443  return obj_;
444 }
445 
446 template<typename Real>
447 const Ptr<Vector<Real>>& Problem<Real>::getPrimalOptimizationVector() {
448  finalize();
449  return xprim_;
450 }
451 
452 template<typename Real>
453 const Ptr<Vector<Real>>& Problem<Real>::getDualOptimizationVector() {
454  finalize();
455  return xdual_;
456 }
457 
458 template<typename Real>
459 const Ptr<BoundConstraint<Real>>& Problem<Real>::getBoundConstraint() {
460  finalize();
461  return bnd_;
462 }
463 
464 template<typename Real>
465 const Ptr<Constraint<Real>>& Problem<Real>::getConstraint() {
466  finalize();
467  return con_;
468 }
469 
470 template<typename Real>
471 const Ptr<Vector<Real>>& Problem<Real>::getMultiplierVector() {
472  finalize();
473  return mul_;
474 }
475 
476 template<typename Real>
477 const Ptr<Vector<Real>>& Problem<Real>::getResidualVector() {
478  finalize();
479  return res_;
480 }
481 
482 template<typename Real>
483 const Ptr<PolyhedralProjection<Real>>& Problem<Real>::getPolyhedralProjection() {
484  finalize();
485  return proj_;
486 }
487 
488 template<typename Real>
490  finalize();
491  return problemType_;
492 }
493 
494 template<typename Real>
495 Real Problem<Real>::checkLinearity(bool printToStream, std::ostream &outStream) const {
496  std::ios_base::fmtflags state(outStream.flags());
497  if (printToStream) {
498  outStream << std::setprecision(8) << std::scientific;
499  outStream << std::endl;
500  outStream << " ROL::Problem::checkLinearity" << std::endl;
501  }
502  const Real one(1), two(2), eps(1e-2*std::sqrt(ROL_EPSILON<Real>()));
503  Real tol(std::sqrt(ROL_EPSILON<Real>())), cnorm(0), err(0), maxerr(0);
504  Ptr<Vector<Real>> x = INPUT_xprim_->clone(); x->randomize(-one,one);
505  Ptr<Vector<Real>> y = INPUT_xprim_->clone(); y->randomize(-one,one);
506  Ptr<Vector<Real>> z = INPUT_xprim_->clone(); z->zero();
507  Ptr<Vector<Real>> xy = INPUT_xprim_->clone();
508  Real alpha = two*static_cast<Real>(rand())/static_cast<Real>(RAND_MAX)-one;
509  xy->set(*x); xy->axpy(alpha,*y);
510  Ptr<Vector<Real>> c1, c2;
511  for (auto it = INPUT_linear_con_.begin(); it != INPUT_linear_con_.end(); ++it) {
512  c1 = it->second.residual->clone();
513  c2 = it->second.residual->clone();
514  it->second.constraint->update(*xy,UpdateType::Temp);
515  it->second.constraint->value(*c1,*xy,tol);
516  cnorm = c1->norm();
517  it->second.constraint->update(*x,UpdateType::Temp);
518  it->second.constraint->value(*c2,*x,tol);
519  c1->axpy(-one,*c2);
520  it->second.constraint->update(*y,UpdateType::Temp);
521  it->second.constraint->value(*c2,*y,tol);
522  c1->axpy(-alpha,*c2);
523  it->second.constraint->update(*z,UpdateType::Temp);
524  it->second.constraint->value(*c2,*z,tol);
525  c1->axpy(alpha,*c2);
526  err = c1->norm();
527  maxerr = std::max(err,maxerr);
528  if (printToStream) {
529  outStream << " Constraint " << it->first;
530  outStream << ": ||c(x+alpha*y) - (c(x)+alpha*(c(y)-c(0)))|| = " << err << std::endl;
531  if (err > eps*cnorm) {
532  outStream << " Constraint " << it->first << " may not be linear!" << std::endl;
533  }
534  }
535  }
536  if (printToStream) {
537  outStream << std::endl;
538  }
539  outStream.flags(state);
540  return maxerr;
541 }
542 
543 template<typename Real>
544 void Problem<Real>::checkVectors(bool printToStream, std::ostream &outStream) const {
545  const Real one(1);
546  Ptr<Vector<Real>> x, y;
547  // Primal optimization space vector
548  x = INPUT_xprim_->clone(); x->randomize(-one,one);
549  y = INPUT_xprim_->clone(); y->randomize(-one,one);
550  if (printToStream) {
551  outStream << std::endl << " Check primal optimization space vector" << std::endl;
552  }
553  INPUT_xprim_->checkVector(*x,*y,printToStream,outStream);
554 
555  // Dual optimization space vector
556  x = INPUT_xdual_->clone(); x->randomize(-one,one);
557  y = INPUT_xdual_->clone(); y->randomize(-one,one);
558  if (printToStream) {
559  outStream << std::endl << " Check dual optimization space vector" << std::endl;
560  }
561  INPUT_xdual_->checkVector(*x,*y,printToStream,outStream);
562 
563  // Check constraint space vectors
564  for (auto it = INPUT_con_.begin(); it != INPUT_con_.end(); ++it) {
565  // Primal constraint space vector
566  x = it->second.residual->clone(); x->randomize(-one,one);
567  y = it->second.residual->clone(); y->randomize(-one,one);
568  if (printToStream) {
569  outStream << std::endl << " " << it->first << ": Check primal constraint space vector" << std::endl;
570  }
571  it->second.residual->checkVector(*x,*y,printToStream,outStream);
572 
573  // Dual optimization space vector
574  x = it->second.multiplier->clone(); x->randomize(-one,one);
575  y = it->second.multiplier->clone(); y->randomize(-one,one);
576  if (printToStream) {
577  outStream << std::endl << " " << it->first << ": Check dual constraint space vector" << std::endl;
578  }
579  it->second.multiplier->checkVector(*x,*y,printToStream,outStream);
580  }
581 
582  // Check constraint space vectors
583  for (auto it = INPUT_linear_con_.begin(); it != INPUT_linear_con_.end(); ++it) {
584  // Primal constraint space vector
585  x = it->second.residual->clone(); x->randomize(-one,one);
586  y = it->second.residual->clone(); y->randomize(-one,one);
587  if (printToStream) {
588  outStream << std::endl << " " << it->first << ": Check primal linear constraint space vector" << std::endl;
589  }
590  it->second.residual->checkVector(*x,*y,printToStream,outStream);
591 
592  // Dual optimization space vector
593  x = it->second.multiplier->clone(); x->randomize(-one,one);
594  y = it->second.multiplier->clone(); y->randomize(-one,one);
595  if (printToStream) {
596  outStream << std::endl << " " << it->first << ": Check dual linear constraint space vector" << std::endl;
597  }
598  it->second.multiplier->checkVector(*x,*y,printToStream,outStream);
599  }
600 }
601 
602 template<typename Real>
603 void Problem<Real>::checkDerivatives(bool printToStream, std::ostream &outStream) const {
604  const Real one(1);
605  Ptr<Vector<Real>> x, d, v, g, c, w;
606  // Objective check
607  x = INPUT_xprim_->clone(); x->randomize(-one,one);
608  d = INPUT_xprim_->clone(); d->randomize(-one,one);
609  v = INPUT_xprim_->clone(); v->randomize(-one,one);
610  g = INPUT_xdual_->clone(); g->randomize(-one,one);
611  if (printToStream) {
612  outStream << std::endl << " Check objective function" << std::endl << std::endl;
613  }
614  INPUT_obj_->checkGradient(*x,*g,*d,printToStream,outStream);
615  INPUT_obj_->checkHessVec(*x,*g,*d,printToStream,outStream);
616  INPUT_obj_->checkHessSym(*x,*g,*d,*v,printToStream,outStream);
617 
618  // Constraint check
619  for (auto it = INPUT_con_.begin(); it != INPUT_con_.end(); ++it) {
620  c = it->second.residual->clone(); c->randomize(-one,one);
621  w = it->second.multiplier->clone(); w->randomize(-one,one);
622  if (printToStream) {
623  outStream << std::endl << " " << it->first << ": Check constraint function" << std::endl << std::endl;
624  }
625  it->second.constraint->checkApplyJacobian(*x,*v,*c,printToStream,outStream);
626  it->second.constraint->checkAdjointConsistencyJacobian(*w,*v,*x,printToStream,outStream);
627  it->second.constraint->checkApplyAdjointHessian(*x,*w,*v,*g,printToStream,outStream);
628  }
629 
630  // Linear constraint check
631  for (auto it = INPUT_linear_con_.begin(); it != INPUT_linear_con_.end(); ++it) {
632  c = it->second.residual->clone(); c->randomize(-one,one);
633  w = it->second.multiplier->clone(); w->randomize(-one,one);
634  if (printToStream) {
635  outStream << std::endl << " " << it->first << ": Check constraint function" << std::endl << std::endl;
636  }
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 
643 template<typename Real>
644 void Problem<Real>::check(bool printToStream, std::ostream &outStream) const {
645  checkVectors(printToStream,outStream);
646  if (hasLinearEquality_ || hasLinearInequality_) {
647  checkLinearity(printToStream,outStream);
648  }
649  checkDerivatives(printToStream,outStream);
650 }
651 
652 template<typename Real>
654  return isFinalized_;
655 }
656 
657 template<typename Real>
659  isFinalized_ = false;
660  rlc_ = nullPtr;
661  proj_ = nullPtr;
662 }
663 
664 template<typename Real>
666  if (rlc_ != nullPtr) {
667  if (!hasInequality_) {
668  rlc_->project(*INPUT_xprim_,*xprim_);
669  INPUT_xprim_->plus(*rlc_->getFeasibleVector());
670  }
671  else {
672  Ptr<Vector<Real>> xprim = dynamic_cast<PartitionedVector<Real>&>(*xprim_).get(0)->clone();
673  xprim->set(*dynamic_cast<PartitionedVector<Real>&>(*xprim_).get(0));
674  rlc_->project(*INPUT_xprim_,*xprim);
675  INPUT_xprim_->plus(*rlc_->getFeasibleVector());
676  }
677  }
678 }
679 
680 } // namespace ROL
681 
682 #endif // ROL_NEWOPTIMIZATIONPROBLEM_DEF_HPP
void removeProxObjective()
Remove an existing prox objective.
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...
Provides the interface to evaluate objective functions.
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.
ROL::Ptr< const Vector< Real > > get(size_type i) const
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.
Defines the linear algebra of vector space on a generic partitioned vector.
const Ptr< BoundConstraint< Real > > & getBoundConstraint() const
void removeLinearConstraint(std::string name)
Remove an existing linear constraint.
const Ptr< BoundConstraint< Real > > & getBoundConstraint()
Get the bound constraint.
const Ptr< Constraint< Real > > & getConstraint() const
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:46
void addProxObjective(const Ptr< ProxObjective< Real >> &prox)
Add a prox objective.
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_
const Ptr< Vector< Real > > & getLinearResidual() const
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...
Provides a wrapper for multiple constraints.
const Ptr< Vector< Real > > & getMultiplier() const
Provides the interface to apply upper and lower bound constraints.
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.
const Ptr< Vector< Real > > & getDualOptVector() const
Ptr< BoundConstraint< Real > > INPUT_bnd_
virtual void edit()
Resume editting optimization problem after finalize has been called.
const Ptr< Vector< Real > > & getOptVector() const
const Ptr< Constraint< Real > > & getLinearConstraint() const
const Ptr< Vector< Real > > & getDualOptimizationVector()
Get the dual optimization space vector.
const Ptr< Vector< Real > > & getLinearMultiplier() const
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_
const Ptr< Vector< Real > > & getResidual() const
Defines the general constraint operator interface.
void setProjectionAlgorithm(ParameterList &parlist)
Set polyhedral projection algorithm.