ROL
example_03.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 "example_03.hpp"
50 
51 typedef double RealT;
52 
53 int main(int argc, char *argv[]) {
54 
55  Teuchos::GlobalMPISession mpiSession(&argc, &argv);
56 
57  // This little trick lets us print to std::cout only if a (dummy) command-line argument is provided.
58  int iprint = argc - 1;
59  Teuchos::RCP<std::ostream> outStream;
60  Teuchos::oblackholestream bhs; // outputs nothing
61  if (iprint > 0)
62  outStream = Teuchos::rcp(&std::cout, false);
63  else
64  outStream = Teuchos::rcp(&bhs, false);
65 
66  int errorFlag = 0;
67 
68  // *** Example body.
69 
70  try {
71  // Initialize full objective function.
72  int nx = 20; // Set spatial discretization.
73  int nt = 20; // Set temporal discretization.
74  RealT T = 1.0; // Set end time.
75  RealT alpha = 0.05; // Set penalty parameter.
76  RealT nu = 1.e-2; // Set viscosity parameter.
77  Objective_BurgersControl<RealT> obj(alpha,nx,nt,T);
78  // Initialize equality constraints
80  // Initialize iteration vectors.
81  Teuchos::RCP<std::vector<RealT> > z_rcp = Teuchos::rcp( new std::vector<RealT> ((nx+2)*(nt+1), 1.0) );
82  Teuchos::RCP<std::vector<RealT> > gz_rcp = Teuchos::rcp( new std::vector<RealT> ((nx+2)*(nt+1), 1.0) );
83  Teuchos::RCP<std::vector<RealT> > yz_rcp = Teuchos::rcp( new std::vector<RealT> ((nx+2)*(nt+1), 1.0) );
84  for (int i=0; i<(nx+2)*(nt+1); i++) {
85  (*z_rcp)[i] = (RealT)rand()/(RealT)RAND_MAX;
86  (*yz_rcp)[i] = (RealT)rand()/(RealT)RAND_MAX;
87  }
88  ROL::StdVector<RealT> z(z_rcp);
89  ROL::StdVector<RealT> gz(gz_rcp);
90  ROL::StdVector<RealT> yz(yz_rcp);
91  Teuchos::RCP<ROL::Vector<RealT> > zp = Teuchos::rcp(&z,false);
92  Teuchos::RCP<ROL::Vector<RealT> > gzp = Teuchos::rcp(&gz,false);
93  Teuchos::RCP<ROL::Vector<RealT> > yzp = Teuchos::rcp(&yz,false);
94 
95  Teuchos::RCP<std::vector<RealT> > u_rcp = Teuchos::rcp( new std::vector<RealT> (nx*nt, 1.0) );
96  Teuchos::RCP<std::vector<RealT> > gu_rcp = Teuchos::rcp( new std::vector<RealT> (nx*nt, 1.0) );
97  Teuchos::RCP<std::vector<RealT> > yu_rcp = Teuchos::rcp( new std::vector<RealT> (nx*nt, 1.0) );
98  for (int i=0; i<nx*nt; i++) {
99  (*u_rcp)[i] = (RealT)rand()/(RealT)RAND_MAX;
100  (*yu_rcp)[i] = (RealT)rand()/(RealT)RAND_MAX;
101  }
102  ROL::StdVector<RealT> u(u_rcp);
103  ROL::StdVector<RealT> gu(gu_rcp);
104  ROL::StdVector<RealT> yu(yu_rcp);
105  Teuchos::RCP<ROL::Vector<RealT> > up = Teuchos::rcp(&u,false);
106  Teuchos::RCP<ROL::Vector<RealT> > gup = Teuchos::rcp(&gu,false);
107  Teuchos::RCP<ROL::Vector<RealT> > yup = Teuchos::rcp(&yu,false);
108 
109  Teuchos::RCP<std::vector<RealT> > c_rcp = Teuchos::rcp( new std::vector<RealT> (nx*nt, 1.0) );
110  Teuchos::RCP<std::vector<RealT> > l_rcp = Teuchos::rcp( new std::vector<RealT> (nx*nt, 1.0) );
111  ROL::StdVector<RealT> c(c_rcp);
112  ROL::StdVector<RealT> l(l_rcp);
113 
114  ROL::Vector_SimOpt<RealT> x(up,zp);
115  ROL::Vector_SimOpt<RealT> g(gup,gzp);
116  ROL::Vector_SimOpt<RealT> y(yup,yzp);
117  // Check derivatives.
118  obj.checkGradient(x,x,y,true,*outStream);
119  obj.checkHessVec(x,x,y,true,*outStream);
120  con.checkApplyJacobian(x,y,c,true,*outStream);
121  con.checkApplyAdjointJacobian(x,yu,c,x,true,*outStream);
122  con.checkApplyAdjointHessian(x,yu,y,x,true,*outStream);
123  // Check consistency of Jacobians and adjoint Jacobians.
124  con.checkAdjointConsistencyJacobian_1(c,yu,u,z,true,*outStream);
125  con.checkAdjointConsistencyJacobian_2(c,yz,u,z,true,*outStream);
126  // Check consistency of solves.
127  con.checkSolve(u,z,c,true,*outStream);
128  con.checkInverseJacobian_1(c,yu,u,z,true,*outStream);
129  con.checkInverseAdjointJacobian_1(yu,c,u,z,true,*outStream);
130 
131  // Initialize reduced objective function.
132  Teuchos::RCP<std::vector<RealT> > p_rcp = Teuchos::rcp( new std::vector<RealT> (nx*nt, 1.0) );
133  ROL::StdVector<RealT> p(p_rcp);
134  Teuchos::RCP<ROL::Vector<RealT> > pp = Teuchos::rcp(&p,false);
135  Teuchos::RCP<ROL::Objective_SimOpt<RealT> > pobj = Teuchos::rcp(&obj,false);
136  Teuchos::RCP<ROL::EqualityConstraint_SimOpt<RealT> > pcon = Teuchos::rcp(&con,false);
137  ROL::Reduced_Objective_SimOpt<RealT> robj(pobj,pcon,up,pp);
138  // Check derivatives.
139  robj.checkGradient(z,z,yz,true,*outStream);
140  robj.checkHessVec(z,z,yz,true,*outStream);
141  // Get input parameter list.
142  std::string filename = "input.xml";
143  Teuchos::RCP<Teuchos::ParameterList> parlist = Teuchos::rcp( new Teuchos::ParameterList() );
144  Teuchos::updateParametersFromXmlFile( filename, parlist.ptr() );
145  parlist->sublist("Status Test").set("Gradient Tolerance",1.e-14);
146  parlist->sublist("Status Test").set("Constraint Tolerance",1.e-14);
147  parlist->sublist("Status Test").set("Step Tolerance",1.e-16);
148  parlist->sublist("Status Test").set("Iteration Limit",100);
149  // Build Algorithm pointer.
150  Teuchos::RCP<ROL::Algorithm<RealT> > algo;
151 
152  // Solve using trust regions.
153  algo = Teuchos::rcp(new ROL::Algorithm<RealT>("Trust Region",*parlist,false));
154  z.zero();
155  std::clock_t timer_tr = std::clock();
156  algo->run(z,robj,true,*outStream);
157  *outStream << "Trust-Region Newton required " << (std::clock()-timer_tr)/(RealT)CLOCKS_PER_SEC
158  << " seconds.\n";
159  Teuchos::RCP<ROL::Vector<RealT> > zTR = z.clone();
160  zTR->set(z);
161 
162  // Solve using a composite step method.
163  algo = Teuchos::rcp(new ROL::Algorithm<RealT>("Composite Step",*parlist,false));
164  x.zero();
165  std::clock_t timer_cs = std::clock();
166  algo->run(x,g,l,c,obj,con,true,*outStream);
167  *outStream << "Composite Step required " << (std::clock()-timer_cs)/(RealT)CLOCKS_PER_SEC
168  << " seconds.\n";
169 
170  // Compute error between solutions
171  Teuchos::RCP<ROL::Vector<RealT> > err = z.clone();
172  err->set(*zTR); err->axpy(-1.,z);
173  errorFlag += (err->norm() > 1.e-8) ? 1 : 0;
174 
175 // std::ofstream control;
176 // control.open("control.txt");
177 // for (int t = 0; t < nt+1; t++) {
178 // for (int n = 0; n < nx+2; n++) {
179 // control << (RealT)t/(RealT)nt << " "
180 // << (RealT)n/((RealT)(nx+1)) << " "
181 // << (*z_rcp)[t*(nx+2)+n] << "\n";
182 // }
183 // }
184 // control.close();
185 //
186 // std::ofstream state;
187 // state.open("state.txt");
188 // for (int t = 0; t < nt; t++) {
189 // for (int n = 0; n < nx; n++) {
190 // state << (RealT)(t+1)/(RealT)nt << " "
191 // << (RealT)(n+1)/((RealT)(nx+1)) << " "
192 // << (*u_rcp)[t*nx+n] << "\n";
193 // }
194 // }
195 // state.close();
196  }
197  catch (std::logic_error err) {
198  *outStream << err.what() << "\n";
199  errorFlag = -1000;
200  }; // end try
201 
202  if (errorFlag != 0)
203  std::cout << "End Result: TEST FAILED\n";
204  else
205  std::cout << "End Result: TEST PASSED\n";
206 
207  return 0;
208 
209 }
210 
virtual Real checkInverseJacobian_1(const Vector< Real > &jv, const Vector< Real > &v, const Vector< Real > &u, const Vector< Real > &z, const bool printToStream=true, std::ostream &outStream=std::cout)
Defines the linear algebra or vector space interface for simulation-based optimization.
virtual std::vector< std::vector< Real > > checkApplyAdjointJacobian(const Vector< Real > &x, const Vector< Real > &v, const Vector< Real > &c, const Vector< Real > &ajv, const bool printToStream=true, std::ostream &outStream=std::cout, const int numSteps=ROL_NUM_CHECKDERIV_STEPS)
Finite-difference check for the application of the adjoint of constraint Jacobian.
virtual void zero()
Set to zero vector.
Definition: ROL_Vector.hpp:157
virtual std::vector< std::vector< Real > > checkGradient(const Vector< Real > &x, const Vector< Real > &d, const bool printToStream=true, std::ostream &outStream=std::cout, const int numSteps=ROL_NUM_CHECKDERIV_STEPS, const int order=1)
Finite-difference gradient check.
virtual std::vector< std::vector< Real > > checkApplyJacobian(const Vector< Real > &x, const Vector< Real > &v, const Vector< Real > &jv, const std::vector< Real > &steps, const bool printToStream=true, std::ostream &outStream=std::cout, const int order=1)
Finite-difference check for the constraint Jacobian application.
Provides the std::vector implementation of the ROL::Vector interface.
virtual Real checkSolve(const ROL::Vector< Real > &u, const ROL::Vector< Real > &z, const ROL::Vector< Real > &c, const bool printToStream=true, std::ostream &outStream=std::cout)
Provides an interface to run optimization algorithms.
int main(int argc, char *argv[])
Definition: example_03.cpp:53
virtual Teuchos::RCP< Vector< Real > > clone() const
Clone to make a new (uninitialized) vector.
virtual std::vector< std::vector< Real > > checkHessVec(const Vector< Real > &x, const Vector< Real > &v, const bool printToStream=true, std::ostream &outStream=std::cout, const int numSteps=ROL_NUM_CHECKDERIV_STEPS, const int order=1)
Finite-difference Hessian-applied-to-vector check.
virtual Real checkAdjointConsistencyJacobian_1(const Vector< Real > &w, const Vector< Real > &v, const Vector< Real > &u, const Vector< Real > &z, const bool printToStream=true, std::ostream &outStream=std::cout)
Check the consistency of the Jacobian and its adjoint. This is the primary interface.
double RealT
virtual Real checkInverseAdjointJacobian_1(const Vector< Real > &jv, const Vector< Real > &v, const Vector< Real > &u, const Vector< Real > &z, const bool printToStream=true, std::ostream &outStream=std::cout)
virtual std::vector< std::vector< Real > > checkApplyAdjointHessian(const Vector< Real > &x, const Vector< Real > &u, const Vector< Real > &v, const Vector< Real > &hv, const std::vector< Real > &step, const bool printToScreen=true, std::ostream &outStream=std::cout, const int order=1)
Finite-difference check for the application of the adjoint of constraint Hessian. ...
double RealT
Definition: example_03.cpp:51
virtual Real checkAdjointConsistencyJacobian_2(const Vector< Real > &w, const Vector< Real > &v, const Vector< Real > &u, const Vector< Real > &z, const bool printToStream=true, std::ostream &outStream=std::cout)
Check the consistency of the Jacobian and its adjoint. This is the primary interface.