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 
56 #include "ROL_ScaledStdVector.hpp"
57 #include "ROL_Objective.hpp"
58 #include "ROL_BoundConstraint.hpp"
59 
60 namespace ROL {
61 namespace 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:
71  uint dim_;
72 
73  public:
74  Objective_BVP(void) : dim_(20) {}
75 
76  Real value( const Vector<Real> &x, Real &tol ) {
77  Teuchos::RCP<const std::vector<Real> > ex
78  = Teuchos::dyn_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  Teuchos::RCP<std::vector<Real> > eg
94  = Teuchos::dyn_cast<StdVector<Real> >(g).getVector();
95  Teuchos::RCP<const std::vector<Real> > ex
96  = Teuchos::dyn_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  Teuchos::RCP<std::vector<Real> > ehv
118  = Teuchos::dyn_cast<StdVector<Real> >(hv).getVector();
119  Teuchos::RCP<const std::vector<Real> > ev
120  = Teuchos::dyn_cast<const StdVector<Real> >(v).getVector();
121  Teuchos::RCP<const std::vector<Real> > ex
122  = Teuchos::dyn_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  void getBVP( Teuchos::RCP<Objective<Real> > &obj,
157  Teuchos::RCP<BoundConstraint<Real> > &con,
158  Teuchos::RCP<Vector<Real> > &x0,
159  Teuchos::RCP<Vector<Real> > &x ) {
160  // Problem dimension
161  int n = 20;
162 
163  // Set scale
164  Teuchos::RCP<std::vector<Real> > scale = Teuchos::rcp(new 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  // Get Initial Guess
187  Teuchos::RCP<std::vector<Real> > x0p = Teuchos::rcp(new std::vector<Real>(n,0.0));
188  Real h = 1.0/((Real)n + 1.0);
189  for ( int i = 0; i < n; i++ ) {
190  (*x0p)[i] = (Real)(i+1)*h*((Real)(i+1)*h - 1.0);
191  }
192  x0 = Teuchos::rcp(new PrimalScaledStdVector<Real>(x0p,scale));
193 
194  // Get Solution
195  Teuchos::RCP<std::vector<Real> > xp = Teuchos::rcp(new std::vector<Real>(n,0.0));
196  (*xp)[0] = 1.2321000000000001e-01;
197  (*xp)[1] = 2.1743122909175336e-01;
198  (*xp)[2] = 2.8625218549543746e-01;
199  (*xp)[3] = 3.3309751851140840e-01;
200  (*xp)[4] = 3.6117201714254760e-01;
201  (*xp)[5] = 3.7342787212179440e-01;
202  (*xp)[6] = 3.7255212003706123e-01;
203  (*xp)[7] = 3.6096984201471016e-01;
204  (*xp)[8] = 3.4085861052124522e-01;
205  (*xp)[9] = 3.1417024791439530e-01;
206  (*xp)[10] = 2.8265678244892922e-01;
207  (*xp)[11] = 2.4789833165179542e-01;
208  (*xp)[12] = 2.1133139591375166e-01;
209  (*xp)[13] = 1.7427666644258599e-01;
210  (*xp)[14] = 1.3796594229036069e-01;
211  (*xp)[15] = 1.0356813245768780e-01;
212  (*xp)[16] = 7.2214621084083663e-02;
213  (*xp)[17] = 4.5024529114833199e-02;
214  (*xp)[18] = 2.3130648161534966e-02;
215  (*xp)[19] = 7.7070870882527927e-03;
216  x = Teuchos::rcp(new PrimalScaledStdVector<Real>(xp,scale));
217 
218  // Instantiate Objective Function
219  obj = Teuchos::rcp(new Objective_BVP<Real>);
220 
221  // Instantiate BoundConstraint
222  Teuchos::RCP<std::vector<Real> > lp = Teuchos::rcp(new std::vector<Real>);
223  Teuchos::RCP<std::vector<Real> > up = Teuchos::rcp(new std::vector<Real>);
224  std::vector<Real> val(n,0.0);
225  val[0] = 0.1*0.2321;
226  val[1] = -0.1*0.4520;
227  val[2] = -0.1*0.6588;
228  val[3] = -0.1*0.8514;
229  val[4] = -0.1*1.0288;
230  val[5] = -0.1*1.1985;
231  val[6] = -0.1*1.3322;
232  val[7] = -0.1*1.4553;
233  val[8] = -0.1*1.5571;
234  val[9] = -0.1*1.6354;
235  val[10] = -0.1*1.6881;
236  val[11] = -0.1*1.7127;
237  val[12] = -0.1*1.7060;
238  val[13] = -0.1*1.6650;
239  val[14] = -0.1*1.5856;
240  val[15] = -0.1*1.4636;
241  val[16] = -0.1*1.2938;
242  val[17] = -0.1*1.0702;
243  val[18] = -0.1*0.7858;
244  val[19] = -0.1*0.4323;
245  for ( int i = 0; i < n; i++ ) {
246  if ( i%2 == 0 ) {
247  lp->push_back(std::max(-0.2*(Real)(n),val[i]+0.1));
248  up->push_back(std::min( 0.2*(Real)(n),val[i]+1.1));
249  }
250  else {
251  lp->push_back(-0.2*(Real)(n));
252  up->push_back( 0.2*(Real)(n));
253  }
254  }
255  Teuchos::RCP<Vector<Real> > l = Teuchos::rcp( new StdVector<Real>(lp) );
256  Teuchos::RCP<Vector<Real> > u = Teuchos::rcp( new StdVector<Real>(up) );
257  con = Teuchos::rcp(new BoundConstraint<Real>(l,u));
258  con->project(*x0);
259  }
260 
261 }// End ZOO Namespace
262 }// End ROL Namespace
263 
264 #endif
Provides the interface to evaluate objective functions.
void gradient(Vector< Real > &g, const Vector< Real > &x, Real &tol)
Compute gradient.
Definition: ROL_BVP.hpp:92
Real value(const Vector< Real > &x, Real &tol)
Compute value.
Definition: ROL_BVP.hpp:76
virtual void hessVec(Vector< Real > &hv, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Apply Hessian approximation to vector.
virtual void zero()
Set to zero vector.
Definition: ROL_Vector.hpp:157
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:74
std::vector< Real >::size_type uint
Definition: ROL_BVP.hpp:68
The discrete boundary value problem.
Definition: ROL_BVP.hpp:66
Provides the interface to apply upper and lower bound constraints.
Provides the std::vector implementation of the ROL::Vector interface that handles scalings in the inn...
void getBVP(Teuchos::RCP< Objective< Real > > &obj, Teuchos::RCP< BoundConstraint< Real > > &con, Teuchos::RCP< Vector< Real > > &x0, Teuchos::RCP< Vector< Real > > &x)
Definition: ROL_BVP.hpp:156