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