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 
44 /*! \file  test_05.cpp
45     \brief Shows how to use the nonlinear least squares interface
46            to find a feasible point for the equality constrained NLP
47            from Nocedal/Wright, 2nd edition, page 574, example 18.2.
48 */
49 
50 #include "ROL_SimpleEqConstrained.hpp"
51 #include "ROL_StdVector.hpp"
52 #include "ROL_NonlinearLeastSquaresObjective.hpp"
53 #include "ROL_Algorithm.hpp"
54 #include "ROL_TrustRegionStep.hpp"
55 #include "ROL_Stream.hpp"
56 #include "Teuchos_GlobalMPISession.hpp"
57 
58 #include <iostream>
59 
60 typedef double RealT;
61 
62 
main(int argc,char * argv[])63 int main(int argc, char *argv[]) {
64 
65   Teuchos::GlobalMPISession mpiSession(&argc, &argv);
66 
67   // This little trick lets us print to std::cout only if a (dummy) command-line argument is provided.
68   int iprint     = argc - 1;
69   ROL::Ptr<std::ostream> outStream;
70   ROL::nullstream bhs; // outputs nothing
71   if (iprint > 0)
72     outStream = ROL::makePtrFromRef(std::cout);
73   else
74     outStream = ROL::makePtrFromRef(bhs);
75 
76   int errorFlag  = 0;
77 
78   // *** Example body.
79 
80   try {
81 
82     ROL::Ptr<ROL::Objective<RealT>> obj;
83     ROL::Ptr<ROL::Constraint<RealT>> constr;
84     ROL::Ptr<ROL::Vector<RealT>> x;
85     ROL::Ptr<ROL::Vector<RealT>> sol;
86 
87     // Retrieve objective, constraint, iteration vector, solution vector.
88     ROL::ZOO::getSimpleEqConstrained<RealT> SEC;
89     obj    = SEC.getObjective();
90     constr = SEC.getEqualityConstraint();
91     x      = SEC.getInitialGuess();
92     sol    = SEC.getSolution();
93 
94     // Inititalize vectors
95     int dim = 5;
96     int nc = 3;
97     RealT left = -1e0, right = 1e0;
98     ROL::Ptr<std::vector<RealT>> xtest_ptr = ROL::makePtr<std::vector<RealT>>(dim, 0.0);
99     ROL::Ptr<std::vector<RealT>> g_ptr = ROL::makePtr<std::vector<RealT>>(dim, 0.0);
100     ROL::Ptr<std::vector<RealT>> d_ptr = ROL::makePtr<std::vector<RealT>>(dim, 0.0);
101     ROL::Ptr<std::vector<RealT>> v_ptr = ROL::makePtr<std::vector<RealT>>(dim, 0.0);
102     ROL::Ptr<std::vector<RealT>> vc_ptr = ROL::makePtr<std::vector<RealT>>(nc, 0.0);
103     ROL::Ptr<std::vector<RealT>> vl_ptr = ROL::makePtr<std::vector<RealT>>(nc, 0.0);
104     ROL::StdVector<RealT> xtest(xtest_ptr);
105     ROL::StdVector<RealT> g(g_ptr);
106     ROL::StdVector<RealT> d(d_ptr);
107     ROL::StdVector<RealT> v(v_ptr);
108     ROL::StdVector<RealT> vc(vc_ptr);
109     ROL::StdVector<RealT> vl(vl_ptr);
110     // set xtest, d, v
111     for (int i=0; i<dim; i++) {
112       (*xtest_ptr)[i] = ( (RealT)rand() / (RealT)RAND_MAX ) * (right - left) + left;
113       (*d_ptr)[i] = ( (RealT)rand() / (RealT)RAND_MAX ) * (right - left) + left;
114       (*v_ptr)[i] = ( (RealT)rand() / (RealT)RAND_MAX ) * (right - left) + left;
115     }
116     // set vc, vl
117     for (int i=0; i<nc; i++) {
118       (*vc_ptr)[i] = ( (RealT)rand() / (RealT)RAND_MAX ) * (right - left) + left;
119       (*vl_ptr)[i] = ( (RealT)rand() / (RealT)RAND_MAX ) * (right - left) + left;
120     }
121 
122     xtest.set(*x);
123 
124     // Initialize nonlinear least squares objectives
125     ROL::NonlinearLeastSquaresObjective<RealT> nlls(constr,*x,vc,false);
126     ROL::NonlinearLeastSquaresObjective<RealT> gnnlls(constr,*x,vc,true);
127 
128     // Check derivatives
129     constr->checkApplyJacobian(xtest, v, vc, true, *outStream);                 *outStream << "\n";
130     constr->checkApplyAdjointJacobian(xtest, vl, vc, xtest, true, *outStream);  *outStream << "\n";
131     constr->checkApplyAdjointHessian(xtest, vl, d, xtest, true, *outStream);    *outStream << "\n";
132     nlls.checkGradient(xtest, d, true, *outStream);                             *outStream << "\n";
133     nlls.checkHessVec(xtest, v, true, *outStream);                              *outStream << "\n";
134     nlls.checkHessSym(xtest, d, v, true, *outStream);                           *outStream << "\n";
135 
136     // Define algorithm.
137     ROL::ParameterList parlist;
138     parlist.sublist("Step").sublist("Trust Region").set("Subproblem Solver","Truncated CG");
139     parlist.sublist("Status Test").set("Gradient Tolerance",1.e-10);
140     parlist.sublist("Status Test").set("Constraint Tolerance",1.e-10);
141     parlist.sublist("Status Test").set("Step Tolerance",1.e-18);
142     parlist.sublist("Status Test").set("Iteration Limit",100);
143     ROL::Ptr<ROL::StatusTest<RealT>> status = ROL::makePtr<ROL::StatusTest<RealT>>(parlist);
144     ROL::Ptr<ROL::Step<RealT>> step = ROL::makePtr<ROL::TrustRegionStep<RealT>>(parlist);
145     ROL::Algorithm<RealT> algo(step,status,false);
146 
147     // Run Algorithm
148     *outStream << "\nSOLVE USING FULL HESSIAN\n";
149     x->set(xtest);
150     algo.run(*x, nlls, true, *outStream);
151     algo.reset();
152     *outStream << "\nSOLVE USING GAUSS-NEWTON HESSIAN\n";
153     x->set(xtest);
154     algo.run(*x, gnnlls, true, *outStream);
155   }
156   catch (std::logic_error& err) {
157     *outStream << err.what() << "\n";
158     errorFlag = -1000;
159   }; // end try
160 
161   if (errorFlag != 0)
162     std::cout << "End Result: TEST FAILED\n";
163   else
164     std::cout << "End Result: TEST PASSED\n";
165 
166   return 0;
167 
168 }
169 
170