ROL
ROL_RiddersProjection_Def.hpp
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
44
45#ifndef ROL_RIDDERSPROJECTION_DEF_H
46#define ROL_RIDDERSPROJECTION_DEF_H
47
48namespace ROL {
49
50template<typename Real>
51RiddersProjection<Real>::RiddersProjection(const Vector<Real> &xprim,
52 const Vector<Real> &xdual,
53 const Ptr<BoundConstraint<Real>> &bnd,
54 const Ptr<Constraint<Real>> &con,
55 const Vector<Real> &mul,
56 const Vector<Real> &res)
57 : PolyhedralProjection<Real>(xprim,xdual,bnd,con,mul,res),
58 DEFAULT_atol_ (std::sqrt(ROL_EPSILON<Real>()*std::sqrt(ROL_EPSILON<Real>()))),
59 DEFAULT_rtol_ (std::sqrt(ROL_EPSILON<Real>())),
60 DEFAULT_ltol_ (ROL_EPSILON<Real>()),
61 DEFAULT_maxit_ (5000),
62 DEFAULT_verbosity_ (0),
63 atol_ (DEFAULT_atol_),
64 rtol_ (DEFAULT_rtol_),
65 ltol_ (DEFAULT_ltol_),
66 maxit_ (DEFAULT_maxit_),
67 verbosity_ (DEFAULT_verbosity_) {
68 dim_ = mul.dimension();
69 ROL_TEST_FOR_EXCEPTION(dim_!=1,std::logic_error,
70 ">>> ROL::RiddersProjection : The range of the linear constraint must be one dimensional!");
71 xnew_ = xprim.clone();
72 mul1_ = static_cast<Real>(0);
73 dlam1_ = static_cast<Real>(2);
74 // con.value(x) = xprim_->dot(x) + b_
75 Real tol(std::sqrt(ROL_EPSILON<Real>()));
76 xprim_->zero();
77 con_->update(*xprim_,UpdateType::Temp);
78 con_->value(*res_,*xprim_,tol);
79 b_ = res_->dot(*res_->basis(0));
80 mul_->setScalar(static_cast<Real>(1));
81 con_->applyAdjointJacobian(*xdual_,*mul_,xprim,tol);
82 xprim_->set(xdual_->dual());
83 cdot_ = xprim_->dot(*xprim_);
84 // Set tolerance
85 //xnew_->zero();
86 //bnd_->project(*xnew_);
87 //Real res0 = std::abs(residual(*xnew_));
88 Real resl = std::abs(residual(*bnd_->getLowerBound()));
89 Real resu = std::abs(residual(*bnd_->getUpperBound()));
90 Real res0 = std::max(resl,resu);
91 if (res0 < atol_) {
92 res0 = static_cast<Real>(1);
93 }
94 ctol_ = std::min(atol_,rtol_*res0);
95}
96
97template<typename Real>
98RiddersProjection<Real>::RiddersProjection(const Vector<Real> &xprim,
99 const Vector<Real> &xdual,
100 const Ptr<BoundConstraint<Real>> &bnd,
101 const Ptr<Constraint<Real>> &con,
102 const Vector<Real> &mul,
103 const Vector<Real> &res,
104 ParameterList &list)
105 : RiddersProjection<Real>(xprim,xdual,bnd,con,mul,res) {
106 atol_ = list.sublist("General").sublist("Polyhedral Projection").get("Absolute Tolerance", DEFAULT_atol_);
107 rtol_ = list.sublist("General").sublist("Polyhedral Projection").get("Relative Tolerance", DEFAULT_rtol_);
108 ltol_ = list.sublist("General").sublist("Polyhedral Projection").get("Multiplier Tolerance", DEFAULT_ltol_);
109 maxit_ = list.sublist("General").sublist("Polyhedral Projection").get("Iteration Limit", DEFAULT_maxit_);
110 verbosity_ = list.sublist("General").get("Output Level", DEFAULT_verbosity_);
111}
112
113template<typename Real>
114void RiddersProjection<Real>::project(Vector<Real> &x, std::ostream &stream) {
115 if (con_ == nullPtr) {
116 bnd_->project(x);
117 }
118 else {
119 mul1_ = -residual(x)/cdot_;
120 //mul1_ = static_cast<Real>(0);
121 dlam1_ = static_cast<Real>(2);
122 //dlam1_ = static_cast<Real>(1)+std::abs(mul1_);
123 project_df(x, mul1_, dlam1_, stream);
124 mul_->setScalar(mul1_);
125 }
126}
127
128template<typename Real>
129Real RiddersProjection<Real>::residual(const Vector<Real> &x) const {
130 return xprim_->dot(x) + b_;
131}
132
133template<typename Real>
134void RiddersProjection<Real>::update_primal(Vector<Real> &y, const Vector<Real> &x, const Real lam) const {
135 y.set(x);
136 y.axpy(lam,*xprim_);
137 bnd_->project(y);
138}
139
140template<typename Real>
141void RiddersProjection<Real>::project_df(Vector<Real> &x, Real &lam, Real &dlam, std::ostream &stream) const {
142 const Real zero(0), one(1), c1(0.1);
143 Real lamLower(0), lamUpper(0), res(0), resLower(0), resUpper(0), s(0);
144 Real rtol = ctol_;
145 int cnt(0);
146 // Compute initial residual
147 update_primal(*xnew_,x,lam);
148 res = residual(*xnew_);
149 if (res == zero) {
150 x.set(*xnew_);
151 return;
152 }
153 std::ios_base::fmtflags streamFlags(stream.flags());
154 if (verbosity_ > 2) {
155 stream << std::scientific << std::setprecision(6);
156 stream << std::endl;
157 stream << " Polyhedral Projection using Ridders' Algorithm" << std::endl;
158 stream << " Bracketing Phase" << std::endl;
159 }
160 // Bracketing phase
161 if ( res < zero ) {
162 lamLower = lam;
163 resLower = res;
164 lam += dlam;
165 update_primal(*xnew_,x,lam);
166 res = residual(*xnew_);
167 if (verbosity_ > 2) {
168 stream << " ";
169 stream << std::setw(6) << std::left << "iter";
170 stream << std::setw(15) << std::left << "lam";
171 stream << std::setw(15) << std::left << "res";
172 stream << std::setw(15) << std::left << "lower lam";
173 stream << std::setw(15) << std::left << "lower res";
174 stream << std::endl;
175 stream << " ";
176 stream << std::setw(6) << std::left << cnt;
177 stream << std::setw(15) << std::left << lam;
178 stream << std::setw(15) << std::left << res;
179 stream << std::setw(15) << std::left << lamLower;
180 stream << std::setw(15) << std::left << resLower;
181 stream << std::endl;
182 }
183 while ( res < zero && std::abs(res) > rtol && cnt < maxit_ ) {
184 s = std::max(resLower/res-one,c1);
185 dlam += dlam/s;
186 lamLower = lam;
187 resLower = res;
188 lam += dlam;
189 update_primal(*xnew_,x,lam);
190 res = residual(*xnew_);
191 cnt++;
192 if (verbosity_ > 2) {
193 stream << " ";
194 stream << std::setw(6) << std::left << cnt;
195 stream << std::setw(15) << std::left << lam;
196 stream << std::setw(15) << std::left << res;
197 stream << std::setw(15) << std::left << lamLower;
198 stream << std::setw(15) << std::left << resLower;
199 stream << std::endl;
200 }
201 }
202 lamUpper = lam;
203 resUpper = res;
204 }
205 else {
206 lamUpper = lam;
207 resUpper = res;
208 lam -= dlam;
209 update_primal(*xnew_,x,lam);
210 res = residual(*xnew_);
211 if (verbosity_ > 2) {
212 stream << " ";
213 stream << std::setw(6) << std::left << "iter";
214 stream << std::setw(15) << std::left << "lam";
215 stream << std::setw(15) << std::left << "res";
216 stream << std::setw(15) << std::left << "upper lam";
217 stream << std::setw(15) << std::left << "upper res";
218 stream << std::endl;
219 stream << " ";
220 stream << std::setw(6) << std::left << cnt;
221 stream << std::setw(15) << std::left << lam;
222 stream << std::setw(15) << std::left << res;
223 stream << std::setw(15) << std::left << lamUpper;
224 stream << std::setw(15) << std::left << resUpper;
225 stream << std::endl;
226 }
227 while ( res > zero && std::abs(res) > rtol && cnt < maxit_ ) {
228 s = std::max(resUpper/res-one,c1);
229 dlam += dlam/s;
230 lamUpper = lam;
231 resUpper = res;
232 lam -= dlam;
233 update_primal(*xnew_,x,lam);
234 res = residual(*xnew_);
235 cnt++;
236 if (verbosity_ > 2) {
237 stream << " ";
238 stream << std::setw(6) << std::left << cnt;
239 stream << std::setw(15) << std::left << lam;
240 stream << std::setw(15) << std::left << res;
241 stream << std::setw(15) << std::left << lamUpper;
242 stream << std::setw(15) << std::left << resUpper;
243 stream << std::endl;
244 }
245 }
246 lamLower = lam;
247 resLower = res;
248 }
249 if (verbosity_ > 2) {
250 stream << " Bracket: ";
251 stream << std::setw(15) << std::left << lamLower;
252 stream << std::setw(15) << std::left << lamUpper;
253 stream << std::endl;
254 }
255
256 // Secant phase
257 //rtol = ctol_*std::max(one,std::min(std::abs(resLower),std::abs(resUpper)));
258 cnt = 0;
259 if (verbosity_ > 2) {
260 stream << std::endl;
261 stream << " Ridders' Phase" << std::endl;
262 stream << " ";
263 stream << std::setw(6) << std::left << "iter";
264 stream << std::setw(15) << std::left << "rtol";
265 stream << std::setw(15) << std::left << "lam";
266 stream << std::setw(15) << std::left << "res";
267 stream << std::setw(15) << std::left << "lam mid";
268 stream << std::setw(15) << std::left << "res mid";
269 stream << std::setw(15) << std::left << "lam low";
270 stream << std::setw(15) << std::left << "res low";
271 stream << std::setw(15) << std::left << "lam up";
272 stream << std::setw(15) << std::left << "res up";
273 stream << std::endl;
274 }
275 const Real half(0.5); //, bsize = std::abs(lamUpper-lamLower);
276 Real lamMid(0), resMid(0);
277 for (cnt = 0; cnt < maxit_; cnt++) {
278 // Exit if residual or bracket length are sufficiently small
279 //if (std::abs(lamUpper-lamLower) < ltol_*bsize) break;
280 if (std::abs(lamUpper-lamLower) < ltol_) break;
281
282 lamMid = half*(lamUpper+lamLower);
283 update_primal(*xnew_,x,lamMid);
284 resMid = residual(*xnew_);
285 if (std::abs(resMid) <= rtol) {
286 res = resMid;
287 lam = lamMid;
288 break;
289 }
290
291 lam = lamMid-(lamMid-lamLower)*resMid/std::sqrt(resMid*resMid-resLower*resUpper);
292 update_primal(*xnew_,x,lam);
293 res = residual(*xnew_);
294 if (std::abs(res) <= rtol) break;
295 else {
296 if (res < -rtol) {
297 if (resMid < -rtol) {
298 resLower = (lam < lamMid ? resMid : res);
299 lamLower = (lam < lamMid ? lamMid : lam);
300 }
301 else {
302 resLower = res;
303 lamLower = lam;
304 resUpper = resMid;
305 lamUpper = lamMid;
306 }
307 }
308 else {
309 if (resMid < -rtol) {
310 resLower = resMid;
311 lamLower = lamMid;
312 resUpper = res;
313 lamUpper = lam;
314 }
315 else {
316 resUpper = (lam < lamMid ? res : resMid);
317 lamUpper = (lam < lamMid ? lam : lamMid);
318 }
319 }
320 }
321
322 if (verbosity_ > 2) {
323 stream << " ";
324 stream << std::setw(6) << std::left << cnt;
325 stream << std::setw(15) << std::left << rtol;
326 stream << std::setw(15) << std::left << lam;
327 stream << std::setw(15) << std::left << res;
328 stream << std::setw(15) << std::left << lamMid;
329 stream << std::setw(15) << std::left << resMid;
330 stream << std::setw(15) << std::left << lamLower;
331 stream << std::setw(15) << std::left << resLower;
332 stream << std::setw(15) << std::left << lamUpper;
333 stream << std::setw(15) << std::left << resUpper;
334 stream << std::endl;
335 }
336 }
337 if (verbosity_ > 2) {
338 if (cnt < maxit_) {
339 stream << " ";
340 stream << std::setw(6) << std::left << cnt;
341 stream << std::setw(15) << std::left << rtol;
342 stream << std::setw(15) << std::left << lam;
343 stream << std::setw(15) << std::left << res;
344 stream << std::setw(15) << std::left << lamMid;
345 stream << std::setw(15) << std::left << resMid;
346 stream << std::setw(15) << std::left << lamLower;
347 stream << std::setw(15) << std::left << resLower;
348 stream << std::setw(15) << std::left << lamUpper;
349 stream << std::setw(15) << std::left << resUpper;
350 stream << std::endl;
351 }
352 stream << std::endl;
353 }
354 // Return projection
355 x.set(*xnew_);
356 if (std::abs(res) > rtol ) {
357 //throw Exception::NotImplemented(">>> ROL::PolyhedralProjection::project : Projection failed!");
358 stream << ">>> ROL::PolyhedralProjection::project : Projection may be inaccurate! rnorm = ";
359 stream << std::abs(res) << " rtol = " << rtol << std::endl;
360 }
361 stream.flags(streamFlags);
362}
363
364} // namespace ROL
365
366#endif
zero
Objective_SerialSimOpt(const Ptr< Obj > &obj, const V &ui) z0_ zero()
Definition: ROL_Objective_SerialSimOpt.hpp:112
ROL::BoundConstraint
Provides the interface to apply upper and lower bound constraints.
Definition: ROL_BoundConstraint.hpp:73
ROL::Constraint
Defines the general constraint operator interface.
Definition: ROL_Constraint.hpp:86
ROL::PolyhedralProjection::con_
const Ptr< Constraint< Real > > con_
Definition: ROL_PolyhedralProjection.hpp:58
ROL::PolyhedralProjection::bnd_
const Ptr< BoundConstraint< Real > > bnd_
Definition: ROL_PolyhedralProjection.hpp:57
ROL::RiddersProjection::update_primal
void update_primal(Vector< Real > &y, const Vector< Real > &x, const Real lam) const
Definition: ROL_RiddersProjection_Def.hpp:134
ROL::RiddersProjection::residual
Real residual(const Vector< Real > &x) const
Definition: ROL_RiddersProjection_Def.hpp:129
ROL::RiddersProjection::project_df
void project_df(Vector< Real > &x, Real &lam, Real &dlam, std::ostream &stream=std::cout) const
Definition: ROL_RiddersProjection_Def.hpp:141
ROL::RiddersProjection::RiddersProjection
RiddersProjection(const Vector< Real > &xprim, const Vector< Real > &xdual, const Ptr< BoundConstraint< Real > > &bnd, const Ptr< Constraint< Real > > &con, const Vector< Real > &mul, const Vector< Real > &res)
Definition: ROL_RiddersProjection_Def.hpp:51
ROL::RiddersProjection::project
void project(Vector< Real > &x, std::ostream &stream=std::cout) override
Definition: ROL_RiddersProjection_Def.hpp:114
ROL::RiddersProjection::xnew_
Ptr< Vector< Real > > xnew_
Definition: ROL_RiddersProjection.hpp:57
ROL::Vector
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:84
ROL::Vector::set
virtual void set(const Vector &x)
Set where .
Definition: ROL_Vector.hpp:209
ROL::Vector::clone
virtual ROL::Ptr< Vector > clone() const =0
Clone to make a new (uninitialized) vector.
ROL::Vector::dimension
virtual int dimension() const
Return dimension of the vector space.
Definition: ROL_Vector.hpp:196
ROL::Vector::axpy
virtual void axpy(const Real alpha, const Vector &x)
Compute where .
Definition: ROL_Vector.hpp:153
ROL::ROL_EPSILON
Real ROL_EPSILON(void)
Platform-dependent machine epsilon.
Definition: ROL_Types.hpp:91
Generated by doxygen 1.9.4