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