ROL
function/test_02.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 
48 #include "ROL_RandomVector.hpp"
49 #include "ROL_StdVector.hpp"
51 
52 #include "Teuchos_oblackholestream.hpp"
53 #include "Teuchos_GlobalMPISession.hpp"
54 #include "Teuchos_ParameterList.hpp"
55 
56 
57 typedef double RealT;
58 
59 int main(int argc, char *argv[]) {
60 
61  using Teuchos::RCP;
62  using Teuchos::rcp;
63 
64  typedef std::vector<RealT> vector;
65  typedef ROL::Vector<RealT> V;
66  typedef ROL::StdVector<RealT> SV;
67 
68  typedef typename vector::size_type uint;
69 
70  Teuchos::GlobalMPISession mpiSession(&argc, &argv);
71 
72  // This little trick lets us print to std::cout only if a (dummy) command-line argument is provided.
73  int iprint = argc - 1;
74  Teuchos::RCP<std::ostream> outStream;
75  Teuchos::oblackholestream bhs; // outputs nothing
76  if (iprint > 0)
77  outStream = Teuchos::rcp(&std::cout, false);
78  else
79  outStream = Teuchos::rcp(&bhs, false);
80 
81  // Save the format state of the original std::cout.
82  Teuchos::oblackholestream oldFormatState;
83  oldFormatState.copyfmt(std::cout);
84 
85 // RealT errtol = std::sqrt(ROL::ROL_THRESHOLD<RealT>());
86 
87  int errorFlag = 0;
88 
89  // *** Test body.
90 
91  try {
92 
93  uint dim = 30;
94  RealT xmin = 0.5;
95  RealT xmax = 2.5;
96 
97  RCP<vector> x_rcp = rcp( new vector(dim,0.0) );
98  RCP<vector> g_rcp = rcp( new vector(dim,0.0) );
99  RCP<vector> v_rcp = rcp( new vector(dim,0.0) );
100  RCP<vector> hv_rcp = rcp( new vector(dim,0.0) );
101 
102  RCP<vector> l_rcp = rcp( new vector(dim,1.0) );
103  RCP<vector> u_rcp = rcp( new vector(dim,2.0) );
104 
105  RCP<vector> xlog0_rcp = rcp( new vector(dim,0.0) );
106  RCP<vector> xlog1_rcp = rcp( new vector(dim,0.0) );
107  RCP<vector> xlog2_rcp = rcp( new vector(dim,0.0) );
108 
109  RCP<vector> xquad0_rcp = rcp( new vector(dim,0.0) );
110  RCP<vector> xquad1_rcp = rcp( new vector(dim,0.0) );
111  RCP<vector> xquad2_rcp = rcp( new vector(dim,0.0) );
112 
113  RCP<vector> xdwell0_rcp = rcp( new vector(dim,0.0) );
114  RCP<vector> xdwell1_rcp = rcp( new vector(dim,0.0) );
115  RCP<vector> xdwell2_rcp = rcp( new vector(dim,0.0) );
116 
117 
118 
119  SV x(x_rcp);
120  SV g(g_rcp);
121  SV v(v_rcp);
122  SV hv(hv_rcp);
123 
124  RCP<SV> xlog0 = rcp( new SV(xlog0_rcp) );
125  RCP<SV> xlog1 = rcp( new SV(xlog1_rcp) );
126  RCP<SV> xlog2 = rcp( new SV(xlog2_rcp) );
127 
128  RCP<SV> xquad0 = rcp( new SV(xquad0_rcp) );
129  RCP<SV> xquad1 = rcp( new SV(xquad1_rcp) );
130  RCP<SV> xquad2 = rcp( new SV(xquad2_rcp) );
131 
132  RCP<SV> xdwell0 = rcp( new SV(xdwell0_rcp) );
133  RCP<SV> xdwell1 = rcp( new SV(xdwell1_rcp) );
134  RCP<SV> xdwell2 = rcp( new SV(xdwell2_rcp) );
135 
136  RCP<V> lo = rcp(new SV(l_rcp) );
137  RCP<V> up = rcp(new SV(u_rcp) );
138 
139  for(uint i=0; i<dim; ++i) {
140  RealT t = static_cast<RealT>(i)/static_cast<RealT>(dim-1);
141  (*x_rcp)[i] = xmin*(1-t) + xmax*t;
142  }
143 
144  // Create bound constraint
145  ROL::BoundConstraint<RealT> bc(lo,up);
146 
147  Teuchos::ParameterList logList;
148  Teuchos::ParameterList quadList;
149  Teuchos::ParameterList dwellList;
150 
151  logList.sublist("Barrier Function").set("Type","Logarithmic");
152  quadList.sublist("Barrier Function").set("Type","Quadratic");
153  dwellList.sublist("Barrier Function").set("Type","Double Well");
154 
155  ROL::ObjectiveFromBoundConstraint<RealT> logObj(bc,logList);
156  ROL::ObjectiveFromBoundConstraint<RealT> quadObj(bc,quadList);
157  ROL::ObjectiveFromBoundConstraint<RealT> dwellObj(bc,dwellList);
158 
159  RealT tol = 0.0;
160 
161 
162  logObj.value(x,tol);
163  xlog0->set(*Teuchos::rcp_static_cast<SV>(logObj.getBarrierVector()));
164 
165  logObj.gradient(g,x,tol);
166  xlog1->set(*Teuchos::rcp_static_cast<SV>(logObj.getBarrierVector()));
167 
168  logObj.hessVec(hv,v,x,tol);
169  xlog2->set(*Teuchos::rcp_static_cast<SV>(logObj.getBarrierVector()));
170 
171 
172  quadObj.value(x,tol);
173  xquad0->set(*Teuchos::rcp_static_cast<SV>(quadObj.getBarrierVector()));
174 
175  quadObj.gradient(g,x,tol);
176  xquad1->set(*Teuchos::rcp_static_cast<SV>(quadObj.getBarrierVector()));
177 
178  quadObj.hessVec(hv,v,x,tol);
179  xquad2->set(*Teuchos::rcp_static_cast<SV>(quadObj.getBarrierVector()));
180 
181 
182  dwellObj.value(x,tol);
183  xdwell0->set(*Teuchos::rcp_static_cast<SV>(dwellObj.getBarrierVector()));
184 
185  dwellObj.gradient(g,x,tol);
186  xdwell1->set(*Teuchos::rcp_static_cast<SV>(dwellObj.getBarrierVector()));
187 
188  dwellObj.hessVec(hv,v,x,tol);
189  xdwell2->set(*Teuchos::rcp_static_cast<SV>(dwellObj.getBarrierVector()));
190 
191 
192  *outStream << std::setw(14) << "x"
193  << std::setw(14) << "log"
194  << std::setw(14) << "D(log)"
195  << std::setw(14) << "D2(log)"
196  << std::setw(14) << "quad"
197  << std::setw(14) << "D(quad)"
198  << std::setw(14) << "D2(quad)"
199  << std::setw(14) << "dwell"
200  << std::setw(14) << "D(dwell)"
201  << std::setw(14) << "D2(dwell)"
202  << std::endl;
203  *outStream << std::string(140,'-') << std::endl;
204 
205  for(uint i=0; i<dim; ++i) {
206  *outStream << std::setw(14) << (*x_rcp)[i]
207  << std::setw(14) << (*xlog0_rcp)[i]
208  << std::setw(14) << (*xlog1_rcp)[i]
209  << std::setw(14) << (*xlog2_rcp)[i]
210  << std::setw(14) << (*xquad0_rcp)[i]
211  << std::setw(14) << (*xquad1_rcp)[i]
212  << std::setw(14) << (*xquad2_rcp)[i]
213  << std::setw(14) << (*xdwell0_rcp)[i]
214  << std::setw(14) << (*xdwell1_rcp)[i]
215  << std::setw(14) << (*xdwell2_rcp)[i]
216  << std::endl;
217  }
218 
219 
220  ROL::RandomizeVector( x, 1.2, 1.8 );
221  ROL::RandomizeVector( v, -0.1, 0.1 );
222 
223  *outStream << "\n\n";
224  *outStream << "Test of logarithmic penalty objective" << std::endl;
225  logObj.checkGradient(x,v,true,*outStream); *outStream << std::endl;
226  logObj.checkHessVec(x,v,true,*outStream); *outStream << std::endl;
227 
228  ROL::RandomizeVector( x, -1.0, 1.0 );
229  ROL::RandomizeVector( v, -1.0, 1.0 );
230 
231  *outStream << "\n\n";
232  *outStream << "Test of piecewise quadratic penalty objective" << std::endl;
233  quadObj.checkGradient(x,v,true,*outStream); *outStream << std::endl;
234  quadObj.checkHessVec(x,v,true,*outStream); *outStream << std::endl;
235 
236 
237  *outStream << "\n\n";
238  *outStream << "Test of double well penalty objective" << std::endl;
239  dwellObj.checkGradient(x,v,true,*outStream); *outStream << std::endl;
240  dwellObj.checkHessVec(x,v,true,*outStream); *outStream << std::endl;
241 
242 
243 
244 
245 
246  }
247  catch (std::logic_error err) {
248  *outStream << err.what() << "\n";
249  errorFlag = -1000;
250  }; // end try
251 
252  if (errorFlag != 0)
253  std::cout << "End Result: TEST FAILED\n";
254  else
255  std::cout << "End Result: TEST PASSED\n";
256 
257  return 0;
258 
259 
260 }
261 
void RandomizeVector(Vector< Real > &x, const Real &lower=0.0, const Real &upper=1.0)
Fill a ROL::Vector with uniformly-distributed random numbers in the interval [lower,upper].
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:74
Teuchos::RCP< Vector< Real > > getBarrierVector(void)
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.
Provides the std::vector implementation of the ROL::Vector interface.
Provides the interface to apply upper and lower bound constraints.
int main(int argc, char *argv[])
Real value(const Vector< Real > &x, Real &tol)
Compute value.
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.
void gradient(Vector< Real > &g, const Vector< Real > &x, Real &tol)
Compute gradient.
double RealT
double RealT
void hessVec(Vector< Real > &hv, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Apply Hessian approximation to vector.
Create a penalty objective from upper and lower bound vectors.