ROL
test_11.hpp
Go to the documentation of this file.
1 // @HEADER
2 // *****************************************************************************
3 // Rapid Optimization Library (ROL) Package
4 //
5 // Copyright 2014 NTESS and the ROL contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 // @HEADER
9 
10 /* \file test_11.hpp
11  \brief Verify that the Coleman-Li Model function produces the
12  correct values for a test problem
13 
14  \f[ \min_x f(x) = \frac{1}{2}(x_1^2+2x_2^2),\quad x_1 \geq 1,\; x_2 <= -1 \f]
15 
16  The gradient is
17  \f[ \nabla f(x) = (x_1,2 x_2 ) \f]
18 
19  and the Hessian is
20  \[f \nabla^2 f(x) = \begin{pmatrix} 1 & 0 \\ 0 & 2 \end{pmatrix} \f]
21 
22  The minimizer is \f$ x^\ast = (1,-1)\f$
23 
24  For feasible \f$x\f$, the Coleman-Li quantities of interest are
25 
26  \f[ v_1 = x_1 - 1,\quad v_2 = -1 - x_2 \f]
27  \f[ D^{-1} = \begin{pmatrix} \sqrt{|x_1-1|} & 0 \\ 0 & \sqrt{|x_2+1|} \end{pmatrix} \f]
28  \f[ J=\begin{pmatrix} 1 & -1 \end{pmatrix} \f]
29  \f[ \hat M_k = \begin{pmatrix} |x_1-1|^2+|x_1| & 0 \\
30  0 & |x_2+1|^2+2|x_2| \end{pmatrix} \f]
31 
32 
33 */
34 
35 #include "ROL_Objective.hpp"
36 #include "ROL_StdVector.hpp"
37 
38 template<class Real>
39 class CLTestObjective : public ROL::Objective<Real> {
40 
41 
42 public:
43 
44  Real value( const ROL::Vector<Real> &x, Real &tol ) {
45  ROL::Ptr<const std::vector<Real> > xp =
46  dynamic_cast<const ROL::StdVector<Real>&>(x).getVector();
47  return 0.5*((*xp)[0]*(*xp)[0] + 2*(*xp)[1]*(*xp)[1]);
48  }
49 
50  void gradient( ROL::Vector<Real> &g, const ROL::Vector<Real> &x, Real &tol ) {
51  ROL::Ptr<std::vector<Real> > gp =
52  dynamic_cast<ROL::StdVector<Real>&>(g).getVector();
53  ROL::Ptr<const std::vector<Real> > xp =
54  dynamic_cast<const ROL::StdVector<Real>&>(x).getVector();
55  (*gp)[0] = (*xp)[0];
56  (*gp)[1] = 2*(*xp)[1];
57  }
58 
60  const ROL::Vector<Real> &v,
61  const ROL::Vector<Real> &x,
62  Real &tol ) {
63  ROL::Ptr<std::vector<Real> > hvp =
64  dynamic_cast<ROL::StdVector<Real>&>(hv).getVector();
65  ROL::Ptr<const std::vector<Real> > vp =
66  dynamic_cast<const ROL::StdVector<Real>&>(v).getVector();
67  (*hvp)[0] = (*vp)[0];
68  (*hvp)[1] = 2*(*vp)[1];
69  }
70 
71 }; // CLTestObjective
72 
73 template<class Real>
74 class CLExactModel : public ROL::Objective<Real> {
75 
76 ROL::Ptr<std::vector<Real> > x_;
77 const ROL::Ptr<const std::vector<Real> > l_;
78 const ROL::Ptr<const std::vector<Real> > u_;
79 ROL::Ptr<std::vector<Real> > g_;
80 ROL::Ptr<std::vector<Real> > di_;
81 ROL::Ptr<std::vector<Real> > j_;
82 ROL::Ptr<ROL::Objective<Real> > obj_;
83 
84 public:
85 
86  CLExactModel( ROL::Ptr<std::vector<Real> > &xp,
87  const ROL::Ptr<const std::vector<Real> > &lp,
88  const ROL::Ptr<const std::vector<Real> > &up ) :
89  x_(xp), l_(lp), u_(up) {
90  g_ = ROL::makePtr<std::vector<double>(x_->size>());
91  di_ = ROL::makePtr<std::vector<double>(x_->size>());
92  j_ = ROL::makePtr<std::vector<double>(x_->size>());
93 
94 
95  obj_ = ROL::makePtr<CLTestObjective<Real>>();
96 
99  Real tol = std::sqrt(ROL::ROL_EPSILON<Real>());
100  obj_->gradient(g,x,tol);
101 
102  std::vector<Real> v(2);
103 
104  for(int i=0; i<2;++i) {
105  (*j_)[i] = 0;
106  // Case (i)
107  if( (*g_)[i]<0 && (*u_)[i] < ROL::ROL_INF<Real>() ) {
108  v[i] = (*u_)[i]-(*x_)[i];
109  (*j_)[i] = -1;
110  }
111  // Case (ii)
112  else if( (*g_)[i]>=0 && (*l_)[i] > ROL::ROL_NINF<Real>() ) {
113  v[i] = (*x_)[i] - (*l_)[i];
114  (*j_)[i] = 1;
115  }
116  // Case (iii)
117  else if( (*g_)[i]<=0 && (*u_)[i] == ROL::ROL_INF<Real>() ) {
118  v[i] = -1;
119  }
120  // Case (iv)
121  else {
122  v[i] = 1;
123  }
124  (*di_)[i] = std::sqrt(std::abs(v[i]));
125  }
126 
127 
128  std::cout << "x[0] = " << (*x_)[0] << std::endl;
129  std::cout << "x[1] = " << (*x_)[1] << std::endl;
130  std::cout << "g[0] = " << (*g_)[0] << std::endl;
131  std::cout << "g[0] = " << (*g_)[1] << std::endl;
132  std::cout << "di[0] = " << (*di_)[0] << std::endl;
133  std::cout << "di[1] = " << (*di_)[1] << std::endl;
134 
135  }
136 
137  void update( const ROL::Vector<Real> &x, bool flag = true, int iter=-1 ) {
138  ROL::Ptr<const std::vector<Real> > xc =
139  dynamic_cast<const ROL::StdVector<Real>&>(x).getVector();
140  (*x_)[0] = (*xc)[0];
141  (*x_)[1] = (*xc)[1];
142 
143  std::vector<Real> v(2);
145  Real tol = std::sqrt(ROL::ROL_EPSILON<Real>());
146  obj_->gradient(g,x,tol);
147 
148  for(int i=0; i<2;++i) {
149  (*j_)[i] = 0;
150  // Case (i)
151  if( (*g_)[i]<0 && (*u_)[i] < ROL::ROL_INF<Real>() ) {
152  v[i] = (*u_)[i]-(*x_)[i];
153  (*j_)[i] = -1;
154  }
155  // Case (ii)
156  else if( (*g_)[i]>=0 && (*l_)[i] > ROL::ROL_NINF<Real>() ) {
157  v[i] = (*x_)[i] - (*l_)[i];
158  (*j_)[i] = 1;
159  }
160  // Case (iii)
161  else if( (*g_)[i]<=0 && (*u_)[i] == ROL::ROL_INF<Real>() ) {
162  v[i] = -1;
163  }
164  // Case (iv)
165  else {
166  v[i] = 1;
167  }
168  (*di_)[i] = std::sqrt(std::abs(v[i]));
169  }
170 
171  std::cout << "x[0] = " << (*x_)[0] << std::endl;
172  std::cout << "x[1] = " << (*x_)[1] << std::endl;
173  std::cout << "g[0] = " << (*g_)[0] << std::endl;
174  std::cout << "g[0] = " << (*g_)[1] << std::endl;
175  std::cout << "di[0] = " << (*di_)[0] << std::endl;
176  std::cout << "di[1] = " << (*di_)[1] << std::endl;
177  }
178 
179  Real value( const ROL::Vector<Real> &s, Real &tol ) {
180  ROL::Ptr<const std::vector<Real> > sp =
181  dynamic_cast<const ROL::StdVector<Real>&>(s).getVector();
182 
183  ROL::Ptr<ROL::Vector<Real> > y = s.clone();
184  hessVec(*y,s,s,tol);
185  Real result = 0.5*y->dot(s);
186  result += (*di_)[0]*(*g_)[0]*(*sp)[0];
187  result += (*di_)[1]*(*g_)[1]*(*sp)[1];
188  return result;
189  }
190 
191  void gradient( ROL::Vector<Real> &g, const ROL::Vector<Real> &s, Real &tol ) {
192  ROL::Ptr<std::vector<Real> > gp =
193  dynamic_cast<ROL::StdVector<Real>&>(g).getVector();
194  hessVec(g,s,s,tol);
195 
196  (*gp)[0] += (*di_)[0]*(*g_)[0];
197  (*gp)[1] += (*di_)[1]*(*g_)[1];
198  }
199 
201  const ROL::Vector<Real> &v,
202  const ROL::Vector<Real> &s,
203  Real &tol ) {
204 
205  ROL::Ptr<std::vector<Real> > hvp =
206  dynamic_cast<ROL::StdVector<Real>&>(hv).getVector();
207  ROL::Ptr<const std::vector<Real> > vp =
208  dynamic_cast<const ROL::StdVector<Real>&>(v).getVector();
209 
210  obj_->hessVec(hv,v,s,tol);
211 
212  for(int i=0; i<2; ++i) {
213  (*hvp)[i] *= (*di_)[i]*(*di_)[i];
214  (*hvp)[i] += (*g_)[i]*(*j_)[i]*(*vp)[i];
215  }
216 
217  }
218 
219 
220 }; // CLExactModel
221 
222 
223 
224 
Real value(const ROL::Vector< Real > &x, Real &tol)
Compute value.
Definition: test_11.hpp:44
Provides the interface to evaluate objective functions.
virtual ROL::Ptr< Vector > clone() const =0
Clone to make a new (uninitialized) vector.
ROL::Ptr< std::vector< Real > > x_
Definition: test_11.hpp:76
void gradient(ROL::Vector< Real > &g, const ROL::Vector< Real > &x, Real &tol)
Compute gradient.
Definition: test_11.hpp:50
const ROL::Ptr< const std::vector< Real > > l_
Definition: test_11.hpp:77
Real value(const ROL::Vector< Real > &s, Real &tol)
Compute value.
Definition: test_11.hpp:179
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:46
ROL::Ptr< std::vector< Real > > di_
Definition: test_11.hpp:80
void hessVec(ROL::Vector< Real > &hv, const ROL::Vector< Real > &v, const ROL::Vector< Real > &s, Real &tol)
Apply Hessian approximation to vector.
Definition: test_11.hpp:200
ROL::Ptr< ROL::Objective< Real > > obj_
Definition: test_11.hpp:82
const ROL::Ptr< const std::vector< Real > > u_
Definition: test_11.hpp:78
ROL::Ptr< std::vector< Real > > j_
Definition: test_11.hpp:81
void gradient(ROL::Vector< Real > &g, const ROL::Vector< Real > &s, Real &tol)
Compute gradient.
Definition: test_11.hpp:191
void update(const ROL::Vector< Real > &x, bool flag=true, int iter=-1)
Update objective function.
Definition: test_11.hpp:137
CLExactModel(ROL::Ptr< std::vector< Real > > &xp, const ROL::Ptr< const std::vector< Real > > &lp, const ROL::Ptr< const std::vector< Real > > &up)
Definition: test_11.hpp:86
void hessVec(ROL::Vector< Real > &hv, const ROL::Vector< Real > &v, const ROL::Vector< Real > &x, Real &tol)
Apply Hessian approximation to vector.
Definition: test_11.hpp:59
ROL::Ptr< std::vector< Real > > g_
Definition: test_11.hpp:79