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_TestProblem.hpp"
57 #include "ROL_StdVector.hpp"
58 #include "Teuchos_SerialDenseVector.hpp"
59 #include "Teuchos_SerialDenseSolver.hpp"
60 
61 namespace ROL {
62 namespace ZOO {
63 
67  template< class Real, class XPrim=StdVector<Real>, class XDual=StdVector<Real> >
68  class Objective_ParaboloidCircle : public Objective<Real> {
69 
70  typedef std::vector<Real> vector;
71  typedef Vector<Real> V;
72 
73  typedef typename vector::size_type uint;
74 
75 
76  private:
77 
78  template<class VectorType>
79  ROL::Ptr<const vector> getVector( const V& x ) {
80 
81  return dynamic_cast<const VectorType&>(x).getVector();
82  }
83 
84  template<class VectorType>
85  ROL::Ptr<vector> getVector( V& x ) {
86 
87  return dynamic_cast<VectorType&>(x).getVector();
88  }
89 
90  public:
91  Objective_ParaboloidCircle() {}
92 
93  Real value( const Vector<Real> &x, Real &tol ) {
94 
95 
96  ROL::Ptr<const vector> xp = getVector<XPrim>(x);
97 
98  uint n = xp->size();
99  ROL_TEST_FOR_EXCEPTION( (n != 2), std::invalid_argument, ">>> ERROR (ROL_ParaboloidCircle, objective value): "
100  "Primal vector x must be of length 2.");
101 
102  Real x1 = (*xp)[0];
103  Real x2 = (*xp)[1];
104 
105  Real val = x1*x1 + x2*x2;
106 
107  return val;
108  }
109 
110  void gradient( Vector<Real> &g, const Vector<Real> &x, Real &tol ) {
111 
112 
113  ROL::Ptr<const vector> xp = getVector<XPrim>(x);
114  ROL::Ptr<vector> gp = getVector<XDual>(g);
115 
116  uint n = xp->size();
117  ROL_TEST_FOR_EXCEPTION( (n != 2), std::invalid_argument, ">>> ERROR (ROL_ParaboloidCircle, objective gradient): "
118  " Primal vector x must be of length 2.");
119 
120  n = gp->size();
121  ROL_TEST_FOR_EXCEPTION( (n != 2), std::invalid_argument, ">>> ERROR (ROL_ParaboloidCircle, objective gradient): "
122  "Gradient vector g must be of length 2.");
123 
124  Real x1 = (*xp)[0];
125  Real x2 = (*xp)[1];
126 
127  Real two(2);
128 
129  (*gp)[0] = two*x1;
130  (*gp)[1] = two*x2;
131  }
132 
133  void hessVec( Vector<Real> &hv, const Vector<Real> &v, const Vector<Real> &x, Real &tol ) {
134 
135 
136  ROL::Ptr<const vector> xp = getVector<XPrim>(x);
137  ROL::Ptr<const vector> vp = getVector<XPrim>(v);
138  ROL::Ptr<vector> hvp = getVector<XDual>(hv);
139 
140  uint n = xp->size();
141  ROL_TEST_FOR_EXCEPTION( (n != 2), std::invalid_argument, ">>> ERROR (ROL_ParaboloidCircle, objective hessVec): "
142  "Primal vector x must be of length 2.");
143 
144  n = vp->size();
145  ROL_TEST_FOR_EXCEPTION( (n != 2), std::invalid_argument, ">>> ERROR (ROL_ParaboloidCircle, objective hessVec): "
146  "Input vector v must be of length 2.");
147 
148  n = hvp->size();
149  ROL_TEST_FOR_EXCEPTION( (n != 2), std::invalid_argument, ">>> ERROR (ROL_ParaboloidCircle, objective hessVec): "
150  "Output vector hv must be of length 2.");
151 
152  Real v1 = (*vp)[0];
153  Real v2 = (*vp)[1];
154 
155  Real two(2);
156 
157  (*hvp)[0] = two*v1;
158  (*hvp)[1] = two*v2;
159  }
160 
161  };
162 
163 
166  template<class Real, class XPrim=StdVector<Real>, class XDual=StdVector<Real>, class CPrim=StdVector<Real>, class CDual=StdVector<Real> >
167  class Constraint_ParaboloidCircle : public Constraint<Real> {
168 
169  typedef std::vector<Real> vector;
170  typedef Vector<Real> V;
171 
172  typedef typename vector::size_type uint;
173 
174  private:
175  template<class VectorType>
176  ROL::Ptr<const vector> getVector( const V& x ) {
177 
178  return dynamic_cast<const VectorType&>(x).getVector();
179  }
180 
181  template<class VectorType>
182  ROL::Ptr<vector> getVector( V& x ) {
183 
184  return dynamic_cast<VectorType&>(x).getVector();
185  }
186 
187  public:
188  Constraint_ParaboloidCircle() {}
189 
190  void value( Vector<Real> &c, const Vector<Real> &x, Real &tol ) {
191 
192 
193  ROL::Ptr<const vector> xp = getVector<XPrim>(x);
194  ROL::Ptr<vector> cp = getVector<CPrim>(c);
195 
196  uint n = xp->size();
197  ROL_TEST_FOR_EXCEPTION( (n != 2), std::invalid_argument, ">>> ERROR (ROL_ParaboloidCircle, constraint value): "
198  "Primal vector x must be of length 2.");
199 
200  uint m = cp->size();
201  ROL_TEST_FOR_EXCEPTION( (m != 1), std::invalid_argument, ">>> ERROR (ROL_ParaboloidCircle, constraint value): "
202  "Constraint vector c must be of length 1.");
203 
204  Real x1 = (*xp)[0];
205  Real x2 = (*xp)[1];
206 
207  Real one(1), two(2);
208 
209  (*cp)[0] = (x1-two)*(x1-two) + x2*x2 - one;
210  }
211 
212  void applyJacobian( Vector<Real> &jv, const Vector<Real> &v, const Vector<Real> &x, Real &tol ) {
213 
214 
215  ROL::Ptr<const vector> xp = getVector<XPrim>(x);
216  ROL::Ptr<const vector> vp = getVector<XPrim>(v);
217  ROL::Ptr<vector> jvp = getVector<CPrim>(jv);
218 
219  uint n = xp->size();
220  ROL_TEST_FOR_EXCEPTION( (n != 2), std::invalid_argument, ">>> ERROR (ROL_ParaboloidCircle, constraint applyJacobian): "
221  "Primal vector x must be of length 2.");
222 
223  uint d = vp->size();
224  ROL_TEST_FOR_EXCEPTION( (d != 2), std::invalid_argument, ">>> ERROR (ROL_ParaboloidCircle, constraint applyJacobian): "
225  "Input vector v must be of length 2.");
226  d = jvp->size();
227  ROL_TEST_FOR_EXCEPTION( (d != 1), std::invalid_argument, ">>> ERROR (ROL_ParaboloidCircle, constraint applyJacobian): "
228  "Output vector jv must be of length 1.");
229 
230  Real x1 = (*xp)[0];
231  Real x2 = (*xp)[1];
232 
233  Real v1 = (*vp)[0];
234  Real v2 = (*vp)[1];
235 
236  Real two(2);
237 
238  (*jvp)[0] = two*(x1-two)*v1 + two*x2*v2;
239  } //applyJacobian
240 
241  void applyAdjointJacobian( Vector<Real> &ajv, const Vector<Real> &v, const Vector<Real> &x, Real &tol ) {
242 
243 
244  ROL::Ptr<const vector> xp = getVector<XPrim>(x);
245  ROL::Ptr<const vector> vp = getVector<CDual>(v);
246  ROL::Ptr<vector> ajvp = getVector<XDual>(ajv);
247 
248  uint n = xp->size();
249  ROL_TEST_FOR_EXCEPTION( (n != 2), std::invalid_argument, ">>> ERROR (ROL_ParaboloidCircle, constraint applyAdjointJacobian): "
250  "Primal vector x must be of length 2.");
251 
252  uint d = vp->size();
253  ROL_TEST_FOR_EXCEPTION( (d != 1), std::invalid_argument, ">>> ERROR (ROL_ParaboloidCircle, constraint applyAdjointJacobian): "
254  "Input vector v must be of length 1.");
255 
256  d = ajvp->size();
257  ROL_TEST_FOR_EXCEPTION( (d != 2), std::invalid_argument, ">>> ERROR (ROL_ParaboloidCircle, constraint applyAdjointJacobian): "
258  "Output vector ajv must be of length 2.");
259 
260  Real x1 = (*xp)[0];
261  Real x2 = (*xp)[1];
262 
263  Real v1 = (*vp)[0];
264 
265  Real two(2);
266 
267  (*ajvp)[0] = two*(x1-two)*v1;
268  (*ajvp)[1] = two*x2*v1;
269 
270  } //applyAdjointJacobian
271 
272  void applyAdjointHessian( Vector<Real> &ahuv, const Vector<Real> &u, const Vector<Real> &v, const Vector<Real> &x, Real &tol ) {
273 
274  bool useFD = true;
275 
276  if (useFD) {
277  Constraint<Real>::applyAdjointHessian( ahuv, u, v, x, tol );
278  }
279  else {
280 
281  ROL::Ptr<const vector> xp = getVector<XPrim>(x);
282  ROL::Ptr<const vector> up = getVector<CDual>(u);
283  ROL::Ptr<const vector> vp = getVector<XPrim>(v);
284  ROL::Ptr<vector> ahuvp = getVector<XDual>(ahuv);
285 
286  uint n = xp->size();
287  ROL_TEST_FOR_EXCEPTION( (n != 2), std::invalid_argument, ">>> ERROR (ROL_ParaboloidCircle, constraint applyAdjointHessian): "
288  "Primal vector x must be of length 2.");
289 
290  n = vp->size();
291  ROL_TEST_FOR_EXCEPTION( (n != 2), std::invalid_argument, ">>> ERROR (ROL_ParaboloidCircle, constraint applyAdjointHessian): "
292  "Direction vector v must be of length 2.");
293 
294  n = ahuvp->size();
295  ROL_TEST_FOR_EXCEPTION( (n != 2), std::invalid_argument, ">>> ERROR (ROL_ParaboloidCircle, constraint applyAdjointHessian): "
296  "Output vector ahuv must be of length 2.");
297  uint d = up->size();
298  ROL_TEST_FOR_EXCEPTION( (d != 1), std::invalid_argument, ">>> ERROR (ROL_ParaboloidCircle, constraint applyAdjointHessian): "
299  "Dual constraint vector u must be of length 1.");
300 
301  Real v1 = (*vp)[0];
302  Real v2 = (*vp)[1];
303 
304  Real u1 = (*up)[0];
305 
306  Real two(2);
307 
308  (*ahuvp)[0] = two*u1*v1;
309  (*ahuvp)[1] = two*u1*v2;
310  }
311  } //applyAdjointHessian
312 
313  };
314 
315 
316  template<class Real, class XPrim=StdVector<Real>, class XDual=StdVector<Real>, class CPrim=StdVector<Real>, class CDual=StdVector<Real> >
317  class getParaboloidCircle : public TestProblem<Real> {
318  typedef std::vector<Real> vector;
319  typedef typename vector::size_type uint;
320  public:
321  getParaboloidCircle(void) {}
322 
323  Ptr<Objective<Real>> getObjective(void) const {
324  // Instantiate objective function.
325  return ROL::makePtr<Objective_ParaboloidCircle<Real,XPrim,XDual>>();
326  }
327 
328  Ptr<Vector<Real>> getInitialGuess(void) const {
329  uint n = 2;
330  // Get initial guess.
331  ROL::Ptr<vector> x0p = makePtr<vector>(n,0.0);
332  (*x0p)[0] = static_cast<Real>(rand())/static_cast<Real>(RAND_MAX);
333  (*x0p)[1] = static_cast<Real>(rand())/static_cast<Real>(RAND_MAX);
334  return makePtr<XPrim>(x0p);
335  }
336 
337  Ptr<Vector<Real>> getSolution(const int i = 0) const {
338  uint n = 2;
339  // Get solution.
340  Real zero(0), one(1);
341  ROL::Ptr<vector> solp = makePtr<vector>(n,0.0);
342  (*solp)[0] = one;
343  (*solp)[1] = zero;
344  return makePtr<XPrim>(solp);
345  }
346 
347  Ptr<Constraint<Real>> getEqualityConstraint(void) const {
348  // Instantiate constraints.
349  return ROL::makePtr<Constraint_ParaboloidCircle<Real,XPrim,XDual,CPrim,CDual>>();
350  }
351 
352  Ptr<Vector<Real>> getEqualityMultiplier(void) const {
353  ROL::Ptr<vector> lp = makePtr<vector>(1,0.0);
354  return makePtr<CDual>(lp);
355  }
356  };
357 
358 } // End ZOO Namespace
359 } // End ROL Namespace
360 
361 #endif
ROL::ZOO::Objective_ParaboloidCircle::getVector
ROL::Ptr< vector > getVector(V &x)
Definition: ROL_ParaboloidCircle.hpp:85
ROL::ZOO::getParaboloidCircle::getObjective
Ptr< Objective< Real > > getObjective(void) const
Definition: ROL_ParaboloidCircle.hpp:323
ROL::ZOO::getParaboloidCircle::getEqualityConstraint
Ptr< Constraint< Real > > getEqualityConstraint(void) const
Definition: ROL_ParaboloidCircle.hpp:347
zero
Objective_SerialSimOpt(const Ptr< Obj > &obj, const V &ui) z0_ zero()
Definition: ROL_Objective_SerialSimOpt.hpp:163
ROL::ZOO::Objective_ParaboloidCircle::vector
std::vector< Real > vector
Definition: ROL_ParaboloidCircle.hpp:70
ROL::ZOO::Constraint_ParaboloidCircle::getVector
ROL::Ptr< vector > getVector(V &x)
Definition: ROL_ParaboloidCircle.hpp:182
ROL::ZOO::Constraint_ParaboloidCircle::value
void value(Vector< Real > &c, const Vector< Real > &x, Real &tol)
Evaluate the constraint operator at .
Definition: ROL_ParaboloidCircle.hpp:190
ROL::ZOO::Constraint_ParaboloidCircle::vector
std::vector< Real > vector
Definition: ROL_ParaboloidCircle.hpp:169
ROL::ZOO::getParaboloidCircle::getEqualityMultiplier
Ptr< Vector< Real > > getEqualityMultiplier(void) const
Definition: ROL_ParaboloidCircle.hpp:352
ROL::ZOO::Constraint_ParaboloidCircle::Constraint_ParaboloidCircle
Constraint_ParaboloidCircle()
Definition: ROL_ParaboloidCircle.hpp:188
ROL::ZOO::Objective_ParaboloidCircle::value
Real value(const Vector< Real > &x, Real &tol)
Compute value.
Definition: ROL_ParaboloidCircle.hpp:93
ROL::Constraint::applyAdjointHessian
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 ,...
Definition: ROL_ConstraintDef.hpp:166
ROL::ZOO::getParaboloidCircle::vector
std::vector< Real > vector
Definition: ROL_ParaboloidCircle.hpp:318
ROL::ZOO::getParaboloidCircle::getParaboloidCircle
getParaboloidCircle(void)
Definition: ROL_ParaboloidCircle.hpp:321
ROL::ZOO::Constraint_ParaboloidCircle
constraint c(x,y) = (x-2)^2 + y^2 - 1.
Definition: ROL_ParaboloidCircle.hpp:167
ROL::ZOO::getParaboloidCircle
Definition: ROL_ParaboloidCircle.hpp:317
ROL::ZOO::Constraint_ParaboloidCircle::V
Vector< Real > V
Definition: ROL_ParaboloidCircle.hpp:170
ROL::ZOO::Objective_ParaboloidCircle::uint
vector::size_type uint
Definition: ROL_ParaboloidCircle.hpp:73
ROL::ZOO::getParaboloidCircle::getInitialGuess
Ptr< Vector< Real > > getInitialGuess(void) const
Definition: ROL_ParaboloidCircle.hpp:328
ROL::ZOO::Objective_ParaboloidCircle
Objective function: f(x,y) = x^2 + y^2.
Definition: ROL_ParaboloidCircle.hpp:68
ROL::Vector
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:80
ROL::ZOO::Objective_ParaboloidCircle::hessVec
void hessVec(Vector< Real > &hv, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Apply Hessian approximation to vector.
Definition: ROL_ParaboloidCircle.hpp:133
ROL::ZOO::Objective_ParaboloidCircle::gradient
void gradient(Vector< Real > &g, const Vector< Real > &x, Real &tol)
Compute gradient.
Definition: ROL_ParaboloidCircle.hpp:110
ROL::Constraint
Defines the general constraint operator interface.
Definition: ROL_Constraint.hpp:85
ROL::TestProblem
Definition: ROL_TestProblem.hpp:57
ROL::ZOO::Constraint_ParaboloidCircle::uint
vector::size_type uint
Definition: ROL_ParaboloidCircle.hpp:172
ROL
Definition: ROL_ElementwiseVector.hpp:61
ROL::ZOO::Constraint_ParaboloidCircle::applyAdjointHessian
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 ,...
Definition: ROL_ParaboloidCircle.hpp:272
ROL::ZOO::Constraint_ParaboloidCircle::applyJacobian
void applyJacobian(Vector< Real > &jv, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Apply the constraint Jacobian at , , to vector .
Definition: ROL_ParaboloidCircle.hpp:212
ROL::ZOO::Objective_ParaboloidCircle::Objective_ParaboloidCircle
Objective_ParaboloidCircle()
Definition: ROL_ParaboloidCircle.hpp:91
ROL::ZOO::getParaboloidCircle::getSolution
Ptr< Vector< Real > > getSolution(const int i=0) const
Definition: ROL_ParaboloidCircle.hpp:337
size_type
typename PV< Real >::size_type size_type
Definition: ROL_Objective_SerialSimOpt.hpp:142
ROL::ZOO::Objective_ParaboloidCircle::getVector
ROL::Ptr< const vector > getVector(const V &x)
Definition: ROL_ParaboloidCircle.hpp:79
ROL::ZOO::getParaboloidCircle::uint
vector::size_type uint
Definition: ROL_ParaboloidCircle.hpp:319
ROL::Objective
Provides the interface to evaluate objective functions.
Definition: ROL_Objective.hpp:77
ROL::ZOO::Objective_ParaboloidCircle::V
Vector< Real > V
Definition: ROL_ParaboloidCircle.hpp:71
ROL_StdVector.hpp
ROL_TestProblem.hpp
Contains definitions of test objective functions.
ROL::ZOO::Constraint_ParaboloidCircle::getVector
ROL::Ptr< const vector > getVector(const V &x)
Definition: ROL_ParaboloidCircle.hpp:176
ROL::ZOO::Constraint_ParaboloidCircle::applyAdjointJacobian
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 .
Definition: ROL_ParaboloidCircle.hpp:241
Generated by   doxygen 1.8.16