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