ROL
function/test_17.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 
15 #include "ROL_Stream.hpp"
16 #include "Teuchos_GlobalMPISession.hpp"
17 
18 #include "ROL_Bounds.hpp"
19 #include "ROL_Constraint.hpp"
21 #include "ROL_StdVector.hpp"
22 #include "ROL_UnaryFunctions.hpp"
23 
24 typedef double RealT;
25 
27  RealT one(1);
28  a.axpy(-one, b);
29  a.applyUnary(ROL::Elementwise::AbsoluteValue<RealT>());
30  return a.reduce(ROL::Elementwise::ReductionMax<RealT>());
31 }
32 
33 int testRandomInputs(int numPoints, RealT tol,
34  ROL::Ptr<std::ostream> outStream) {
35 
36  int errorFlag = 0;
37  RealT errorInftyNorm;
38 
39  // Generate standard vectors that hold data.
40  ROL::Ptr<std::vector<RealT>> vsv
41  = ROL::makePtr<std::vector<RealT>>(numPoints, 0.0);
42  ROL::Ptr<std::vector<RealT>> xsv
43  = ROL::makePtr<std::vector<RealT>>(numPoints, 0.0);
44  ROL::Ptr<std::vector<RealT>> gsv
45  = ROL::makePtr<std::vector<RealT>>(numPoints, 0.0);
46  ROL::Ptr<std::vector<RealT>> lsv
47  = ROL::makePtr<std::vector<RealT>>(numPoints, 0.0);
48  ROL::Ptr<std::vector<RealT>> usv
49  = ROL::makePtr<std::vector<RealT>>(numPoints, 0.0);
50  // Include space for storing results.
51  ROL::Ptr<std::vector<RealT>> out1sv
52  = ROL::makePtr<std::vector<RealT>>(numPoints, 0.0);
53  ROL::Ptr<std::vector<RealT>> out2sv
54  = ROL::makePtr<std::vector<RealT>>(numPoints, 0.0);
55 
56  // Use these standard vectors to define ROL::StdVectors (or, in the case of lp
57  // and up, pointers to ROL::Vectors).
58  ROL::StdVector<RealT> v(vsv);
59  ROL::StdVector<RealT> x(xsv);
60  ROL::StdVector<RealT> g(gsv);
61  ROL::Ptr<ROL::Vector<RealT>> lp = ROL::makePtr<ROL::StdVector<RealT>>(lsv);
62  ROL::Ptr<ROL::Vector<RealT>> up = ROL::makePtr<ROL::StdVector<RealT>>(usv);
63  ROL::StdVector<RealT> out1(out1sv);
64  ROL::StdVector<RealT> out2(out2sv);
65 
66  // Initialize.
67  lp->randomize(-10.0, 10.0);
68  up->randomize( 0.0, 10.0);
69  up->plus(*lp);
70  x.randomize(-20.0, 20.0);
71  v.randomize(- 3.0, 3.0);
72  g.randomize(-20.0, 20.0);
73 
74  ROL::Bounds<RealT> boundsBC( lp , up );
75  ROL::StdBoundConstraint<RealT> stdvecBC(*lsv,*usv);
76  boundsBC.projectInterior(x);
77 
78  // Test 1 - Check that the Elementwise applyInverseScalingFunction does
79  // indeed apply a scale factor to the vector v.
80 
81  boundsBC.applyInverseScalingFunction(out1, v, x, g);
82  out1.applyUnary(ROL::Elementwise::Reciprocal<RealT>());
83  boundsBC.applyInverseScalingFunction(out2, out1, x, g);
84  out2.applyUnary(ROL::Elementwise::Reciprocal<RealT>());
85  errorInftyNorm = calcError(out2, v);
86  errorFlag += errorInftyNorm > tol;
87 
88  *outStream << std::endl;
89  *outStream << "Scaling Check at " << numPoints
90  << " Randomly Sampled Points -- " << std::endl
91  << " Infinity Norm of | v - 1/f(1/f(v)) | = "
92  << errorInftyNorm << std::endl;
93  *outStream << std::endl;
94 
95  // Test 2 - Use finite differences to check that the Elementwise
96  // applyScalingFunctionJacobian and applyInverseScalingFunction are
97  // consistent with each other.
98  // This test is meant to be visually inspected; it cannot cause the
99  // cpp file to fail.
100 
101  class TestWrapper : public ROL::Constraint<RealT> {
102  private:
103  ROL::Bounds<RealT> boundsBC_;
106 
107  public:
108  TestWrapper(ROL::Bounds<RealT>& boundsBC, ROL::StdVector<RealT>& g)
109  : boundsBC_(boundsBC), g_(g), v_(g) {
110  RealT one(1);
111  v_.setScalar(one);
112  }
113 
114  void value(ROL::Vector<RealT>& c, const ROL::Vector<RealT>& x,
115  RealT& tol) override {
116  boundsBC_.applyInverseScalingFunction(c, v_, x, g_);
117  c.applyUnary( ROL::Elementwise::Reciprocal<RealT>());
118  c.applyBinary(ROL::Elementwise::Multiply<RealT>(), g_);
119  }
120 
122  const ROL::Vector<RealT>& x, RealT& tol) override {
123  boundsBC_.applyScalingFunctionJacobian(jv, v, x, g_);
124  }
125  } testWrapper(boundsBC, g);
126 
127  // Use out1 and and out2 as working arrays to build a point comfortably
128  // within our bounds (so that the finite difference increments don't all step
129  // out). Larger values of gamma => larger separation between this point
130  // and our bounds.
131  RealT gamma = 1e-8;
132  out1.randomize(-1.0, 1.0);
133  out2.set(*up);
134  out2.axpy(-1,*lp);
135  out2.scale(1 - gamma);
136  out1.applyBinary(ROL::Elementwise::Multiply<RealT>(), out2);
137  out1.plus(*lp);
138  out1.plus(*up);
139  out1.scale(0.5); // the point at which we check the Jacobian
140 
141  *outStream << "Elementwise Jacobian Check:" << std::endl;
142  testWrapper.checkApplyJacobian(out1, v, out2, true, *outStream, 15);
143  *outStream << std::endl;
144 
145  // Test 3 - Check that applyInverseScalingFunction and
146  // applyScalingFunctionJacobian agree between the Elementwise and
147  // StdVector implementations.
148 
149  boundsBC.applyInverseScalingFunction(out1, v, x, g);
150  stdvecBC.applyInverseScalingFunction(out2, v, x, g);
151  errorInftyNorm = 100;
152  errorInftyNorm = calcError(out1, out2);
153  errorFlag += errorInftyNorm > tol;
154 
155  *outStream << "Consistency Check at " << numPoints
156  << " Randomly Sampled Points -- " << std::endl
157  << " Infinity Norm of | StdBoundConstraint - Elementwise |:"
158  << std::endl
159  << " Inverse = " << errorInftyNorm << std::endl;
160 
161  boundsBC.applyScalingFunctionJacobian(out1, v, x, g);
162  stdvecBC.applyScalingFunctionJacobian(out2, v, x, g);
163  errorInftyNorm = calcError(out1, out2);
164  errorFlag += errorInftyNorm > tol;
165 
166  *outStream << " Jacobian = " << errorInftyNorm << std::endl;
167  *outStream << std::endl;
168 
169  return errorFlag;
170 }
171 
172 int testCases(RealT tol, ROL::Ptr<std::ostream> outStream) {
173 
174  // Test 4 - Check the Elementwise and StdVector implementations of
175  // applyInverseScalingFunction and applyScalingFunctionJacobian on
176  // specific test cases.
177 
178  int numCases = 3;
179 
180  std::vector<RealT> ewErrors, svErrors;
181 
182  // Generate standard vectors that hold data.
183  ROL::Ptr<std::vector<RealT>> vsv
184  = ROL::makePtr<std::vector<RealT>>(numCases, 0.0);
185  ROL::Ptr<std::vector<RealT>> xsv
186  = ROL::makePtr<std::vector<RealT>>(numCases, 0.0);
187  ROL::Ptr<std::vector<RealT>> gsv
188  = ROL::makePtr<std::vector<RealT>>(numCases, 0.0);
189  ROL::Ptr<std::vector<RealT>> lsv
190  = ROL::makePtr<std::vector<RealT>>(numCases, 0.0);
191  ROL::Ptr<std::vector<RealT>> usv
192  = ROL::makePtr<std::vector<RealT>>(numCases, 0.0);
193  // Include space for storing results.
194  ROL::Ptr<std::vector<RealT>> resultp
195  = ROL::makePtr<std::vector<RealT>>(numCases, 0.0);
196  ROL::Ptr<std::vector<RealT>> targetp
197  = ROL::makePtr<std::vector<RealT>>(numCases, 0.0);
198 
199  // Use the standard vectors above to define ROL::StdVectors (or, in the
200  // case of lp and up, pointers to ROL::Vectors).
201  ROL::StdVector<RealT> v(vsv);
202  ROL::StdVector<RealT> x(xsv);
203  ROL::StdVector<RealT> g(gsv);
204  ROL::Ptr<ROL::Vector<RealT>> lp = ROL::makePtr<ROL::StdVector<RealT>>(lsv);
205  ROL::Ptr<ROL::Vector<RealT>> up = ROL::makePtr<ROL::StdVector<RealT>>(usv);
206  ROL::StdVector<RealT> result(resultp);
207  ROL::StdVector<RealT> target(targetp);
208 
209  // Problem 1
210  (*vsv)[0] = 4.0;
211  (*xsv)[0] = 1.9;
212  (*gsv)[0] = 0.5;
213  (*lsv)[0] = 0.0;
214  (*usv)[0] = 2.0;
215 
216  // Problem 2
217  (*vsv)[1] = -1.0;
218  (*xsv)[1] = 10.0;
219  (*gsv)[1] = 0.0002;
220  (*lsv)[1] = ROL::ROL_NINF<RealT>();
221  (*usv)[1] = ROL::ROL_INF<RealT>();
222 
223  // Problem 3
224  (*vsv)[2] = -0.0002;
225  (*xsv)[2] = 1.0;
226  (*gsv)[2] = 0.5;
227  (*lsv)[2] = 0.0;
228  (*usv)[2] = ROL::ROL_INF<RealT>();
229 
230  ROL::Bounds<RealT> boundsBC( lp, up );
231  ROL::StdBoundConstraint<RealT> stdvecBC(*lsv,*usv);
232 
233  // Expected results when applying the scaling function to v.
234  (*targetp)[0] = (*vsv)[0]*(*gsv)[0];
235  (*targetp)[1] = (*vsv)[1]*1.0;
236  (*targetp)[2] = (*vsv)[2]*1.0;
237 
238  target.applyBinary(ROL::Elementwise::DivideAndInvert<RealT>(), v);
239  target.applyBinary(ROL::Elementwise::Multiply<RealT>(), v);
240  boundsBC.applyInverseScalingFunction(result, v, x, g);
241  ewErrors.push_back(calcError(result, target));
242  stdvecBC.applyInverseScalingFunction(result, v, x, g);
243  svErrors.push_back(calcError(result, target));
244 
245  // Expected results when applying the Jacobian to v.
246  (*targetp)[0] = (*vsv)[0]*(*gsv)[0];
247  (*targetp)[1] = 0.0;
248  (*targetp)[2] = 0.0;
249  boundsBC.applyScalingFunctionJacobian(result, v, x, g);
250  ewErrors.push_back(calcError(result, target));
251  stdvecBC.applyScalingFunctionJacobian(result, v, x, g);
252  svErrors.push_back(calcError(result, target));
253 
254  *outStream << "Elementwise Test Case Errors (Infinity Norm):" << std::endl
255  << " Inverse = " << ewErrors[1] << std::endl
256  << " Jacobian = " << ewErrors[2] << std::endl;
257  *outStream << "StdBoundConstraint Test Case Errors (Infinity Norm):" << std::endl
258  << " Inverse = " << svErrors[1] << std::endl
259  << " Jacobian = " << svErrors[2] << std::endl;
260  *outStream << std::endl;
261 
262  RealT maxError = std::max(*std::max_element(svErrors.begin(), svErrors.end()),
263  *std::max_element(ewErrors.begin(), ewErrors.end()));
264  return maxError > tol;
265 }
266 
267 int main(int argc, char *argv[]) {
268 
269  Teuchos::GlobalMPISession mpiSession(&argc, &argv);
270 
271  // This little trick lets us print to std::cout only if a
272  // (dummy) command-line argument is provided.
273  int iprint = argc - 1;
274  ROL::Ptr<std::ostream> outStream;
275  ROL::nullstream bhs; // outputs nothing
276  if (iprint > 0)
277  outStream = ROL::makePtrFromRef(std::cout);
278  else
279  outStream = ROL::makePtrFromRef(bhs);
280 
281  int errorFlag = 0;
282 
283  // *** Test body.
284 
285  try {
286  RealT tol = 1e-8; // tolerance
287  int n = 1e+3; // number of random test points
288  errorFlag += testRandomInputs(n, tol, outStream);
289  errorFlag += testCases(tol, outStream);
290  }
291 
292  catch (std::logic_error& err) {
293  *outStream << err.what() << "\n";
294  errorFlag = -1000;
295  }; // end try
296 
297  if (errorFlag != 0)
298  std::cout << "End Result: TEST FAILED\n";
299  else
300  std::cout << "End Result: TEST PASSED\n";
301 
302  return 0;
303 }
Contains definitions for std::vector bound constraints.
virtual void applyJacobian(Vector< Real > &jv, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Apply the constraint Jacobian at , , to vector .
void axpy(const Real alpha, const Vector< Real > &x)
Compute where .
void scale(const Real alpha)
Compute where .
int testCases(RealT tol, ROL::Ptr< std::ostream > outStream)
virtual void axpy(const Real alpha, const Vector &x)
Compute where .
Definition: ROL_Vector.hpp:119
RealT calcError(ROL::Vector< RealT > &a, const ROL::Vector< RealT > &b)
void randomize(const Real l=0.0, const Real u=1.0)
Set vector to be uniform random between [l,u].
virtual Real reduce(const Elementwise::ReductionOp< Real > &r) const
Definition: ROL_Vector.hpp:221
virtual void applyBinary(const Elementwise::BinaryFunction< Real > &f, const Vector &x)
Definition: ROL_Vector.hpp:214
void applyBinary(const Elementwise::BinaryFunction< Real > &f, const Vector< Real > &x)
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:46
Defines a no-output stream class ROL::NullStream and a function makeStreamPtr which either wraps a re...
virtual void value(Vector< Real > &c, const Vector< Real > &x, Real &tol)=0
Evaluate the constraint operator at .
Provides the ROL::Vector interface for scalar values, to be used, for example, with scalar constraint...
void applyInverseScalingFunction(Vector< Real > &dv, const Vector< Real > &v, const Vector< Real > &x, const Vector< Real > &g) const override
Apply inverse scaling function.
Provides the elementwise interface to apply upper and lower bound constraints.
Definition: ROL_Bounds.hpp:25
virtual void applyUnary(const Elementwise::UnaryFunction< Real > &f)
Definition: ROL_Vector.hpp:208
void set(const Vector< Real > &x)
Set where .
void applyScalingFunctionJacobian(Vector< Real > &dv, const Vector< Real > &v, const Vector< Real > &x, const Vector< Real > &g) const override
Apply scaling function Jacobian.
basic_nullstream< char, char_traits< char >> nullstream
Definition: ROL_Stream.hpp:38
int main(int argc, char *argv[])
void applyUnary(const Elementwise::UnaryFunction< Real > &f)
void applyInverseScalingFunction(Vector< Real > &dv, const Vector< Real > &v, const Vector< Real > &x, const Vector< Real > &g) const override
Apply inverse scaling function.
void projectInterior(Vector< Real > &x) override
Project optimization variables into the interior of the feasible set.
int testRandomInputs(int numPoints, RealT tol, ROL::Ptr< std::ostream > outStream)
void applyScalingFunctionJacobian(Vector< Real > &dv, const Vector< Real > &v, const Vector< Real > &x, const Vector< Real > &g) const override
Apply scaling function Jacobian.
Defines the general constraint operator interface.