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