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 "ROL_ParameterList.hpp"
45
46 #include "ROL_Stream.hpp"
47 #include "Teuchos_GlobalMPISession.hpp"
48
49 #include "ROL_ScaledStdVector.hpp"
50 #include "ROL_StdBoundConstraint.hpp"
51 #include "ROL_StdConstraint.hpp"
52 #include "ROL_Types.hpp"
53 #include "ROL_Algorithm.hpp"
54
55 #include "ROL_OptimizationProblem.hpp"
56 #include "ROL_OptimizationSolver.hpp"
57 #include "ROL_StdObjective.hpp"
58 #include "ROL_BatchManager.hpp"
59 #include "ROL_UserInputGenerator.hpp"
60
61 typedef double RealT;
62
63 template<class Real>
64 class ObjectiveEx12 : public ROL::StdObjective<Real> {
65 private:
66 const Real alpha_;
67 public:
ObjectiveEx12(const Real alpha=1e-4)68 ObjectiveEx12(const Real alpha = 1e-4) : alpha_(alpha) {}
69
value(const std::vector<Real> & x,Real & tol)70 Real value( const std::vector<Real> &x, Real &tol ) {
71 unsigned size = x.size();
72 Real val(0);
73 for ( unsigned i = 0; i < size; i++ ) {
74 val += static_cast<Real>(0.5)*alpha_*x[i]*x[i] + x[i];
75 }
76 return val;
77 }
78
gradient(std::vector<Real> & g,const std::vector<Real> & x,Real & tol)79 void gradient( std::vector<Real> &g, const std::vector<Real> &x, Real &tol ) {
80 unsigned size = x.size();
81 for ( unsigned i = 0; i < size; i++ ) {
82 g[i] = alpha_*x[i] + static_cast<Real>(1);
83 }
84 }
85
hessVec(std::vector<Real> & hv,const std::vector<Real> & v,const std::vector<Real> & x,Real & tol)86 void hessVec( std::vector<Real> &hv, const std::vector<Real> &v, const std::vector<Real> &x, Real &tol ) {
87 unsigned size = x.size();
88 for ( unsigned i = 0; i < size; i++ ) {
89 hv[i] = alpha_*v[i];
90 }
91 }
92
precond(std::vector<Real> & pv,const std::vector<Real> & v,const std::vector<Real> & x,Real & tol)93 void precond( std::vector<Real> &pv, const std::vector<Real> &v, const std::vector<Real> &x, Real &tol ) {
94 unsigned size = x.size();
95 for ( unsigned i = 0; i < size; i++ ) {
96 pv[i] = v[i]/alpha_;
97 }
98 }
99 };
100
101 template<class Real>
102 class ConstraintEx12 : public ROL::StdConstraint<Real> {
103 private:
104 const std::vector<Real> zeta_;
105
106 public:
ConstraintEx12(const std::vector<Real> & zeta)107 ConstraintEx12(const std::vector<Real> &zeta) : zeta_(zeta) {}
108
value(std::vector<Real> & c,const std::vector<Real> & x,Real & tol)109 void value( std::vector<Real> &c,
110 const std::vector<Real> &x,
111 Real &tol ) {
112 unsigned size = x.size();
113 const std::vector<Real> param = ROL::Constraint<Real>::getParameter();
114 for ( unsigned i = 0; i < size; ++i ) {
115 c[i] = std::exp(param[7])/zeta_[i] - std::exp(param[i])*x[i];
116 }
117 }
118
applyJacobian(std::vector<Real> & jv,const std::vector<Real> & v,const std::vector<Real> & x,Real & tol)119 void applyJacobian( std::vector<Real> &jv,
120 const std::vector<Real> &v,
121 const std::vector<Real> &x,
122 Real &tol ) {
123 unsigned size = x.size();
124 const std::vector<Real> param = ROL::Constraint<Real>::getParameter();
125 for ( unsigned i = 0; i < size; ++i ) {
126 jv[i] = -std::exp(param[i])*v[i];
127 }
128 }
129
applyAdjointJacobian(std::vector<Real> & ajv,const std::vector<Real> & v,const std::vector<Real> & x,Real & tol)130 void applyAdjointJacobian( std::vector<Real> &ajv,
131 const std::vector<Real> &v,
132 const std::vector<Real> &x,
133 Real &tol ) {
134 unsigned size = x.size();
135 const std::vector<Real> param = ROL::Constraint<Real>::getParameter();
136 for ( unsigned i = 0; i < size; ++i ) {
137 ajv[i] = -std::exp(param[i])*v[i];
138 }
139 }
140
applyAdjointHessian(std::vector<Real> & ahuv,const std::vector<Real> & u,const std::vector<Real> & v,const std::vector<Real> & x,Real & tol)141 void applyAdjointHessian( std::vector<Real> &ahuv,
142 const std::vector<Real> &u,
143 const std::vector<Real> &v,
144 const std::vector<Real> &x,
145 Real &tol ) {
146 unsigned size = x.size();
147 for ( unsigned i = 0; i < size; ++i ) {
148 ahuv[i] = static_cast<Real>(0);
149 }
150 }
151
applyPreconditioner(std::vector<Real> & pv,const std::vector<Real> & v,const std::vector<Real> & x,const std::vector<Real> & g,Real & tol)152 void applyPreconditioner( std::vector<Real> &pv,
153 const std::vector<Real> &v,
154 const std::vector<Real> &x,
155 const std::vector<Real> &g,
156 Real &tol ) {
157 unsigned size = x.size();
158 const std::vector<Real> param = ROL::Constraint<Real>::getParameter();
159 for ( unsigned i = 0; i < size; ++i ) {
160 pv[i] = v[i]/std::pow(std::exp(param[i]),2);
161 }
162 }
163 };
164
165 template<class Real>
setUpAndSolve(ROL::ParameterList & list,ROL::Ptr<ROL::Objective<Real>> & obj,ROL::Ptr<ROL::Vector<Real>> & x,ROL::Ptr<ROL::BoundConstraint<Real>> & bnd,ROL::Ptr<ROL::Constraint<Real>> & con,ROL::Ptr<ROL::Vector<Real>> & mul,ROL::Ptr<ROL::BoundConstraint<Real>> & ibnd,ROL::Ptr<ROL::SampleGenerator<Real>> & sampler,std::ostream & outStream)166 Real setUpAndSolve(ROL::ParameterList &list,
167 ROL::Ptr<ROL::Objective<Real> > &obj,
168 ROL::Ptr<ROL::Vector<Real> > &x,
169 ROL::Ptr<ROL::BoundConstraint<Real> > &bnd,
170 ROL::Ptr<ROL::Constraint<Real> > &con,
171 ROL::Ptr<ROL::Vector<Real> > &mul,
172 ROL::Ptr<ROL::BoundConstraint<Real> > &ibnd,
173 ROL::Ptr<ROL::SampleGenerator<Real> > &sampler,
174 std::ostream & outStream) {
175 ROL::OptimizationProblem<Real> optProblem(obj,x,bnd,con,mul,ibnd);
176 optProblem.setRiskNeutralInequality(sampler,sampler->getBatchManager());
177 optProblem.check(outStream);
178 // Run ROL algorithm
179 ROL::OptimizationSolver<Real> optSolver(optProblem, list);
180 optSolver.solve(outStream);
181 ROL::Ptr<ROL::Objective<Real> > robj = optProblem.getObjective();
182 Real tol(1.e-8);
183 return robj->value(*(optProblem.getSolutionVector()),tol);
184 }
185
186 template<class Real>
printSolution(const std::vector<Real> & x,std::ostream & outStream)187 void printSolution(const std::vector<Real> &x,
188 std::ostream & outStream) {
189 unsigned dim = x.size();
190 outStream << "x = (";
191 for ( unsigned i = 0; i < dim-1; i++ ) {
192 outStream << x[i] << ", ";
193 }
194 outStream << x[dim-1] << ")\n";
195 }
196
main(int argc,char * argv[])197 int main(int argc, char* argv[]) {
198
199 Teuchos::GlobalMPISession mpiSession(&argc, &argv);
200
201 // This little trick lets us print to std::cout only if a (dummy) command-line argument is provided.
202 int iprint = argc - 1;
203 ROL::Ptr<std::ostream> outStream;
204 ROL::nullstream bhs; // outputs nothing
205 if (iprint > 0)
206 outStream = ROL::makePtrFromRef(std::cout);
207 else
208 outStream = ROL::makePtrFromRef(bhs);
209
210 int errorFlag = 0;
211
212 try {
213 /**********************************************************************************************/
214 /************************* CONSTRUCT ROL ALGORITHM ********************************************/
215 /**********************************************************************************************/
216 // Get ROL parameterlist
217 std::string filename = "input_12.xml";
218
219 auto parlist = ROL::getParametersFromXmlFile( filename );
220 ROL::ParameterList list = *parlist;
221 /**********************************************************************************************/
222 /************************* CONSTRUCT SOL COMPONENTS *******************************************/
223 /**********************************************************************************************/
224 // Build vectors
225 const RealT zero(0), half(0.5), one(1), two(2);
226 unsigned dim = 7;
227 ROL::Ptr<std::vector<RealT> > x_ptr;
228 x_ptr = ROL::makePtr<std::vector<RealT>>(dim,zero);
229 ROL::Ptr<ROL::Vector<RealT> > x;
230 x = ROL::makePtr<ROL::StdVector<RealT>>(x_ptr);
231 // Build samplers
232 int nSamp = 50;
233 unsigned sdim = dim + 1;
234 ROL::Ptr<ROL::BatchManager<RealT> > bman
235 = ROL::makePtr<ROL::BatchManager<RealT>>();
236 ROL::Ptr<ROL::SampleGenerator<RealT> > sampler
237 = ROL::makePtr<ROL::UserInputGenerator<RealT>>("points.txt","weights.txt",nSamp,sdim,bman);
238 // Build objective function
239 ROL::Ptr<ROL::Objective<RealT> > obj
240 = ROL::makePtr<ObjectiveEx12<RealT>>();
241 // Build bound constraints
242 std::vector<RealT> lx(dim,half);
243 std::vector<RealT> ux(dim,two);
244 ROL::Ptr<ROL::BoundConstraint<RealT> > bnd
245 = ROL::makePtr<ROL::StdBoundConstraint<RealT>>(lx,ux);
246 // Build inequality constraint
247 std::vector<RealT> zeta(dim,one/static_cast<RealT>(std::sqrt(3.)));
248 zeta[0] *= half;
249 zeta[1] *= half;
250 ROL::Ptr<ROL::Constraint<RealT> > con
251 = ROL::makePtr<ConstraintEx12<RealT>>(zeta);
252 // Build inequality bound constraint
253 std::vector<RealT> lc(dim,ROL::ROL_NINF<RealT>());
254 std::vector<RealT> uc(dim,zero);
255 ROL::Ptr<ROL::BoundConstraint<RealT> > ibnd
256 = ROL::makePtr<ROL::StdBoundConstraint<RealT>>(lc,uc);
257 ibnd->deactivateLower();
258 // Build multipliers
259 std::vector<RealT> mean(sdim,zero), gmean(sdim,zero);
260 for (int i = 0; i < sampler->numMySamples(); ++i) {
261 std::vector<RealT> pt = sampler->getMyPoint(i);
262 RealT wt = sampler->getMyWeight(i);
263 for (unsigned j = 0; j < sdim; ++j) {
264 mean[j] += wt*std::exp(pt[j]);
265 }
266 }
267 sampler->sumAll(&mean[0],&gmean[0],sdim);
268 ROL::Ptr<std::vector<RealT> > scaling_vec
269 = ROL::makePtr<std::vector<RealT>>(dim,zero);
270 for (unsigned i = 0; i < dim; ++i) {
271 RealT cl = std::abs(gmean[dim]/zeta[i] - gmean[i]*lx[i]);
272 RealT cu = std::abs(gmean[dim]/zeta[i] - gmean[i]*ux[i]);
273 RealT scale = static_cast<RealT>(1e2)/std::pow(std::max(cl,cu),two);
274 (*scaling_vec)[i] = (scale > std::sqrt(ROL::ROL_EPSILON<RealT>()))
275 ? scale : one;
276 }
277 ROL::Ptr<std::vector<RealT> > l_ptr
278 = ROL::makePtr<std::vector<RealT>>(dim,one);
279 ROL::Ptr<ROL::Vector<RealT> > l
280 = ROL::makePtr<ROL::DualScaledStdVector<RealT>>(l_ptr,scaling_vec);
281
282 /**********************************************************************************************/
283 /************************* SUPER QUANTILE QUADRANGLE ******************************************/
284 /**********************************************************************************************/
285 RealT val = setUpAndSolve(list,obj,x,bnd,con,l,ibnd,sampler,*outStream);
286 *outStream << "Computed Solution" << std::endl;
287 printSolution<RealT>(*x_ptr,*outStream);
288
289 // Compute exact solution
290 ROL::Ptr<std::vector<RealT> > ximean, gximean, xsol;
291 ximean = ROL::makePtr<std::vector<RealT>>(dim+1,0);
292 gximean = ROL::makePtr<std::vector<RealT>>(dim+1);
293 xsol = ROL::makePtr<std::vector<RealT>>(dim);
294 for (int i = 0; i < sampler->numMySamples(); ++i) {
295 std::vector<RealT> pt = sampler->getMyPoint(i);
296 RealT wt = sampler->getMyWeight(i);
297 for (unsigned j = 0; j < dim+1; ++j) {
298 (*ximean)[j] += wt*std::exp(pt[j]);
299 }
300 }
301 bman->reduceAll(&(*ximean)[0],&(*gximean)[0],dim+1,ROL::Elementwise::ReductionMax<RealT>());
302 for (unsigned i = 0; i < dim; ++i) {
303 (*xsol)[i] = std::max(half,((*gximean)[dim]/(*gximean)[i])/zeta[i]);
304 }
305 ROL::Ptr<ROL::Vector<RealT> > xtrue
306 = ROL::makePtr<ROL::StdVector<RealT>>(xsol);
307 *outStream << "True Solution" << std::endl;
308 printSolution<RealT>(*xsol,*outStream);
309 xtrue->axpy(-one,*x);
310
311 RealT error = xtrue->norm();
312 *outStream << std::right
313 << std::setw(20) << "obj val"
314 << std::setw(20) << "error"
315 << std::endl;
316 *outStream << std::fixed << std::setprecision(0) << std::right
317 << std::scientific << std::setprecision(11) << std::right
318 << std::setw(20) << val
319 << std::setw(20) << error
320 << std::endl;
321
322 errorFlag = (error < std::sqrt(ROL::ROL_EPSILON<RealT>())) ? 0 : 1;
323 }
324 catch (std::logic_error& err) {
325 *outStream << err.what() << "\n";
326 errorFlag = -1000;
327 }; // end try
328
329 if (errorFlag != 0)
330 std::cout << "End Result: TEST FAILED\n";
331 else
332 std::cout << "End Result: TEST PASSED\n";
333
334 return 0;
335 }
336