ROL
ROL_ObjectiveFromBoundConstraint.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_OBJECTIVE_FROM_BOUND_CONSTRAINT_H
11 #define ROL_OBJECTIVE_FROM_BOUND_CONSTRAINT_H
12 
13 #include "ROL_Objective.hpp"
14 #include "ROL_BoundConstraint.hpp"
15 
16 #include "ROL_ParameterList.hpp"
17 
18 namespace ROL {
19 
20 
26 template <class Real>
28 
29  typedef Vector<Real> V;
30 
31  typedef Elementwise::Fill<Real> Fill;
32  typedef Elementwise::Reciprocal<Real> Reciprocal;
33  typedef Elementwise::Power<Real> Power;
34  typedef Elementwise::Logarithm<Real> Logarithm;
35  typedef Elementwise::Multiply<Real> Multiply;
36  typedef Elementwise::Heaviside<Real> Heaviside;
37  typedef Elementwise::ThresholdUpper<Real> ThresholdUpper;
38  typedef Elementwise::ThresholdLower<Real> ThresholdLower;
39  typedef Elementwise::ReductionSum<Real> Sum;
40  typedef Elementwise::UnaryFunction<Real> UnaryFunction;
41 
42 
43 
44  enum EBarrierType {
49  } eBarrierType_;
50 
51  inline std::string EBarrierToString( EBarrierType type ) {
52  std::string retString;
53  switch(type) {
54  case BARRIER_LOGARITHM:
55  retString = "Logarithmic";
56  break;
57  case BARRIER_QUADRATIC:
58  retString = "Quadratic";
59  break;
60  case BARRIER_DOUBLEWELL:
61  retString = "Double Well";
62  break;
63  case BARRIER_LAST:
64  retString = "Last Type (Dummy)";
65  break;
66  default:
67  retString = "Invalid EBarrierType";
68  break;
69  }
70  return retString;
71  }
72 
73  inline EBarrierType StringToEBarrierType( std::string s ) {
74  s = removeStringFormat(s);
76  for( int to = BARRIER_LOGARITHM; to != BARRIER_LAST; ++to ) {
77  type = static_cast<EBarrierType>(to);
78  if( !s.compare(removeStringFormat(EBarrierToString(type))) ) {
79  return type;
80  }
81  }
82  return type;
83  }
84 
85 private:
86  const ROL::Ptr<const V> lo_;
87  const ROL::Ptr<const V> up_;
88  ROL::Ptr<V> a_; // scratch vector
89  ROL::Ptr<V> b_; // scratch vector
93 
94 public:
95 
97  ROL::ParameterList &parlist ) :
98  lo_( bc.getLowerBound() ),
99  up_( bc.getUpperBound() ) {
100 
103 
104  a_ = lo_->clone();
105  b_ = up_->clone();
106 
107  std::string bfstring = parlist.sublist("Barrier Function").get("Type","Logarithmic");
108  btype_ = StringToEBarrierType(bfstring);
109  }
110 
112  lo_( bc.getLowerBound() ),
113  up_( bc.getUpperBound() ),
115 
118 
119  a_ = lo_->clone();
120  b_ = up_->clone();
121  }
122 
123 
124  Real value( const Vector<Real> &x, Real &tol ) {
125  const Real zero(0), one(1), two(2);
126 
127  ROL::Ptr<UnaryFunction> func;
128 
129  a_->zero(); b_->zero();
130  switch(btype_) {
131  case BARRIER_LOGARITHM:
132 
133  if ( isLowerActivated_ ) {
134  a_->set(x); // a = x
135  a_->axpy(-one,*lo_); // a = x-l
136  a_->applyUnary(Logarithm()); // a = log(x-l)
137  }
138 
139  if ( isUpperActivated_ ) {
140  b_->set(*up_); // b = u
141  b_->axpy(-one,x); // b = u-x
142  b_->applyUnary(Logarithm()); // b = log(u-x)
143  }
144 
145  b_->plus(*a_); // b = log(x-l)+log(u-x)
146  b_->scale(-one); // b = -log(x-l)-log(u-x)
147 
148  break;
149 
150  case BARRIER_QUADRATIC:
151 
152  if ( isLowerActivated_ ) {
153  a_->set(x); // a = x
154  a_->axpy(-one,*lo_); // a = x-l
155  a_->applyUnary(ThresholdLower(zero)); // a = min(x-l,0)
156  a_->applyUnary(Power(two)); // a = min(x-l,0)^2
157  }
158 
159  if ( isUpperActivated_ ) {
160  b_->set(*up_); // b = u
161  b_->axpy(-one,x); // b = u-x
162  b_->applyUnary(ThresholdUpper(zero)); // b = max(x-u,0)
163  b_->applyUnary(Power(two)); // b = max(x-u,0)^2
164  }
165 
166  b_->plus(*a_); // b = min(x-l,0)^2 + max(x-u,0)^2
167 
168  break;
169 
170  case BARRIER_DOUBLEWELL:
171 
172  if ( isLowerActivated_ ) {
173  a_->set(x); // a = x
174  a_->axpy(-one,*lo_); // a = x-l
175  a_->applyUnary(Power(two)); // a = (x-l)^2
176  }
177  else {
178  a_->applyUnary(Fill(one));
179  }
180 
181  if ( isUpperActivated_ ) {
182  b_->set(*up_); // b = u
183  b_->axpy(-one,x); // b = u-x
184  b_->applyUnary(Power(two)); // b = (u-x)^2
185  }
186  else {
187  b_->applyUnary(Fill(one));
188  }
189 
190  b_->applyBinary(Multiply(),*a_); // b = (x-l)^2*(u-x)^2
191 
192  break;
193 
194  default:
195 
196  ROL_TEST_FOR_EXCEPTION(true,std::invalid_argument,
197  ">>>(ObjectiveFromBoundConstraint::value): Undefined barrier function type!");
198 
199  }
200 
201  Real result = b_->reduce(Sum());
202  return result;
203 
204  }
205 
206  void gradient( Vector<Real> &g, const Vector<Real> &x, Real &tol ) {
207  const Real zero(0), one(1), two(2);
208 
209  a_->zero(); b_->zero();
210  switch(btype_) {
211  case BARRIER_LOGARITHM:
212 
213  if ( isLowerActivated_ ) {
214  a_->set(*lo_); // a = l
215  a_->axpy(-one,x); // a = l-x
216  a_->applyUnary(Reciprocal()); // a = -1/(x-l)
217  }
218 
219  if ( isUpperActivated_ ) {
220  b_->set(*up_); // b = u
221  b_->axpy(-one,x); // b = u-x
222  b_->applyUnary(Reciprocal()); // b = 1/(u-x)
223  }
224 
225  b_->plus(*a_); // b = -1/(x-l)+1/(u-x)
226 
227  break;
228 
229  case BARRIER_QUADRATIC:
230 
231  if ( isLowerActivated_ ) {
232  a_->set(x); // a = x
233  a_->axpy(-one,*lo_); // a = x-l
234  a_->applyUnary(ThresholdLower(zero)); // a = min(x-l,0)
235  }
236 
237  if ( isUpperActivated_ ) {
238  b_->set(*up_); // b = u
239  b_->axpy(-one,x); // b = u-x
240  b_->applyUnary(ThresholdUpper(zero)); // b = max(x-u,0)
241  }
242 
243  b_->plus(*a_); // b = max(x-u,0) + min(x-l,0)
244  b_->scale(two); // b = 2*max(x-u,0) + 2*min(x-l,0)
245  break;
246 
247  case BARRIER_DOUBLEWELL:
248 
249  if ( isLowerActivated_ ) {
250  a_->set(x); // a = l
251  a_->axpy(-one,*lo_); // a = x-l
252  }
253  else {
254  a_->applyUnary(Fill(one));
255  }
256 
257  if ( isUpperActivated_ ) {
258  b_->set(*up_); // b = x
259  b_->axpy(-one,x); // b = u-x
260  }
261  else {
262  b_->applyUnary(Fill(one));
263  }
264 
265  b_->applyBinary(Multiply(),*a_); // b = (x-l)*(u-x)
266  b_->scale(two); // b = 2*(x-l)*(u-x)
267 
269  a_->set(*up_); // a = u
270  a_->axpy(-two,x); // a = (u-x)-x
271  a_->plus(*lo_); // a = (u-x)-(x-l)
272  b_->applyBinary(Multiply(),*b_); // b = 2*(x-l)*(u-x)*[(u-x)-(x-l)]
273  }
274 
275  break;
276 
277  default:
278 
279  ROL_TEST_FOR_EXCEPTION(true,std::invalid_argument,
280  ">>>(ObjectiveFromBoundConstraint::gradient): Undefined barrier function type!");
281 
282  }
283 
284  g.set(*b_);
285 
286  }
287 
288  void hessVec( Vector<Real> &hv, const Vector<Real> &v, const Vector<Real> &x, Real &tol ) {
289  const Real one(1), two(2), eight(8);
290 
291  switch(btype_) {
292  case BARRIER_LOGARITHM:
293 
294  if ( isLowerActivated_ ) {
295  a_->set(x); // a = x
296  a_->axpy(-one,*lo_); // a = x-l
297  a_->applyUnary(Reciprocal()); // a = 1/(x-l)
298  a_->applyUnary(Power(two)); // a = 1/(x-l)^2
299  }
300 
301  if ( isUpperActivated_ ) {
302  b_->set(*up_); // b = u
303  b_->axpy(-one,x); // b = u-x
304  b_->applyUnary(Reciprocal()); // b = 1/(u-x)
305  b_->applyUnary(Power(two)); // b = 1/(u-x)^2
306  }
307 
308  b_->plus(*a_); // b = 1/(x-l)^2 + 1/(u-x)^2
309 
310  break;
311 
312  case BARRIER_QUADRATIC:
313 
314  if ( isLowerActivated_ ) {
315  a_->set(*lo_); // a = l
316  a_->axpy(-one,x); // a = l-x
317  a_->applyUnary(Heaviside()); // a = theta(l-x)
318  }
319 
320  if ( isUpperActivated_ ) {
321  b_->set(x); // b = x
322  b_->axpy(-one,*up_); // b = x-u
323  b_->applyUnary(Heaviside()); // b = theta(x-u)
324  }
325 
326  b_->plus(*a_); // b = theta(l-x) + theta(x-u)
327  b_->scale(two); // b = 2*theta(l-x) + 2*theta(x-u)
328 
329  break;
330 
331  case BARRIER_DOUBLEWELL:
332 
334  a_->set(x); // a = x
335  a_->axpy(-one,*lo_); // a = x-l
336 
337  b_->set(*up_); // b = u
338  b_->axpy(-one,x); // b = u-x
339 
340  b_->applyBinary(Multiply(),*a_); // b = (u-x)*(x-l)
341  b_->scale(-eight); // b = -8*(u-x)*(x-l)
342 
343  a_->applyUnary(Power(two)); // a = (x-l)^2
344  a_->scale(two); // a = 2*(x-l)^2
345 
346  b_->plus(*a_); // b = 2*(x-l)^2-8*(u-x)*(x-l)
347 
348  a_->set(*up_); // a = u
349  a_->axpy(-one,x); // a = u-x
350  a_->applyUnary(Power(two)); // a = (u-x)^2
351  a_->scale(two); // a = 2*(u-x)^2
352 
353  b_->plus(*a_); // b = 2*(u-x)^2-8*(u-x)*(x-l)+2*(x-l)^2
354  }
355  else {
356  b_->applyUnary(Fill(two));
357  }
358 
359  break;
360 
361  default:
362 
363  ROL_TEST_FOR_EXCEPTION(true,std::invalid_argument,
364  ">>>(ObjectiveFromBoundConstraint::hessVec): Undefined barrier function type!");
365 
366  }
367 
368  hv.set(v);
369  hv.applyBinary(Multiply(),*b_);
370 
371  }
372 
373  // For testing purposes
374  ROL::Ptr<Vector<Real> > getBarrierVector(void) {
375  return b_;
376  }
377 
378 
379 }; // class ObjectiveFromBoundConstraint
380 
381 }
382 
383 
384 #endif // ROL_OBJECTIVE_FROM_BOUND_CONSTRAINT_H
385 
Provides the interface to evaluate objective functions.
virtual void applyBinary(const Elementwise::BinaryFunction< Real > &f, const Vector &x)
Definition: ROL_Vector.hpp:214
std::string removeStringFormat(std::string s)
Definition: ROL_Types.hpp:215
Elementwise::ThresholdUpper< Real > ThresholdUpper
ROL::Ptr< Vector< Real > > getBarrierVector(void)
ObjectiveFromBoundConstraint(const BoundConstraint< Real > &bc)
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:46
ObjectiveFromBoundConstraint(const BoundConstraint< Real > &bc, ROL::ParameterList &parlist)
Objective_SerialSimOpt(const Ptr< Obj > &obj, const V &ui) z0_ zero()
enum ROL::ObjectiveFromBoundConstraint::EBarrierType eBarrierType_
Elementwise::UnaryFunction< Real > UnaryFunction
Provides the interface to apply upper and lower bound constraints.
Real value(const Vector< Real > &x, Real &tol)
Compute value.
void gradient(Vector< Real > &g, const Vector< Real > &x, Real &tol)
Compute gradient.
virtual void set(const Vector &x)
Set where .
Definition: ROL_Vector.hpp:175
bool isLowerActivated(void) const
Check if lower bound are on.
void hessVec(Vector< Real > &hv, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Apply Hessian approximation to vector.
Elementwise::ThresholdLower< Real > ThresholdLower
Create a penalty objective from upper and lower bound vectors.
bool isUpperActivated(void) const
Check if upper bound are on.