ROL
example_06.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
49#include "ROL_Algorithm.hpp"
51#include "ROL_StatusTest.hpp"
52
54//#include "ROL_HMCRObjective.hpp"
55#include "ROL_RiskVector.hpp"
57#include "ROL_ParameterList.hpp"
58
60
61#include "ROL_Stream.hpp"
62#include "Teuchos_GlobalMPISession.hpp"
63#include "Teuchos_Comm.hpp"
64#include "Teuchos_DefaultComm.hpp"
65#include "Teuchos_CommHelpers.hpp"
66
67#include <iostream>
68#include <algorithm>
69
70#include "example_06.hpp"
71
72typedef double RealT;
79
80int main(int argc, char *argv[]) {
81
82 Teuchos::GlobalMPISession mpiSession(&argc, &argv);
83
84 auto comm = ROL::toPtr(Teuchos::DefaultComm<int>::getComm());
85
86 // This little trick lets us print to std::cout only if a (dummy) command-line argument is provided.
87 int iprint = argc - 1;
88 bool print = (iprint>0);
89 ROL::Ptr<std::ostream> outStream;
90 ROL::nullstream bhs; // outputs nothing
91 if (print)
92 outStream = ROL::makePtrFromRef(std::cout);
93 else
94 outStream = ROL::makePtrFromRef(bhs);
95
96 bool print0 = print && !comm->getRank();
97 ROL::Ptr<std::ostream> outStream0;
98 if (print0)
99 outStream0 = ROL::makePtrFromRef(std::cout);
100 else
101 outStream0 = ROL::makePtrFromRef(bhs);
102
103 int errorFlag = 0;
104
105 // *** Example body.
106
107 try {
108 /*************************************************************************/
109 /************* INITIALIZE BURGERS FEM CLASS ******************************/
110 /*************************************************************************/
111 int nx = 256; // Set spatial discretization.
112 RealT alpha = 1.e-3; // Set penalty parameter.
113 RealT nl = 1.0; // Nonlinearity parameter (1 = Burgers, 0 = linear).
114 RealT cH1 = 1.0; // Scale for derivative term in H1 norm.
115 RealT cL2 = 0.0; // Scale for mass term in H1 norm.
116 ROL::Ptr<BurgersFEM<RealT> > fem
117 = ROL::makePtr<BurgersFEM<RealT>>(nx,nl,cH1,cL2);
118 fem->test_inverse_mass(*outStream0);
119 fem->test_inverse_H1(*outStream0);
120 /*************************************************************************/
121 /************* INITIALIZE SIMOPT OBJECTIVE FUNCTION **********************/
122 /*************************************************************************/
123 ROL::Ptr<std::vector<RealT> > ud_ptr
124 = ROL::makePtr<std::vector<RealT>>(nx, 1.0);
125 ROL::Ptr<ROL::Vector<RealT> > ud
126 = ROL::makePtr<L2VectorPrimal<RealT>>(ud_ptr,fem);
127 ROL::Ptr<ROL::Objective_SimOpt<RealT> > pobj
128 = ROL::makePtr<Objective_BurgersControl<RealT>>(fem,ud,alpha);
129 /*************************************************************************/
130 /************* INITIALIZE SIMOPT EQUALITY CONSTRAINT *********************/
131 /*************************************************************************/
132 bool hess = true;
133 ROL::Ptr<ROL::Constraint_SimOpt<RealT> > pcon
134 = ROL::makePtr<Constraint_BurgersControl<RealT>>(fem,hess);
135 /*************************************************************************/
136 /************* INITIALIZE VECTOR STORAGE *********************************/
137 /*************************************************************************/
138 // INITIALIZE CONTROL VECTORS
139 ROL::Ptr<std::vector<RealT> > z_ptr
140 = ROL::makePtr<std::vector<RealT>>(nx+2, 1.0);
141 ROL::Ptr<std::vector<RealT> > gz_ptr
142 = ROL::makePtr<std::vector<RealT>>(nx+2, 1.0);
143 ROL::Ptr<std::vector<RealT> > yz_ptr
144 = ROL::makePtr<std::vector<RealT>>(nx+2, 1.0);
145 for (int i=0; i<nx+2; i++) {
146 (*yz_ptr)[i] = 2.0*random<RealT>(comm)-1.0;
147 }
148 ROL::Ptr<ROL::Vector<RealT> > zp
149 = ROL::makePtr<PrimalControlVector>(z_ptr,fem);
150 ROL::Ptr<ROL::Vector<RealT> > gzp
151 = ROL::makePtr<DualControlVector>(gz_ptr,fem);
152 ROL::Ptr<ROL::Vector<RealT> > yzp
153 = ROL::makePtr<PrimalControlVector>(yz_ptr,fem);
154 RealT zvar = 0.0*random<RealT>(comm);
155 RealT gvar = random<RealT>(comm);
156 RealT yvar = random<RealT>(comm);
157 ROL::Ptr<ROL::ParameterList> hmcrlist = ROL::makePtr<ROL::ParameterList>();
158 hmcrlist->sublist("SOL").sublist("Risk Measure").set("Name","HMCR");
159 ROL::RiskVector<RealT> z(hmcrlist,zp,zvar), g(hmcrlist,gzp,gvar), y(hmcrlist,yzp,yvar);
160 // INITIALIZE STATE VECTORS
161 ROL::Ptr<std::vector<RealT> > u_ptr
162 = ROL::makePtr<std::vector<RealT>>(nx, 1.0);
163 ROL::Ptr<std::vector<RealT> > gu_ptr
164 = ROL::makePtr<std::vector<RealT>>(nx, 1.0);
165 ROL::Ptr<ROL::Vector<RealT> > up
166 = ROL::makePtr<PrimalStateVector>(u_ptr,fem);
167 ROL::Ptr<ROL::Vector<RealT> > gup
168 = ROL::makePtr<DualStateVector>(gu_ptr,fem);
169 // INITIALIZE CONSTRAINT VECTORS
170 ROL::Ptr<std::vector<RealT> > c_ptr
171 = ROL::makePtr<std::vector<RealT>>(nx, 1.0);
172 ROL::Ptr<std::vector<RealT> > l_ptr
173 = ROL::makePtr<std::vector<RealT>>(nx, 1.0);
174 for (int i=0; i<nx; i++) {
175 (*l_ptr)[i] = random<RealT>(comm);
176 }
177 ROL::Ptr<ROL::Vector<RealT> > cp
178 = ROL::makePtr<PrimalConstraintVector>(c_ptr,fem);
179 ROL::Ptr<ROL::Vector<RealT> > lp
180 = ROL::makePtr<DualConstraintVector>(l_ptr,fem);
181 /*************************************************************************/
182 /************* INITIALIZE SAMPLE GENERATOR *******************************/
183 /*************************************************************************/
184 int dim = 4, nSamp = 1000;
185 std::vector<RealT> tmp(2,0.0); tmp[0] = -1.0; tmp[1] = 1.0;
186 std::vector<std::vector<RealT> > bounds(dim,tmp);
187 ROL::Ptr<ROL::BatchManager<RealT> > bman
188 = ROL::makePtr<L2VectorBatchManager<RealT,int>>(comm);
189 ROL::Ptr<ROL::SampleGenerator<RealT> > sampler
190 = ROL::makePtr<ROL::MonteCarloGenerator<RealT>>(
191 nSamp,bounds,bman,false,false,100);
192 /*************************************************************************/
193 /************* INITIALIZE RISK-AVERSE OBJECTIVE FUNCTION *****************/
194 /*************************************************************************/
195 bool storage = true, fdhess = false;
196 ROL::Ptr<ROL::Objective<RealT> > robj
197 = ROL::makePtr<ROL::Reduced_Objective_SimOpt<RealT>>(
198 pobj,pcon,up,zp,lp,gup,gzp,cp,storage,fdhess);
199 //RealT order = 2.0, prob = 0.95;
200 //ROL::Ptr<ROL::Objective<RealT> > obj
201 // = ROL::makePtr<ROL::HMCRObjective<RealT>>(
202 // robj,order,prob,sampler,storage);
203 hmcrlist->sublist("SOL").sublist("Risk Measure").sublist("HMCR").set("Order",2);
204 hmcrlist->sublist("SOL").sublist("Risk Measure").sublist("HMCR").set("Confidence Level",0.95);
205 hmcrlist->sublist("SOL").sublist("Risk Measure").sublist("HMCR").set("Convex Combination Parameter",0.0);
206 ROL::Ptr<ROL::Objective<RealT> > obj
207 = ROL::makePtr<ROL::StochasticObjective<RealT> >(robj,*hmcrlist,sampler);
208 /*************************************************************************/
209 /************* CHECK DERIVATIVES AND CONSISTENCY *************************/
210 /*************************************************************************/
211 // CHECK OBJECTIVE DERIVATIVES
212 bool derivcheck = false;
213 if (derivcheck) {
214 int nranks = sampler->numBatches();
215 for (int pid = 0; pid < nranks; pid++) {
216 if ( pid == sampler->batchID() ) {
217 for (int i = sampler->start(); i < sampler->numMySamples(); i++) {
218 *outStream << "Sample " << i << " Rank " << sampler->batchID() << "\n";
219 *outStream << "(" << sampler->getMyPoint(i)[0] << ", "
220 << sampler->getMyPoint(i)[1] << ", "
221 << sampler->getMyPoint(i)[2] << ", "
222 << sampler->getMyPoint(i)[3] << ")\n";
223 pcon->setParameter(sampler->getMyPoint(i));
224 pcon->checkSolve(*up,*zp,*cp,print,*outStream);
225 robj->setParameter(sampler->getMyPoint(i));
226 *outStream << "\n";
227 robj->checkGradient(*zp,*gzp,*yzp,print,*outStream);
228 robj->checkHessVec(*zp,*gzp,*yzp,print,*outStream);
229 *outStream << "\n\n";
230 }
231 }
232 comm->barrier();
233 }
234 }
235 obj->checkGradient(z,g,y,print0,*outStream0);
236 obj->checkHessVec(z,g,y,print0,*outStream0);
237 /*************************************************************************/
238 /************* RUN OPTIMIZATION ******************************************/
239 /*************************************************************************/
240 // READ IN XML INPUT
241 std::string filename = "input.xml";
242 auto parlist = ROL::getParametersFromXmlFile( filename );
243 // DEFINE ALGORITHM
244 ROL::Ptr<ROL::Step<RealT>>
245 step = ROL::makePtr<ROL::TrustRegionStep<RealT>>(*parlist);
246 ROL::Ptr<ROL::StatusTest<RealT>>
247 status = ROL::makePtr<ROL::StatusTest<RealT>>(*parlist);
248 ROL::Algorithm<RealT> algo(step,status,false);
249 // RUN OPTIMIZATION
250 z.zero();
251 algo.run(z, g, *obj, print0, *outStream0);
252 /*************************************************************************/
253 /************* PRINT CONTROL AND STATE TO SCREEN *************************/
254 /*************************************************************************/
255 *outStream0 << "\n";
256 for ( int i = 0; i < nx+2; i++ ) {
257 *outStream0 << std::scientific << std::setprecision(10);
258 *outStream0 << std::setw(20) << std::left << (RealT)i/((RealT)nx+1.0);
259 *outStream0 << std::setw(20) << std::left << (*z_ptr)[i];
260 *outStream0 << "\n";
261 }
262 *outStream0 << "\n";
263 *outStream0 << "Scalar Parameter: " << z.getStatistic(0) << "\n";
264 }
265 catch (std::logic_error& err) {
266 *outStream << err.what() << "\n";
267 errorFlag = -1000;
268 }; // end try
269
270 comm->barrier();
271 if (errorFlag != 0)
272 std::cout << "End Result: TEST FAILED\n";
273 else
274 std::cout << "End Result: TEST PASSED\n";
275
276 return 0;
277}
Defines a no-output stream class ROL::NullStream and a function makeStreamPtr which either wraps a re...
Provides an interface to run optimization algorithms.
virtual std::vector< std::string > run(Vector< Real > &x, Objective< Real > &obj, bool print=false, std::ostream &outStream=std::cout, bool printVectors=false, std::ostream &vectorStream=std::cout)
Run algorithm on unconstrained problems (Type-U). This is the primary Type-U interface.
int main(int argc, char *argv[])
L2VectorPrimal< RealT > PrimalControlVector
H1VectorPrimal< RealT > DualConstraintVector
H1VectorPrimal< RealT > PrimalStateVector
H1VectorDual< RealT > DualStateVector
L2VectorDual< RealT > DualControlVector
H1VectorDual< RealT > PrimalConstraintVector
double RealT
constexpr auto dim