ROL
function/test_14.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 
49 #include "ROL_HS41.hpp"
50 #include "ROL_HS53.hpp"
51 #include "ROL_HS55.hpp"
52 #include "ROL_Bounds.hpp"
54 
55 #include "ROL_Stream.hpp"
56 #include "Teuchos_GlobalMPISession.hpp"
57 
58 typedef double RealT;
59 
60 int main(int argc, char *argv[]) {
61 
62  Teuchos::GlobalMPISession mpiSession(&argc, &argv);
63 
64  // This little trick lets us print to std::cout only if a
65  // (dummy) command-line argument is provided.
66  int iprint = argc - 1;
67  ROL::Ptr<std::ostream> outStream;
68  ROL::nullstream bhs; // outputs nothing
69  if (iprint > 0)
70  outStream = ROL::makePtrFromRef(std::cout);
71  else
72  outStream = ROL::makePtrFromRef(bhs);
73 
74  int errorFlag = 0;
75 
76  try {
77  RealT tol = std::sqrt(ROL::ROL_EPSILON<RealT>());
78  RealT cnorm(0);
79  ROL::Ptr<ROL::Vector<RealT>> sol, mul, x, lam, l, u, c;
80  ROL::Ptr<ROL::Objective<RealT>> obj;
81  ROL::Ptr<ROL::Constraint<RealT>> con;
82  ROL::Ptr<ROL::BoundConstraint<RealT>> bnd;
83  ROL::Ptr<ROL::PolyhedralProjection<RealT>> proj;
84  ROL::ParameterList list;
85  list.sublist("General").set("Output Level",2);
86  std::vector<RealT> data;
87 
88  *outStream << std::endl << "Hock and Schittkowski Problem #41" << std::endl << std::endl;
90  obj = HS41.getObjective();
91  sol = HS41.getInitialGuess();
92  con = HS41.getEqualityConstraint();
93  mul = HS41.getEqualityMultiplier();
94  bnd = HS41.getBoundConstraint();
95 
96  lam = mul->clone(); lam->set(*mul);
97  x = sol->clone(); x->set(*sol);
98  l = sol->clone(); l->zero();
99  u = sol->clone(); u->setScalar(static_cast<RealT>(1));
100  c = mul->dual().clone();
101 
102  list.sublist("General").sublist("Polyhedral Projection").set("Type","Dai-Fletcher");
103  proj = ROL::PolyhedralProjectionFactory<RealT>(*sol,sol->dual(),bnd,con,*lam,*c,list);
104  proj->project(*x,*outStream);
105 
106  con->value(*c,*x,tol);
107  cnorm = c->norm();
108 
109  data = *ROL::staticPtrCast<ROL::StdVector<RealT>>(sol)->getVector();
110  *outStream << " Initial: x1 = " << data[0] << " x2 = " << data[1]
111  << " x3 = " << data[2] << std::endl;
112  data = *ROL::staticPtrCast<ROL::StdVector<RealT>>(x)->getVector();
113  *outStream << " Result: x1 = " << data[0] << " x2 = " << data[1]
114  << " x3 = " << data[2] << std::endl;
115  data = *ROL::staticPtrCast<ROL::StdVector<RealT>>(lam)->getVector();
116  *outStream << " Multiplier: l1 = " << data[0] << std::endl;
117 
118  *outStream << std::endl;
119  *outStream << " is equality feasible = " << (cnorm<=tol) << std::endl
120  << " are bounds feasible = " << bnd->isFeasible(*x) << std::endl;
121 
122  errorFlag += !bnd->isFeasible(*x);
123  errorFlag += (cnorm > tol);
124 
125  *outStream << std::endl << "Hock and Schittkowski Problem #41" << std::endl << std::endl;
127  obj = HS41a.getObjective();
128  sol = HS41a.getInitialGuess();
129  con = HS41a.getEqualityConstraint();
130  mul = HS41a.getEqualityMultiplier();
131  bnd = HS41a.getBoundConstraint();
132 
133  lam = mul->clone(); lam->set(*mul);
134  x = sol->clone(); x->set(*sol);
135  l = sol->clone(); l->zero();
136  u = sol->clone(); u->setScalar(static_cast<RealT>(1));
137  c = mul->dual().clone();
138 
139  list.sublist("General").sublist("Polyhedral Projection").set("Type","Ridders");
140  proj = ROL::PolyhedralProjectionFactory<RealT>(*sol,sol->dual(),bnd,con,*lam,*c,list);
141  proj->project(*x,*outStream);
142 
143  con->value(*c,*x,tol);
144  cnorm = c->norm();
145 
146  data = *ROL::staticPtrCast<ROL::StdVector<RealT>>(sol)->getVector();
147  *outStream << " Initial: x1 = " << data[0] << " x2 = " << data[1]
148  << " x3 = " << data[2] << std::endl;
149  data = *ROL::staticPtrCast<ROL::StdVector<RealT>>(x)->getVector();
150  *outStream << " Result: x1 = " << data[0] << " x2 = " << data[1]
151  << " x3 = " << data[2] << std::endl;
152  data = *ROL::staticPtrCast<ROL::StdVector<RealT>>(lam)->getVector();
153  *outStream << " Multiplier: l1 = " << data[0] << std::endl;
154 
155  *outStream << std::endl;
156  *outStream << " is equality feasible = " << (cnorm<=tol) << std::endl
157  << " are bounds feasible = " << bnd->isFeasible(*x) << std::endl;
158 
159  errorFlag += !bnd->isFeasible(*x);
160  errorFlag += (cnorm > tol);
161 
162  *outStream << std::endl << "Hock and Schittkowski Problem #41" << std::endl << std::endl;
164  obj = HS41b.getObjective();
165  sol = HS41b.getInitialGuess();
166  con = HS41b.getEqualityConstraint();
167  mul = HS41b.getEqualityMultiplier();
168  bnd = HS41b.getBoundConstraint();
169 
170  lam = mul->clone(); lam->set(*mul);
171  x = sol->clone(); x->set(*sol);
172  l = sol->clone(); l->zero();
173  u = sol->clone(); u->setScalar(static_cast<RealT>(1));
174  c = mul->dual().clone();
175 
176  list.sublist("General").sublist("Polyhedral Projection").set("Type","Brents");
177  proj = ROL::PolyhedralProjectionFactory<RealT>(*sol,sol->dual(),bnd,con,*lam,*c,list);
178  proj->project(*x,*outStream);
179 
180  con->value(*c,*x,tol);
181  cnorm = c->norm();
182 
183  data = *ROL::staticPtrCast<ROL::StdVector<RealT>>(sol)->getVector();
184  *outStream << " Initial: x1 = " << data[0] << " x2 = " << data[1]
185  << " x3 = " << data[2] << std::endl;
186  data = *ROL::staticPtrCast<ROL::StdVector<RealT>>(x)->getVector();
187  *outStream << " Result: x1 = " << data[0] << " x2 = " << data[1]
188  << " x3 = " << data[2] << std::endl;
189  data = *ROL::staticPtrCast<ROL::StdVector<RealT>>(lam)->getVector();
190  *outStream << " Multiplier: l1 = " << data[0] << std::endl;
191 
192  *outStream << std::endl;
193  *outStream << " is equality feasible = " << (cnorm<=tol) << std::endl
194  << " are bounds feasible = " << bnd->isFeasible(*x) << std::endl;
195 
196  errorFlag += !bnd->isFeasible(*x);
197  errorFlag += (cnorm > tol);
198 
199  *outStream << std::endl << "Hock and Schittkowski Problem #53" << std::endl << std::endl;
201  obj = HS53.getObjective();
202  sol = HS53.getInitialGuess();
203  con = HS53.getEqualityConstraint();
204  mul = HS53.getEqualityMultiplier();
205  bnd = HS53.getBoundConstraint();
206 
207  lam = mul->clone(); lam->set(*mul);
208  x = sol->clone(); x->set(*sol);
209  l = sol->clone(); l->zero();
210  u = sol->clone(); u->setScalar(static_cast<RealT>(1));
211  c = mul->dual().clone();
212 
213  list.sublist("General").sublist("Polyhedral Projection").set("Type","Dykstra");
214  proj = ROL::PolyhedralProjectionFactory<RealT>(*sol,sol->dual(),bnd,con,*lam,*c,list);
215  proj->project(*x,*outStream);
216 
217  con->value(*c,*x,tol);
218  cnorm = c->norm();
219 
220  data = *ROL::staticPtrCast<ROL::StdVector<RealT>>(sol)->getVector();
221  *outStream << " Initial: x1 = " << data[0] << " x2 = " << data[1]
222  << " x3 = " << data[2] << " x4 = " << data[3]
223  << " x5 = " << data[4] << std::endl;
224  data = *ROL::staticPtrCast<ROL::StdVector<RealT>>(x)->getVector();
225  *outStream << " Result: x1 = " << data[0] << " x2 = " << data[1]
226  << " x3 = " << data[2] << " x4 = " << data[3]
227  << " x5 = " << data[4] << std::endl;
228  data = *ROL::staticPtrCast<ROL::StdVector<RealT>>(lam)->getVector();
229  *outStream << " Multiplier: l1 = " << data[0] << " l2 = " << data[1]
230  << " l3 = " << data[2] << std::endl;
231 
232  *outStream << std::endl;
233  *outStream << " is equality feasible = " << (cnorm<=tol) << std::endl
234  << " are bounds feasible = " << bnd->isFeasible(*x) << std::endl;
235 
236  errorFlag += !bnd->isFeasible(*x);
237  errorFlag += (cnorm > tol);
238 
239  *outStream << std::endl << "Hock and Schittkowski Problem #53" << std::endl << std::endl;
241  obj = HS53a.getObjective();
242  sol = HS53a.getInitialGuess();
243  con = HS53a.getEqualityConstraint();
244  mul = HS53a.getEqualityMultiplier();
245  bnd = HS53a.getBoundConstraint();
246 
247  lam = mul->clone(); lam->set(*mul);
248  x = sol->clone(); x->set(*sol);
249  l = sol->clone(); l->zero();
250  u = sol->clone(); u->setScalar(static_cast<RealT>(1));
251  c = mul->dual().clone();
252 
253  list.sublist("General").sublist("Polyhedral Projection").set("Type","Douglas-Rachford");
254  proj = ROL::PolyhedralProjectionFactory<RealT>(*sol,sol->dual(),bnd,con,*lam,*c,list);
255  proj->project(*x,*outStream);
256 
257  con->value(*c,*x,tol);
258  cnorm = c->norm();
259 
260  data = *ROL::staticPtrCast<ROL::StdVector<RealT>>(sol)->getVector();
261  *outStream << " Initial: x1 = " << data[0] << " x2 = " << data[1]
262  << " x3 = " << data[2] << " x4 = " << data[3]
263  << " x5 = " << data[4] << std::endl;
264  data = *ROL::staticPtrCast<ROL::StdVector<RealT>>(x)->getVector();
265  *outStream << " Result: x1 = " << data[0] << " x2 = " << data[1]
266  << " x3 = " << data[2] << " x4 = " << data[3]
267  << " x5 = " << data[4] << std::endl;
268  data = *ROL::staticPtrCast<ROL::StdVector<RealT>>(lam)->getVector();
269  *outStream << " Multiplier: l1 = " << data[0] << " l2 = " << data[1]
270  << " l3 = " << data[2] << std::endl;
271 
272  *outStream << std::endl;
273  *outStream << " is equality feasible = " << (cnorm<=tol) << std::endl
274  << " are bounds feasible = " << bnd->isFeasible(*x) << std::endl;
275 
276  errorFlag += !bnd->isFeasible(*x);
277  errorFlag += (cnorm > tol);
278 
279  *outStream << std::endl << "Hock and Schittkowski Problem #55" << std::endl << std::endl;
281  obj = HS55.getObjective();
282  sol = HS55.getInitialGuess();
283  con = HS55.getEqualityConstraint();
284  mul = HS55.getEqualityMultiplier();
285  bnd = HS55.getBoundConstraint();
286 
287  //ROL::Ptr<ROL::OptimizationProblem<RealT>> problem;
288  //ROL::Ptr<ROL::Vector<RealT>> xt;
289  //std::vector<ROL::Ptr<ROL::Vector<RealT>>> xv;
290  //HS55.get(problem,xt,xv);
291  //problem->check(*outStream);
292 
293  lam = mul->clone(); lam->set(*mul);
294  x = sol->clone(); x->set(*sol);
295  l = sol->clone(); l->zero();
296  u = sol->clone(); u->setScalar(static_cast<RealT>(1));
297  c = mul->dual().clone();
298 
299  list.sublist("General").sublist("Polyhedral Projection").set("Type","Semismooth Newton");
300  proj = ROL::PolyhedralProjectionFactory<RealT>(*sol,sol->dual(),bnd,con,*lam,*c,list);
301  proj->project(*x,*outStream);
302 
303  con->value(*c,*x,tol);
304  cnorm = c->norm();
305 
306  data = *ROL::staticPtrCast<ROL::StdVector<RealT>>(sol)->getVector();
307  *outStream << " Initial: x1 = " << data[0] << " x2 = " << data[1]
308  << " x3 = " << data[2] << " x4 = " << data[3]
309  << " x5 = " << data[4] << " x6 = " << data[5] << std::endl;
310  data = *ROL::staticPtrCast<ROL::StdVector<RealT>>(x)->getVector();
311  *outStream << " Result: x1 = " << data[0] << " x2 = " << data[1]
312  << " x3 = " << data[2] << " x4 = " << data[3]
313  << " x5 = " << data[4] << " x6 = " << data[5] << std::endl;
314  data = *ROL::staticPtrCast<ROL::StdVector<RealT>>(lam)->getVector();
315  *outStream << " Multiplier: l1 = " << data[0] << " l2 = " << data[1]
316  << " l3 = " << data[2] << " l4 = " << data[3]
317  << " l5 = " << data[4] << " l6 = " << data[5] << std::endl;
318 
319  *outStream << std::endl;
320  *outStream << " is equality feasible = " << (cnorm<=tol) << std::endl
321  << " are bounds feasible = " << bnd->isFeasible(*x) << std::endl;
322  *outStream << std::endl;
323 
324  errorFlag += !bnd->isFeasible(*x);
325  errorFlag += (cnorm > tol);
326  }
327 
328  catch (std::logic_error& err) {
329  *outStream << err.what() << "\n";
330  errorFlag = -1000;
331  }; // end try
332 
333  if (errorFlag != 0)
334  std::cout << "End Result: TEST FAILED\n";
335  else
336  std::cout << "End Result: TEST PASSED\n";
337 
338  return 0;
339 }
340 
Ptr< Vector< Real > > getInitialGuess(void) const
Definition: ROL_HS53.hpp:141
Contains definitions for W. Hock and K. Schittkowski 55th test function.
Ptr< BoundConstraint< Real > > getBoundConstraint(void) const
Definition: ROL_HS55.hpp:184
Ptr< Constraint< Real > > getEqualityConstraint(void) const
Definition: ROL_HS41.hpp:144
Contains definitions for W. Hock and K. Schittkowski 41th test function.
Contains definitions for W. Hock and K. Schittkowski 53th test function.
Defines a no-output stream class ROL::NullStream and a function makeStreamPtr which either wraps a re...
Ptr< Constraint< Real > > getEqualityConstraint(void) const
Definition: ROL_HS53.hpp:157
Ptr< BoundConstraint< Real > > getBoundConstraint(void) const
Definition: ROL_HS41.hpp:153
Ptr< Objective< Real > > getObjective(void) const
Definition: ROL_HS53.hpp:137
Ptr< Vector< Real > > getEqualityMultiplier(void) const
Definition: ROL_HS55.hpp:179
Ptr< Vector< Real > > getEqualityMultiplier(void) const
Definition: ROL_HS53.hpp:161
Ptr< Vector< Real > > getInitialGuess(void) const
Definition: ROL_HS55.hpp:151
Ptr< Objective< Real > > getObjective(void) const
Definition: ROL_HS41.hpp:127
basic_nullstream< char, char_traits< char >> nullstream
Definition: ROL_Stream.hpp:72
int main(int argc, char *argv[])
Ptr< Objective< Real > > getObjective(void) const
Definition: ROL_HS55.hpp:147
Ptr< BoundConstraint< Real > > getBoundConstraint(void) const
Definition: ROL_HS53.hpp:166
Ptr< Vector< Real > > getInitialGuess(void) const
Definition: ROL_HS41.hpp:131
Ptr< Vector< Real > > getEqualityMultiplier(void) const
Definition: ROL_HS41.hpp:148
Ptr< Constraint< Real > > getEqualityConstraint(void) const
Definition: ROL_HS55.hpp:175