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
61namespace ROL {
62namespace ZOO {
63
67 template< class Real, class XPrim=StdVector<Real>, class XDual=StdVector<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:
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> >
168
169 typedef std::vector<Real> vector;
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:
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:
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
Objective_SerialSimOpt(const Ptr< Obj > &obj, const V &ui) z0_ zero()
Contains definitions of test objective functions.
Defines the general constraint operator interface.
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 ,...
Provides the interface to evaluate objective functions.
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:84
constraint c(x,y) = (x-2)^2 + y^2 - 1.
void value(Vector< Real > &c, const Vector< Real > &x, Real &tol)
Evaluate the constraint operator at .
void applyJacobian(Vector< Real > &jv, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Apply the constraint Jacobian at , , to vector .
ROL::Ptr< const vector > getVector(const V &x)
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 .
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.
Real value(const Vector< Real > &x, Real &tol)
Compute value.
ROL::Ptr< const vector > getVector(const V &x)
Ptr< Vector< Real > > getSolution(const int i=0) const
Ptr< Constraint< Real > > getEqualityConstraint(void) const
Ptr< Objective< Real > > getObjective(void) const
Ptr< Vector< Real > > getEqualityMultiplier(void) const
Ptr< Vector< Real > > getInitialGuess(void) const