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