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 
69 template<class Real>
70 void 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
100 template<class Real>
101 void 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
151 template<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 
207 typedef double RealT;
208 
209 int 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 
Provides the interface to evaluate objective functions.
PartitionedVector< Real > PV
typename PV< Real >::size_type size_type
Defines the linear algebra of vector space on a generic partitioned vector.
ROL::Ptr< Vector< Real > > CreatePartitionedVector(const ROL::Ptr< Vector< Real >> &a)
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,upper].
Symmetrized form of the KKT operator for the Type-EB problem with equality and bound multipliers...
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:80
void applyJacobian(ROL::Vector< Real > &jv, const ROL::Vector< Real > &v, const ROL::Vector< Real > &sol)
Vector< Real > V
void printVector(const ROL::Vector< Real > &x, std::ostream &outStream)
Provides the interface to evaluate the Interior Pointy log barrier penalty function with upper and lo...
Provides the interface to apply upper and lower bound constraints.
void value(ROL::Vector< Real > &c, const ROL::Vector< Real > &sol, const Real &mu)
basic_nullstream< char, char_traits< char >> nullstream
Definition: ROL_Stream.hpp:72
int main(int argc, char *argv[])
Defines the general constraint operator interface.