ROL
ROL_ObjectiveFromBoundConstraint.hpp
Go to the documentation of this file.
1 // Redistribution and use in source and binary forms, with or without
2 // modification, are permitted provided that the following conditions are
3 // met:
4 //
5 // 1. Redistributions of source code must retain the above copyright
6 // notice, this list of conditions and the following disclaimer.
7 //
8 // 2. Redistributions in binary form must reproduce the above copyright
9 // notice, this list of conditions and the following disclaimer in the
10 // documentation and/or other materials provided with the distribution.
11 //
12 // 3. Neither the name of the Corporation nor the names of the
13 // contributors may be used to endorse or promote products derived from
14 // this software without specific prior written permission.
15 //
16 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
17 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
18 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
19 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
20 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
21 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
22 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
23 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
24 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
25 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
26 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
27 //
28 // Questions? Contact lead developers:
29 // Drew Kouri (dpkouri@sandia.gov) and
30 // Denis Ridzal (dridzal@sandia.gov)
31 //
32 // ************************************************************************
33 // @HEADER
34 
35 #ifndef ROL_OBJECTIVE_FROM_BOUND_CONSTRAINT_H
36 #define ROL_OBJECTIVE_FROM_BOUND_CONSTRAINT_H
37 
38 #include "ROL_Objective.hpp"
39 #include "ROL_BoundConstraint.hpp"
40 
41 #include "ROL_ParameterList.hpp"
42 
43 namespace ROL {
44 
45 
51 template <class Real>
53 
54  typedef Vector<Real> V;
55 
56  typedef Elementwise::Fill<Real> Fill;
57  typedef Elementwise::Reciprocal<Real> Reciprocal;
58  typedef Elementwise::Power<Real> Power;
59  typedef Elementwise::Logarithm<Real> Logarithm;
60  typedef Elementwise::Multiply<Real> Multiply;
61  typedef Elementwise::Heaviside<Real> Heaviside;
62  typedef Elementwise::ThresholdUpper<Real> ThresholdUpper;
63  typedef Elementwise::ThresholdLower<Real> ThresholdLower;
64  typedef Elementwise::ReductionSum<Real> Sum;
65  typedef Elementwise::UnaryFunction<Real> UnaryFunction;
66 
67 
68 
69  enum EBarrierType {
74  } eBarrierType_;
75 
76  inline std::string EBarrierToString( EBarrierType type ) {
77  std::string retString;
78  switch(type) {
79  case BARRIER_LOGARITHM:
80  retString = "Logarithmic";
81  break;
82  case BARRIER_QUADRATIC:
83  retString = "Quadratic";
84  break;
85  case BARRIER_DOUBLEWELL:
86  retString = "Double Well";
87  break;
88  case BARRIER_LAST:
89  retString = "Last Type (Dummy)";
90  break;
91  default:
92  retString = "Invalid EBarrierType";
93  break;
94  }
95  return retString;
96  }
97 
98  inline EBarrierType StringToEBarrierType( std::string s ) {
99  s = removeStringFormat(s);
101  for( int to = BARRIER_LOGARITHM; to != BARRIER_LAST; ++to ) {
102  type = static_cast<EBarrierType>(to);
103  if( !s.compare(removeStringFormat(EBarrierToString(type))) ) {
104  return type;
105  }
106  }
107  return type;
108  }
109 
110 private:
111  const ROL::Ptr<const V> lo_;
112  const ROL::Ptr<const V> up_;
113  ROL::Ptr<V> a_; // scratch vector
114  ROL::Ptr<V> b_; // scratch vector
118 
119 public:
120 
122  ROL::ParameterList &parlist ) :
123  lo_( bc.getLowerBound() ),
124  up_( bc.getUpperBound() ) {
125 
128 
129  a_ = lo_->clone();
130  b_ = up_->clone();
131 
132  std::string bfstring = parlist.sublist("Barrier Function").get("Type","Logarithmic");
133  btype_ = StringToEBarrierType(bfstring);
134  }
135 
137  lo_( bc.getLowerBound() ),
138  up_( bc.getUpperBound() ),
140 
143 
144  a_ = lo_->clone();
145  b_ = up_->clone();
146  }
147 
148 
149  Real value( const Vector<Real> &x, Real &tol ) {
150  const Real zero(0), one(1), two(2);
151 
152  ROL::Ptr<UnaryFunction> func;
153 
154  a_->zero(); b_->zero();
155  switch(btype_) {
156  case BARRIER_LOGARITHM:
157 
158  if ( isLowerActivated_ ) {
159  a_->set(x); // a = x
160  a_->axpy(-one,*lo_); // a = x-l
161  a_->applyUnary(Logarithm()); // a = log(x-l)
162  }
163 
164  if ( isUpperActivated_ ) {
165  b_->set(*up_); // b = u
166  b_->axpy(-one,x); // b = u-x
167  b_->applyUnary(Logarithm()); // b = log(u-x)
168  }
169 
170  b_->plus(*a_); // b = log(x-l)+log(u-x)
171  b_->scale(-one); // b = -log(x-l)-log(u-x)
172 
173  break;
174 
175  case BARRIER_QUADRATIC:
176 
177  if ( isLowerActivated_ ) {
178  a_->set(x); // a = x
179  a_->axpy(-one,*lo_); // a = x-l
180  a_->applyUnary(ThresholdLower(zero)); // a = min(x-l,0)
181  a_->applyUnary(Power(two)); // a = min(x-l,0)^2
182  }
183 
184  if ( isUpperActivated_ ) {
185  b_->set(*up_); // b = u
186  b_->axpy(-one,x); // b = u-x
187  b_->applyUnary(ThresholdUpper(zero)); // b = max(x-u,0)
188  b_->applyUnary(Power(two)); // b = max(x-u,0)^2
189  }
190 
191  b_->plus(*a_); // b = min(x-l,0)^2 + max(x-u,0)^2
192 
193  break;
194 
195  case BARRIER_DOUBLEWELL:
196 
197  if ( isLowerActivated_ ) {
198  a_->set(x); // a = x
199  a_->axpy(-one,*lo_); // a = x-l
200  a_->applyUnary(Power(two)); // a = (x-l)^2
201  }
202  else {
203  a_->applyUnary(Fill(one));
204  }
205 
206  if ( isUpperActivated_ ) {
207  b_->set(*up_); // b = u
208  b_->axpy(-one,x); // b = u-x
209  b_->applyUnary(Power(two)); // b = (u-x)^2
210  }
211  else {
212  b_->applyUnary(Fill(one));
213  }
214 
215  b_->applyBinary(Multiply(),*a_); // b = (x-l)^2*(u-x)^2
216 
217  break;
218 
219  default:
220  ROL_TEST_FOR_EXCEPTION(true,std::invalid_argument,
221  ">>>(ObjectiveFromBoundConstraint::value): Undefined barrier function type!");
222 
223  break;
224  }
225 
226  Real result = b_->reduce(Sum());
227  return result;
228 
229  }
230 
231  void gradient( Vector<Real> &g, const Vector<Real> &x, Real &tol ) {
232  const Real zero(0), one(1), two(2);
233 
234  a_->zero(); b_->zero();
235  switch(btype_) {
236  case BARRIER_LOGARITHM:
237 
238  if ( isLowerActivated_ ) {
239  a_->set(*lo_); // a = l
240  a_->axpy(-one,x); // a = l-x
241  a_->applyUnary(Reciprocal()); // a = -1/(x-l)
242  }
243 
244  if ( isUpperActivated_ ) {
245  b_->set(*up_); // b = u
246  b_->axpy(-one,x); // b = u-x
247  b_->applyUnary(Reciprocal()); // b = 1/(u-x)
248  }
249 
250  b_->plus(*a_); // b = -1/(x-l)+1/(u-x)
251 
252  break;
253 
254  case BARRIER_QUADRATIC:
255 
256  if ( isLowerActivated_ ) {
257  a_->set(x); // a = x
258  a_->axpy(-one,*lo_); // a = x-l
259  a_->applyUnary(ThresholdLower(zero)); // a = min(x-l,0)
260  }
261 
262  if ( isUpperActivated_ ) {
263  b_->set(*up_); // b = u
264  b_->axpy(-one,x); // b = u-x
265  b_->applyUnary(ThresholdUpper(zero)); // b = max(x-u,0)
266  }
267 
268  b_->plus(*a_); // b = max(x-u,0) + min(x-l,0)
269  b_->scale(two); // b = 2*max(x-u,0) + 2*min(x-l,0)
270  break;
271 
272  case BARRIER_DOUBLEWELL:
273 
274  if ( isLowerActivated_ ) {
275  a_->set(x); // a = l
276  a_->axpy(-one,*lo_); // a = x-l
277  }
278  else {
279  a_->applyUnary(Fill(one));
280  }
281 
282  if ( isUpperActivated_ ) {
283  b_->set(*up_); // b = x
284  b_->axpy(-one,x); // b = u-x
285  }
286  else {
287  b_->applyUnary(Fill(one));
288  }
289 
290  b_->applyBinary(Multiply(),*a_); // b = (x-l)*(u-x)
291  b_->scale(two); // b = 2*(x-l)*(u-x)
292 
294  a_->set(*up_); // a = u
295  a_->axpy(-two,x); // a = (u-x)-x
296  a_->plus(*lo_); // a = (u-x)-(x-l)
297  b_->applyBinary(Multiply(),*b_); // b = 2*(x-l)*(u-x)*[(u-x)-(x-l)]
298  }
299 
300  break;
301 
302  default:
303  ROL_TEST_FOR_EXCEPTION(true,std::invalid_argument,
304  ">>>(ObjectiveFromBoundConstraint::gradient): Undefined barrier function type!");
305 
306  break;
307  }
308 
309  g.set(*b_);
310 
311  }
312 
313  void hessVec( Vector<Real> &hv, const Vector<Real> &v, const Vector<Real> &x, Real &tol ) {
314  const Real one(1), two(2), eight(8);
315 
316  switch(btype_) {
317  case BARRIER_LOGARITHM:
318 
319  if ( isLowerActivated_ ) {
320  a_->set(x); // a = x
321  a_->axpy(-one,*lo_); // a = x-l
322  a_->applyUnary(Reciprocal()); // a = 1/(x-l)
323  a_->applyUnary(Power(two)); // a = 1/(x-l)^2
324  }
325 
326  if ( isUpperActivated_ ) {
327  b_->set(*up_); // b = u
328  b_->axpy(-one,x); // b = u-x
329  b_->applyUnary(Reciprocal()); // b = 1/(u-x)
330  b_->applyUnary(Power(two)); // b = 1/(u-x)^2
331  }
332 
333  b_->plus(*a_); // b = 1/(x-l)^2 + 1/(u-x)^2
334 
335  break;
336 
337  case BARRIER_QUADRATIC:
338 
339  if ( isLowerActivated_ ) {
340  a_->set(*lo_); // a = l
341  a_->axpy(-one,x); // a = l-x
342  a_->applyUnary(Heaviside()); // a = theta(l-x)
343  }
344 
345  if ( isUpperActivated_ ) {
346  b_->set(x); // b = x
347  b_->axpy(-one,*up_); // b = x-u
348  b_->applyUnary(Heaviside()); // b = theta(x-u)
349  }
350 
351  b_->plus(*a_); // b = theta(l-x) + theta(x-u)
352  b_->scale(two); // b = 2*theta(l-x) + 2*theta(x-u)
353 
354  break;
355 
356  case BARRIER_DOUBLEWELL:
357 
359  a_->set(x); // a = x
360  a_->axpy(-one,*lo_); // a = x-l
361 
362  b_->set(*up_); // b = u
363  b_->axpy(-one,x); // b = u-x
364 
365  b_->applyBinary(Multiply(),*a_); // b = (u-x)*(x-l)
366  b_->scale(-eight); // b = -8*(u-x)*(x-l)
367 
368  a_->applyUnary(Power(two)); // a = (x-l)^2
369  a_->scale(two); // a = 2*(x-l)^2
370 
371  b_->plus(*a_); // b = 2*(x-l)^2-8*(u-x)*(x-l)
372 
373  a_->set(*up_); // a = u
374  a_->axpy(-one,x); // a = u-x
375  a_->applyUnary(Power(two)); // a = (u-x)^2
376  a_->scale(two); // a = 2*(u-x)^2
377 
378  b_->plus(*a_); // b = 2*(u-x)^2-8*(u-x)*(x-l)+2*(x-l)^2
379  }
380  else {
381  b_->applyUnary(Fill(two));
382  }
383 
384  break;
385 
386  default:
387  ROL_TEST_FOR_EXCEPTION(true,std::invalid_argument,
388  ">>>(ObjectiveFromBoundConstraint::hessVec): Undefined barrier function type!");
389 
390  break;
391  }
392 
393  hv.set(v);
394  hv.applyBinary(Multiply(),*b_);
395 
396  }
397 
398  // For testing purposes
399  ROL::Ptr<Vector<Real> > getBarrierVector(void) {
400  return b_;
401  }
402 
403 
404 }; // class ObjectiveFromBoundConstraint
405 
406 }
407 
408 
409 #endif // ROL_OBJECTIVE_FROM_BOUND_CONSTRAINT_H
410 
Provides the interface to evaluate objective functions.
virtual void applyBinary(const Elementwise::BinaryFunction< Real > &f, const Vector &x)
Definition: ROL_Vector.hpp:236
std::string removeStringFormat(std::string s)
Definition: ROL_Types.hpp:247
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:80
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:209
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.