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 
38 typedef double RealT;
39 typedef H1VectorPrimal<RealT> PrimalStateVector;
40 typedef H1VectorDual<RealT> DualStateVector;
41 typedef L2VectorPrimal<RealT> PrimalControlVector;
42 typedef L2VectorDual<RealT> DualControlVector;
43 typedef H1VectorDual<RealT> PrimalConstraintVector;
44 typedef H1VectorPrimal<RealT> DualConstraintVector;
45 
46 int 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_StochasticObjective.hpp
ROL_Algorithm.hpp
ROL_Reduced_Objective_SimOpt.hpp
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
PrimalControlVector
L2VectorPrimal< RealT > PrimalControlVector
Definition: function/test_04.cpp:30
ROL_Stream.hpp
Defines a no-output stream class ROL::NullStream and a function makeStreamPtr which either wraps a re...
ROL::details::nullstream
basic_nullstream< char, std::char_traits< char >> nullstream
Definition: ROL_Stream.hpp:36
DualStateVector
H1VectorDual< RealT > DualStateVector
Definition: function/test_04.cpp:29
ROL_StatusTest.hpp
RealT
double RealT
Definition: function/constraint/test_01.cpp:29
ROL_RiskVector.hpp
DualControlVector
L2VectorDual< RealT > DualControlVector
Definition: function/test_04.cpp:31
ROL::Algorithm
Provides an interface to run optimization algorithms.
Definition: ROL_Algorithm.hpp:29
PrimalConstraintVector
H1VectorDual< RealT > PrimalConstraintVector
Definition: function/test_04.cpp:32
ROL_TrustRegionStep.hpp
example_06.hpp
main
int main(int argc, char *argv[])
Definition: function/constraint/test_01.cpp:31
PrimalStateVector
H1VectorPrimal< RealT > PrimalStateVector
Definition: function/test_04.cpp:28
ROL_MonteCarloGenerator.hpp
dim
constexpr auto dim
Definition: vector/test_11.cpp:23
DualConstraintVector
H1VectorPrimal< RealT > DualConstraintVector
Definition: function/test_04.cpp:33
Generated by   doxygen 1.8.5