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"
52 #include "Teuchos_oblackholestream.hpp"
53 #include "Teuchos_GlobalMPISession.hpp"
54 
55 
56 typedef double RealT;
57 
58 int main(int argc, char *argv[]) {
59 
60  using Teuchos::RCP;
61  using Teuchos::rcp;
62 
63  typedef std::vector<RealT> vector;
64  typedef ROL::Vector<RealT> V;
65  typedef ROL::StdVector<RealT> SV;
66 
67  Teuchos::GlobalMPISession mpiSession(&argc, &argv);
68 
69  // This little trick lets us print to std::cout only if a (dummy) command-line argument is provided.
70  int iprint = argc - 1;
71  Teuchos::RCP<std::ostream> outStream;
72  Teuchos::oblackholestream bhs; // outputs nothing
73  if (iprint > 0)
74  outStream = Teuchos::rcp(&std::cout, false);
75  else
76  outStream = Teuchos::rcp(&bhs, false);
77 
78  // Save the format state of the original std::cout.
79  Teuchos::oblackholestream oldFormatState;
80  oldFormatState.copyfmt(std::cout);
81 
82  RealT errtol = std::sqrt(ROL::ROL_THRESHOLD);
83 
84  int errorFlag = 0;
85 
86  // Specify interval on which to generate uniform random numbers.
87  RealT left = 0.1, right = 0.9;
88 
89  // *** Test body.
90 
91  try {
92 
93  int dim = 1;
94  RCP<vector> x_rcp = rcp( new vector(dim,0.0) );
95  RCP<vector> l_rcp = rcp( new vector(dim,0.0) ); // Lower bound
96  RCP<vector> u_rcp = rcp( new vector(dim,1.0) ); // Upper bound
97  RCP<vector> y_rcp = rcp( new vector(dim,0.0) );
98  RCP<vector> v_rcp = rcp( new vector(dim,0.0) );
99  RCP<vector> d_rcp = rcp( new vector(dim,0.0) );
100  RCP<vector> gx_rcp = rcp( new vector(dim,0.0) );
101  RCP<vector> gy_rcp = rcp( new vector(dim,0.0) );
102  RCP<vector> hv_rcp = rcp( new vector(dim,0.0) );
103 
104  SV x( x_rcp);
105  SV y( y_rcp);
106  SV v( v_rcp);
107  SV d( d_rcp);
108  SV gx(gx_rcp);
109  SV gy(gy_rcp);
110  SV hv(hv_rcp);
111 
112  RandomizeVector(x,left,right);
113  RandomizeVector(v,left,right);
114  RandomizeVector(d,left,right);
115 
116  RCP<V> l = rcp( new SV(l_rcp) );
117  RCP<V> u = rcp( new SV(u_rcp) );
118 
120 
122 
123  // Fixed difference step size
124  RealT delta = 1.e-8;
125 
126  y.set(x); // y = x
127  y.axpy(delta,d); // y = x+delta*d
128 
130 
131  // Do step size sweep
132  *outStream << "Test of single logarithmic penalty objective" << std::endl;
133  log_obj.checkGradient(x, d, true, *outStream); *outStream << "\n";
134  log_obj.checkHessVec(x, v, true, *outStream); *outStream << "\n";
135 
136  *outStream << "Test of bound constraint as logarithmic penalty objective" << std::endl;
137  bc_obj.checkGradient(x, d, true, *outStream); *outStream << "\n";
138  bc_obj.checkHessVec(x, v, true, *outStream); *outStream << "\n";
139 
140 
141 
142 
143  RealT tol = 0;
144 
145  // Compute objective at x and y
146  RealT fx = log_obj.value(x,tol);
147  RealT fy = log_obj.value(y,tol);
148 
149  // Compute gradient at x and y
150  log_obj.gradient(gx,x,tol);
151  log_obj.gradient(gy,y,tol);
152 
153  // Compute action of Hessian on v at x
154  log_obj.hessVec(hv,v,x,tol);
155 
156  // FD gradient error
157  RealT graderr = (fy - fx)/delta - gx.dot(d);
158 
159  // FD Hessian error
160  RCP<V> dg = gx.clone();
161  dg->set(gy);
162  dg->axpy(-1.0,gx);
163 
164  RealT hesserr = ( dg->dot(v) )/delta - hv.dot(d);
165 
166 
167 
168 
169 
170 
171  if( std::abs(graderr) > errtol ) {
172  ++errorFlag;
173  }
174 
175  if( std::abs(hesserr) > errtol ) {
176  ++errorFlag;
177  }
178 
179  }
180  catch (std::logic_error err) {
181  *outStream << err.what() << "\n";
182  errorFlag = -1000;
183  }; // end try
184 
185  if (errorFlag != 0)
186  std::cout << "End Result: TEST FAILED\n";
187  else
188  std::cout << "End Result: TEST PASSED\n";
189 
190  return 0;
191 
192 
193 }
194 
void gradient(Vector< Real > &g, const Vector< Real > &x, Real &tol)
Compute gradient.
void hessVec(Vector< Real > &hv, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Apply Hessian approximation to vector.
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].
static const double ROL_THRESHOLD
Tolerance for various equality tests.
Definition: ROL_Types.hpp:122
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:74
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.
Real value(const Vector< Real > &x, Real &tol)
Compute value.
Provides the interface to apply upper and lower bound constraints.
int main(int argc, char *argv[])
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.
Log barrier objective for interior point methods.
double RealT
double RealT
Create a logarithmic penalty objective from upper and lower bound vectors.