ROL
ROL_DaiFletcherProjection_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_DAIFLETCHERPROJECTION_DEF_H
11#define ROL_DAIFLETCHERPROJECTION_DEF_H
12
13namespace ROL {
14
15template<typename Real>
16DaiFletcherProjection<Real>::DaiFletcherProjection(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>
37DaiFletcherProjection<Real>::DaiFletcherProjection(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 DaiFletcherProjection<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::DaiFletcherProjection : 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 DaiFletcherProjection<Real>::project(Vector<Real> &x, std::ostream &stream) {
101 if (con_ == nullPtr) {
102 bnd_->project(x);
103 }
104 else {
105 Px_->set(x); bnd_->project(*Px_);
106 mul1_ = -residual(*Px_)/cdot_;
107 //mul1_ = -residual(x)/cdot_;
108 //mul1_ = static_cast<Real>(0);
109 dlam1_ = static_cast<Real>(2);
110 //dlam1_ = static_cast<Real>(1)+std::abs(mul1_);
111 project_df(x, mul1_, dlam1_, stream);
112 mul_->setScalar(mul1_);
113 }
114}
115
116template<typename Real>
117Real DaiFletcherProjection<Real>::residual(const Vector<Real> &x) const {
118 return xprim_->dot(x) + b_;
119}
120
121template<typename Real>
122void DaiFletcherProjection<Real>::update_primal(Vector<Real> &y, const Vector<Real> &x, const Real lam) const {
123 y.set(x);
124 y.axpy(lam,*xprim_);
125 bnd_->project(y);
126}
127
128template<typename Real>
129void DaiFletcherProjection<Real>::project_df(Vector<Real> &x, Real &lam, Real &dlam, std::ostream &stream) const {
130 const Real zero(0), one(1), two(2), c1(0.1), c2(0.75), c3(0.25);
131 Real lamLower(0), lamUpper(0), lamNew(0), res(0), resLower(0), resUpper(0), s(0);
132 Real rtol = ctol_;
133 int cnt(0);
134 // Compute initial residual
135 update_primal(*xnew_,x,lam);
136 res = residual(*xnew_);
137 if (res == zero) {
138 x.set(*xnew_);
139 return;
140 }
141 std::ios_base::fmtflags streamFlags(stream.flags());
142 if (verbosity_ > 2) {
143 stream << std::scientific << std::setprecision(6);
144 stream << std::endl;
145 stream << " Polyhedral Projection using the Dai-Fletcher Algorithm" << std::endl;
146 stream << " Bracketing Phase" << std::endl;
147 }
148 // Bracketing phase
149 if ( res < zero ) {
150 lamLower = lam;
151 resLower = res;
152 lam += dlam;
153 update_primal(*xnew_,x,lam);
154 res = residual(*xnew_);
155 if (verbosity_ > 2) {
156 stream << " ";
157 stream << std::setw(6) << std::left << "iter";
158 stream << std::setw(15) << std::left << "lam";
159 stream << std::setw(15) << std::left << "res";
160 stream << std::setw(15) << std::left << "lower lam";
161 stream << std::setw(15) << std::left << "lower res";
162 stream << std::endl;
163 stream << " ";
164 stream << std::setw(6) << std::left << cnt;
165 stream << std::setw(15) << std::left << lam;
166 stream << std::setw(15) << std::left << res;
167 stream << std::setw(15) << std::left << lamLower;
168 stream << std::setw(15) << std::left << resLower;
169 stream << std::endl;
170 }
171 while ( res < zero && std::abs(res) > rtol && cnt < maxit_ ) {
172 s = std::max(resLower/res-one,c1);
173 dlam += dlam/s;
174 lamLower = lam;
175 resLower = res;
176 lam += dlam;
177 update_primal(*xnew_,x,lam);
178 res = residual(*xnew_);
179 cnt++;
180 if (verbosity_ > 2) {
181 stream << " ";
182 stream << std::setw(6) << std::left << cnt;
183 stream << std::setw(15) << std::left << lam;
184 stream << std::setw(15) << std::left << res;
185 stream << std::setw(15) << std::left << lamLower;
186 stream << std::setw(15) << std::left << resLower;
187 stream << std::endl;
188 }
189 }
190 lamUpper = lam;
191 resUpper = res;
192 }
193 else {
194 lamUpper = lam;
195 resUpper = res;
196 lam -= dlam;
197 update_primal(*xnew_,x,lam);
198 res = residual(*xnew_);
199 if (verbosity_ > 2) {
200 stream << " ";
201 stream << std::setw(6) << std::left << "iter";
202 stream << std::setw(15) << std::left << "lam";
203 stream << std::setw(15) << std::left << "res";
204 stream << std::setw(15) << std::left << "upper lam";
205 stream << std::setw(15) << std::left << "upper res";
206 stream << std::endl;
207 stream << " ";
208 stream << std::setw(6) << std::left << cnt;
209 stream << std::setw(15) << std::left << lam;
210 stream << std::setw(15) << std::left << res;
211 stream << std::setw(15) << std::left << lamUpper;
212 stream << std::setw(15) << std::left << resUpper;
213 stream << std::endl;
214 }
215 while ( res > zero && std::abs(res) > rtol && cnt < maxit_ ) {
216 s = std::max(resUpper/res-one,c1);
217 dlam += dlam/s;
218 lamUpper = lam;
219 resUpper = res;
220 lam -= dlam;
221 update_primal(*xnew_,x,lam);
222 res = residual(*xnew_);
223 cnt++;
224 if (verbosity_ > 2) {
225 stream << " ";
226 stream << std::setw(6) << std::left << cnt;
227 stream << std::setw(15) << std::left << lam;
228 stream << std::setw(15) << std::left << res;
229 stream << std::setw(15) << std::left << lamUpper;
230 stream << std::setw(15) << std::left << resUpper;
231 stream << std::endl;
232 }
233 }
234 lamLower = lam;
235 resLower = res;
236 }
237 if (verbosity_ > 2) {
238 stream << " Bracket: ";
239 stream << std::setw(15) << std::left << lamLower;
240 stream << std::setw(15) << std::left << lamUpper;
241 stream << std::endl;
242 }
243
244 // Secant phase
245 rtol = ctol_*std::max(one,std::min(std::abs(resLower),std::abs(resUpper)));
246 //s = one - resLower / resUpper;
247 //dlam = (lamUpper - lamLower) / s;
248 //lam = lamUpper - dlam;
249 s = (resUpper - resLower) / resUpper;
250 lam = (resUpper * lamLower - resLower * lamUpper) / (resUpper - resLower);
251 dlam = lamUpper - lam;
252 update_primal(*xnew_,x,lam);
253 res = residual(*xnew_);
254 cnt = 0;
255 if (verbosity_ > 2) {
256 stream << std::endl;
257 stream << " Secant Phase" << std::endl;
258 stream << " ";
259 stream << std::setw(6) << std::left << "iter";
260 stream << std::setw(15) << std::left << "lam";
261 stream << std::setw(15) << std::left << "res";
262 stream << std::setw(15) << std::left << "stepsize";
263 stream << std::setw(15) << std::left << "rtol";
264 stream << std::setw(15) << std::left << "lbnd";
265 stream << std::setw(15) << std::left << "lres";
266 stream << std::setw(15) << std::left << "ubnd";
267 stream << std::setw(15) << std::left << "ures";
268 stream << std::endl;
269 stream << " ";
270 stream << std::setw(6) << std::left << cnt;
271 stream << std::setw(15) << std::left << lam;
272 stream << std::setw(15) << std::left << res;
273 stream << std::setw(15) << std::left << dlam;
274 stream << std::setw(15) << std::left << rtol;
275 stream << std::setw(15) << std::left << lamLower;
276 stream << std::setw(15) << std::left << resLower;
277 stream << std::setw(15) << std::left << lamUpper;
278 stream << std::setw(15) << std::left << resUpper;
279 stream << std::endl;
280 }
281 for (cnt = 1; cnt < maxit_; cnt++) {
282 // Exit if residual or bracket length are sufficiently small
283 if ( std::abs(res) <= rtol ||
284 std::abs(lamUpper-lamLower) < ltol_*std::max(std::abs(lamUpper),std::abs(lamLower)) ) {
285 break;
286 }
287
288 if ( res > zero ) {
289 if ( s <= two ) {
290 lamUpper = lam;
291 resUpper = res;
292 //s = one - resLower / resUpper;
293 //dlam = (lamUpper - lamLower) / s;
294 //lam = lamUpper - dlam;
295 s = (resUpper - resLower) / resUpper;
296 lam = (lamLower * resUpper - lamUpper * resLower) / (resUpper - resLower);
297 dlam = lamUpper - lam;
298 }
299 else {
300 //s = std::max(resUpper / res - one, c1);
301 //dlam = (lamUpper - lam) / s;
302 //lamNew = std::max(lam - dlam, c2*lamLower + c3*lam);
303 if (resUpper <= (c1+one)*res) {
304 dlam = (lamUpper - lam) / c1;
305 lamNew = std::max(lam - dlam, c2*lamLower + c3*lam);
306 }
307 else {
308 lamNew = std::max((lam * resUpper - lamUpper * res) / (resUpper - res),
309 c2*lamLower + c3*lam);
310 dlam = lam - lamNew;
311 }
312 lamUpper = lam;
313 resUpper = res;
314 lam = lamNew;
315 s = (lamUpper - lamLower) / (lamUpper - lam);
316 }
317 }
318 else {
319 if ( s >= two ) {
320 lamLower = lam;
321 resLower = res;
322 //s = one - resLower / resUpper;
323 //dlam = (lamUpper - lamLower) / s;
324 //lam = lamUpper - dlam;
325 s = (resUpper - resLower) / resUpper;
326 lam = (lamLower * resUpper - lamUpper * resLower) / (resUpper - resLower);
327 dlam = lamUpper - lam;
328 }
329 else {
330 //s = std::max(resLower / res - one, c1);
331 //dlam = (lam + lamLower) / s;
332 //lamNew = std::min(lam + dlam, c2*lamUpper + c3*lam);
333 if (resLower >= (c1+one)*res) {
334 dlam = (lam - lamLower) / c1;
335 lamNew = std::max(lam + dlam, c2*lamUpper + c3*lam);
336 }
337 else {
338 lamNew = std::max((lamLower * res - lam * resLower) / (res - resLower),
339 c2*lamUpper + c3*lam);
340 dlam = lamNew - lamLower;
341 }
342 lamLower = lam;
343 resLower = res;
344 lam = lamNew;
345 s = (lamUpper - lamLower) / (lamUpper - lam);
346 }
347 }
348 update_primal(*xnew_,x,lam);
349 res = residual(*xnew_);
350
351 if (verbosity_ > 2) {
352 stream << " ";
353 stream << std::setw(6) << std::left << cnt;
354 stream << std::setw(15) << std::left << lam;
355 stream << std::setw(15) << std::left << res;
356 stream << std::setw(15) << std::left << dlam;
357 stream << std::setw(15) << std::left << rtol;
358 stream << std::setw(15) << std::left << lamLower;
359 stream << std::setw(15) << std::left << resLower;
360 stream << std::setw(15) << std::left << lamUpper;
361 stream << std::setw(15) << std::left << resUpper;
362 stream << std::endl;
363 }
364 }
365 if (verbosity_ > 2) {
366 stream << std::endl;
367 }
368 // Return projection
369 x.set(*xnew_);
370 if (std::abs(res) > rtol ) {
371 //throw Exception::NotImplemented(">>> ROL::PolyhedralProjection::project : Projection failed!");
372 stream << ">>> ROL::PolyhedralProjection::project : Projection may be inaccurate! rnorm = ";
373 stream << std::abs(res) << " rtol = " << rtol << std::endl;
374 }
375 stream.flags(streamFlags);
376}
377
378} // namespace ROL
379
380#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::DaiFletcherProjection::project_df
void project_df(Vector< Real > &x, Real &lam, Real &dlam, std::ostream &stream=std::cout) const
Definition ROL_DaiFletcherProjection_Def.hpp:129
ROL::DaiFletcherProjection::project
void project(Vector< Real > &x, std::ostream &stream=std::cout) override
Definition ROL_DaiFletcherProjection_Def.hpp:100
ROL::DaiFletcherProjection::DaiFletcherProjection
DaiFletcherProjection(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_DaiFletcherProjection_Def.hpp:16
ROL::DaiFletcherProjection::residual
Real residual(const Vector< Real > &x) const
Definition ROL_DaiFletcherProjection_Def.hpp:117
ROL::DaiFletcherProjection::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_DaiFletcherProjection_Def.hpp:64
ROL::DaiFletcherProjection::update_primal
void update_primal(Vector< Real > &y, const Vector< Real > &x, const Real lam) const
Definition ROL_DaiFletcherProjection_Def.hpp:122
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