ROL
example_08.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"
50
54#include "ROL_ParameterList.hpp"
55#include "ROL_Stream.hpp"
56
57#include "Teuchos_GlobalMPISession.hpp"
58#include "Teuchos_Comm.hpp"
59#include "Teuchos_DefaultComm.hpp"
60#include "Teuchos_CommHelpers.hpp"
61
62#include <iostream>
63#include <fstream>
64#include <algorithm>
65
66#include "example_08.hpp"
67
68typedef double RealT;
75
76int main(int argc, char *argv[]) {
77
78 Teuchos::GlobalMPISession mpiSession(&argc, &argv);
79 ROL::Ptr<const Teuchos::Comm<int> > comm
80 = ROL::toPtr(Teuchos::DefaultComm<int>::getComm());
81
82 // This little trick lets us print to std::cout only if a (dummy) command-line argument is provided.
83 int iprint = argc - 1;
84 bool print = (iprint>0); // && !(comm->getRank());
85 ROL::Ptr<std::ostream> outStream;
86 ROL::nullstream bhs; // outputs nothing
87 if (print)
88 outStream = ROL::makePtrFromRef(std::cout);
89 else
90 outStream = ROL::makePtrFromRef(bhs);
91
92 bool print0 = print && !comm->getRank();
93 ROL::Ptr<std::ostream> outStream0;
94 if (print0)
95 outStream0 = ROL::makePtrFromRef(std::cout);
96 else
97 outStream0 = ROL::makePtrFromRef(bhs);
98
99 int errorFlag = 0;
100
101 // *** Example body.
102
103 try {
104 /*************************************************************************/
105 /************* INITIALIZE BURGERS FEM CLASS ******************************/
106 /*************************************************************************/
107 int nx = 512; // Set spatial discretization.
108 RealT alpha = 1.e-3; // Set penalty parameter.
109 RealT nl = 1.0; // Nonlinearity parameter (1 = Burgers, 0 = linear).
110 RealT cH1 = 1.0; // Scale for derivative term in H1 norm.
111 RealT cL2 = 0.0; // Scale for mass term in H1 norm.
112 ROL::Ptr<BurgersFEM<RealT> > fem
113 = ROL::makePtr<BurgersFEM<RealT>>(nx,nl,cH1,cL2);
114 fem->test_inverse_mass(*outStream0);
115 fem->test_inverse_H1(*outStream0);
116 /*************************************************************************/
117 /************* INITIALIZE SIMOPT OBJECTIVE FUNCTION **********************/
118 /*************************************************************************/
119 ROL::Ptr<std::vector<RealT> > ud_ptr
120 = ROL::makePtr<std::vector<RealT>>(nx, 1.0);
121 ROL::Ptr<ROL::Vector<RealT> > ud
122 = ROL::makePtr<L2VectorPrimal<RealT>>(ud_ptr,fem);
123 ROL::Ptr<ROL::Objective_SimOpt<RealT> > pobj
124 = ROL::makePtr<Objective_BurgersControl<RealT>>(fem,ud,alpha);
125 /*************************************************************************/
126 /************* INITIALIZE SIMOPT EQUALITY CONSTRAINT *********************/
127 /*************************************************************************/
128 bool hess = true;
129 ROL::Ptr<ROL::Constraint_SimOpt<RealT> > pcon
130 = ROL::makePtr<Constraint_BurgersControl<RealT>>(fem,hess);
131 /*************************************************************************/
132 /************* INITIALIZE VECTOR STORAGE *********************************/
133 /*************************************************************************/
134 // INITIALIZE CONTROL VECTORS
135 ROL::Ptr<std::vector<RealT> > z_ptr
136 = ROL::makePtr<std::vector<RealT>>(nx+2, 1.0);
137 ROL::Ptr<std::vector<RealT> > gz_ptr
138 = ROL::makePtr<std::vector<RealT>>(nx+2, 1.0);
139 ROL::Ptr<std::vector<RealT> > yz_ptr
140 = ROL::makePtr<std::vector<RealT>>(nx+2, 1.0);
141 for (int i=0; i<nx+2; i++) {
142 (*yz_ptr)[i] = 2.0*random<RealT>(comm)-1.0;
143 }
144 ROL::Ptr<ROL::Vector<RealT> > zp
145 = ROL::makePtr<PrimalControlVector>(z_ptr,fem);
146 ROL::Ptr<ROL::Vector<RealT> > gzp
147 = ROL::makePtr<DualControlVector>(gz_ptr,fem);
148 ROL::Ptr<ROL::Vector<RealT> > yzp
149 = ROL::makePtr<PrimalControlVector>(yz_ptr,fem);
150 // INITIALIZE STATE VECTORS
151 ROL::Ptr<std::vector<RealT> > u_ptr
152 = ROL::makePtr<std::vector<RealT>>(nx, 1.0);
153 ROL::Ptr<std::vector<RealT> > gu_ptr
154 = ROL::makePtr<std::vector<RealT>>(nx, 1.0);
155 ROL::Ptr<ROL::Vector<RealT> > up
156 = ROL::makePtr<PrimalStateVector>(u_ptr,fem);
157 ROL::Ptr<ROL::Vector<RealT> > gup
158 = ROL::makePtr<DualStateVector>(gu_ptr,fem);
159 // INITIALIZE CONSTRAINT VECTORS
160 ROL::Ptr<std::vector<RealT> > c_ptr
161 = ROL::makePtr<std::vector<RealT>>(nx, 1.0);
162 ROL::Ptr<std::vector<RealT> > l_ptr
163 = ROL::makePtr<std::vector<RealT>>(nx, 1.0);
164 for (int i=0; i<nx; i++) {
165 (*l_ptr)[i] = random<RealT>(comm);
166 }
167 ROL::Ptr<ROL::Vector<RealT> > cp
168 = ROL::makePtr<PrimalConstraintVector>(c_ptr,fem);
169 ROL::Ptr<ROL::Vector<RealT> > lp
170 = ROL::makePtr<DualConstraintVector>(l_ptr,fem);
171 /*************************************************************************/
172 /************* INITIALIZE SAMPLE GENERATOR *******************************/
173 /*************************************************************************/
174 int dim = 4, nSamp = 1000;
175 std::vector<RealT> tmp(2,0.0); tmp[0] = -1.0; tmp[1] = 1.0;
176 std::vector<std::vector<RealT> > bounds(dim,tmp);
177 ROL::Ptr<ROL::BatchManager<RealT> > bman
178 = ROL::makePtr<L2VectorBatchManager<RealT,int>>(comm);
179 ROL::Ptr<ROL::SampleGenerator<RealT> > sampler
180 = ROL::makePtr<ROL::MonteCarloGenerator<RealT>>(
181 nSamp,bounds,bman,false,false,100);
182 /*************************************************************************/
183 /************* INITIALIZE REDUCED OBJECTIVE FUNCTION *********************/
184 /*************************************************************************/
185 bool storage = true, fdhess = false;
186 ROL::Ptr<ROL::Objective<RealT> > robj
187 = ROL::makePtr<ROL::Reduced_Objective_SimOpt<RealT>>(
188 pobj,pcon,up,zp,lp,gup,gzp,cp,storage,fdhess);
189 /*************************************************************************/
190 /************* INITIALIZE BOUND CONSTRAINTS ******************************/
191 /*************************************************************************/
192 std::vector<RealT> Zlo(nx+2,0.0), Zhi(nx+2,10.0);
193 for (int i = 0; i < nx+2; i++) {
194 if ( i < (int)((nx+2)/3) ) {
195 Zlo[i] = -1.0;
196 Zhi[i] = 1.0;
197 }
198 if ( i >= (int)((nx+2)/3) && i < (int)(2*(nx+2)/3) ) {
199 Zlo[i] = 1.0;
200 Zhi[i] = 5.0;
201 }
202 if ( i >= (int)(2*(nx+2)/3) ) {
203 Zlo[i] = 5.0;
204 Zhi[i] = 10.0;
205 }
206 }
207 ROL::Ptr<ROL::BoundConstraint<RealT> > Zbnd
208 = ROL::makePtr<L2BoundConstraint<RealT>>(Zlo,Zhi,fem);
209 /*************************************************************************/
210 /************* INITIALIZE OPTIMIZATION PROBLEM ***************************/
211 /*************************************************************************/
212 ROL::ParameterList SOLlist;
213 SOLlist.sublist("SOL").set("Type","Risk Averse");
214 SOLlist.sublist("SOL").set("Store Sampled Value and Gradient",storage);
215 SOLlist.sublist("SOL").sublist("Risk Measure").set("Name","KL Divergence");
216 SOLlist.sublist("SOL").sublist("Risk Measure").sublist("KL Divergence").set("Threshold",1.e-2);
217 ROL::OptimizationProblem<RealT> optProb(robj,zp,Zbnd);
218 optProb.setStochasticObjective(SOLlist,sampler);
219 /*************************************************************************/
220 /************* CHECK DERIVATIVES AND CONSISTENCY *************************/
221 /*************************************************************************/
222 // CHECK OBJECTIVE DERIVATIVES
223 bool derivcheck = false;
224 if (derivcheck) {
225 int nranks = sampler->numBatches();
226 for (int pid = 0; pid < nranks; pid++) {
227 if ( pid == sampler->batchID() ) {
228 for (int i = sampler->start(); i < sampler->numMySamples(); i++) {
229 *outStream << "Sample " << i << " Rank " << sampler->batchID() << "\n";
230 *outStream << "(" << sampler->getMyPoint(i)[0] << ", "
231 << sampler->getMyPoint(i)[1] << ", "
232 << sampler->getMyPoint(i)[2] << ", "
233 << sampler->getMyPoint(i)[3] << ")\n";
234 pcon->setParameter(sampler->getMyPoint(i));
235 pcon->checkSolve(*up,*zp,*cp,print,*outStream);
236 robj->setParameter(sampler->getMyPoint(i));
237 *outStream << "\n";
238 robj->checkGradient(*zp,*gzp,*yzp,print,*outStream);
239 robj->checkHessVec(*zp,*gzp,*yzp,print,*outStream);
240 *outStream << "\n\n";
241 }
242 }
243 comm->barrier();
244 }
245 }
246 optProb.check(*outStream0);
247 /*************************************************************************/
248 /************* RUN OPTIMIZATION ******************************************/
249 /*************************************************************************/
250 // READ IN XML INPUT
251 std::string filename = "input.xml";
252 auto parlist = ROL::getParametersFromXmlFile( filename );
253
254 // RUN OPTIMIZATION
255 parlist->sublist("Step").set("Type","Trust Region");
256 ROL::OptimizationSolver<RealT> solver(optProb,*parlist);
257 zp->zero();
258 solver.solve(*outStream0);
259 /*************************************************************************/
260 /************* PRINT CONTROL AND STATE TO SCREEN *************************/
261 /*************************************************************************/
262 if ( print0 ) {
263 std::ofstream ofs;
264 ofs.open("output_example_08.txt",std::ofstream::out);
265 for ( int i = 0; i < nx+2; i++ ) {
266 ofs << std::scientific << std::setprecision(10);
267 ofs << std::setw(20) << std::left << (RealT)i/((RealT)nx+1.0);
268 ofs << std::setw(20) << std::left << (*z_ptr)[i];
269 ofs << "\n";
270 }
271 ofs.close();
272 }
273 *outStream0 << "Scalar Parameter: " << optProb.getSolutionStatistic() << "\n\n";
274 }
275 catch (std::logic_error& err) {
276 *outStream << err.what() << "\n";
277 errorFlag = -1000;
278 }; // end try
279
280 comm->barrier();
281 if (errorFlag != 0)
282 std::cout << "End Result: TEST FAILED\n";
283 else
284 std::cout << "End Result: TEST PASSED\n";
285
286 return 0;
287}
Defines a no-output stream class ROL::NullStream and a function makeStreamPtr which either wraps a re...
void check(std::ostream &outStream=std::cout, const int numSteps=ROL_NUM_CHECKDERIV_STEPS, const int order=1)
Real getSolutionStatistic(int comp=0, int index=0)
Returns the statistic from the soluton vector.
void setStochasticObjective(ParameterList &parlist, const Ptr< SampleGenerator< Real > > &vsampler, const Ptr< SampleGenerator< Real > > &gsampler=nullPtr, const Ptr< SampleGenerator< Real > > &hsampler=nullPtr)
Provides a simplified interface for solving a wide range of optimization problems.
int solve(const ROL::Ptr< StatusTest< Real > > &status=ROL::nullPtr, const bool combineStatus=true)
Solve optimization problem with no iteration output.
int main(int argc, char *argv[])
Definition: example_08.cpp:76
L2VectorPrimal< RealT > PrimalControlVector
Definition: example_08.cpp:71
H1VectorPrimal< RealT > DualConstraintVector
Definition: example_08.cpp:74
H1VectorPrimal< RealT > PrimalStateVector
Definition: example_08.cpp:69
H1VectorDual< RealT > DualStateVector
Definition: example_08.cpp:70
L2VectorDual< RealT > DualControlVector
Definition: example_08.cpp:72
H1VectorDual< RealT > PrimalConstraintVector
Definition: example_08.cpp:73
double RealT
Definition: example_08.cpp:68
constexpr auto dim