ROL
ROL_BVP.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
49#ifndef USE_HESSVEC
50#define USE_HESSVEC 1
51#endif
52
53#ifndef ROL_BVP_HPP
54#define ROL_BVP_HPP
55
57#include "ROL_TestProblem.hpp"
58#include "ROL_Bounds.hpp"
59
60namespace ROL {
61namespace ZOO {
62
65 template<class Real>
66 class Objective_BVP : public Objective<Real> {
67
68 typedef typename std::vector<Real>::size_type uint;
69
70 private:
72
73 public:
74 Objective_BVP(void) : dim_(20) {}
75
76 Real value( const Vector<Real> &x, Real &tol ) {
77 Ptr<const std::vector<Real> > ex
78 = dynamic_cast<const StdVector<Real>&>(x).getVector();
79
80 Real val = 0.0;
81 Real f = 0.0;
82 Real h = 1.0/((Real)(dim_) + 1.0);
83 for ( uint i = 0; i < dim_; i++ ) {
84 f = 2.0*(*ex)[i] + h*h*std::pow((*ex)[i] + (Real)(i+1)*h + 1.0,3.0)/2.0;
85 if ( i < (dim_-1) ) { f -= (*ex)[i+1]; }
86 if ( i > 0 ) { f -= (*ex)[i-1]; }
87 val += f*f;
88 }
89 return val;
90 }
91
92 void gradient( Vector<Real> &g, const Vector<Real> &x, Real &tol ) {
93 Ptr<std::vector<Real> > eg
94 = dynamic_cast<StdVector<Real>&>(g).getVector();
95 Ptr<const std::vector<Real> > ex
96 = dynamic_cast<const StdVector<Real>&>(x).getVector();
97
98 g.zero();
99 Real h = 1.0/((Real)(dim_) + 1.0);
100 std::vector<Real> f(dim_,0.0);
101
102 for ( uint i = 0; i < dim_; i++ ) {
103 f[i] = 2.0*(*ex)[i] + h*h*std::pow((*ex)[i] + (Real)(i+1)*h + 1.0,3.0)/2.0;
104 if ( i < (dim_-1) ) { f[i] -= (*ex)[i+1]; }
105 if ( i > 0) { f[i] -= (*ex)[i-1]; }
106 }
107 Real df = 0.0;
108 for ( uint i = 0; i < dim_; i++ ) {
109 df = (2.0 + 3.0*h*h*std::pow((*ex)[i] + (Real)(i+1)*h + 1.0,2.0)/2.0)*f[i];
110 if ( i < (dim_-1) ) { df -= f[i+1]; }
111 if ( i > 0 ) { df -= f[i-1]; }
112 (*eg)[i] += 2.0*df;
113 }
114 }
115#if USE_HESSVEC
116 void hessVec( Vector<Real> &hv, const Vector<Real> &v, const Vector<Real> &x, Real &tol ) {
117 Ptr<std::vector<Real> > ehv
118 = dynamic_cast<StdVector<Real>&>(hv).getVector();
119 Ptr<const std::vector<Real> > ev
120 = dynamic_cast<const StdVector<Real>&>(v).getVector();
121 Ptr<const std::vector<Real> > ex
122 = dynamic_cast<const StdVector<Real>&>(x).getVector();
123
124 hv.zero();
125 Real h = 1.0/((Real)(dim_) + 1.0);
126 Real f = 0.0, df = 0.0, dfn = 0.0, hf = 0.0;
127 for ( uint i = 0; i < dim_; i++ ) {
128 f = 2.0*(*ex)[i] + h*h*std::pow((*ex)[i] + (Real)(i+1)*h + 1.0,3.0)/2.0;
129 df = 2.0 + 3.0/2.0 * h*h * std::pow((*ex)[i] + (Real)(i+1)*h + 1.0,2.0);
130 hf = 3.0 * h*h * ((*ex)[i] + (Real)(i+1)*h + 1.0);
131 if ( i < (dim_-2) ) {
132 (*ehv)[i] += 2.0*(*ev)[i+2];
133 }
134 if ( i < (dim_-1) ) {
135 f -= (*ex)[i+1];
136 dfn = 2.0 + 3.0/2.0 * h*h * std::pow((*ex)[i+1] + (Real)(i+2)*h + 1.0,2.0);
137 (*ehv)[i] -= 2.0*(df + dfn)*(*ev)[i+1];
138 (*ehv)[i] += 2.0*(*ev)[i];
139 }
140 if ( i > 0 ) {
141 f -= (*ex)[i-1];
142 dfn = 2.0 + 3.0/2.0 * h*h * std::pow((*ex)[i-1] + (Real)(i)*h + 1.0,2.0);
143 (*ehv)[i] -= 2.0*(df + dfn)*(*ev)[i-1];
144 (*ehv)[i] += 2.0*(*ev)[i];
145 }
146 if ( i > 1 ) {
147 (*ehv)[i] += 2.0*(*ev)[i-2];
148 }
149 (*ehv)[i] += 2.0*(hf*f + df*df)*(*ev)[i];
150 }
151 }
152#endif
153 };
154
155 template<class Real>
156 class getBVP : public TestProblem<Real> {
157 private:
158 int n_;
159 Ptr<std::vector<Real> > scale_;
160
161 public:
162 getBVP(void) {
163 n_ = 20;
164 scale_ = makePtr<std::vector<Real>>(n_,0.0);
165 (*scale_)[0] = 1.e2;
166 (*scale_)[1] = 1.e2;
167 (*scale_)[2] = 1.e2;
168 (*scale_)[3] = 1.e2;
169 (*scale_)[4] = 1.e2;
170 (*scale_)[5] = 1.e2;
171 (*scale_)[6] = 1.e2;
172 (*scale_)[7] = 1.e2;
173 (*scale_)[8] = 1.e2;
174 (*scale_)[9] = 1.e2;
175 (*scale_)[10] = 1.e2;
176 (*scale_)[11] = 1.e2;
177 (*scale_)[12] = 1.e2;
178 (*scale_)[13] = 1.e2;
179 (*scale_)[14] = 1.e2;
180 (*scale_)[15] = 1.e2;
181 (*scale_)[16] = 1.e4;
182 (*scale_)[17] = 1.e4;
183 (*scale_)[18] = 1.e4;
184 (*scale_)[19] = 1.e6;
185 }
186
187 Ptr<Objective<Real>> getObjective(void) const {
188 return makePtr<Objective_BVP<Real>>();
189 }
190
191 Ptr<Vector<Real>> getInitialGuess(void) const {
192 // Get Initial Guess
193 Ptr<std::vector<Real> > x0p = makePtr<std::vector<Real>>(n_,0.0);
194 Real h = 1.0/((Real)n_ + 1.0);
195 for ( int i = 0; i < n_; i++ ) {
196 (*x0p)[i] = (Real)(i+1)*h*((Real)(i+1)*h - 1.0);
197 }
198 return makePtr<PrimalScaledStdVector<Real>>(x0p,scale_);
199 }
200
201 Ptr<Vector<Real>> getSolution(const int i = 0) const {
202 // Get Solution
203 Ptr<std::vector<Real> > xp = makePtr<std::vector<Real>>(n_,0.0);
204 (*xp)[0] = 1.2321000000000001e-01;
205 (*xp)[1] = 2.1743122909175336e-01;
206 (*xp)[2] = 2.8625218549543746e-01;
207 (*xp)[3] = 3.3309751851140840e-01;
208 (*xp)[4] = 3.6117201714254760e-01;
209 (*xp)[5] = 3.7342787212179440e-01;
210 (*xp)[6] = 3.7255212003706123e-01;
211 (*xp)[7] = 3.6096984201471016e-01;
212 (*xp)[8] = 3.4085861052124522e-01;
213 (*xp)[9] = 3.1417024791439530e-01;
214 (*xp)[10] = 2.8265678244892922e-01;
215 (*xp)[11] = 2.4789833165179542e-01;
216 (*xp)[12] = 2.1133139591375166e-01;
217 (*xp)[13] = 1.7427666644258599e-01;
218 (*xp)[14] = 1.3796594229036069e-01;
219 (*xp)[15] = 1.0356813245768780e-01;
220 (*xp)[16] = 7.2214621084083663e-02;
221 (*xp)[17] = 4.5024529114833199e-02;
222 (*xp)[18] = 2.3130648161534966e-02;
223 (*xp)[19] = 7.7070870882527927e-03;
224 return makePtr<PrimalScaledStdVector<Real>>(xp,scale_);
225 }
226
227 Ptr<BoundConstraint<Real>> getBoundConstraint(void) const {
228 // Instantiate BoundConstraint
229 Ptr<std::vector<Real> > lp = makePtr<std::vector<Real>>();
230 Ptr<std::vector<Real> > up = makePtr<std::vector<Real>>();
231 std::vector<Real> val(n_,0.0);
232 val[0] = 0.1*0.2321;
233 val[1] = -0.1*0.4520;
234 val[2] = -0.1*0.6588;
235 val[3] = -0.1*0.8514;
236 val[4] = -0.1*1.0288;
237 val[5] = -0.1*1.1985;
238 val[6] = -0.1*1.3322;
239 val[7] = -0.1*1.4553;
240 val[8] = -0.1*1.5571;
241 val[9] = -0.1*1.6354;
242 val[10] = -0.1*1.6881;
243 val[11] = -0.1*1.7127;
244 val[12] = -0.1*1.7060;
245 val[13] = -0.1*1.6650;
246 val[14] = -0.1*1.5856;
247 val[15] = -0.1*1.4636;
248 val[16] = -0.1*1.2938;
249 val[17] = -0.1*1.0702;
250 val[18] = -0.1*0.7858;
251 val[19] = -0.1*0.4323;
252 for ( int i = 0; i < n_; i++ ) {
253 if ( i%2 == 0 ) {
254 lp->push_back(std::max(-0.2*(Real)(n_),val[i]+0.1));
255 up->push_back(std::min( 0.2*(Real)(n_),val[i]+1.1));
256 }
257 else {
258 lp->push_back(-0.2*(Real)(n_));
259 up->push_back( 0.2*(Real)(n_));
260 }
261 }
262 Ptr<Vector<Real> > l = makePtr<StdVector<Real>>(lp);
263 Ptr<Vector<Real> > u = makePtr<StdVector<Real>>(up);
264 return makePtr<Bounds<Real>>(l,u);
265 }
266 };
267
268}// End ZOO Namespace
269}// End ROL Namespace
270
271#endif
Contains definitions of test objective functions.
Provides the interface to evaluate objective functions.
virtual void hessVec(Vector< Real > &hv, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Apply Hessian approximation to vector.
Provides the ROL::Vector interface for scalar values, to be used, for example, with scalar constraint...
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:84
virtual void zero()
Set to zero vector.
Definition: ROL_Vector.hpp:167
The discrete boundary value problem.
Definition: ROL_BVP.hpp:66
Real value(const Vector< Real > &x, Real &tol)
Compute value.
Definition: ROL_BVP.hpp:76
void gradient(Vector< Real > &g, const Vector< Real > &x, Real &tol)
Compute gradient.
Definition: ROL_BVP.hpp:92
std::vector< Real >::size_type uint
Definition: ROL_BVP.hpp:68
Ptr< BoundConstraint< Real > > getBoundConstraint(void) const
Definition: ROL_BVP.hpp:227
Ptr< std::vector< Real > > scale_
Definition: ROL_BVP.hpp:159
Ptr< Objective< Real > > getObjective(void) const
Definition: ROL_BVP.hpp:187
Ptr< Vector< Real > > getSolution(const int i=0) const
Definition: ROL_BVP.hpp:201
Ptr< Vector< Real > > getInitialGuess(void) const
Definition: ROL_BVP.hpp:191