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