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 "Teuchos_ParameterList.hpp"
42 
43 namespace ROL {
44 
45 
51 template <class Real>
53 
54  typedef Vector<Real> V;
55 
56  typedef Elementwise::Axpy<Real> Axpy;
57  typedef Elementwise::Aypx<Real> Aypx;
58  typedef Elementwise::Scale<Real> Scale;
59  typedef Elementwise::Reciprocal<Real> Reciprocal;
60  typedef Elementwise::Power<Real> Power;
61  typedef Elementwise::Logarithm<Real> Logarithm;
62  typedef Elementwise::Multiply<Real> Multiply;
63  typedef Elementwise::Heaviside<Real> Heaviside;
64  typedef Elementwise::ThresholdUpper<Real> ThresholdUpper;
65  typedef Elementwise::ThresholdLower<Real> ThresholdLower;
66  typedef Elementwise::ReductionSum<Real> Sum;
67  typedef Elementwise::UnaryFunction<Real> UnaryFunction;
68 
69 
70 
71  enum EBarrierType {
76  } eBarrierType_;
77 
78  inline std::string EBarrierToString( EBarrierType type ) {
79  std::string retString;
80  switch(type) {
81  case BARRIER_LOGARITHM:
82  retString = "Logarithmic";
83  break;
84  case BARRIER_QUADRATIC:
85  retString = "Quadratic";
86  break;
87  case BARRIER_DOUBLEWELL:
88  retString = "Double Well";
89  break;
90  case BARRIER_LAST:
91  retString = "Last Type (Dummy)";
92  break;
93  default:
94  retString = "Invalid EBarrierType";
95  break;
96  }
97  return retString;
98  }
99 
100  inline EBarrierType StringToEBarrierType( std::string s ) {
101  s = removeStringFormat(s);
103  for( int to = BARRIER_LOGARITHM; to != BARRIER_LAST; ++to ) {
104  type = static_cast<EBarrierType>(to);
105  if( !s.compare(removeStringFormat(EBarrierToString(type))) ) {
106  return type;
107  }
108  }
109  return type;
110  }
111 
112 
113 
114 
115 private:
116  const Teuchos::RCP<const V> lo_;
117  const Teuchos::RCP<const V> up_;
118  Teuchos::RCP<V> a_; // scratch vector
119  Teuchos::RCP<V> b_; // scratch vector
121 
122 public:
123 
125  Teuchos::ParameterList &parlist ) :
126  lo_( bc.getLowerVectorRCP() ),
127  up_( bc.getUpperVectorRCP() ) {
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.getLowerVectorRCP() ),
138  up_( bc.getUpperVectorRCP() )
139  {
140  a_ = lo_->clone();
141  b_ = up_->clone();
142  }
143 
144 
145  Real value( const Vector<Real> &x, Real &tol ) {
146 
147  Teuchos::RCP<UnaryFunction> func;
148 
149  switch(btype_) {
150  case BARRIER_LOGARITHM:
151 
152  a_->set(*lo_); // a = l
153  a_->applyBinary(Aypx(-1.0),x); // a = x-l
154  a_->applyUnary(Logarithm()); // a = log(x-l)
155  a_->applyUnary(Scale(-1.0)); // a = -log(x-l)
156 
157  b_->set(x); // b = x
158  b_->applyBinary(Aypx(-1.0),*up_); // b = u-x
159  b_->applyUnary(Logarithm()); // b = log(u-x)
160  b_->applyUnary(Scale(-1.0)); // b = -log(u-x)
161 
162  b_->plus(*a_); // b = -log(x-l)-log(u-x)
163 
164  break;
165 
166  case BARRIER_QUADRATIC:
167 
168  a_->set(*lo_); // a = l
169  a_->applyBinary(Aypx(-1.0),x); // a = x-l
170  a_->applyUnary(ThresholdLower(0.0)); // a = min(x-l,0)
171  a_->applyUnary(Power(2.0)); // a = min(x-l,0)^2
172 
173  b_->set(*up_); // b = u
174  b_->applyBinary(Aypx(-1.0),x); // b = x-u
175  b_->applyUnary(ThresholdUpper(0.0)); // b = max(x-u,0)
176  b_->applyUnary(Power(2.0)); // b = max(x-u,0)^2
177 
178  b_->plus(*a_); // b = min(x-l,0)^2 + max(x-u,0)^2
179 
180  break;
181 
182  case BARRIER_DOUBLEWELL:
183 
184  a_->set(*lo_); // a = l
185  a_->applyBinary(Aypx(-1.0),x); // a = x-l
186  a_->applyUnary(Power(2.0)); // a = (x-l)^2
187 
188  b_->set(x); // b = x
189  b_->applyBinary(Aypx(-1.0),*up_); // b = u-x
190  b_->applyUnary(Power(2.0)); // b = (u-x)^2
191 
192  b_->applyBinary(Multiply(),*a_); // b = (x-l)^2*(u-x)^2
193 
194  break;
195 
196  default:
197  TEUCHOS_TEST_FOR_EXCEPTION(true,std::invalid_argument,
198  ">>>(ObjectiveFromBoundConstraint::value): Undefined barrier function type!");
199 
200  break;
201  }
202 
203  Real result = b_->reduce(Sum());
204  return result;
205 
206  }
207 
208  void gradient( Vector<Real> &g, const Vector<Real> &x, Real &tol ) {
209 
210  switch(btype_) {
211  case BARRIER_LOGARITHM:
212 
213  a_->set(*lo_); // a = l
214  a_->applyBinary(Aypx(-1.0),x); // a = x-l
215  a_->applyUnary(Reciprocal()); // a = 1/(x-l)
216  a_->applyUnary(Scale(-1.0)); // a = -1/(x-l)
217 
218  b_->set(x); // b = x
219  b_->applyBinary(Aypx(-1.0),*up_); // b = u-x
220  b_->applyUnary(Reciprocal()); // b = 1/(u-x)
221 
222  b_->plus(*a_); // b = -1/(x-l)+1/(u-x)
223 
224  break;
225 
226  case BARRIER_QUADRATIC:
227 
228  a_->set(*lo_); // a = l
229  a_->applyBinary(Aypx(-1.0),x); // a = x-l
230  a_->applyUnary(ThresholdLower(0.0)); // a = min(x-l,0)
231 
232  b_->set(*up_); // b = u
233  b_->applyBinary(Aypx(-1.0),x); // b = x-u
234  b_->applyUnary(ThresholdUpper(0.0)); // b = max(x-u,0)
235 
236  b_->plus(*a_); // b = max(x-u,0) + min(x-l,0)
237  b_->scale(2.0); // b = 2*max(x-u,0) + 2*min(x-l,0)
238  break;
239 
240  case BARRIER_DOUBLEWELL:
241 
242  a_->set(*lo_); // a = l
243  a_->applyBinary(Aypx(-1.0),x); // a = x-l
244 
245  b_->set(x); // b = x
246  b_->applyBinary(Aypx(-1.0),*up_); // b = u-x
247 
248  a_->applyBinary(Aypx(-1.0),*b_); // a = (u-x)-(x-l)
249  a_->applyUnary(Scale(2.0)); // a = 2*[(u-x)-(x-l)]
250 
251  b_->set(x); // b = x
252  b_->applyBinary(Aypx(-1.0),*up_); // b = u-x
253 
254  a_->applyBinary(Multiply(),*b_); // a = [(u-x)-(x-l)]*(u-x);
255 
256 
257  b_->set(*lo_); // b = l
258  b_->applyBinary(Aypx(-1.0),x); // b = x-l
259  b_->applyBinary(Multiply(),*a_); // b = [(u-x)-(x-l)]*(u-x)*(x-l)
260 
261  break;
262 
263  default:
264  TEUCHOS_TEST_FOR_EXCEPTION(true,std::invalid_argument,
265  ">>>(ObjectiveFromBoundConstraint::gradient): Undefined barrier function type!");
266 
267  break;
268  }
269 
270  g.set(*b_);
271 
272  }
273 
274  void hessVec( Vector<Real> &hv, const Vector<Real> &v, const Vector<Real> &x, Real &tol ) {
275 
276  switch(btype_) {
277  case BARRIER_LOGARITHM:
278 
279  a_->set(*lo_); // a = l
280  a_->applyBinary(Axpy(-1.0),x); // a = x-l
281  a_->applyUnary(Reciprocal()); // a = 1/(x-l)
282  a_->applyUnary(Power(2.0)); // a = 1/(x-l)^2
283 
284  b_->set(x); // b = x
285  b_->applyBinary(Axpy(-1.0),*up_); // b = u-x
286  b_->applyUnary(Reciprocal()); // b = 1/(u-x)
287  b_->applyUnary(Power(2.0)); // b = 1/(u-x)^2
288 
289  b_->plus(*a_); // b = 1/(x-l)^2 + 1/(u-x)^2
290 
291  break;
292 
293  case BARRIER_QUADRATIC:
294 
295  a_->set(*lo_); // a = l
296  a_->applyBinary(Axpy(-1.0),x); // a = x-l
297  a_->scale(-1.0); // a = l-x
298  a_->applyUnary(Heaviside()); // a = theta(l-x)
299 
300  b_->set(*up_); // b = u
301  b_->applyBinary(Axpy(-1.0),x); // b = x-u
302  b_->applyUnary(Heaviside()); // b = theta(x-u)
303  b_->plus(*a_); // b = theta(l-x) + theta(x-u)
304  b_->scale(2.0); // b = 2*theta(l-x) + 2*theta(x-u)
305 
306  break;
307 
308  case BARRIER_DOUBLEWELL:
309 
310  a_->set(*lo_); // a = l
311  a_->applyBinary(Axpy(-1.0),x); // a = x-l
312 
313  b_->set(x); // b = x
314  b_->applyBinary(Axpy(-1.0),*up_); // b = u-x
315 
316  b_->applyBinary(Multiply(),*a_); // b = (u-x)*(x-l)
317  b_->applyUnary(Scale(-8.0)); // b = -8*(u-x)*(x-l)
318 
319  a_->applyUnary(Power(2.0)); // a = (x-l)^2
320  a_->applyUnary(Scale(2.0)); // a = 2*(x-l)^2
321 
322  b_->applyBinary(Axpy(1.0),*a_); // b = 2*(x-l)^2-8*(u-x)*(x-l)
323 
324  a_->set(x); // a = x
325  a_->applyBinary(Axpy(-1.0),*up_); // a = u-x
326  a_->applyUnary(Power(2.0)); // a = (u-x)^2
327  a_->applyUnary(Scale(2.0)); // a = 2*(u-x)^2
328 
329  b_->plus(*a_); // b = 2*(u-x)^2-8*(u-x)*(x-l)+2*(x-l)
330 
331 
332  break;
333 
334  default:
335  TEUCHOS_TEST_FOR_EXCEPTION(true,std::invalid_argument,
336  ">>>(ObjectiveFromBoundConstraint::hessVec): Undefined barrier function type!");
337 
338  break;
339  }
340 
341  hv.set(v);
342  hv.applyBinary(Multiply(),*b_);
343 
344  }
345 
346  // For testing purposes
347  Teuchos::RCP<Vector<Real> > getBarrierVector(void) {
348  return b_;
349  }
350 
351 
352 }; // class ObjectiveFromBoundConstraint
353 
354 }
355 
356 
357 #endif // ROL_OBJECTIVE_FROM_BOUND_CONSTRAINT_H
358 
Provides the interface to evaluate objective functions.
virtual void applyBinary(const Elementwise::BinaryFunction< Real > &f, const Vector &x)
Definition: ROL_Vector.hpp:222
std::string removeStringFormat(std::string s)
Definition: ROL_Types.hpp:174
Elementwise::ThresholdUpper< Real > ThresholdUpper
ObjectiveFromBoundConstraint(const BoundConstraint< Real > &bc)
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:74
Teuchos::RCP< Vector< Real > > getBarrierVector(void)
enum ROL::ObjectiveFromBoundConstraint::EBarrierType eBarrierType_
ObjectiveFromBoundConstraint(const BoundConstraint< Real > &bc, Teuchos::ParameterList &parlist)
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:196
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.