ROL
step/interiorpoint/test_02.cpp
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
44#define OPTIMIZATION_PROBLEM_REFACTOR
45
46#include "Teuchos_GlobalMPISession.hpp"
47
48#include "ROL_RandomVector.hpp"
49#include "ROL_StdVector.hpp"
50#include "ROL_NonlinearProgram.hpp"
54
55#include "HS_Problem_041.hpp"
56
57#include <iomanip>
58
69template<class Real>
70void printVector( const ROL::Vector<Real> &x, std::ostream &outStream ) {
71
72 try {
73 ROL::Ptr<const std::vector<Real> > xp =
74 dynamic_cast<const ROL::StdVector<Real>&>(x).getVector();
75
76 outStream << "Standard Vector" << std::endl;
77 for( size_t i=0; i<xp->size(); ++i ) {
78 outStream << (*xp)[i] << std::endl;
79 }
80 }
81 catch( const std::bad_cast& e ) {
82 outStream << "Partitioned Vector" << std::endl;
83
85 typedef typename PV::size_type size_type;
86
87 const PV &xpv = dynamic_cast<const PV&>(x);
88
89 for( size_type i=0; i<xpv.numVectors(); ++i ) {
90 outStream << "--------------------" << std::endl;
91 printVector( *(xpv.get(i)), outStream );
92 }
93 outStream << "--------------------" << std::endl;
94 }
95
96}
97
98
99// Exact residual for H&S Problem 41
100template<class Real>
101void value( ROL::Vector<Real> &c, const ROL::Vector<Real> &sol, const Real &mu ) {
102
103 typedef std::vector<Real> vector;
104 typedef ROL::StdVector<Real> SV;
106
107 typedef typename PV::size_type size_type;
108
109
110
111 using ROL::dynamicPtrCast;
112
113 const size_type OPT = 0;
114 const size_type EQUAL = 1;
115 const size_type LOWER = 2;
116 const size_type UPPER = 3;
117
118 const PV &sol_pv = dynamic_cast<const PV&>(sol);
119 const vector &x = *(ROL::dynamicPtrCast<const SV>(sol_pv.get(OPT))->getVector());
120 const vector &l = *(ROL::dynamicPtrCast<const SV>(sol_pv.get(EQUAL))->getVector());
121 const vector &zl = *(ROL::dynamicPtrCast<const SV>(sol_pv.get(LOWER))->getVector());
122 const vector &zu = *(ROL::dynamicPtrCast<const SV>(sol_pv.get(UPPER))->getVector());
123
124 PV &c_pv = dynamic_cast<PV&>(c);
125 vector &cx = *(ROL::dynamicPtrCast<SV>(c_pv.get(OPT))->getVector());
126 vector &cl = *(ROL::dynamicPtrCast<SV>(c_pv.get(EQUAL))->getVector());
127 vector &czl = *(ROL::dynamicPtrCast<SV>(c_pv.get(LOWER))->getVector());
128 vector &czu = *(ROL::dynamicPtrCast<SV>(c_pv.get(UPPER))->getVector());
129
130 cx[0] = -x[1]*x[2] + l[0] - zl[0] + zu[0];
131 cx[1] = -x[0]*x[1] + 2*l[0] - zl[1] + zu[1];
132 cx[2] = -x[0]*x[1] + 2*l[0] - zl[2] + zu[2];
133 cx[3] = - l[0] - zl[3] + zu[3];
134
135 cl[0] = x[0] + 2*x[1] + 2*x[2] - x[3];
136
137 czl[0] = x[0]*zl[0] - mu;
138 czl[1] = x[1]*zl[1] - mu;
139 czl[2] = x[2]*zl[2] - mu;
140 czl[3] = x[3]*zl[3] - mu;
141
142 czu[0] = (1.0-x[0])*zu[0] - mu;
143 czu[1] = (1.0-x[1])*zu[1] - mu;
144 czu[2] = (1.0-x[2])*zl[2] - mu;
145 czu[3] = (2.0-x[3])*zl[3] - mu;
146}
147
148
149
150// Exact residual for H&S Problem 41
151template<class Real>
153
154 typedef std::vector<Real> vector;
155 typedef ROL::StdVector<Real> SV;
157
158 typedef typename PV::size_type size_type;
159
160
161
162 using ROL::dynamicPtrCast;
163
164 const size_type OPT = 0;
165 const size_type EQUAL = 1;
166 const size_type LOWER = 2;
167 const size_type UPPER = 3;
168
169 const PV &sol_pv = dynamic_cast<const PV&>(sol);
170 const vector &x = *(ROL::dynamicPtrCast<const SV>(sol_pv.get(OPT))->getVector());
171//const vector &l = *(ROL::dynamicPtrCast<const SV>(sol_pv.get(EQUAL))->getVector());
172 const vector &zl = *(ROL::dynamicPtrCast<const SV>(sol_pv.get(LOWER))->getVector());
173 const vector &zu = *(ROL::dynamicPtrCast<const SV>(sol_pv.get(UPPER))->getVector());
174
175 const PV &v_pv = dynamic_cast<const PV&>(v);
176 const vector &vx = *(ROL::dynamicPtrCast<const SV>(v_pv.get(OPT))->getVector());
177 const vector &vl = *(ROL::dynamicPtrCast<const SV>(v_pv.get(EQUAL))->getVector());
178 const vector &vzl = *(ROL::dynamicPtrCast<const SV>(v_pv.get(LOWER))->getVector());
179 const vector &vzu = *(ROL::dynamicPtrCast<const SV>(v_pv.get(UPPER))->getVector());
180
181 PV &jv_pv = dynamic_cast<PV&>(jv);
182 vector &jvx = *(ROL::dynamicPtrCast<SV>(jv_pv.get(OPT))->getVector());
183 vector &jvl = *(ROL::dynamicPtrCast<SV>(jv_pv.get(EQUAL))->getVector());
184 vector &jvzl = *(ROL::dynamicPtrCast<SV>(jv_pv.get(LOWER))->getVector());
185 vector &jvzu = *(ROL::dynamicPtrCast<SV>(jv_pv.get(UPPER))->getVector());
186
187 jvx[0] = -x[1]*vx[2] - x[2]*vx[1] + vl[0] - vzl[0] + vzu[0];
188 jvx[1] = -x[0]*vx[2] - x[2]*vx[0] + 2*vl[0] - vzl[1] + vzu[1];
189 jvx[2] = -x[0]*vx[1] - x[1]*vx[0] + 2*vl[0] - vzl[2] + vzu[2];
190 jvx[3] = - vl[0] - vzl[3] + vzu[3];
191
192 jvl[0] = vx[0] + 2*vx[1] + 2*vx[2] - vx[3];
193
194 jvzl[0] = zl[0]*vx[0] + vzl[0]*x[0];
195 jvzl[1] = zl[1]*vx[1] + vzl[1]*x[1];
196 jvzl[2] = zl[2]*vx[2] + vzl[2]*x[2];
197 jvzl[3] = zl[3]*vx[3] + vzl[3]*x[3];
198
199 jvzu[0] = -zu[0]*vx[0] + vzu[0]*(1.0-x[0]);
200 jvzu[1] = -zu[1]*vx[1] + vzu[1]*(1.0-x[1]);
201 jvzu[2] = -zu[2]*vx[2] + vzu[2]*(1.0-x[2]);
202 jvzu[3] = -zu[3]*vx[3] + vzu[3]*(2.0-x[3]);
203
204}
205
206
207typedef double RealT;
208
209int main(int argc, char *argv[]) {
210
211// typedef std::vector<RealT> vector;
212
213 typedef ROL::ParameterList PL;
214
215 typedef ROL::Vector<RealT> V;
217 typedef ROL::Objective<RealT> OBJ;
218 typedef ROL::Constraint<RealT> CON;
219 typedef ROL::BoundConstraint<RealT> BND;
221 typedef ROL::NonlinearProgram<RealT> NLP;
222
223 typedef ROL::InteriorPointPenalty<RealT> PENALTY;
225
226
227
228 Teuchos::GlobalMPISession mpiSession(&argc, &argv);
229
230 int iprint = argc - 1;
231 ROL::Ptr<std::ostream> outStream;
232 ROL::nullstream bhs;
233 if( iprint > 0 )
234 outStream = ROL::makePtrFromRef(std::cout);
235 else
236 outStream = ROL::makePtrFromRef(bhs);
237
238 int errorFlag = 0;
239
240 try {
241
242 RealT mu = 0.1;
243
244 RealT tol = std::sqrt(ROL::ROL_EPSILON<RealT>());
245
246 PL parlist;
247
248 PL &iplist = parlist.sublist("Step").sublist("Primal Dual Interior Point");
249 PL &lblist = iplist.sublist("Barrier Objective");
250
251 iplist.set("Symmetrize Primal Dual System",false);
252
253 lblist.set("Use Linear Damping", true);
254 lblist.set("Linear Damping Coefficient",1.e-4);
255 lblist.set("Initial Barrier Parameter", mu);
256
257 // Need an example problem that satisfies the following criteria:
258 // 1) Has an equality constraint
259 // 2) Has a bound constraint where all variables have finite upper and lower bounds
260
261 ROL::Ptr<NLP> nlp = ROL::makePtr<HS::Problem_041<RealT>>();
262 ROL::Ptr<OPT> opt = nlp->getOptimizationProblem();
263
264 ROL::Ptr<V> x = opt->getSolutionVector();
265 ROL::Ptr<V> l = opt->getMultiplierVector();
266 ROL::Ptr<V> zl = x->clone();
267 ROL::Ptr<V> zu = x->clone();
268
269 ROL::Ptr<V> scratch = x->clone();
270
271 ROL::Ptr<PV> x_pv = ROL::dynamicPtrCast<PV>(x);
272
273 ROL::Ptr<V> sol = CreatePartitionedVector(x,l,zl,zu);
274
275 ROL::Ptr<V> c = sol->clone();
276 ROL::Ptr<V> v = sol->clone();
277 ROL::Ptr<V> jv = sol->clone();
278
279 ROL::Ptr<V> c_exact = c->clone();
280 ROL::Ptr<V> jv_exact = jv->clone();
281
282 ROL::RandomizeVector(*l, -1.0, 1.0);
283 ROL::RandomizeVector(*v, 0.0, 1.0);
284
285
286 ROL::Ptr<OBJ> obj = opt->getObjective();
287 ROL::Ptr<CON> con = opt->getConstraint();
288 ROL::Ptr<BND> bnd = opt->getBoundConstraint();
289
290 PENALTY penalty(obj,bnd,parlist);
291
292 ROL::Ptr<const V> maskL = penalty.getLowerMask();
293 ROL::Ptr<const V> maskU = penalty.getUpperMask();
294
295 zl->set(*maskL);
296 zu->set(*maskU);
297
298 ROL::Ptr<CON> res = ROL::makePtr<RESIDUAL>(obj,con,bnd,*sol,maskL,maskU,scratch,mu,false);
299
300
301 *outStream << "\n[x|lambda|zl|zu]" << std::endl;
302 printVector(*sol,*outStream);
303
304 res->value(*c,*sol,tol);
305 value(*c_exact,*sol,mu);
306
307 *outStream << "\n[cx|clambda|czl|czu]" << std::endl;
308
309 printVector(*c,*outStream);
310
311 c->axpy(-1.0,*c_exact);
312
313 RealT cerror = c->norm();
314
315 if( cerror>tol ) {
316 ++errorFlag;
317 }
318
319 *outStream << "\n\n||c-c_exact|| = " << cerror << std::endl;
320
321 res->applyJacobian(*jv,*v,*sol,tol);
322 applyJacobian(*jv_exact,*v,*sol);
323
324 *outStream << "\n[vx|vlambda|vzl|vzu]" << std::endl;
325
326 printVector(*v,*outStream);
327
328 *outStream << "\n[jvx|jvlambda|jvzl|jvzu]" << std::endl;
329
330 printVector(*jv,*outStream);
331
332 jv->axpy(-1.0,*jv_exact);
333
334 RealT jverror = jv->norm();
335
336 if( jverror > tol ) {
337 ++errorFlag;
338 }
339
340 *outStream << "\n\n||jv-jv_exact|| = " << jverror << std::endl;
341
342 *outStream << "Residual Jacobian Finite Difference Check" << std::endl;
343 res->checkApplyJacobian(*sol,*v,*v);
344
345
346
347 }
348 catch (std::logic_error& err) {
349 *outStream << err.what() << std::endl;
350 errorFlag = -1000;
351 }
352
353 if (errorFlag != 0)
354 std::cout << "End Result: TEST FAILED" << std::endl;
355 else
356 std::cout << "End Result: TEST PASSED" << std::endl;
357
358 return 0;
359}
360
Vector< Real > V
PartitionedVector< Real > PV
typename PV< Real >::size_type size_type
Provides the interface to apply upper and lower bound constraints.
Defines the general constraint operator interface.
Provides the interface to evaluate the Interior Pointy log barrier penalty function with upper and lo...
Provides the interface to evaluate objective functions.
Defines the linear algebra of vector space on a generic partitioned vector.
Symmetrized form of the KKT operator for the Type-EB problem with equality and bound multipliers.
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
void RandomizeVector(Vector< Real > &x, const Real &lower=0.0, const Real &upper=1.0)
Fill a ROL::Vector with uniformly-distributed random numbers in the interval [lower,...
int main(int argc, char *argv[])
void applyJacobian(ROL::Vector< Real > &jv, const ROL::Vector< Real > &v, const ROL::Vector< Real > &sol)
void value(ROL::Vector< Real > &c, const ROL::Vector< Real > &sol, const Real &mu)
void printVector(const ROL::Vector< Real > &x, std::ostream &outStream)