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 
62 template<class Real>
63 class NullObjective : public ROL::Objective<Real> {
65 public:
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 
77 template<class Real>
78 void 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 
90 typedef double RealT;
91 
92 int 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;
99  typedef ROL::BoundConstraint<RealT> BND;
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 }
Provides the interface to evaluate objective functions.
void gradient(V &g, const V &x, Real &tol)
Compute gradient.
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,upper].
virtual void zero()
Set to zero vector.
Definition: ROL_Vector.hpp:167
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:80
Vector< Real > V
Real value(const V &x, Real &tol)
Compute value.
void printVector(const ROL::Vector< Real > &x, std::ostream &outStream)
Provides the interface to evaluate the Interior Pointy log barrier penalty function with upper and lo...
Provides the interface to apply upper and lower bound constraints.
basic_nullstream< char, char_traits< char >> nullstream
Definition: ROL_Stream.hpp:72
void hessVec(V &hv, const V &v, const V &x, Real &tol)
Apply Hessian approximation to vector.
int main(int argc, char *argv[])
ROL::Vector< Real > V