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 
56 #include "ROL_Stream.hpp"
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 
68 typedef double RealT;
75 
76 int 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("Stochastic Component 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  ROL::Algorithm<RealT> algo("Trust Region",*parlist,false);
256  zp->zero();
257  algo.run(optProb,print0,*outStream0);
258  /*************************************************************************/
259  /************* PRINT CONTROL AND STATE TO SCREEN *************************/
260  /*************************************************************************/
261  if ( print0 ) {
262  std::ofstream ofs;
263  ofs.open("output_example_08.txt",std::ofstream::out);
264  for ( int i = 0; i < nx+2; i++ ) {
265  ofs << std::scientific << std::setprecision(10);
266  ofs << std::setw(20) << std::left << (RealT)i/((RealT)nx+1.0);
267  ofs << std::setw(20) << std::left << (*z_ptr)[i];
268  ofs << "\n";
269  }
270  ofs.close();
271  }
272  *outStream0 << "Scalar Parameter: " << optProb.getSolutionStatistic() << "\n\n";
273  }
274  catch (std::logic_error err) {
275  *outStream << err.what() << "\n";
276  errorFlag = -1000;
277  }; // end try
278 
279  comm->barrier();
280  if (errorFlag != 0)
281  std::cout << "End Result: TEST FAILED\n";
282  else
283  std::cout << "End Result: TEST PASSED\n";
284 
285  return 0;
286 }
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.
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
Provides an interface to run optimization algorithms.
H1VectorDual< RealT > PrimalConstraintVector
basic_nullstream< char, char_traits< char >> nullstream
Definition: ROL_Stream.hpp:72
int main(int argc, char *argv[])
void check(std::ostream &outStream=std::cout, const int numSteps=ROL_NUM_CHECKDERIV_STEPS, const int order=1)
H1VectorPrimal< RealT > PrimalStateVector
H1VectorPrimal< RealT > DualConstraintVector