ROL
dual-spaces/rosenbrock-1/example_02.cpp
Go to the documentation of this file.
1 // @HEADER
2 // *****************************************************************************
3 // Rapid Optimization Library (ROL) Package
4 //
5 // Copyright 2014 NTESS and the ROL contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 // @HEADER
9 
14 #define USE_HESSVEC 1
15 
16 #include "ROL_Rosenbrock.hpp"
17 #include "ROL_Solver.hpp"
18 #include "ROL_Bounds.hpp"
19 #include "ROL_Stream.hpp"
20 #include "Teuchos_GlobalMPISession.hpp"
21 #include <iostream>
22 
23 typedef double RealT;
24 
25 
26 /*** Declare two vector spaces. ***/
27 
28 // Forward declarations:
29 
30 template <class Real, class Element=Real>
31 class OptStdVector; // Optimization space.
32 
33 template <class Real, class Element=Real>
34 class OptDualStdVector; // Dual optimization space.
35 
36 
37 // Vector space definitions:
38 
39 // Optimization space.
40 template <class Real, class Element>
41 class OptStdVector : public ROL::Vector<Real> {
42 
43  typedef std::vector<Element> vector;
45 
46  typedef typename vector::size_type uint;
47 
48 private:
49 ROL::Ptr<std::vector<Element> > std_vec_;
50 mutable ROL::Ptr<OptDualStdVector<Real> > dual_vec_;
51 
52 public:
53 
54 OptStdVector(const ROL::Ptr<std::vector<Element> > & std_vec) : std_vec_(std_vec), dual_vec_(ROL::nullPtr) {}
55 
56 void plus( const ROL::Vector<Real> &x ) {
57 
58  ROL::Ptr<const vector> xvalptr = dynamic_cast<const OptStdVector&>(x).getVector();
59  uint dimension = std_vec_->size();
60  for (uint i=0; i<dimension; i++) {
61  (*std_vec_)[i] += (*xvalptr)[i];
62  }
63 }
64 
65 void scale( const Real alpha ) {
66  uint dimension = std_vec_->size();
67  for (uint i=0; i<dimension; i++) {
68  (*std_vec_)[i] *= alpha;
69  }
70 }
71 
72 Real dot( const ROL::Vector<Real> &x ) const {
73  Real val = 0;
74 
75  ROL::Ptr<const vector> xvalptr = dynamic_cast<const OptStdVector&>(x).getVector();
76  uint dimension = std_vec_->size();
77  for (uint i=0; i<dimension; i++) {
78  val += (*std_vec_)[i]*(*xvalptr)[i];
79  }
80  return val;
81 }
82 
83 Real norm() const {
84  Real val = 0;
85  val = std::sqrt( dot(*this) );
86  return val;
87 }
88 
89 ROL::Ptr<ROL::Vector<Real> > clone() const {
90  return ROL::makePtr<OptStdVector>( ROL::makePtr<std::vector<Element>>(std_vec_->size()) );
91 }
92 
93 ROL::Ptr<const std::vector<Element> > getVector() const {
94  return std_vec_;
95 }
96 
97 ROL::Ptr<std::vector<Element> > getVector() {
98  return std_vec_;
99 }
100 
101 ROL::Ptr<ROL::Vector<Real> > basis( const int i ) const {
102 
103  ROL::Ptr<vector> e_ptr = ROL::makePtr<vector>(std_vec_->size(),0.0);
104  (*e_ptr)[i] = 1.0;
105 
106  ROL::Ptr<V> e = ROL::makePtr<OptStdVector>( e_ptr );
107 
108  return e;
109 }
110 
111 int dimension() const {return static_cast<int>(std_vec_->size());}
112 
113 const ROL::Vector<Real> & dual() const {
114  dual_vec_ = ROL::makePtr<OptDualStdVector<Real>>( ROL::makePtr<std::vector<Element>>(*std_vec_) );
115  return *dual_vec_;
116 }
117 
118 Real apply( const ROL::Vector<Real> &x ) const {
119  Real val = 0;
120 
121  ROL::Ptr<const vector> xvalptr = dynamic_cast<const OptDualStdVector<Real,Element>&>(x).getVector();
122  uint dimension = std_vec_->size();
123  for (uint i=0; i<dimension; i++) {
124  val += (*std_vec_)[i]*(*xvalptr)[i];
125  }
126  return val;
127 }
128 
129 void applyUnary( const ROL::Elementwise::UnaryFunction<Real> &f ) {
130  for( auto& e : *std_vec_ ) e = f.apply(e);
131 }
132 
133 void applyBinary( const ROL::Elementwise::BinaryFunction<Real> &f, const ROL::Vector<Real> &x ) {
134  ROL::Ptr<const vector> xvalptr = dynamic_cast<const OptStdVector&>(x).getVector();
135  uint dim = std_vec_->size();
136  for (uint i=0; i<dim; i++) {
137  (*std_vec_)[i] = f.apply((*std_vec_)[i],(*xvalptr)[i]);
138  }
139 }
140 
141 Real reduce( const ROL::Elementwise::ReductionOp<Real> &r ) const {
142  Real result = r.initialValue();
143  uint dim = std_vec_->size();
144  for(uint i=0; i<dim; ++i) {
145  r.reduce((*std_vec_)[i],result);
146  }
147  return result;
148 }
149 
150 }; // class OptStdVector
151 
152 
153 // Dual optimization space.
154 template <class Real, class Element>
155 class OptDualStdVector : public ROL::Vector<Real> {
156 
157  typedef std::vector<Element> vector;
159 
160  typedef typename vector::size_type uint;
161 
162 private:
163 ROL::Ptr<std::vector<Element> > std_vec_;
164 mutable ROL::Ptr<OptStdVector<Real> > dual_vec_;
165 
166 public:
167 
168 OptDualStdVector(const ROL::Ptr<std::vector<Element> > & std_vec) : std_vec_(std_vec), dual_vec_(ROL::nullPtr) {}
169 
170 void plus( const ROL::Vector<Real> &x ) {
171 
172  ROL::Ptr<const vector> xvalptr = dynamic_cast<const OptDualStdVector&>(x).getVector();
173 
174  uint dimension = std_vec_->size();
175  for (uint i=0; i<dimension; i++) {
176  (*std_vec_)[i] += (*xvalptr)[i];
177  }
178 }
179 
180 void scale( const Real alpha ) {
181  uint dimension = std_vec_->size();
182  for (uint i=0; i<dimension; i++) {
183  (*std_vec_)[i] *= alpha;
184  }
185 }
186 
187 Real dot( const ROL::Vector<Real> &x ) const {
188  Real val = 0;
189 
190  ROL::Ptr<const vector> xvalptr = dynamic_cast<const OptDualStdVector&>(x).getVector();
191  uint dimension = std_vec_->size();
192  for (uint i=0; i<dimension; i++) {
193  val += (*std_vec_)[i]*(*xvalptr)[i];
194  }
195  return val;
196 }
197 
198 Real norm() const {
199  Real val = 0;
200  val = std::sqrt( dot(*this) );
201  return val;
202 }
203 
204 ROL::Ptr<ROL::Vector<Real> > clone() const {
205  return ROL::makePtr<OptDualStdVector>( ROL::makePtr<std::vector<Element>>(std_vec_->size()) );
206 }
207 
208 ROL::Ptr<const std::vector<Element> > getVector() const {
209  return std_vec_;
210 }
211 
212 ROL::Ptr<std::vector<Element> > getVector() {
213  return std_vec_;
214 }
215 
216 ROL::Ptr<ROL::Vector<Real> > basis( const int i ) const {
217 
218  ROL::Ptr<vector> e_ptr = ROL::makePtr<vector>(std_vec_->size(),0.0);
219  (*e_ptr)[i] = 1.0;
220 
221  ROL::Ptr<V> e = ROL::makePtr<OptDualStdVector>( e_ptr );
222  return e;
223 }
224 
225 int dimension() const {return std_vec_->size();}
226 
227 const ROL::Vector<Real> & dual() const {
228  dual_vec_ = ROL::makePtr<OptStdVector<Real>>( ROL::makePtr<std::vector<Element>>(*std_vec_) );
229  return *dual_vec_;
230 }
231 
232 Real apply( const ROL::Vector<Real> &x ) const {
233  Real val = 0;
234 
235  ROL::Ptr<const vector> xvalptr = dynamic_cast<const OptStdVector<Real,Element>&>(x).getVector();
236  uint dimension = std_vec_->size();
237  for (uint i=0; i<dimension; i++) {
238  val += (*std_vec_)[i]*(*xvalptr)[i];
239  }
240  return val;
241 }
242 
243 void applyUnary( const ROL::Elementwise::UnaryFunction<Real> &f ) {
244  for( auto& e : *std_vec_ ) e = f.apply(e);
245 }
246 
247 void applyBinary( const ROL::Elementwise::BinaryFunction<Real> &f, const ROL::Vector<Real> &x ) {
248  ROL::Ptr<const vector> xvalptr = dynamic_cast<const OptDualStdVector&>(x).getVector();
249  uint dim = std_vec_->size();
250  for (uint i=0; i<dim; i++) {
251  (*std_vec_)[i] = f.apply((*std_vec_)[i],(*xvalptr)[i]);
252  }
253 }
254 
255 Real reduce( const ROL::Elementwise::ReductionOp<Real> &r ) const {
256  Real result = r.initialValue();
257  uint dim = std_vec_->size();
258  for(uint i=0; i<dim; ++i) {
259  r.reduce((*std_vec_)[i],result);
260  }
261  return result;
262 }
263 
264 }; // class OptDualStdVector
265 
266 
267 /*** End of declaration of two vector spaces. ***/
268 
269 
270 
271 
272 
273 
274 int main(int argc, char *argv[]) {
275 
276  Teuchos::GlobalMPISession mpiSession(&argc, &argv);
277 
278  // This little trick lets us print to std::cout only if a (dummy) command-line argument is provided.
279  int iprint = argc - 1;
280  ROL::Ptr<std::ostream> outStream;
281  ROL::nullstream bhs; // outputs nothing
282  if (iprint > 0)
283  outStream = ROL::makePtrFromRef(std::cout);
284  else
285  outStream = ROL::makePtrFromRef(bhs);
286 
287  int errorFlag = 0;
288 
289  // *** Example body.
290 
291  try {
292 
293  // Define objective function.
294  ROL::Ptr<ROL::ZOO::Objective_Rosenbrock<RealT, OptStdVector<RealT>, OptDualStdVector<RealT>>> obj
295  = ROL::makePtr<ROL::ZOO::Objective_Rosenbrock<RealT, OptStdVector<RealT>, OptDualStdVector<RealT>>>();
296  int dim = 100; // Set problem dimension. Must be even.
297 
298  // Iteration Vector
299  ROL::Ptr<std::vector<RealT>> x_ptr = ROL::makePtr<std::vector<RealT>>(dim, 0.0);
300  ROL::Ptr<std::vector<RealT>> g_ptr = ROL::makePtr<std::vector<RealT>>(dim, 0.0);
301  ROL::Ptr<std::vector<RealT>> l_ptr = ROL::makePtr<std::vector<RealT>>(dim, 0.0);
302  // Set Initial Guess
303  for (int i=0; i<dim/2; i++) {
304  (*x_ptr)[2*i] = -1.2;
305  (*x_ptr)[2*i+1] = 1.0;
306  (*g_ptr)[2*i] = 0;
307  (*g_ptr)[2*i+1] = 0;
308  (*l_ptr)[2*i] = ROL::ROL_NINF<RealT>();
309  (*l_ptr)[2*i+1] = -1.5;
310  }
311  ROL::Ptr<OptStdVector<RealT>> x = ROL::makePtr<OptStdVector<RealT>>(x_ptr); // Iteration Vector
312  ROL::Ptr<OptDualStdVector<RealT>> g = ROL::makePtr<OptDualStdVector<RealT>>(g_ptr); // zeroed gradient vector in dual space
313  ROL::Ptr<OptStdVector<RealT>> l = ROL::makePtr<OptStdVector<RealT>>(l_ptr); // Lower bound Vector
314 
315  // Bound constraint
316  ROL::Ptr<ROL::Bounds<RealT>> bnd = ROL::makePtr<ROL::Bounds<RealT>>(*l);
317 
318  // Check vector.
319  ROL::Ptr<std::vector<RealT> > aa_ptr = ROL::makePtr<std::vector<RealT>>(1, 1.0);
320  OptDualStdVector<RealT> av(aa_ptr);
321  ROL::Ptr<std::vector<RealT> > bb_ptr = ROL::makePtr<std::vector<RealT>>(1, 2.0);
322  OptDualStdVector<RealT> bv(bb_ptr);
323  ROL::Ptr<std::vector<RealT> > cc_ptr = ROL::makePtr<std::vector<RealT>>(1, 3.0);
324  OptDualStdVector<RealT> cv(cc_ptr);
325  std::vector<RealT> std_vec_err = av.checkVector(bv,cv,true,*outStream);
326 
327  // Build optimization problem.
328  ROL::Ptr<ROL::Problem<RealT>> problem
329  = ROL::makePtr<ROL::Problem<RealT>>(obj,x,g);
330  problem->addBoundConstraint(bnd);
331  problem->finalize(false,true,*outStream);
332 
333  // Define algorithm.
334  ROL::ParameterList parlist;
335  std::string stepname = "Trust Region";
336  parlist.sublist("Step").set("Type",stepname);
337  //parlist.sublist("Step").sublist(stepname).sublist("Descent Method").set("Type","Newton-Krylov");
338  parlist.sublist("Step").sublist(stepname).set("Subproblem Solver", "Truncated CG");
339  parlist.sublist("Step").sublist(stepname).set("Subproblem Model", "Lin-More");
340  parlist.sublist("General").set("Output Level",1);
341  parlist.sublist("General").sublist("Krylov").set("Relative Tolerance",1e-2);
342  parlist.sublist("General").sublist("Krylov").set("Iteration Limit",10);
343  parlist.sublist("General").sublist("Krylov").set("Absolute Tolerance",1e-4);
344  parlist.sublist("General").sublist("Secant").set("Type","Limited-Memory BFGS");
345  parlist.sublist("General").sublist("Secant").set("Use as Hessian",false);
346  parlist.sublist("Status Test").set("Gradient Tolerance",1.e-12);
347  parlist.sublist("Status Test").set("Step Tolerance",1.e-14);
348  parlist.sublist("Status Test").set("Iteration Limit",100);
349  ROL::Ptr<ROL::Solver<RealT>> solver
350  = ROL::makePtr<ROL::Solver<RealT>>(problem,parlist);
351 
352  // Run Algorithm
353  solver->solve(*outStream);
354 
355  // Get True Solution
356  ROL::Ptr<std::vector<RealT> > xtrue_ptr = ROL::makePtr<std::vector<RealT>>(dim, 1.0);
357  OptStdVector<RealT> xtrue(xtrue_ptr);
358 
359  // Compute Errors
360  x->axpy(-1.0, xtrue);
361  RealT abserr = x->norm();
362  RealT relerr = abserr/xtrue.norm();
363  *outStream << std::scientific << "\n Absolute solution error: " << abserr;
364  *outStream << std::scientific << "\n Relative solution error: " << relerr;
365  if ( relerr > sqrt(ROL::ROL_EPSILON<RealT>()) ) {
366  errorFlag += 1;
367  }
368  ROL::Ptr<std::vector<RealT> > vec_err_ptr = ROL::makePtr<std::vector<RealT>>(std_vec_err);
369  ROL::StdVector<RealT> vec_err(vec_err_ptr);
370  *outStream << std::scientific << "\n Linear algebra error: " << vec_err.norm() << std::endl;
371  if ( vec_err.norm() > 1e2*ROL::ROL_EPSILON<RealT>() ) {
372  errorFlag += 1;
373  }
374  }
375  catch (std::logic_error& err) {
376  *outStream << err.what() << "\n";
377  errorFlag = -1000;
378  }; // end try
379 
380  if (errorFlag != 0)
381  std::cout << "End Result: TEST FAILED\n";
382  else
383  std::cout << "End Result: TEST PASSED\n";
384 
385  return 0;
386 
387 }
388 
typename PV< Real >::size_type size_type
void scale(const Real alpha)
Compute where .
const ROL::Vector< Real > & dual() const
Return dual representation of , for example, the result of applying a Riesz map, or change of basis...
ROL::Ptr< std::vector< Element > > getVector()
OptDualStdVector(const ROL::Ptr< std::vector< Element > > &std_vec)
virtual std::vector< Real > checkVector(const Vector< Real > &x, const Vector< Real > &y, const bool printToStream=true, std::ostream &outStream=std::cout) const
Verify vector-space methods.
Definition: ROL_Vector.hpp:291
Contains definitions for Rosenbrock&#39;s function.
ROL::Ptr< ROL::Vector< Real > > basis(const int i) const
Return i-th basis vector.
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:46
Defines a no-output stream class ROL::NullStream and a function makeStreamPtr which either wraps a re...
void applyUnary(const ROL::Elementwise::UnaryFunction< Real > &f)
Real apply(const ROL::Vector< Real > &x) const
Apply to a dual vector. This is equivalent to the call .
ROL::Ptr< ROL::Vector< Real > > clone() const
Clone to make a new (uninitialized) vector.
ROL::Ptr< OptStdVector< Real > > dual_vec_
ROL::Ptr< std::vector< Element > > std_vec_
ROL::Ptr< std::vector< Element > > getVector()
const ROL::Vector< Real > & dual() const
Return dual representation of , for example, the result of applying a Riesz map, or change of basis...
ROL::Ptr< ROL::Vector< Real > > clone() const
Clone to make a new (uninitialized) vector.
void applyUnary(const ROL::Elementwise::UnaryFunction< Real > &f)
Real reduce(const ROL::Elementwise::ReductionOp< Real > &r) const
Provides the ROL::Vector interface for scalar values, to be used, for example, with scalar constraint...
void scale(const Real alpha)
Compute where .
ROL::Ptr< const std::vector< Element > > getVector() const
Real dot(const ROL::Vector< Real > &x) const
Compute where .
int dimension() const
Return dimension of the vector space.
Real reduce(const ROL::Elementwise::ReductionOp< Real > &r) const
ROL::Ptr< OptDualStdVector< Real > > dual_vec_
Real apply(const ROL::Vector< Real > &x) const
Apply to a dual vector. This is equivalent to the call .
void applyBinary(const ROL::Elementwise::BinaryFunction< Real > &f, const ROL::Vector< Real > &x)
Real norm() const
Returns where .
int dimension() const
Return dimension of the vector space.
ROL::Ptr< std::vector< Element > > std_vec_
basic_nullstream< char, char_traits< char >> nullstream
Definition: ROL_Stream.hpp:38
int main(int argc, char *argv[])
void applyBinary(const ROL::Elementwise::BinaryFunction< Real > &f, const ROL::Vector< Real > &x)
ROL::Ptr< const std::vector< Element > > getVector() const
Real dot(const ROL::Vector< Real > &x) const
Compute where .
constexpr auto dim
void plus(const ROL::Vector< Real > &x)
Compute , where .
void plus(const ROL::Vector< Real > &x)
Compute , where .
Real norm() const
Returns where .
OptStdVector(const ROL::Ptr< std::vector< Element > > &std_vec)
ROL::Ptr< ROL::Vector< Real > > basis(const int i) const
Return i-th basis vector.