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 test_11.hpp
45     \brief Verify that the Coleman-Li Model function produces the
46            correct values for a test problem
47 
48     \f[ \min_x f(x) = \frac{1}{2}(x_1^2+2x_2^2),\quad x_1 \geq 1,\; x_2 <= -1 \f]
49 
50     The gradient is
51     \f[  \nabla f(x) = (x_1,2 x_2 ) \f]
52 
53     and the Hessian is
54     \[f \nabla^2 f(x) = \begin{pmatrix} 1 & 0 \\ 0 & 2 \end{pmatrix} \f]
55 
56     The minimizer is \f$ x^\ast = (1,-1)\f$
57 
58     For feasible \f$x\f$, the Coleman-Li quantities of interest are
59 
60     \f[ v_1 = x_1 - 1,\quad v_2 = -1 - x_2 \f]
61     \f[ D^{-1} = \begin{pmatrix} \sqrt{|x_1-1|} & 0 \\ 0 & \sqrt{|x_2+1|} \end{pmatrix} \f]
62     \f[ J=\begin{pmatrix} 1 & -1 \end{pmatrix} \f]
63     \f[ \hat M_k = \begin{pmatrix} |x_1-1|^2+|x_1| & 0  \\
64                                    0 & |x_2+1|^2+2|x_2| \end{pmatrix} \f]
65 
66 
67 */
68 
69 #include "ROL_Objective.hpp"
70 #include "ROL_StdVector.hpp"
71 
72 template<class Real>
73 class CLTestObjective : public ROL::Objective<Real> {
74 
75 
76 public:
77 
value(const ROL::Vector<Real> & x,Real & tol)78   Real value( const ROL::Vector<Real> &x, Real &tol ) {
79     ROL::Ptr<const std::vector<Real> > xp =
80       dynamic_cast<const ROL::StdVector<Real>&>(x).getVector();
81       return 0.5*((*xp)[0]*(*xp)[0] + 2*(*xp)[1]*(*xp)[1]);
82   }
83 
gradient(ROL::Vector<Real> & g,const ROL::Vector<Real> & x,Real & tol)84   void gradient( ROL::Vector<Real> &g, const ROL::Vector<Real> &x, Real &tol ) {
85     ROL::Ptr<std::vector<Real> > gp =
86       dynamic_cast<ROL::StdVector<Real>&>(g).getVector();
87     ROL::Ptr<const std::vector<Real> > xp =
88       dynamic_cast<const ROL::StdVector<Real>&>(x).getVector();
89     (*gp)[0] =   (*xp)[0];
90     (*gp)[1] = 2*(*xp)[1];
91   }
92 
hessVec(ROL::Vector<Real> & hv,const ROL::Vector<Real> & v,const ROL::Vector<Real> & x,Real & tol)93   void hessVec( ROL::Vector<Real> &hv,
94                 const ROL::Vector<Real> &v,
95                 const ROL::Vector<Real> &x,
96                 Real &tol ) {
97     ROL::Ptr<std::vector<Real> > hvp =
98       dynamic_cast<ROL::StdVector<Real>&>(hv).getVector();
99     ROL::Ptr<const std::vector<Real> > vp =
100       dynamic_cast<const ROL::StdVector<Real>&>(v).getVector();
101     (*hvp)[0] =   (*vp)[0];
102     (*hvp)[1] = 2*(*vp)[1];
103   }
104 
105 }; // CLTestObjective
106 
107 template<class Real>
108 class CLExactModel : public ROL::Objective<Real> {
109 
110 ROL::Ptr<std::vector<Real> > x_;
111 const ROL::Ptr<const std::vector<Real> > l_;
112 const ROL::Ptr<const std::vector<Real> > u_;
113 ROL::Ptr<std::vector<Real> > g_;
114 ROL::Ptr<std::vector<Real> > di_;
115 ROL::Ptr<std::vector<Real> > j_;
116 ROL::Ptr<ROL::Objective<Real> > obj_;
117 
118 public:
119 
CLExactModel(ROL::Ptr<std::vector<Real>> & xp,const ROL::Ptr<const std::vector<Real>> & lp,const ROL::Ptr<const std::vector<Real>> & up)120   CLExactModel( ROL::Ptr<std::vector<Real> > &xp,
121                 const ROL::Ptr<const std::vector<Real> > &lp,
122                 const ROL::Ptr<const std::vector<Real> > &up ) :
123                 x_(xp), l_(lp), u_(up) {
124     g_  = ROL::makePtr<std::vector<double>(x_->size>());
125     di_ = ROL::makePtr<std::vector<double>(x_->size>());
126     j_  = ROL::makePtr<std::vector<double>(x_->size>());
127 
128 
129     obj_ = ROL::makePtr<CLTestObjective<Real>>();
130 
131     ROL::StdVector<Real> g(g_);
132     ROL::StdVector<Real> x(x_);
133     Real tol = std::sqrt(ROL::ROL_EPSILON<Real>());
134     obj_->gradient(g,x,tol);
135 
136     std::vector<Real> v(2);
137 
138     for(int i=0; i<2;++i) {
139       (*j_)[i] = 0;
140       // Case (i)
141       if( (*g_)[i]<0 && (*u_)[i] < ROL::ROL_INF<Real>() ) {
142         v[i] = (*u_)[i]-(*x_)[i];
143         (*j_)[i] = -1;
144       }
145       // Case (ii)
146       else if( (*g_)[i]>=0 && (*l_)[i] > ROL::ROL_NINF<Real>() ) {
147          v[i] = (*x_)[i] - (*l_)[i];
148          (*j_)[i] = 1;
149       }
150       // Case (iii)
151       else if( (*g_)[i]<=0 && (*u_)[i] == ROL::ROL_INF<Real>() ) {
152         v[i] = -1;
153       }
154       // Case (iv)
155       else {
156         v[i] = 1;
157       }
158       (*di_)[i] = std::sqrt(std::abs(v[i]));
159     }
160 
161 
162     std::cout << "x[0]  = " << (*x_)[0] << std::endl;
163     std::cout << "x[1]  = " << (*x_)[1] << std::endl;
164     std::cout << "g[0]  = " << (*g_)[0] << std::endl;
165     std::cout << "g[0]  = " << (*g_)[1] << std::endl;
166     std::cout << "di[0] = " << (*di_)[0] << std::endl;
167     std::cout << "di[1] = " << (*di_)[1] << std::endl;
168 
169   }
170 
update(const ROL::Vector<Real> & x,bool flag=true,int iter=-1)171   void update( const ROL::Vector<Real> &x, bool flag = true, int iter=-1 ) {
172     ROL::Ptr<const std::vector<Real> > xc =
173       dynamic_cast<const ROL::StdVector<Real>&>(x).getVector();
174     (*x_)[0]  = (*xc)[0];
175     (*x_)[1]  = (*xc)[1];
176 
177     std::vector<Real> v(2);
178     ROL::StdVector<Real> g(g_);
179     Real tol = std::sqrt(ROL::ROL_EPSILON<Real>());
180     obj_->gradient(g,x,tol);
181 
182     for(int i=0; i<2;++i) {
183       (*j_)[i] = 0;
184       // Case (i)
185       if( (*g_)[i]<0 && (*u_)[i] < ROL::ROL_INF<Real>() ) {
186         v[i] = (*u_)[i]-(*x_)[i];
187         (*j_)[i] = -1;
188       }
189       // Case (ii)
190       else if( (*g_)[i]>=0 && (*l_)[i] > ROL::ROL_NINF<Real>() ) {
191          v[i] = (*x_)[i] - (*l_)[i];
192          (*j_)[i] = 1;
193       }
194       // Case (iii)
195       else if( (*g_)[i]<=0 && (*u_)[i] == ROL::ROL_INF<Real>() ) {
196         v[i] = -1;
197       }
198       // Case (iv)
199       else {
200         v[i] = 1;
201       }
202       (*di_)[i] = std::sqrt(std::abs(v[i]));
203     }
204 
205     std::cout << "x[0]  = " << (*x_)[0] << std::endl;
206     std::cout << "x[1]  = " << (*x_)[1] << std::endl;
207     std::cout << "g[0]  = " << (*g_)[0] << std::endl;
208     std::cout << "g[0]  = " << (*g_)[1] << std::endl;
209     std::cout << "di[0] = " << (*di_)[0] << std::endl;
210     std::cout << "di[1] = " << (*di_)[1] << std::endl;
211   }
212 
value(const ROL::Vector<Real> & s,Real & tol)213   Real value( const ROL::Vector<Real> &s, Real &tol ) {
214     ROL::Ptr<const std::vector<Real> > sp =
215       dynamic_cast<const ROL::StdVector<Real>&>(s).getVector();
216 
217     ROL::Ptr<ROL::Vector<Real> > y = s.clone();
218     hessVec(*y,s,s,tol);
219     Real result = 0.5*y->dot(s);
220     result += (*di_)[0]*(*g_)[0]*(*sp)[0];
221     result += (*di_)[1]*(*g_)[1]*(*sp)[1];
222     return result;
223   }
224 
gradient(ROL::Vector<Real> & g,const ROL::Vector<Real> & s,Real & tol)225   void gradient( ROL::Vector<Real> &g, const ROL::Vector<Real> &s, Real &tol ) {
226     ROL::Ptr<std::vector<Real> > gp =
227       dynamic_cast<ROL::StdVector<Real>&>(g).getVector();
228     hessVec(g,s,s,tol);
229 
230     (*gp)[0] += (*di_)[0]*(*g_)[0];
231     (*gp)[1] += (*di_)[1]*(*g_)[1];
232   }
233 
hessVec(ROL::Vector<Real> & hv,const ROL::Vector<Real> & v,const ROL::Vector<Real> & s,Real & tol)234   void hessVec( ROL::Vector<Real> &hv,
235                 const ROL::Vector<Real> &v,
236                 const ROL::Vector<Real> &s,
237                 Real &tol ) {
238 
239     ROL::Ptr<std::vector<Real> > hvp =
240       dynamic_cast<ROL::StdVector<Real>&>(hv).getVector();
241     ROL::Ptr<const std::vector<Real> > vp =
242       dynamic_cast<const ROL::StdVector<Real>&>(v).getVector();
243 
244     obj_->hessVec(hv,v,s,tol);
245 
246     for(int i=0; i<2; ++i) {
247       (*hvp)[i] *=  (*di_)[i]*(*di_)[i];
248       (*hvp)[i] += (*g_)[i]*(*j_)[i]*(*vp)[i];
249     }
250 
251   }
252 
253 
254 }; // CLExactModel
255 
256 
257 
258 
259