ROL
function/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
48#include "ROL_RandomVector.hpp"
49#include "ROL_StdVector.hpp"
50#include "ROL_Bounds.hpp"
52
53#include "ROL_Stream.hpp"
54#include "Teuchos_GlobalMPISession.hpp"
55#include "ROL_ParameterList.hpp"
56
57
58typedef double RealT;
59
60int main(int argc, char *argv[]) {
61
62 typedef std::vector<RealT> vector;
63 typedef ROL::Vector<RealT> V;
64 typedef ROL::StdVector<RealT> SV;
65
66 typedef typename vector::size_type uint;
67
68 Teuchos::GlobalMPISession mpiSession(&argc, &argv);
69
70 // This little trick lets us print to std::cout only if a (dummy) command-line argument is provided.
71 int iprint = argc - 1;
72 ROL::Ptr<std::ostream> outStream;
73 ROL::nullstream bhs; // outputs nothing
74 if (iprint > 0)
75 outStream = ROL::makePtrFromRef(std::cout);
76 else
77 outStream = ROL::makePtrFromRef(bhs);
78
79 // Save the format state of the original std::cout.
80 ROL::nullstream oldFormatState;
81 oldFormatState.copyfmt(std::cout);
82
83// RealT errtol = std::sqrt(ROL::ROL_THRESHOLD<RealT>());
84
85 int errorFlag = 0;
86
87 // *** Test body.
88
89 try {
90
91 uint dim = 30;
92 RealT xmin = 0.5;
93 RealT xmax = 2.5;
94
95 ROL::Ptr<vector> x_ptr = ROL::makePtr<vector>(dim,0.0);
96 ROL::Ptr<vector> g_ptr = ROL::makePtr<vector>(dim,0.0);
97 ROL::Ptr<vector> v_ptr = ROL::makePtr<vector>(dim,0.0);
98 ROL::Ptr<vector> hv_ptr = ROL::makePtr<vector>(dim,0.0);
99
100 ROL::Ptr<vector> l_ptr = ROL::makePtr<vector>(dim,1.0);
101 ROL::Ptr<vector> u_ptr = ROL::makePtr<vector>(dim,2.0);
102
103 ROL::Ptr<vector> xlog0_ptr = ROL::makePtr<vector>(dim,0.0);
104 ROL::Ptr<vector> xlog1_ptr = ROL::makePtr<vector>(dim,0.0);
105 ROL::Ptr<vector> xlog2_ptr = ROL::makePtr<vector>(dim,0.0);
106
107 ROL::Ptr<vector> xquad0_ptr = ROL::makePtr<vector>(dim,0.0);
108 ROL::Ptr<vector> xquad1_ptr = ROL::makePtr<vector>(dim,0.0);
109 ROL::Ptr<vector> xquad2_ptr = ROL::makePtr<vector>(dim,0.0);
110
111 ROL::Ptr<vector> xdwell0_ptr = ROL::makePtr<vector>(dim,0.0);
112 ROL::Ptr<vector> xdwell1_ptr = ROL::makePtr<vector>(dim,0.0);
113 ROL::Ptr<vector> xdwell2_ptr = ROL::makePtr<vector>(dim,0.0);
114
115
116
117 SV x(x_ptr);
118 SV g(g_ptr);
119 SV v(v_ptr);
120 SV hv(hv_ptr);
121
122 ROL::Ptr<SV> xlog0 = ROL::makePtr<SV>(xlog0_ptr);
123 ROL::Ptr<SV> xlog1 = ROL::makePtr<SV>(xlog1_ptr);
124 ROL::Ptr<SV> xlog2 = ROL::makePtr<SV>(xlog2_ptr);
125
126 ROL::Ptr<SV> xquad0 = ROL::makePtr<SV>(xquad0_ptr);
127 ROL::Ptr<SV> xquad1 = ROL::makePtr<SV>(xquad1_ptr);
128 ROL::Ptr<SV> xquad2 = ROL::makePtr<SV>(xquad2_ptr);
129
130 ROL::Ptr<SV> xdwell0 = ROL::makePtr<SV>(xdwell0_ptr);
131 ROL::Ptr<SV> xdwell1 = ROL::makePtr<SV>(xdwell1_ptr);
132 ROL::Ptr<SV> xdwell2 = ROL::makePtr<SV>(xdwell2_ptr);
133
134 ROL::Ptr<V> lo = ROL::makePtr<SV>(l_ptr);
135 ROL::Ptr<V> up = ROL::makePtr<SV>(u_ptr);
136
137 for(uint i=0; i<dim; ++i) {
138 RealT t = static_cast<RealT>(i)/static_cast<RealT>(dim-1);
139 (*x_ptr)[i] = xmin*(1-t) + xmax*t;
140 }
141
142 // Create bound constraint
143 ROL::Bounds<RealT> bc(lo,up);
144
145 ROL::ParameterList logList;
146 ROL::ParameterList quadList;
147 ROL::ParameterList dwellList;
148
149 logList.sublist("Barrier Function").set("Type","Logarithmic");
150 quadList.sublist("Barrier Function").set("Type","Quadratic");
151 dwellList.sublist("Barrier Function").set("Type","Double Well");
152
154 ROL::ObjectiveFromBoundConstraint<RealT> quadObj(bc,quadList);
155 ROL::ObjectiveFromBoundConstraint<RealT> dwellObj(bc,dwellList);
156
157 RealT tol = 0.0;
158
159
160 logObj.value(x,tol);
161 xlog0->set(*ROL::staticPtrCast<SV>(logObj.getBarrierVector()));
162
163 logObj.gradient(g,x,tol);
164 xlog1->set(*ROL::staticPtrCast<SV>(logObj.getBarrierVector()));
165
166 logObj.hessVec(hv,v,x,tol);
167 xlog2->set(*ROL::staticPtrCast<SV>(logObj.getBarrierVector()));
168
169
170 quadObj.value(x,tol);
171 xquad0->set(*ROL::staticPtrCast<SV>(quadObj.getBarrierVector()));
172
173 quadObj.gradient(g,x,tol);
174 xquad1->set(*ROL::staticPtrCast<SV>(quadObj.getBarrierVector()));
175
176 quadObj.hessVec(hv,v,x,tol);
177 xquad2->set(*ROL::staticPtrCast<SV>(quadObj.getBarrierVector()));
178
179
180 dwellObj.value(x,tol);
181 xdwell0->set(*ROL::staticPtrCast<SV>(dwellObj.getBarrierVector()));
182
183 dwellObj.gradient(g,x,tol);
184 xdwell1->set(*ROL::staticPtrCast<SV>(dwellObj.getBarrierVector()));
185
186 dwellObj.hessVec(hv,v,x,tol);
187 xdwell2->set(*ROL::staticPtrCast<SV>(dwellObj.getBarrierVector()));
188
189
190 *outStream << std::setw(14) << "x"
191 << std::setw(14) << "log"
192 << std::setw(14) << "D(log)"
193 << std::setw(14) << "D2(log)"
194 << std::setw(14) << "quad"
195 << std::setw(14) << "D(quad)"
196 << std::setw(14) << "D2(quad)"
197 << std::setw(14) << "dwell"
198 << std::setw(14) << "D(dwell)"
199 << std::setw(14) << "D2(dwell)"
200 << std::endl;
201 *outStream << std::string(140,'-') << std::endl;
202
203 for(uint i=0; i<dim; ++i) {
204 *outStream << std::setw(14) << (*x_ptr)[i]
205 << std::setw(14) << (*xlog0_ptr)[i]
206 << std::setw(14) << (*xlog1_ptr)[i]
207 << std::setw(14) << (*xlog2_ptr)[i]
208 << std::setw(14) << (*xquad0_ptr)[i]
209 << std::setw(14) << (*xquad1_ptr)[i]
210 << std::setw(14) << (*xquad2_ptr)[i]
211 << std::setw(14) << (*xdwell0_ptr)[i]
212 << std::setw(14) << (*xdwell1_ptr)[i]
213 << std::setw(14) << (*xdwell2_ptr)[i]
214 << std::endl;
215 }
216
217
218 ROL::RandomizeVector( x, 1.2, 1.8 );
219 ROL::RandomizeVector( v, -0.1, 0.1 );
220
221 *outStream << "\n\n";
222 *outStream << "Test of logarithmic penalty objective" << std::endl;
223 logObj.checkGradient(x,v,true,*outStream); *outStream << std::endl;
224 logObj.checkHessVec(x,v,true,*outStream); *outStream << std::endl;
225
226 ROL::RandomizeVector( x, -1.0, 1.0 );
227 ROL::RandomizeVector( v, -1.0, 1.0 );
228
229 *outStream << "\n\n";
230 *outStream << "Test of piecewise quadratic penalty objective" << std::endl;
231 quadObj.checkGradient(x,v,true,*outStream); *outStream << std::endl;
232 quadObj.checkHessVec(x,v,true,*outStream); *outStream << std::endl;
233
234
235 *outStream << "\n\n";
236 *outStream << "Test of double well penalty objective" << std::endl;
237 dwellObj.checkGradient(x,v,true,*outStream); *outStream << std::endl;
238 dwellObj.checkHessVec(x,v,true,*outStream); *outStream << std::endl;
239
240
241
242
243
244 }
245 catch (std::logic_error& err) {
246 *outStream << err.what() << "\n";
247 errorFlag = -1000;
248 }; // end try
249
250 if (errorFlag != 0)
251 std::cout << "End Result: TEST FAILED\n";
252 else
253 std::cout << "End Result: TEST PASSED\n";
254
255 return 0;
256
257
258}
259
Vector< Real > V
Defines a no-output stream class ROL::NullStream and a function makeStreamPtr which either wraps a re...
Provides the elementwise interface to apply upper and lower bound constraints.
Create a penalty objective from upper and lower bound vectors.
Real value(const Vector< Real > &x, Real &tol)
Compute value.
void gradient(Vector< Real > &g, const Vector< Real > &x, Real &tol)
Compute gradient.
void hessVec(Vector< Real > &hv, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Apply Hessian approximation to vector.
virtual std::vector< std::vector< Real > > checkGradient(const Vector< Real > &x, const Vector< Real > &d, const bool printToStream=true, std::ostream &outStream=std::cout, const int numSteps=ROL_NUM_CHECKDERIV_STEPS, const int order=1)
Finite-difference gradient check.
virtual std::vector< std::vector< Real > > checkHessVec(const Vector< Real > &x, const Vector< Real > &v, const bool printToStream=true, std::ostream &outStream=std::cout, const int numSteps=ROL_NUM_CHECKDERIV_STEPS, const int order=1)
Finite-difference Hessian-applied-to-vector check.
Provides the ROL::Vector interface for scalar values, to be used, for example, with scalar constraint...
Defines the linear algebra or vector space interface.
int main(int argc, char *argv[])
double RealT
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,...
constexpr auto dim