ROL
ROL_ParaboloidCircle.hpp
Go to the documentation of this file.
1 // @HEADER
2 // ************************************************************************
3 //
4 // Rapid Optimization Library (ROL) Package
5 // Copyright (2014) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact lead developers:
38 // Drew Kouri (dpkouri@sandia.gov) and
39 // Denis Ridzal (dridzal@sandia.gov)
40 //
41 // ************************************************************************
42 // @HEADER
43 
53 #ifndef ROL_PARABOLOIDCIRCLE_HPP
54 #define ROL_PARABOLOIDCIRCLE_HPP
55 
56 #include "ROL_Objective.hpp"
57 #include "ROL_StdVector.hpp"
59 #include "Teuchos_SerialDenseVector.hpp"
60 #include "Teuchos_SerialDenseSolver.hpp"
61 
62 namespace ROL {
63 namespace ZOO {
64 
68  template< class Real, class XPrim=StdVector<Real>, class XDual=StdVector<Real> >
69  class Objective_ParaboloidCircle : public Objective<Real> {
70 
71  typedef std::vector<Real> vector;
72  typedef Vector<Real> V;
73 
74  typedef typename vector::size_type uint;
75 
76 
77  private:
78 
79  template<class VectorType>
80  Teuchos::RCP<const vector> getVector( const V& x ) {
81  using Teuchos::dyn_cast;
82  return dyn_cast<const VectorType>(x).getVector();
83  }
84 
85  template<class VectorType>
86  Teuchos::RCP<vector> getVector( V& x ) {
87  using Teuchos::dyn_cast;
88  return dyn_cast<VectorType>(x).getVector();
89  }
90 
91  public:
93 
94  Real value( const Vector<Real> &x, Real &tol ) {
95 
96  using Teuchos::RCP;
97  RCP<const vector> xp = getVector<XPrim>(x);
98 
99  uint n = xp->size();
100  TEUCHOS_TEST_FOR_EXCEPTION( (n != 2), std::invalid_argument, ">>> ERROR (ROL_ParaboloidCircle, objective value): "
101  "Primal vector x must be of length 2.");
102 
103  Real x1 = (*xp)[0];
104  Real x2 = (*xp)[1];
105 
106  Real val = x1*x1 + x2*x2;
107 
108  return val;
109  }
110 
111  void gradient( Vector<Real> &g, const Vector<Real> &x, Real &tol ) {
112 
113  using Teuchos::RCP;
114  RCP<const vector> xp = getVector<XPrim>(x);
115  RCP<vector> gp = getVector<XDual>(g);
116 
117  uint n = xp->size();
118  TEUCHOS_TEST_FOR_EXCEPTION( (n != 2), std::invalid_argument, ">>> ERROR (ROL_ParaboloidCircle, objective gradient): "
119  " Primal vector x must be of length 2.");
120 
121  n = gp->size();
122  TEUCHOS_TEST_FOR_EXCEPTION( (n != 2), std::invalid_argument, ">>> ERROR (ROL_ParaboloidCircle, objective gradient): "
123  "Gradient vector g must be of length 2.");
124 
125  Real x1 = (*xp)[0];
126  Real x2 = (*xp)[1];
127 
128  Real two(2);
129 
130  (*gp)[0] = two*x1;
131  (*gp)[1] = two*x2;
132  }
133 
134  void hessVec( Vector<Real> &hv, const Vector<Real> &v, const Vector<Real> &x, Real &tol ) {
135 
136  using Teuchos::RCP;
137  RCP<const vector> xp = getVector<XPrim>(x);
138  RCP<const vector> vp = getVector<XPrim>(v);
139  RCP<vector> hvp = getVector<XDual>(hv);
140 
141  uint n = xp->size();
142  TEUCHOS_TEST_FOR_EXCEPTION( (n != 2), std::invalid_argument, ">>> ERROR (ROL_ParaboloidCircle, objective hessVec): "
143  "Primal vector x must be of length 2.");
144 
145  n = vp->size();
146  TEUCHOS_TEST_FOR_EXCEPTION( (n != 2), std::invalid_argument, ">>> ERROR (ROL_ParaboloidCircle, objective hessVec): "
147  "Input vector v must be of length 2.");
148 
149  n = hvp->size();
150  TEUCHOS_TEST_FOR_EXCEPTION( (n != 2), std::invalid_argument, ">>> ERROR (ROL_ParaboloidCircle, objective hessVec): "
151  "Output vector hv must be of length 2.");
152 
153  Real v1 = (*vp)[0];
154  Real v2 = (*vp)[1];
155 
156  Real two(2);
157 
158  (*hvp)[0] = two*v1;
159  (*hvp)[1] = two*v2;
160  }
161 
162  };
163 
164 
167  template<class Real, class XPrim=StdVector<Real>, class XDual=StdVector<Real>, class CPrim=StdVector<Real>, class CDual=StdVector<Real> >
169 
170  typedef std::vector<Real> vector;
171  typedef Vector<Real> V;
172 
173  typedef typename vector::size_type uint;
174 
175  private:
176  template<class VectorType>
177  Teuchos::RCP<const vector> getVector( const V& x ) {
178  using Teuchos::dyn_cast;
179  return dyn_cast<const VectorType>(x).getVector();
180  }
181 
182  template<class VectorType>
183  Teuchos::RCP<vector> getVector( V& x ) {
184  using Teuchos::dyn_cast;
185  return dyn_cast<VectorType>(x).getVector();
186  }
187 
188  public:
190 
191  void value( Vector<Real> &c, const Vector<Real> &x, Real &tol ) {
192 
193  using Teuchos::RCP;
194  RCP<const vector> xp = getVector<XPrim>(x);
195  RCP<vector> cp = getVector<CPrim>(c);
196 
197  uint n = xp->size();
198  TEUCHOS_TEST_FOR_EXCEPTION( (n != 2), std::invalid_argument, ">>> ERROR (ROL_ParaboloidCircle, constraint value): "
199  "Primal vector x must be of length 2.");
200 
201  uint m = cp->size();
202  TEUCHOS_TEST_FOR_EXCEPTION( (m != 1), std::invalid_argument, ">>> ERROR (ROL_ParaboloidCircle, constraint value): "
203  "Constraint vector c must be of length 1.");
204 
205  Real x1 = (*xp)[0];
206  Real x2 = (*xp)[1];
207 
208  Real one(1), two(2);
209 
210  (*cp)[0] = (x1-two)*(x1-two) + x2*x2 - one;
211  }
212 
213  void applyJacobian( Vector<Real> &jv, const Vector<Real> &v, const Vector<Real> &x, Real &tol ) {
214 
215  using Teuchos::RCP;
216  RCP<const vector> xp = getVector<XPrim>(x);
217  RCP<const vector> vp = getVector<XPrim>(v);
218  RCP<vector> jvp = getVector<CPrim>(jv);
219 
220  uint n = xp->size();
221  TEUCHOS_TEST_FOR_EXCEPTION( (n != 2), std::invalid_argument, ">>> ERROR (ROL_ParaboloidCircle, constraint applyJacobian): "
222  "Primal vector x must be of length 2.");
223 
224  uint d = vp->size();
225  TEUCHOS_TEST_FOR_EXCEPTION( (d != 2), std::invalid_argument, ">>> ERROR (ROL_ParaboloidCircle, constraint applyJacobian): "
226  "Input vector v must be of length 2.");
227  d = jvp->size();
228  TEUCHOS_TEST_FOR_EXCEPTION( (d != 1), std::invalid_argument, ">>> ERROR (ROL_ParaboloidCircle, constraint applyJacobian): "
229  "Output vector jv must be of length 1.");
230 
231  Real x1 = (*xp)[0];
232  Real x2 = (*xp)[1];
233 
234  Real v1 = (*vp)[0];
235  Real v2 = (*vp)[1];
236 
237  Real two(2);
238 
239  (*jvp)[0] = two*(x1-two)*v1 + two*x2*v2;
240  } //applyJacobian
241 
242  void applyAdjointJacobian( Vector<Real> &ajv, const Vector<Real> &v, const Vector<Real> &x, Real &tol ) {
243 
244  using Teuchos::RCP;
245  RCP<const vector> xp = getVector<XPrim>(x);
246  RCP<const vector> vp = getVector<CDual>(v);
247  RCP<vector> ajvp = getVector<XDual>(ajv);
248 
249  uint n = xp->size();
250  TEUCHOS_TEST_FOR_EXCEPTION( (n != 2), std::invalid_argument, ">>> ERROR (ROL_ParaboloidCircle, constraint applyAdjointJacobian): "
251  "Primal vector x must be of length 2.");
252 
253  uint d = vp->size();
254  TEUCHOS_TEST_FOR_EXCEPTION( (d != 1), std::invalid_argument, ">>> ERROR (ROL_ParaboloidCircle, constraint applyAdjointJacobian): "
255  "Input vector v must be of length 1.");
256 
257  d = ajvp->size();
258  TEUCHOS_TEST_FOR_EXCEPTION( (d != 2), std::invalid_argument, ">>> ERROR (ROL_ParaboloidCircle, constraint applyAdjointJacobian): "
259  "Output vector ajv must be of length 2.");
260 
261  Real x1 = (*xp)[0];
262  Real x2 = (*xp)[1];
263 
264  Real v1 = (*vp)[0];
265 
266  Real two(2);
267 
268  (*ajvp)[0] = two*(x1-two)*v1;
269  (*ajvp)[1] = two*x2*v1;
270 
271  } //applyAdjointJacobian
272 
273  void applyAdjointHessian( Vector<Real> &ahuv, const Vector<Real> &u, const Vector<Real> &v, const Vector<Real> &x, Real &tol ) {
274 
275  bool useFD = true;
276 
277  if (useFD) {
278  EqualityConstraint<Real>::applyAdjointHessian( ahuv, u, v, x, tol );
279  }
280  else {
281  using Teuchos::RCP;
282  RCP<const vector> xp = getVector<XPrim>(x);
283  RCP<const vector> up = getVector<CDual>(u);
284  RCP<const vector> vp = getVector<XPrim>(v);
285  RCP<vector> ahuvp = getVector<XDual>(ahuv);
286 
287  uint n = xp->size();
288  TEUCHOS_TEST_FOR_EXCEPTION( (n != 2), std::invalid_argument, ">>> ERROR (ROL_ParaboloidCircle, constraint applyAdjointHessian): "
289  "Primal vector x must be of length 2.");
290 
291  n = vp->size();
292  TEUCHOS_TEST_FOR_EXCEPTION( (n != 2), std::invalid_argument, ">>> ERROR (ROL_ParaboloidCircle, constraint applyAdjointHessian): "
293  "Direction vector v must be of length 2.");
294 
295  n = ahuvp->size();
296  TEUCHOS_TEST_FOR_EXCEPTION( (n != 2), std::invalid_argument, ">>> ERROR (ROL_ParaboloidCircle, constraint applyAdjointHessian): "
297  "Output vector ahuv must be of length 2.");
298  uint d = up->size();
299  TEUCHOS_TEST_FOR_EXCEPTION( (d != 1), std::invalid_argument, ">>> ERROR (ROL_ParaboloidCircle, constraint applyAdjointHessian): "
300  "Dual constraint vector u must be of length 1.");
301 
302  Real v1 = (*vp)[0];
303  Real v2 = (*vp)[1];
304 
305  Real u1 = (*up)[0];
306 
307  Real two(2);
308 
309  (*ahuvp)[0] = two*u1*v1;
310  (*ahuvp)[1] = two*u1*v2;
311  }
312  } //applyAdjointHessian
313 
314  };
315 
316 
317  template<class Real, class XPrim, class XDual, class CPrim, class CDual>
318  void getParaboloidCircle( Teuchos::RCP<Objective<Real> > &obj,
319  Teuchos::RCP<EqualityConstraint<Real> > &constr,
320  Vector<Real> &x0,
321  Vector<Real> &sol ) {
322 
323  typedef std::vector<Real> vector;
324 
325  typedef typename vector::size_type uint;
326 
327  using Teuchos::RCP; using Teuchos::rcp;
328  using Teuchos::dyn_cast;
329 
330  // Cast initial guess and solution vectors.
331  RCP<vector> x0p = dyn_cast<XPrim>(x0).getVector();
332  RCP<vector> solp = dyn_cast<XPrim>(sol).getVector();
333 
334  uint n = 2;
335 
336  // Resize vectors.
337  x0p->resize(n);
338  solp->resize(n);
339  // Instantiate objective function.
340  obj = Teuchos::rcp( new Objective_ParaboloidCircle<Real, XPrim, XDual> );
341  // Instantiate constraints.
343  // later we will bundle equality constraints into constraints ...
344  //std::vector<Teuchos::RCP<EqualityConstraint<Real> > > eqc( 1, Teuchos::rcp( new EqualityConstraint_ParaboloidCircle<Real> ) );
345  //constr = Teuchos::rcp( new Constraints<Real>(eqc) );
346 
347  // Get initial guess.
348  Real zero(0), one(1);
349  (*x0p)[0] = static_cast<Real>(rand())/static_cast<Real>(RAND_MAX);
350  (*x0p)[1] = static_cast<Real>(rand())/static_cast<Real>(RAND_MAX);
351  // Get solution.
352  (*solp)[0] = one;
353  (*solp)[1] = zero;
354  }
355 
356 } // End ZOO Namespace
357 } // End ROL Namespace
358 
359 #endif
Provides the interface to evaluate objective functions.
Equality constraint c(x,y) = (x-2)^2 + y^2 - 1.
void applyJacobian(Vector< Real > &jv, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Apply the constraint Jacobian at , , to vector .
virtual 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 ...
Teuchos::RCP< vector > getVector(V &x)
Teuchos::RCP< const vector > getVector(const V &x)
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:74
void getParaboloidCircle(Teuchos::RCP< Objective< Real > > &obj, Teuchos::RCP< EqualityConstraint< Real > > &constr, Vector< Real > &x0, Vector< Real > &sol)
Defines the equality constraint operator interface.
void value(Vector< Real > &c, const Vector< Real > &x, Real &tol)
Evaluate the constraint operator at .
Real value(const Vector< Real > &x, Real &tol)
Compute value.
Teuchos::RCP< const vector > getVector(const V &x)
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 ...
Objective function: f(x,y) = x^2 + y^2.
void hessVec(Vector< Real > &hv, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Apply Hessian approximation to vector.
void gradient(Vector< Real > &g, const Vector< Real > &x, Real &tol)
Compute gradient.
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 .