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