ROL
function/test_16.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 
50 #include "ROL_Bounds.hpp"
51 #include "ROL_ScaledStdVector.hpp"
52 #include "ROL_StdConstraint.hpp"
53 
54 #include "ROL_Stream.hpp"
55 #include "Teuchos_GlobalMPISession.hpp"
56 
57 template<typename Real>
58 class con2d : public ROL::StdConstraint<Real> {
59 public:
60  void value(std::vector<Real> &c, const std::vector<Real> &x, Real &tol) {
61  c[0] = x[0]+x[1]-static_cast<Real>(1);
62  }
63  void applyJacobian(std::vector<Real> &jv, const std::vector<Real> &v, const std::vector<Real> &x, Real &tol) {
64  jv[0] = v[0]+v[1];
65  }
66  void applyAdjointJacobian(std::vector<Real> &ajv, const std::vector<Real> &v, const std::vector<Real> &x, Real &tol) {
67  ajv[0] = v[0];
68  ajv[1] = v[0];
69  }
70 };
71 
72 typedef double RealT;
73 
74 int main(int argc, char *argv[]) {
75 
76  Teuchos::GlobalMPISession mpiSession(&argc, &argv);
77 
78  // This little trick lets us print to std::cout only if a
79  // (dummy) command-line argument is provided.
80  int iprint = argc - 1;
81  ROL::Ptr<std::ostream> outStream;
82  ROL::nullstream bhs; // outputs nothing
83  if (iprint > 0)
84  outStream = ROL::makePtrFromRef(std::cout);
85  else
86  outStream = ROL::makePtrFromRef(bhs);
87 
88  int errorFlag = 0;
89 
90  try {
91  const RealT zero(0), half(0.5), one(1);
92  RealT tol = std::sqrt(ROL::ROL_EPSILON<RealT>());
93  RealT err(0);
94  ROL::Ptr<con2d<RealT>> con = ROL::makePtr<con2d<RealT>>();
96  ROL::ParameterList list;
97  list.sublist("General").set("Output Level",2);
98  list.sublist("General").sublist("Polyhedral Projection").set("Type","Dai-Fletcher");
99  //list.sublist("General").sublist("Polyhedral Projection").set("Type","Ridders");
100  //list.sublist("General").sublist("Polyhedral Projection").set("Type","Dykstra");
101  //list.sublist("General").sublist("Polyhedral Projection").set("Type","Semismooth Newton");
102 
103  ROL::Ptr<std::vector<RealT>> yptr = ROL::makePtr<std::vector<RealT>>(2);
104  (*yptr)[0] = static_cast<RealT>(10)*(static_cast<RealT>(rand())/static_cast<RealT>(RAND_MAX)-half);
105  (*yptr)[1] = static_cast<RealT>(rand())/static_cast<RealT>(RAND_MAX);
106  //(*yptr)[1] = static_cast<RealT>(10)*(static_cast<RealT>(rand())/static_cast<RealT>(RAND_MAX)-half);
107  ROL::StdVector<RealT> y(yptr);
108 
109  ROL::Ptr<std::vector<RealT>> xptr = ROL::makePtr<std::vector<RealT>>(2);
110  (*xptr)[0] = (*yptr)[0];
111  (*xptr)[1] = (*yptr)[1];
112  ROL::StdVector<RealT> x(xptr);
113 
114  ROL::Ptr<std::vector<RealT>> Pxptr = ROL::makePtr<std::vector<RealT>>(2,0.0);
115  (*Pxptr)[0] = (*yptr)[0];
116  (*Pxptr)[1] = (*yptr)[1];
117  ROL::StdVector<RealT> Px(Pxptr);
118 
119  ROL::Ptr<ROL::Vector<RealT>> l0 = x.clone(); l0->setScalar(static_cast<RealT>(0));
120  ROL::Ptr<ROL::Vector<RealT>> u0 = x.clone(); u0->setScalar(static_cast<RealT>(1));
121  ROL::Ptr<ROL::Bounds<RealT>> bnd0 = ROL::makePtr<ROL::Bounds<RealT>>(l0,u0);
122 
123  ROL::Ptr<ROL::PolyhedralProjection<RealT>> pp0 = ROL::PolyhedralProjectionFactory<RealT>(x,x.dual(),bnd0,con,r,r.dual(),list);
124  pp0->project(Px,*outStream);
125 
126  ROL::Ptr<std::vector<RealT>> x0ptr = ROL::makePtr<std::vector<RealT>>(2);
127  RealT k0 = std::max(zero,std::min(one,half*(one+(*yptr)[0]-(*yptr)[1])));
128  (*x0ptr)[0] = k0;
129  (*x0ptr)[1] = one-k0;
130  ROL::StdVector<RealT> x0(x0ptr);
131 
132  ROL::StdVector<RealT> e0(2);
133 
134  *outStream << std::setprecision(6) << std::scientific << std::endl;
135  *outStream << " x[0] = " << (*xptr)[0] << " x[1] = " << (*xptr)[1] << std::endl;
136  *outStream << " Px[0] = " << (*Pxptr)[0] << " Px[1] = " << (*Pxptr)[1] << std::endl;
137  *outStream << " x*[0] = " << (*x0ptr)[0] << " x*[1] = " << (*x0ptr)[1] << std::endl;
138 
139  e0.set(x0); e0.axpy(static_cast<RealT>(-1),Px);
140  err = e0.norm();
141  *outStream << " Error in Euclidean Projection: " << err << std::endl;
142 
143  e0.set(x); e0.axpy(static_cast<RealT>(-1),x0);
144  *outStream << " ||x*-x||^2 = " << e0.norm() << std::endl;
145 
146  e0.set(x); e0.axpy(static_cast<RealT>(-1),Px);
147  *outStream << " ||Px-x||^2 = " << e0.norm() << std::endl << std::endl;
148 
149  errorFlag += (err > tol);
150 
151  ROL::Ptr<std::vector<RealT>> dptr = ROL::makePtr<std::vector<RealT>>(2);
152  (*dptr)[0] = static_cast<RealT>(1)+static_cast<RealT>(2)*static_cast<RealT>(rand())/static_cast<RealT>(RAND_MAX);
153  (*dptr)[1] = static_cast<RealT>(1)+static_cast<RealT>(5)*static_cast<RealT>(rand())/static_cast<RealT>(RAND_MAX);
154 
155  ROL::Ptr<std::vector<RealT>> x1ptr = ROL::makePtr<std::vector<RealT>>(2);
156  RealT k1 = std::max(zero,std::min(one,((*dptr)[1]*(one-(*yptr)[1])+(*dptr)[0]*(*yptr)[0])/((*dptr)[0]+(*dptr)[1])));
157  (*x1ptr)[0] = k1;
158  (*x1ptr)[1] = one-k1;
159  ROL::PrimalScaledStdVector<RealT> x1(x1ptr,dptr);
160 
161  ROL::Ptr<std::vector<RealT>> zptr = ROL::makePtr<std::vector<RealT>>(2);
162  (*zptr)[0] = (*yptr)[0];
163  (*zptr)[1] = (*yptr)[1];
165 
166  ROL::Ptr<std::vector<RealT>> Pzptr = ROL::makePtr<std::vector<RealT>>(2,0.0);
167  (*Pzptr)[0] = (*yptr)[0];
168  (*Pzptr)[1] = (*yptr)[1];
169  ROL::PrimalScaledStdVector<RealT> Pz(Pzptr,dptr);
170 
171  ROL::Ptr<ROL::Vector<RealT>> l1 = z.clone(); l1->setScalar(static_cast<RealT>(0));
172  ROL::Ptr<ROL::Vector<RealT>> u1 = z.clone(); u1->setScalar(static_cast<RealT>(1));
173  ROL::Ptr<ROL::Bounds<RealT>> bnd1 = ROL::makePtr<ROL::Bounds<RealT>>(l1,u1);
174 
175  ROL::Ptr<ROL::PolyhedralProjection<RealT>> pp1 = ROL::PolyhedralProjectionFactory<RealT>(z,z.dual(),bnd1,con,r,r.dual(),list);
176  pp1->project(Pz,*outStream);
177 
178  ROL::Ptr<std::vector<RealT>> e1ptr = ROL::makePtr<std::vector<RealT>>(2);
179  ROL::PrimalScaledStdVector<RealT> e1(e1ptr,dptr);
180 
181  *outStream << std::endl;
182  *outStream << " x[0] = " << (*zptr)[0] << " x[1] = " << (*zptr)[1] << std::endl;
183  *outStream << " Px[0] = " << (*Pzptr)[0] << " Px[1] = " << (*Pzptr)[1] << std::endl;
184  *outStream << " x*[0] = " << (*x1ptr)[0] << " x*[1] = " << (*x1ptr)[1] << std::endl;
185 
186  e1.set(x1); e1.axpy(static_cast<RealT>(-1),Pz);
187  err = e1.norm();
188  *outStream << " Error in Scaled Projection: " << err << std::endl;
189 
190  e1.set(z); e1.axpy(static_cast<RealT>(-1),x1);
191  *outStream << " ||x*-x||^2 = " << e1.norm() << std::endl;
192 
193  e1.set(z); e1.axpy(static_cast<RealT>(-1),Pz);
194  *outStream << " ||Px-x||^2 = " << e1.norm() << std::endl << std::endl;
195 
196  errorFlag += (err > tol);
197  }
198 
199  catch (std::logic_error& err) {
200  *outStream << err.what() << "\n";
201  errorFlag = -1000;
202  }; // end try
203 
204  if (errorFlag != 0)
205  std::cout << "End Result: TEST FAILED\n";
206  else
207  std::cout << "End Result: TEST PASSED\n";
208 
209  return 0;
210 }
211 
void axpy(const Real alpha, const Vector< Real > &x)
Compute where .
basic_nullstream< char, char_traits< char > > nullstream
Definition: ROL_Stream.hpp:72
double RealT
int main(int argc, char *argv[])
Defines the equality constraint operator interface for StdVectors.
Defines a no-output stream class ROL::NullStream and a function makeStreamPtr which either wraps a re...
Objective_SerialSimOpt(const Ptr< Obj > &obj, const V &ui) z0_ zero()
virtual const Vector & dual() const
Return dual representation of , for example, the result of applying a Riesz map, or change of basis...
Definition: ROL_Vector.hpp:226
Provides the ROL::Vector interface for scalar values, to be used, for example, with scalar constraint...
void applyJacobian(std::vector< Real > &jv, const std::vector< Real > &v, const std::vector< Real > &x, Real &tol)
Provides the std::vector implementation of the ROL::Vector interface that handles scalings in the inn...
void set(const Vector< Real > &x)
Set where .
virtual Ptr< Vector< Real > > clone() const
Clone to make a new (uninitialized) vector.
Ptr< Vector< Real > > clone() const
Clone to make a new (uninitialized) vector.
const Vector< Real > & dual() const
Return dual representation of , for example, the result of applying a Riesz map, or change of basis...
void applyAdjointJacobian(std::vector< Real > &ajv, const std::vector< Real > &v, const std::vector< Real > &x, Real &tol)
void value(std::vector< Real > &c, const std::vector< Real > &x, Real &tol)
Real norm() const
Returns where .