ROL
step/interiorpoint/test_01.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
44#include "Teuchos_GlobalMPISession.hpp"
45
46#include "Teuchos_GlobalMPISession.hpp"
47
49#include "ROL_RandomVector.hpp"
50#include "ROL_StdVector.hpp"
51#include "ROL_Bounds.hpp"
52#include "ROL_ParameterList.hpp"
53
54#include <iomanip>
55
56
62template<class Real>
63class NullObjective : public ROL::Objective<Real> {
65public:
66 Real value( const V &x, Real &tol ) {
67 return Real(0.0);
68 }
69 void gradient( V &g, const V &x, Real &tol ) {
70 g.zero();
71 }
72 void hessVec( V &hv, const V &v, const V &x, Real &tol ) {
73 hv.zero();
74 }
75};
76
77template<class Real>
78void printVector( const ROL::Vector<Real> &x, std::ostream &outStream ) {
79 ROL::Ptr<const std::vector<Real> > xp =
80 dynamic_cast<const ROL::StdVector<Real>&>(x).getVector();
81
82 for( size_t i=0; i<xp->size(); ++i ) {
83 outStream << (*xp)[i] << std::endl;
84 }
85}
86
87
88
89
90typedef double RealT;
91
92int main(int argc, char *argv[]) {
93
94 typedef std::vector<RealT> vector;
95
96 typedef ROL::Vector<RealT> V;
97 typedef ROL::StdVector<RealT> SV;
98 typedef ROL::Objective<RealT> OBJ;
100
101 typedef ROL::ParameterList PL;
102
103 Teuchos::GlobalMPISession mpiSession(&argc, &argv);
104
105 int iprint = argc - 1;
106 ROL::Ptr<std::ostream> outStream;
107 ROL::nullstream bhs;
108 if( iprint > 0 )
109 outStream = ROL::makePtrFromRef(std::cout);
110 else
111 outStream = ROL::makePtrFromRef(bhs);
112
113 int errorFlag = 0;
114
115 try {
116
117 PL parlist;
118 PL &iplist = parlist.sublist("Step").sublist("Primal Dual Interior Point");
119 PL &lblist = iplist.sublist("Barrier Objective");
120
121 RealT mu = 1.e2;
122 RealT kappaD = 1.e-4;
123 bool useLinearDamping = true;
124
125 lblist.set("Use Linear Damping", useLinearDamping);
126 lblist.set("Linear Damping Coefficient",kappaD);
127 lblist.set("Initial Barrier Parameter",mu);
128
129 RealT ninf = ROL::ROL_NINF<RealT>();
130 RealT inf = ROL::ROL_INF<RealT>();
131
132 int dim = 4;
133 int numTestVectors = 19;
134
135 ROL::Ptr<vector> x_ptr = ROL::makePtr<vector>(dim, 0.0);
136 ROL::Ptr<vector> d_ptr = ROL::makePtr<vector>(dim, 0.0);
137 ROL::Ptr<vector> v_ptr = ROL::makePtr<vector>(dim, 0.0);
138 ROL::Ptr<vector> l_ptr = ROL::makePtr<vector>(dim, 0.0);
139 ROL::Ptr<vector> u_ptr = ROL::makePtr<vector>(dim, 0.0);
140 ROL::Ptr<vector> e0_ptr = ROL::makePtr<vector>(dim, 0.0); // First canonical vector
141
142 (*e0_ptr)[0] = 1.0;
143
144 SV e0(e0_ptr);
145
146 // Lower Bounds // Upper Bounds
147 (*l_ptr)[0] = ninf; (*u_ptr)[0] = 5.0;
148 (*l_ptr)[1] = ninf; (*u_ptr)[1] = inf;
149 (*l_ptr)[2] = -5.0; (*u_ptr)[2] = inf;
150 (*l_ptr)[3] = -5.0; (*u_ptr)[3] = 5.0;
151
152 RealT left = -1.0; RealT right = 1.0;
153
154 RealT xmax = 4.99;
155
156 ROL::Ptr<V> x = ROL::makePtr<SV>( x_ptr );
157 ROL::Ptr<V> d = ROL::makePtr<SV>( d_ptr );
158 ROL::Ptr<V> v = ROL::makePtr<SV>( v_ptr );
159 ROL::Ptr<V> l = ROL::makePtr<SV>( l_ptr );
160 ROL::Ptr<V> u = ROL::makePtr<SV>( u_ptr );
161
162 ROL::Ptr<const V> maskL, maskU;
163
164 ROL::RandomizeVector(*d,left,right);
165 ROL::RandomizeVector(*v,left,right);
166
167 std::vector<RealT> values(numTestVectors); // Computed objective value for each
168 std::vector<RealT> exact_values(numTestVectors);
169
170 std::vector<ROL::Ptr<V> > x_test;
171
172 for(int i=0; i<numTestVectors; ++i) {
173 x_test.push_back(x->clone());
174 RealT t = static_cast<RealT>(i)/static_cast<RealT>(numTestVectors-1);
175 RealT fillValue = xmax*(2.0*t-1.0);
176 x_test[i]->applyUnary(ROL::Elementwise::Fill<RealT>(fillValue));
177 }
178
179 ROL::Ptr<OBJ> obj = ROL::makePtr<NullObjective<RealT>>();
180 ROL::Ptr<BND> bnd = ROL::makePtr<ROL::Bounds<RealT>>(l,u);
181
182 ROL::InteriorPointPenalty<RealT> ipobj(obj,bnd,parlist);
183
184 maskL = ipobj.getLowerMask();
185 maskU = ipobj.getUpperMask();
186
187 ROL::Ptr<const std::vector<RealT> > maskL_ptr = dynamic_cast<const SV&>(*maskL).getVector();
188 ROL::Ptr<const std::vector<RealT> > maskU_ptr = dynamic_cast<const SV&>(*maskU).getVector();
189
190 *outStream << "\nLower bound vector" << std::endl;
191 printVector(*l,*outStream);
192
193 *outStream << "\nUpper bound vector" << std::endl;
194 printVector(*u,*outStream);
195
196 *outStream << "\nLower mask vector" << std::endl;
197 printVector(*maskL, *outStream);
198
199 *outStream << "\nUpper mask vector" << std::endl;
200 printVector(*maskU, *outStream);
201
202 *outStream << "\nChecking Objective value" << std::endl;
203
204 RealT tol = std::sqrt(ROL::ROL_EPSILON<RealT>());
205 *outStream << std::setw(16) << "x[i], i=0,1,2,3"
206 << std::setw(20) << "Computed Objective"
207 << std::setw(20) << "Exact Objective" << std::endl;
208
209 RealT valueError(0.0);
210
211 for(int i=0; i<numTestVectors; ++i) {
212 values[i] = ipobj.value(*(x_test[i]),tol);
213
214 exact_values[i] = 0;
215
216 // Extract the value from the test vector that is in every element
217 RealT xval = x_test[i]->dot(e0);
218
219
220 for(int j=0; j<dim; ++j) {
221 if( (*maskL_ptr)[j] ) {
222 RealT diff = xval-(*l_ptr)[j];
223 exact_values[i] -= mu*std::log(diff);
224
225 if( useLinearDamping && !(*maskU_ptr)[j] ) {
226 exact_values[i] += mu*kappaD*diff;
227 }
228
229 }
230 if( (*maskU_ptr)[j] ) {
231 RealT diff = (*u_ptr)[j]-xval;
232 exact_values[i] -= mu*std::log(diff);
233
234 if(useLinearDamping && !(*maskL_ptr)[j] ) {
235 exact_values[i] += mu*kappaD*diff;
236 }
237
238 }
239 } // end loop over elements
240
241 *outStream << std::setw(16) << xval
242 << std::setw(20) << values[i]
243 << std::setw(20) << exact_values[i] << std::endl;
244 RealT valDiff = exact_values[i] - values[i];
245 valueError += valDiff*valDiff;
246 } // end loop over vectors
247
248 if(valueError>ROL::ROL_EPSILON<RealT>()) {
249 errorFlag++;
250 }
251
252 *outStream << "\nPerforming finite difference checks" << std::endl;
253
254 ipobj.checkGradient(*x,*v,true,*outStream); *outStream << std::endl;
255 ipobj.checkHessVec(*x,*d,true,*outStream); *outStream << std::endl;
256 ipobj.checkHessSym(*x,*d,*v,true,*outStream); *outStream << std::endl;
257
258 }
259 catch (std::logic_error& err) {
260 *outStream << err.what() << std::endl;
261 errorFlag = -1000;
262 }
263
264 if (errorFlag != 0)
265 std::cout << "End Result: TEST FAILED" << std::endl;
266 else
267 std::cout << "End Result: TEST PASSED" << std::endl;
268
269 return 0;
270}
Vector< Real > V
void gradient(V &g, const V &x, Real &tol)
Compute gradient.
void hessVec(V &hv, const V &v, const V &x, Real &tol)
Apply Hessian approximation to vector.
Real value(const V &x, Real &tol)
Compute value.
Provides the interface to apply upper and lower bound constraints.
Provides the interface to evaluate the Interior Pointy log barrier penalty function with upper and lo...
Provides the interface to evaluate objective functions.
Provides the ROL::Vector interface for scalar values, to be used, for example, with scalar constraint...
Defines the linear algebra or vector space interface.
virtual void zero()
Set to zero vector.
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,...
int main(int argc, char *argv[])
void printVector(const ROL::Vector< Real > &x, std::ostream &outStream)
constexpr auto dim