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 
44 /** \file
45     \brief  Contains definitions for W. Hock and K. Schittkowski 39th test function.
46     \author Created by G. von Winckel
47  */
48 
49 
50 #ifndef ROL_HS39_HPP
51 #define ROL_HS39_HPP
52 
53 #include "ROL_StdVector.hpp"
54 #include "ROL_TestProblem.hpp"
55 #include "ROL_Bounds.hpp"
56 #include "ROL_Constraint_Partitioned.hpp"
57 #include "ROL_Types.hpp"
58 
59 namespace ROL {
60 namespace ZOO {
61 
62 /** \brief W. Hock and K. Schittkowski 39th test function.
63  *
64  *  Exact solution x* = (1, 0, 0, 0)
65  *  f(x*) = -1
66  */
67 
68 template<class Real>
69 class Objective_HS39 : public Objective<Real> {
70 
71   typedef std::vector<Real>    vector;
72 
73   typedef StdVector<Real>   SV;
74 
75 
76 
77 public:
78 
value(const Vector<Real> & x,Real & tol)79   Real value( const Vector<Real> &x, Real &tol ) {
80     ROL::Ptr<const vector> xp = dynamic_cast<const SV&>(x).getVector();
81     return -(*xp)[0];
82   }
83 
gradient(Vector<Real> & g,const Vector<Real> & x,Real & tol)84   void gradient( Vector<Real> &g, const Vector<Real> &x, Real &tol ) {
85     ROL::Ptr<const vector> xp = dynamic_cast<const SV&>(x).getVector();
86     ROL::Ptr<vector> gp = dynamic_cast<SV&>(g).getVector();
87 
88     (*gp)[0] = -1.0;
89     (*gp)[1] =  0.0;
90     (*gp)[2] =  0.0;
91     (*gp)[3] =  0.0;
92 
93   }
94 
hessVec(Vector<Real> & hv,const Vector<Real> & v,const Vector<Real> & x,Real & tol)95   void hessVec( Vector<Real> &hv, const Vector<Real> &v, const Vector<Real> &x, Real &tol ) {
96     hv.zero();
97   }
98 };
99 
100 // First of two equality constraints
101 template<class Real>
102 class Constraint_HS39a : public Constraint<Real> {
103 
104   typedef std::vector<Real> vector;
105   typedef StdVector<Real>   SV;
106 
107 public:
Constraint_HS39a(void)108   Constraint_HS39a(void) {}
109 
value(Vector<Real> & c,const Vector<Real> & x,Real & tol)110   void value( Vector<Real> &c, const Vector<Real> &x, Real &tol ) {
111 
112     ROL::Ptr<vector> cp = dynamic_cast<SV&>(c).getVector();
113     ROL::Ptr<const vector> xp = dynamic_cast<const SV&>(x).getVector();
114 
115     (*cp)[0] = (*xp)[1]-std::pow((*xp)[0],3)-std::pow((*xp)[2],2);
116   }
117 
applyJacobian(Vector<Real> & jv,const Vector<Real> & v,const Vector<Real> & x,Real & tol)118   void applyJacobian(Vector<Real> &jv, const Vector<Real> &v,
119                      const Vector<Real> &x, Real &tol) {
120 
121     ROL::Ptr<vector> jvp = dynamic_cast<SV&>(jv).getVector();
122     ROL::Ptr<const vector> vp = dynamic_cast<const SV&>(v).getVector();
123     ROL::Ptr<const vector> xp = dynamic_cast<const SV&>(x).getVector();
124 
125     (*jvp)[0] = (*vp)[1] - 3.0*(*xp)[0]*(*xp)[0]*(*vp)[0] - 2.0*(*xp)[2]*(*vp)[2];
126 
127   }
128 
applyAdjointJacobian(Vector<Real> & ajv,const Vector<Real> & v,const Vector<Real> & x,Real & tol)129   void applyAdjointJacobian( Vector<Real> &ajv, const Vector<Real> &v,
130                              const Vector<Real> &x, Real &tol ) {
131 
132     ROL::Ptr<vector> ajvp = dynamic_cast<SV&>(ajv).getVector();
133     ROL::Ptr<const vector> vp = dynamic_cast<const SV&>(v).getVector();
134     ROL::Ptr<const vector> xp = dynamic_cast<const SV&>(x).getVector();
135 
136     (*ajvp)[0] = -3.0*(*xp)[0]*(*xp)[0]*(*vp)[0];
137     (*ajvp)[1] = (*vp)[0];
138     (*ajvp)[2] = -2.0*(*xp)[2]*(*vp)[0];
139     (*ajvp)[3] = 0.0;
140   }
141 
applyAdjointHessian(Vector<Real> & ahuv,const Vector<Real> & u,const Vector<Real> & v,const Vector<Real> & x,Real & tol)142   void applyAdjointHessian(Vector<Real> &ahuv, const Vector<Real> &u,
143                            const Vector<Real> &v, const Vector<Real> &x,
144                            Real &tol) {
145 
146     ROL::Ptr<vector> ahuvp = dynamic_cast<SV&>(ahuv).getVector();
147     ROL::Ptr<const vector> up = dynamic_cast<const SV&>(u).getVector();
148     ROL::Ptr<const vector> vp = dynamic_cast<const SV&>(v).getVector();
149     ROL::Ptr<const vector> xp = dynamic_cast<const SV&>(x).getVector();
150 
151     (*ahuvp)[0] = -6.0*(*up)[0]*(*xp)[0]*(*vp)[0];
152     (*ahuvp)[1] = 0.0;
153     (*ahuvp)[2] = -2.0*(*up)[0]*(*vp)[2];
154     (*ahuvp)[3] = 0.0;
155 
156   }
157 
158 
159 };
160 
161 // Second of two equality constraints
162 template<class Real>
163 class Constraint_HS39b : public Constraint<Real> {
164 
165   typedef std::vector<Real> vector;
166   typedef StdVector<Real>   SV;
167 
168 public:
Constraint_HS39b(void)169   Constraint_HS39b(void) {}
170 
value(Vector<Real> & c,const Vector<Real> & x,Real & tol)171   void value( Vector<Real> &c, const Vector<Real> &x, Real &tol ) {
172     ROL::Ptr<vector> cp = dynamic_cast<SV&>(c).getVector();
173     ROL::Ptr<const vector> xp = dynamic_cast<const SV&>(x).getVector();
174 
175     (*cp)[0] = std::pow((*xp)[0],2)-(*xp)[1]-std::pow((*xp)[3],2);
176   }
177 
applyJacobian(Vector<Real> & jv,const Vector<Real> & v,const Vector<Real> & x,Real & tol)178   void applyJacobian(Vector<Real> &jv, const Vector<Real> &v,
179                      const Vector<Real> &x, Real &tol) {
180 
181     ROL::Ptr<vector> jvp = dynamic_cast<SV&>(jv).getVector();
182     ROL::Ptr<const vector> vp = dynamic_cast<const SV&>(v).getVector();
183     ROL::Ptr<const vector> xp = dynamic_cast<const SV&>(x).getVector();
184 
185     (*jvp)[0] = 2.0*(*xp)[0]*(*vp)[0] - (*vp)[1] - 2.0*(*xp)[3]*(*vp)[3];
186 
187   }
188 
applyAdjointJacobian(Vector<Real> & ajv,const Vector<Real> & v,const Vector<Real> & x,Real & tol)189   void applyAdjointJacobian( Vector<Real> &ajv, const Vector<Real> &v,
190                              const Vector<Real> &x, Real &tol ) {
191 
192     ROL::Ptr<vector> ajvp = dynamic_cast<SV&>(ajv).getVector();
193     ROL::Ptr<const vector> vp = dynamic_cast<const SV&>(v).getVector();
194     ROL::Ptr<const vector> xp = dynamic_cast<const SV&>(x).getVector();
195 
196     (*ajvp)[0] = 2.0*(*xp)[0]*(*vp)[0];
197     (*ajvp)[1] = -(*vp)[0];
198     (*ajvp)[2] = 0.0;
199     (*ajvp)[3] = -2.0*(*vp)[0]*(*xp)[3];
200   }
201 
applyAdjointHessian(Vector<Real> & ahuv,const Vector<Real> & u,const Vector<Real> & v,const Vector<Real> & x,Real & tol)202   void applyAdjointHessian(Vector<Real> &ahuv, const Vector<Real> &u,
203                            const Vector<Real> &v, const Vector<Real> &x,
204                            Real &tol) {
205 
206     ROL::Ptr<vector> ahuvp = dynamic_cast<SV&>(ahuv).getVector();
207     ROL::Ptr<const vector> up = dynamic_cast<const SV&>(u).getVector();
208     ROL::Ptr<const vector> vp = dynamic_cast<const SV&>(v).getVector();
209     ROL::Ptr<const vector> xp = dynamic_cast<const SV&>(x).getVector();
210 
211   //  (*cp)[0] = std::pow((*xp)[0],2)-(*xp)[1]-std::pow((*xp)[3],2);
212 
213     (*ahuvp)[0] = 2.0*(*up)[0]*(*vp)[0];
214     (*ahuvp)[1] = 0.0;
215     (*ahuvp)[2] = 0.0;
216     (*ahuvp)[3] = -2.0*(*up)[0]*(*vp)[3];
217   }
218 
219 
220 };
221 
222 template<class Real>
223 class getHS39 : public TestProblem<Real> {
224 public:
getHS39(void)225   getHS39(void) {}
226 
getObjective(void) const227   Ptr<Objective<Real>> getObjective(void) const {
228     // Instantiate Objective Function
229     return ROL::makePtr<Objective_HS39<Real>>();
230   }
231 
getInitialGuess(void) const232   Ptr<Vector<Real>> getInitialGuess(void) const {
233     // Problem dimension
234     int n = 4;
235     // Get Initial Guess
236     ROL::Ptr<std::vector<Real> > x0p = ROL::makePtr<std::vector<Real>>(n,2.0);
237     return ROL::makePtr<StdVector<Real>>(x0p);
238   }
239 
getSolution(const int i=0) const240   Ptr<Vector<Real>> getSolution(const int i = 0) const {
241     // Problem dimension
242     int n = 4;
243     // Get Solution
244     ROL::Ptr<std::vector<Real> > xp = ROL::makePtr<std::vector<Real>>(n,0.0);
245     (*xp)[0] = static_cast<Real>(1);
246     (*xp)[1] = static_cast<Real>(1);
247     (*xp)[2] = static_cast<Real>(0);
248     (*xp)[3] = static_cast<Real>(0);
249     return ROL::makePtr<StdVector<Real>>(xp);
250   }
251 
getEqualityConstraint(void) const252   Ptr<Constraint<Real>> getEqualityConstraint(void) const {
253     std::vector<Ptr<Constraint<Real>>> cvec(2);
254     cvec[0] = makePtr<Constraint_HS39a<Real>>();
255     cvec[1] = makePtr<Constraint_HS39b<Real>>();
256     return ROL::makePtr<Constraint_Partitioned<Real>>(cvec);
257   }
258 
getEqualityMultiplier(void) const259   Ptr<Vector<Real>> getEqualityMultiplier(void) const {
260     std::vector<Ptr<Vector<Real>>> lvec(2);
261     lvec[0] = makePtr<StdVector<Real>>(makePtr<std::vector<Real>>(1,0.0));
262     lvec[1] = makePtr<StdVector<Real>>(makePtr<std::vector<Real>>(1,0.0));
263     return ROL::makePtr<PartitionedVector<Real>>(lvec);
264   }
265 };
266 
267 
268 
269 } // End ZOO Namespace
270 } // End ROL Namespace
271 
272 #endif
273