ROL
dual-spaces/rosenbrock-1/example_01.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_Stream.hpp"
19 #include "Teuchos_GlobalMPISession.hpp"
20 #include <iostream>
21 
22 typedef double RealT;
23 
24 
25 /*** Declare two vector spaces. ***/
26 
27 // Forward declarations:
28 
29 template <class Real, class Element=Real>
30 class OptStdVector; // Optimization space.
31 
32 template <class Real, class Element=Real>
33 class OptDualStdVector; // Dual optimization space.
34 
35 
36 // Vector space definitions:
37 
38 // Optimization space.
39 template <class Real, class Element>
40 class OptStdVector : public ROL::Vector<Real> {
41 
42  typedef std::vector<Element> vector;
44 
45  typedef typename vector::size_type uint;
46 
47 private:
48 ROL::Ptr<std::vector<Element> > std_vec_;
49 mutable ROL::Ptr<OptDualStdVector<Real> > dual_vec_;
50 
51 public:
52 
53 OptStdVector(const ROL::Ptr<std::vector<Element> > & std_vec) : std_vec_(std_vec), dual_vec_(ROL::nullPtr) {}
54 
55 void plus( const ROL::Vector<Real> &x ) {
56 
57  ROL::Ptr<const vector> xvalptr = dynamic_cast<const OptStdVector&>(x).getVector();
58  uint dimension = std_vec_->size();
59  for (uint i=0; i<dimension; i++) {
60  (*std_vec_)[i] += (*xvalptr)[i];
61  }
62 }
63 
64 void scale( const Real alpha ) {
65  uint dimension = std_vec_->size();
66  for (uint i=0; i<dimension; i++) {
67  (*std_vec_)[i] *= alpha;
68  }
69 }
70 
71 Real dot( const ROL::Vector<Real> &x ) const {
72  Real val = 0;
73 
74  ROL::Ptr<const vector> xvalptr = dynamic_cast<const OptStdVector&>(x).getVector();
75  uint dimension = std_vec_->size();
76  for (uint i=0; i<dimension; i++) {
77  val += (*std_vec_)[i]*(*xvalptr)[i];
78  }
79  return val;
80 }
81 
82 Real norm() const {
83  Real val = 0;
84  val = std::sqrt( dot(*this) );
85  return val;
86 }
87 
88 ROL::Ptr<ROL::Vector<Real> > clone() const {
89  return ROL::makePtr<OptStdVector>( ROL::makePtr<std::vector<Element>>(std_vec_->size()) );
90 }
91 
92 ROL::Ptr<const std::vector<Element> > getVector() const {
93  return std_vec_;
94 }
95 
96 ROL::Ptr<std::vector<Element> > getVector() {
97  return std_vec_;
98 }
99 
100 ROL::Ptr<ROL::Vector<Real> > basis( const int i ) const {
101 
102  ROL::Ptr<vector> e_ptr = ROL::makePtr<vector>(std_vec_->size(),0.0);
103  (*e_ptr)[i] = 1.0;
104 
105  ROL::Ptr<V> e = ROL::makePtr<OptStdVector>( e_ptr );
106 
107  return e;
108 }
109 
110 int dimension() const {return static_cast<int>(std_vec_->size());}
111 
112 const ROL::Vector<Real> & dual() const {
113  dual_vec_ = ROL::makePtr<OptDualStdVector<Real>>( ROL::makePtr<std::vector<Element>>(*std_vec_) );
114  return *dual_vec_;
115 }
116 
117 Real apply( const ROL::Vector<Real> &x ) const {
118  Real val = 0;
119 
120  ROL::Ptr<const vector> xvalptr = dynamic_cast<const OptDualStdVector<Real,Element>&>(x).getVector();
121  uint dimension = std_vec_->size();
122  for (uint i=0; i<dimension; i++) {
123  val += (*std_vec_)[i]*(*xvalptr)[i];
124  }
125  return val;
126 }
127 
128 }; // class OptStdVector
129 
130 
131 // Dual optimization space.
132 template <class Real, class Element>
133 class OptDualStdVector : public ROL::Vector<Real> {
134 
135  typedef std::vector<Element> vector;
137 
138  typedef typename vector::size_type uint;
139 
140 private:
141 ROL::Ptr<std::vector<Element> > std_vec_;
142 mutable ROL::Ptr<OptStdVector<Real> > dual_vec_;
143 
144 public:
145 
146 OptDualStdVector(const ROL::Ptr<std::vector<Element> > & std_vec) : std_vec_(std_vec), dual_vec_(ROL::nullPtr) {}
147 
148 void plus( const ROL::Vector<Real> &x ) {
149 
150  ROL::Ptr<const vector> xvalptr = dynamic_cast<const OptDualStdVector&>(x).getVector();
151 
152  uint dimension = std_vec_->size();
153  for (uint i=0; i<dimension; i++) {
154  (*std_vec_)[i] += (*xvalptr)[i];
155  }
156 }
157 
158 void scale( const Real alpha ) {
159  uint dimension = std_vec_->size();
160  for (uint i=0; i<dimension; i++) {
161  (*std_vec_)[i] *= alpha;
162  }
163 }
164 
165 Real dot( const ROL::Vector<Real> &x ) const {
166  Real val = 0;
167 
168  ROL::Ptr<const vector> xvalptr = dynamic_cast<const OptDualStdVector&>(x).getVector();
169  uint dimension = std_vec_->size();
170  for (uint i=0; i<dimension; i++) {
171  val += (*std_vec_)[i]*(*xvalptr)[i];
172  }
173  return val;
174 }
175 
176 Real norm() const {
177  Real val = 0;
178  val = std::sqrt( dot(*this) );
179  return val;
180 }
181 
182 ROL::Ptr<ROL::Vector<Real> > clone() const {
183  return ROL::makePtr<OptDualStdVector>( ROL::makePtr<std::vector<Element>>(std_vec_->size()) );
184 }
185 
186 ROL::Ptr<const std::vector<Element> > getVector() const {
187  return std_vec_;
188 }
189 
190 ROL::Ptr<std::vector<Element> > getVector() {
191  return std_vec_;
192 }
193 
194 ROL::Ptr<ROL::Vector<Real> > basis( const int i ) const {
195 
196  ROL::Ptr<vector> e_ptr = ROL::makePtr<vector>(std_vec_->size(),0.0);
197  (*e_ptr)[i] = 1.0;
198 
199  ROL::Ptr<V> e = ROL::makePtr<OptDualStdVector>( e_ptr );
200  return e;
201 }
202 
203 int dimension() const {return std_vec_->size();}
204 
205 const ROL::Vector<Real> & dual() const {
206  dual_vec_ = ROL::makePtr<OptStdVector<Real>>( ROL::makePtr<std::vector<Element>>(*std_vec_) );
207  return *dual_vec_;
208 }
209 
210 Real apply( const ROL::Vector<Real> &x ) const {
211  Real val = 0;
212 
213  ROL::Ptr<const vector> xvalptr = dynamic_cast<const OptStdVector<Real,Element>&>(x).getVector();
214  uint dimension = std_vec_->size();
215  for (uint i=0; i<dimension; i++) {
216  val += (*std_vec_)[i]*(*xvalptr)[i];
217  }
218  return val;
219 }
220 
221 }; // class OptDualStdVector
222 
223 
224 /*** End of declaration of two vector spaces. ***/
225 
226 
227 
228 
229 
230 
231 int main(int argc, char *argv[]) {
232 
233  Teuchos::GlobalMPISession mpiSession(&argc, &argv);
234 
235  // This little trick lets us print to std::cout only if a (dummy) command-line argument is provided.
236  int iprint = argc - 1;
237  ROL::Ptr<std::ostream> outStream;
238  ROL::nullstream bhs; // outputs nothing
239  if (iprint > 0)
240  outStream = ROL::makePtrFromRef(std::cout);
241  else
242  outStream = ROL::makePtrFromRef(bhs);
243 
244  int errorFlag = 0;
245 
246  // *** Example body.
247 
248  try {
249 
250  // Define objective function.
251  ROL::Ptr<ROL::ZOO::Objective_Rosenbrock<RealT, OptStdVector<RealT>, OptDualStdVector<RealT>>> obj
252  = ROL::makePtr<ROL::ZOO::Objective_Rosenbrock<RealT, OptStdVector<RealT>, OptDualStdVector<RealT>>>();
253  int dim = 100; // Set problem dimension. Must be even.
254 
255  // Iteration Vector
256  ROL::Ptr<std::vector<RealT>> x_ptr = ROL::makePtr<std::vector<RealT>>(dim, 0.0);
257  ROL::Ptr<std::vector<RealT>> g_ptr = ROL::makePtr<std::vector<RealT>>(dim, 0.0);
258  // Set Initial Guess
259  for (int i=0; i<dim/2; i++) {
260  (*x_ptr)[2*i] = -1.2;
261  (*x_ptr)[2*i+1] = 1.0;
262  (*g_ptr)[2*i] = 0;
263  (*g_ptr)[2*i+1] = 0;
264  }
265  ROL::Ptr<OptStdVector<RealT>> x = ROL::makePtr<OptStdVector<RealT>>(x_ptr); // Iteration Vector
266  ROL::Ptr<OptDualStdVector<RealT>> g = ROL::makePtr<OptDualStdVector<RealT>>(g_ptr); // zeroed gradient vector in dual space
267 
268  // Check vector.
269  ROL::Ptr<std::vector<RealT> > aa_ptr = ROL::makePtr<std::vector<RealT>>(1, 1.0);
270  OptDualStdVector<RealT> av(aa_ptr);
271  ROL::Ptr<std::vector<RealT> > bb_ptr = ROL::makePtr<std::vector<RealT>>(1, 2.0);
272  OptDualStdVector<RealT> bv(bb_ptr);
273  ROL::Ptr<std::vector<RealT> > cc_ptr = ROL::makePtr<std::vector<RealT>>(1, 3.0);
274  OptDualStdVector<RealT> cv(cc_ptr);
275  std::vector<RealT> std_vec_err = av.checkVector(bv,cv,true,*outStream);
276 
277  // Build optimization problem.
278  ROL::Ptr<ROL::Problem<RealT>> problem
279  = ROL::makePtr<ROL::Problem<RealT>>(obj,x,g);
280  problem->finalize(false,true,*outStream);
281 
282  // Define algorithm.
283  ROL::ParameterList parlist;
284  std::string stepname = "Trust Region";
285  parlist.sublist("Step").set("Type",stepname);
286  //parlist.sublist("Step").sublist(stepname).sublist("Descent Method").set("Type","Quasi-Newton Method");
287  parlist.sublist("Step").sublist(stepname).set("Subproblem Solver", "Truncated CG");
288  parlist.sublist("General").set("Output Level",1);
289  parlist.sublist("General").sublist("Krylov").set("Relative Tolerance",1e-2);
290  parlist.sublist("General").sublist("Krylov").set("Iteration Limit",10);
291  parlist.sublist("General").sublist("Krylov").set("Absolute Tolerance",1e-4);
292  parlist.sublist("General").sublist("Secant").set("Type","Limited-Memory BFGS");
293  parlist.sublist("General").sublist("Secant").set("Use as Hessian",true);
294  parlist.sublist("Status Test").set("Gradient Tolerance",1.e-12);
295  parlist.sublist("Status Test").set("Step Tolerance",1.e-14);
296  parlist.sublist("Status Test").set("Iteration Limit",100);
297  ROL::Ptr<ROL::Solver<RealT>> solver
298  = ROL::makePtr<ROL::Solver<RealT>>(problem,parlist);
299 
300  // Run Algorithm
301  solver->solve(*outStream);
302 
303  // Get True Solution
304  ROL::Ptr<std::vector<RealT> > xtrue_ptr = ROL::makePtr<std::vector<RealT>>(dim, 1.0);
305  OptStdVector<RealT> xtrue(xtrue_ptr);
306 
307  // Compute Errors
308  x->axpy(-1.0, xtrue);
309  RealT abserr = x->norm();
310  RealT relerr = abserr/xtrue.norm();
311  *outStream << std::scientific << "\n Absolute solution error: " << abserr;
312  *outStream << std::scientific << "\n Relative solution error: " << relerr;
313  if ( relerr > sqrt(ROL::ROL_EPSILON<RealT>()) ) {
314  errorFlag += 1;
315  }
316  ROL::Ptr<std::vector<RealT> > vec_err_ptr = ROL::makePtr<std::vector<RealT>>(std_vec_err);
317  ROL::StdVector<RealT> vec_err(vec_err_ptr);
318  *outStream << std::scientific << "\n Linear algebra error: " << vec_err.norm() << std::endl;
319  if ( vec_err.norm() > 1e2*ROL::ROL_EPSILON<RealT>() ) {
320  errorFlag += 1;
321  }
322  }
323  catch (std::logic_error& err) {
324  *outStream << err.what() << "\n";
325  errorFlag = -1000;
326  }; // end try
327 
328  if (errorFlag != 0)
329  std::cout << "End Result: TEST FAILED\n";
330  else
331  std::cout << "End Result: TEST PASSED\n";
332 
333  return 0;
334 
335 }
336 
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...
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.
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.
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 .
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[])
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.