ROL
ROL_InteriorPoint.hpp
Go to the documentation of this file.
1 // Rapid Optimization Library (ROL) Package
2 // Copyright (2014) Sandia Corporation
3 //
4 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
5 // license for use of this work by or on behalf of the U.S. Government.
6 //
7 // Redistribution and use in source and binary forms, with or without
8 // modification, are permitted provided that the following conditions are
9 // met:
10 //
11 // 1. Redistributions of source code must retain the above copyright
12 // notice, this list of conditions and the following disclaimer.
13 //
14 // 2. Redistributions in binary form must reproduce the above copyright
15 // notice, this list of conditions and the following disclaimer in the
16 // documentation and/or other materials provided with the distribution.
17 //
18 // 3. Neither the name of the Corporation nor the names of the
19 // contributors may be used to endorse or promote products derived from
20 // this software without specific prior written permission.
21 //
22 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
23 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
24 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
25 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
26 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
27 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
28 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
29 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
30 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
31 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
32 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
33 //
34 // Questions? Contact lead developers:
35 // Drew Kouri (dpkouri@sandia.gov) and
36 // Denis Ridzal (dridzal@sandia.gov)
37 //
38 // ************************************************************************
39 // @HEADER
40 
41 #ifndef ROL_INTERIORPOINT_H
42 #define ROL_INTERIORPOINT_H
43 
45 #include "ROL_Objective.hpp"
47 
48 namespace ROL {
49 
50 namespace InteriorPoint {
51 
56 template <class Real>
57 class PenalizedObjective : public ROL::Objective<Real> {
58 private:
59 
60  typedef Vector<Real> V;
62  typedef typename PV::size_type size_type;
63 
64  const static size_type OPT = 0;
65  const static size_type SLACK = 1;
66 
67  Teuchos::RCP<Objective<Real> > obj_;
68  Teuchos::RCP<Objective<Real> > slack_barrier_;
69  Teuchos::RCP<Objective<Real> > bc_barrier_;
70  Teuchos::RCP<PV> x_;
71  Teuchos::RCP<PV> g_;
72  Teuchos::RCP<V> scratch_;
73 
74 
75  Real mu_;
76  int nfval_;
77  int ngval_;
78  Real fval_;
79  Real gnorm_;
81 
82 public:
83 
84  // Constructor without BoundConstraint
85  PenalizedObjective( const Teuchos::RCP<Objective<Real> > &obj,
86  const Teuchos::RCP<Objective<Real> > &slack_barrier,
87  const Vector<Real> &x,
88  Real mu ) :
89  obj_(obj), slack_barrier_(slack_barrier), bc_barrier_(Teuchos::null),
90  x_(Teuchos::null), g_(Teuchos::null), scratch_(Teuchos::null),
91  mu_(mu), nfval_(0), ngval_(0), fval_(0.0), gnorm_(0.0),
92  hasBoundConstraint_(false) {
93 
94  const PV &xpv = Teuchos::dyn_cast<const PV>(x);
95 
96  x_ = Teuchos::rcp_static_cast<PV>(xpv.clone());
97  g_ = Teuchos::rcp_static_cast<PV>(xpv.dual().clone());
98  }
99 
100 
101  // Constructor with BoundConstraint
102  PenalizedObjective( const Teuchos::RCP<Objective<Real> > &obj,
103  const Teuchos::RCP<Objective<Real> > &slack_barrier,
104  const Teuchos::RCP<Objective<Real> > &bc_barrier,
105  const Vector<Real> &x,
106  Real mu ) :
107  obj_(obj), slack_barrier_(slack_barrier), bc_barrier_(bc_barrier),
108  x_(Teuchos::null), g_(Teuchos::null), scratch_(Teuchos::null),
109  mu_(mu), nfval_(0), ngval_(0), fval_(0.0), gnorm_(0.0),
110  hasBoundConstraint_(true) {
111 
112  const PV &xpv = Teuchos::dyn_cast<const PV>(x);
113 
114  x_ = Teuchos::rcp_static_cast<PV>(xpv.clone());
115  g_ = Teuchos::rcp_static_cast<PV>(xpv.dual().clone());
116 
117  scratch_ = g_->get(OPT)->clone();
118  }
119 
120 
121 
122 
123  void updatePenalty( Real mu ) {
124  mu_ = mu;
125  }
126 
128  return nfval_;
129  }
130 
132  return ngval_;
133  }
134 
135  void reset(void) {
136  nfval_ = 0.; nfval_ = 0.;
137  }
138 
147  void update( const Vector<Real> &x, bool flag = true, int iter = -1 ) {
148 
149  const PV &xpv = Teuchos::dyn_cast<const PV>(x);
150 
151  Teuchos::RCP<const V> xo = xpv.get(OPT);
152  Teuchos::RCP<const V> xs = xpv.get(SLACK);
153 
154  obj_->update(*xo,flag,iter);
155  slack_barrier_->update(*xs,flag,iter);
156 
157  if(hasBoundConstraint_) {
158  bc_barrier_->update(*xo,flag,iter);
159  }
160 
161  }
162 
169  Real value( const Vector<Real> &x, Real &tol ) {
170 
171  const PV &xpv = Teuchos::dyn_cast<const PV>(x);
172 
173  Teuchos::RCP<const V> xo = xpv.get(OPT);
174  Teuchos::RCP<const V> xs = xpv.get(SLACK);
175 
176  Real val = 0;
177 
178  // Compute objective function value
179  fval_ = obj_->value(*xo,tol);
180  Real pval = slack_barrier_->value(*xs,tol);
181 
182  val = fval_ + mu_*pval;
183 
184  if( hasBoundConstraint_ ) {
185  Real bval = bc_barrier_->value(*xo,tol);
186  val += mu_*bval;
187  }
188 
189  ++nfval_;
190 
191  return val;
192  }
193 
195  return fval_;
196  }
197 
198 
206  void gradient( Vector<Real> &g, const Vector<Real> &x, Real &tol ) {
207 
208  // Compute objective function gradient
209  const PV &xpv = Teuchos::dyn_cast<const PV>(x);
210  PV &gpv = Teuchos::dyn_cast<PV>(g);
211 
212  Teuchos::RCP<const V> xo = xpv.get(OPT);
213  Teuchos::RCP<const V> xs = xpv.get(SLACK);
214 
215  Teuchos::RCP<V> go = gpv.get(OPT);
216  Teuchos::RCP<V> gs = gpv.get(SLACK);
217 
218  obj_->gradient(*go,*xo,tol);
219 
220  if( hasBoundConstraint_ ) {
221  bc_barrier_->gradient(*scratch_,*xo,tol);
222  scratch_->scale(mu_);
223  go->plus(*scratch_);
224  }
225 
226  slack_barrier_->gradient(*gs,*xs,tol);
227  gs->scale(mu_);
228 
229  g_->set(g);
230  g_->zero(SLACK);
231 
232  gnorm_ = g.norm();
233 
234  ++ngval_;
235 
236  }
237 
239 
240  }
241 
243  return gnorm_;
244  }
245 
255  void hessVec( Vector<Real> &hv, const Vector<Real> &v,
256  const Vector<Real> &x, Real &tol ) {
257 
258  using Teuchos::RCP; using Teuchos::dyn_cast;
259 
260  const PV &xpv = dyn_cast<const PV>(x);
261  const PV &vpv = dyn_cast<const PV>(v);
262  PV &hvpv = dyn_cast<PV>(hv);
263 
264  RCP<const V> xo = xpv.get(OPT);
265  RCP<const V> xs = xpv.get(SLACK);
266 
267  RCP<const V> vo = vpv.get(OPT);
268  RCP<const V> vs = vpv.get(SLACK);
269 
270  RCP<V> hvo = hvpv.get(OPT);
271  RCP<V> hvs = hvpv.get(SLACK);
272 
273  obj_->hessVec(*hvo, *vo, *xo, tol);
274 
275  if( hasBoundConstraint_ ) {
276  bc_barrier_->hessVec(*scratch_,*vo,*xo,tol);
277  scratch_->scale(mu_);
278  hvo->plus(*scratch_);
279 
280  }
281 
282  slack_barrier_->hessVec(*hvs, *vs, *xs, tol);
283  hvs->scale(mu_);
284 
285  }
286 
287 }; // class InteriorPointObjective
288 
289 
290 
297 template<class Real>
299 private:
300 
301  typedef Vector<Real> V;
303  typedef typename PV::size_type size_type;
304 
305  const static size_type OPT = 0;
306  const static size_type SLACK = 1;
307 
308  const static size_type INEQ = 0;
309  const static size_type EQUAL = 1;
310 
311  Teuchos::RCP<InequalityConstraint<Real> > incon_;
312  Teuchos::RCP<EqualityConstraint<Real> > eqcon_;
313 
314  bool hasEquality_; // True if an equality constraint is present
315  int ncval_; // Number of constraint evaluations
316 
317 
318 public:
319 
320  // Constructor with inequality and equality constraints
321  CompositeConstraint( const Teuchos::RCP<InequalityConstraint<Real> > &incon,
322  const Teuchos::RCP<EqualityConstraint<Real> > &eqcon ) :
323  incon_(incon), eqcon_(eqcon),
324  hasEquality_(true), ncval_(0) { }
325 
326  // Constructor with inequality constraint only
327  CompositeConstraint( const Teuchos::RCP<InequalityConstraint<Real> > &incon ) :
328  incon_(incon), eqcon_(Teuchos::null),
329  hasEquality_(false), ncval_(0) { }
330 
331 
333  return ncval_;
334  }
335 
336  void update( const Vector<Real> &x, bool flag = true, int iter = -1 ) {
337 
338  const PV &xpv = Teuchos::dyn_cast<const PV>(x);
339 
340  Teuchos::RCP<const V> xo = xpv.get(OPT);
341  Teuchos::RCP<const V> xs = xpv.get(SLACK);
342 
343  incon_->update(*xo,flag,iter);
344 
345  if( hasEquality_ ) {
346  eqcon_->update(*xo,flag,iter);
347  }
348 
349  }
350 
351  void value( Vector<Real> &c, const Vector<Real> &x, Real &tol ) {
352 
353  PV &cpv = Teuchos::dyn_cast<PV>(c);
354  const PV &xpv = Teuchos::dyn_cast<const PV>(x);
355 
356  Teuchos::RCP<const V> xo = xpv.get(OPT);
357  Teuchos::RCP<const V> xs = xpv.get(SLACK);
358 
359  Teuchos::RCP<V> ci = cpv.get(INEQ);
360  Teuchos::RCP<V> ce;
361 
362  incon_->value(*ci, *xo, tol);
363  ci->axpy(-1.0,*xs);
364 
365  if(hasEquality_) {
366  ce = cpv.get(EQUAL);
367  eqcon_->value(*ce, *xo, tol);
368  }
369 
370  ++ncval_;
371 
372  }
373 
375  const Vector<Real> &v,
376  const Vector<Real> &x,
377  Real &tol ) {
378 
379  using Teuchos::RCP; using Teuchos::dyn_cast;
380 
381  // Partition vectors and extract subvectors
382  const PV &xpv = dyn_cast<const PV>(x);
383  const PV &vpv = dyn_cast<const PV>(v);
384 
385  RCP<const V> xo = xpv.get(OPT);
386  RCP<const V> xs = xpv.get(SLACK);
387 
388  RCP<const V> vo = vpv.get(OPT);
389  RCP<const V> vs = vpv.get(SLACK);
390 
391  PV &jvpv = dyn_cast<PV>(jv);
392 
393  RCP<V> jvi = jvpv.get(INEQ);
394  incon_->applyJacobian(*jvi, *vo, *xo, tol);
395  jvi->axpy(-1.0,*vs);
396 
397  if(hasEquality_) {
398  RCP<V> jve = jvpv.get(EQUAL);
399  eqcon_->applyJacobian(*jve, *vo, *xo, tol);
400  }
401 
402  }
403 
405  const Vector<Real> &v,
406  const Vector<Real> &x,
407  Real &tol ) {
408 
409  using Teuchos::RCP; using Teuchos::dyn_cast;
410 
411  // Partition vectors and extract subvectors
412  const PV &xpv = dyn_cast<const PV>(x);
413  PV &ajvpv = dyn_cast<PV>(ajv);
414 
415  RCP<const V> xo = xpv.get(OPT);
416  RCP<const V> xs = xpv.get(SLACK);
417 
418  RCP<V> ajvo = ajvpv.get(OPT);
419  RCP<V> ajvs = ajvpv.get(SLACK);
420 
421  const PV &vpv = dyn_cast<const PV>(v);
422 
423  RCP<const V> vi = vpv.get(INEQ);
424 
425  incon_->applyAdjointJacobian(*ajvo,*vi,*xo,tol);
426 
427  ajvs->set(*vi);
428  ajvs->scale(-1.0);
429 
430  if(hasEquality_) {
431 
432  RCP<const V> ve = vpv.get(EQUAL);
433  RCP<V> temp = ajvo->clone();
434  eqcon_->applyAdjointJacobian(*temp,*ve,*xo,tol);
435  ajvo->plus(*temp);
436 
437  }
438 
439  }
440 
442  const Vector<Real> &u,
443  const Vector<Real> &v,
444  const Vector<Real> &x,
445  Real &tol ) {
446 
447  using Teuchos::RCP; using Teuchos::dyn_cast;
448 
449  const PV &xpv = dyn_cast<const PV>(x);
450  const PV &vpv = dyn_cast<const PV>(v);
451  PV &ahuvpv = dyn_cast<PV>(ahuv);
452 
453  RCP<const V> xo = xpv.get(OPT);
454  RCP<const V> xs = xpv.get(SLACK);
455 
456  RCP<const V> vo = vpv.get(OPT);
457 
458  RCP<V> ahuvo = ahuvpv.get(OPT);
459  RCP<V> ahuvs = ahuvpv.get(SLACK);
460 
461  RCP<V> temp = ahuvo->clone();
462 
463  const PV &upv = dyn_cast<const PV>(u);
464 
465  RCP<const V> ui = upv.get(INEQ);
466 
467  incon_->applyAdjointHessian(*ahuvo,*ui,*vo,*xo,tol);
468  ahuvs->zero();
469 
470  if(hasEquality_) {
471  RCP<const V> ue = upv.get(EQUAL);
472  eqcon_->applyAdjointHessian(*temp,*ue,*vo,*xo,tol);
473  ahuvo->plus(*temp);
474  }
475 
476  }
477 
478 }; // class CompositeConstraint
479 
480 } // namespace InteriorPoint
481 } // namespace ROL
482 
483 #endif
Provides the interface to evaluate objective functions.
void applyAdjointHessian(Vector< Real > &ahuv, const Vector< Real > &u, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Apply the derivative of the adjoint of the constraint Jacobian at to vector in direction ...
void gradient(Vector< Real > &g, const Vector< Real > &x, Real &tol)
Compute gradient.
CompositeConstraint(const Teuchos::RCP< InequalityConstraint< Real > > &incon, const Teuchos::RCP< EqualityConstraint< Real > > &eqcon)
void applyAdjointJacobian(Vector< Real > &ajv, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Apply the adjoint of the the constraint Jacobian at , , to vector .
Defines the linear algebra of vector space on a generic partitioned vector.
Has both inequality and equality constraints. Treat inequality constraint as equality with slack vari...
const V & dual(void) const
Return dual representation of , for example, the result of applying a Riesz map, or change of basis...
virtual Teuchos::RCP< Vector > clone() const =0
Clone to make a new (uninitialized) vector.
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:74
PenalizedObjective(const Teuchos::RCP< Objective< Real > > &obj, const Teuchos::RCP< Objective< Real > > &slack_barrier, const Vector< Real > &x, Real mu)
Real value(const Vector< Real > &x, Real &tol)
Compute value.
Defines the equality constraint operator interface.
Teuchos::RCP< const Vector< Real > > get(size_type i) const
PenalizedObjective(const Teuchos::RCP< Objective< Real > > &obj, const Teuchos::RCP< Objective< Real > > &slack_barrier, const Teuchos::RCP< Objective< Real > > &bc_barrier, const Vector< Real > &x, Real mu)
CompositeConstraint(const Teuchos::RCP< InequalityConstraint< Real > > &incon)
void update(const Vector< Real > &x, bool flag=true, int iter=-1)
Update barrier penalized objective function.
Teuchos::RCP< InequalityConstraint< Real > > incon_
Teuchos::RCP< Objective< Real > > slack_barrier_
void applyJacobian(Vector< Real > &jv, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Apply the constraint Jacobian at , , to vector .
void hessVec(Vector< Real > &hv, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Apply Hessian approximation to vector.
RCPV clone() const
Clone to make a new (uninitialized) vector.
std::vector< PV >::size_type size_type
Teuchos::RCP< Objective< Real > > obj_
Provides a unique argument for inequality constraints, which otherwise behave exactly as equality con...
virtual Real norm() const =0
Returns where .
Teuchos::RCP< Objective< Real > > bc_barrier_
void update(const Vector< Real > &x, bool flag=true, int iter=-1)
Update constraint functions. x is the optimization variable, flag = true if optimization variable is ...
void value(Vector< Real > &c, const Vector< Real > &x, Real &tol)
Evaluate the constraint operator at .
Teuchos::RCP< EqualityConstraint< Real > > eqcon_