ROL
ROL_Bounds_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#ifndef ROL_BOUNDS_DEF_H
45#define ROL_BOUNDS_DEF_H
46
48
56namespace ROL {
57
58template<typename Real>
59Bounds<Real>::Bounds(const Vector<Real> &x, bool isLower, Real scale, Real feasTol)
60 : scale_(scale), feasTol_(feasTol), mask_(x.clone()), min_diff_(ROL_INF<Real>()) {
61 lower_ = x.clone();
62 upper_ = x.clone();
63 if (isLower) {
64 lower_->set(x);
65 upper_->applyUnary(Elementwise::Fill<Real>( BoundConstraint<Real>::computeInf(x)));
67 }
68 else {
69 lower_->applyUnary(Elementwise::Fill<Real>(-BoundConstraint<Real>::computeInf(x)));
70 upper_->set(x);
72 }
73}
74
75template<typename Real>
76Bounds<Real>::Bounds(const Ptr<Vector<Real>> &x_lo, const Ptr<Vector<Real>> &x_up,
77 const Real scale, const Real feasTol)
78 : scale_(scale), feasTol_(feasTol), mask_(x_lo->clone()) {
79 lower_ = x_lo;
80 upper_ = x_up;
81 const Real half(0.5), one(1);
82 // Compute difference between upper and lower bounds
83 mask_->set(*upper_);
84 mask_->axpy(-one,*lower_);
85 // Compute minimum difference
86 min_diff_ = mask_->reduce(minimum_);
87 min_diff_ *= half;
88}
89
90template<typename Real>
92 struct Lesser : public Elementwise::BinaryFunction<Real> {
93 Real apply(const Real &xc, const Real &yc) const { return xc<yc ? xc : yc; }
94 } lesser;
95
96 struct Greater : public Elementwise::BinaryFunction<Real> {
97 Real apply(const Real &xc, const Real &yc) const { return xc>yc ? xc : yc; }
98 } greater;
99
101 x.applyBinary(lesser, *upper_); // Set x to the elementwise minimum of x and upper_
102 }
104 x.applyBinary(greater,*lower_); // Set x to the elementwise maximum of x and lower_
105 }
106}
107
108template<typename Real>
110 // Make vector strictly feasible
111 // Lower feasibility
113 class LowerFeasible : public Elementwise::BinaryFunction<Real> {
114 private:
115 const Real eps_;
116 const Real diff_;
117 public:
118 LowerFeasible(const Real eps, const Real diff)
119 : eps_(eps), diff_(diff) {}
120 Real apply( const Real &x, const Real &y ) const {
121 const Real tol = static_cast<Real>(100)*ROL_EPSILON<Real>();
122 const Real one(1);
123 Real val = ((y <-tol) ? y*(one-eps_)
124 : ((y > tol) ? y*(one+eps_)
125 : y+eps_));
126 val = std::min(y+eps_*diff_, val);
127 return x < val ? val : x;
128 }
129 };
130 x.applyBinary(LowerFeasible(feasTol_,min_diff_), *lower_);
131 }
132 // Upper feasibility
134 class UpperFeasible : public Elementwise::BinaryFunction<Real> {
135 private:
136 const Real eps_;
137 const Real diff_;
138 public:
139 UpperFeasible(const Real eps, const Real diff)
140 : eps_(eps), diff_(diff) {}
141 Real apply( const Real &x, const Real &y ) const {
142 const Real tol = static_cast<Real>(100)*ROL_EPSILON<Real>();
143 const Real one(1);
144 Real val = ((y <-tol) ? y*(one+eps_)
145 : ((y > tol) ? y*(one-eps_)
146 : y-eps_));
147 val = std::max(y-eps_*diff_, val);
148 return x > val ? val : x;
149 }
150 };
151 x.applyBinary(UpperFeasible(feasTol_,min_diff_), *upper_);
152 }
153}
154
155template<typename Real>
158 Real one(1), epsn(std::min(scale_*eps,static_cast<Real>(0.1)*min_diff_));
159
160 mask_->set(*upper_);
161 mask_->axpy(-one,x);
162
163 Active op(epsn);
164 v.applyBinary(op,*mask_);
165 }
166}
167
168template<typename Real>
169void Bounds<Real>::pruneUpperActive( Vector<Real> &v, const Vector<Real> &g, const Vector<Real> &x, Real xeps, Real geps ) {
171 Real one(1), epsn(std::min(scale_*xeps,static_cast<Real>(0.1)*min_diff_));
172
173 mask_->set(*upper_);
174 mask_->axpy(-one,x);
175
176 UpperBinding op(epsn,geps);
177 mask_->applyBinary(op,g);
178
179 v.applyBinary(prune_,*mask_);
180 }
181}
182
183template<typename Real>
186 Real one(1), epsn(std::min(scale_*eps,static_cast<Real>(0.1)*min_diff_));
187
188 mask_->set(x);
189 mask_->axpy(-one,*lower_);
190
191 Active op(epsn);
192 v.applyBinary(op,*mask_);
193 }
194}
195
196template<typename Real>
197void Bounds<Real>::pruneLowerActive( Vector<Real> &v, const Vector<Real> &g, const Vector<Real> &x, Real xeps, Real geps ) {
199 Real one(1), epsn(std::min(scale_*xeps,static_cast<Real>(0.1)*min_diff_));
200
201 mask_->set(x);
202 mask_->axpy(-one,*lower_);
203
204 LowerBinding op(epsn,geps);
205 mask_->applyBinary(op,g);
206
207 v.applyBinary(prune_,*mask_);
208 }
209}
210
211template<typename Real>
213 const Real one(1);
214 bool flagU = false, flagL = false;
216 mask_->set(*upper_);
217 mask_->axpy(-one,v);
218 Real uminusv = mask_->reduce(minimum_);
219 flagU = ((uminusv<0) ? true : false);
220 }
222 mask_->set(v);
223 mask_->axpy(-one,*lower_);
224 Real vminusl = mask_->reduce(minimum_);
225
226 flagL = ((vminusl<0) ? true : false);
227 }
228 return ((flagU || flagL) ? false : true);
229}
230
231template<typename Real>
233 // TODO: Cache values?
234
235 const Real zero(0), one(1);
236
237 // This implementation handles -l and/or u = infinity properly.
238
239 /*
240 The first if statement defining d_II (Equation 4.4 of [Ulbrich et al. 1999])
241 is equivalent to x - l <= min(u - x, -g) = - max(g, x - u).
242 * Our x, l, u represent u, a, b in the paper.
243 * Indeed, if x - l = -g < u - x, then min{|g|,c} = min{x - l, u - x, c}.
244 */
245
246 d.set(x);
247 d.axpy(-one,*upper_);
248 d.applyBinary(Elementwise::Max<Real>(),g);
249 d.plus(x);
250 d.axpy(-one,*lower_); // = x - l + max(g, x - u)
251
252 mask_->set(x);
253 mask_->axpy(-one,*lower_);
254 mask_->applyBinary(Elementwise::Min<Real>(),g);
255 mask_->plus(x);
256 mask_->axpy(-one,*upper_);
257 mask_->scale(-one); // = u - x - min(g, x - l)
258
259 mask_->applyBinary(Elementwise::Min<Real>(),d);
260
261 d.setScalar(-one);
262 d.applyBinary(Active(zero),*mask_);
263 mask_->setScalar(one);
264 d.plus(*mask_);
265 // d[i] = 1 when one of the if conditions in (4.4) are true else 0
266
267 mask_->set(g);
268 mask_->applyUnary(Elementwise::AbsoluteValue<Real>());
269 d.applyBinary(Elementwise::Multiply<Real>(),*mask_);
270 // d[i] = |g[i]| when one of the if conditions in (4.4) are true else 0
271
272 // Compute min(x - l, u - x).
273 // * We use the identity min(p, q) = min(p + r, q + r) - r to handle the case
274 // where -l or u = infinity.
275 mask_->set(x);
276 mask_->axpy(-one,*lower_);
277 mask_->plus(x);
278 mask_->applyBinary(Elementwise::Min<Real>(),*upper_);
279 mask_->axpy(-one,x); // = min(x - l, u - x)
280
281 // When one of the if conditions in (4.4) are true |g|[i] >= (*mask_)[i].
282 // Thus by taking
283 d.applyBinary(Elementwise::Max<Real>(),*mask_);
284 // d_II follows as min(d, c), i.e.,
285 mask_->set(*upper_);
286 mask_->axpy(-one,*lower_);
287 mask_->applyUnary(buildC_);
288 d.applyBinary(Elementwise::Min<Real>(),*mask_);
289}
290
291template<typename Real>
294 dv.applyBinary(Elementwise::DivideAndInvert<Real>(),v);
295}
296
297template<typename Real>
299 const Real one(1), two(2), three(3);
300
301 // This implementation builds Equation (5.1) of [Ulbrich et al. 1999].
302
303 Bounds<Real>::buildScalingFunction(dv, x, g); // dv = d_II
304
305 mask_->set(*upper_);
306 mask_->axpy(-one,*lower_);
307 mask_->applyUnary(buildC_); // = c
308 mask_->axpy(-one,dv); // = c - d_II
309 dv.setScalar(one);
310 dv.applyBinary(Active(0),*mask_); // = \chi
311
312 mask_->setScalar(three);
313 dv.applyBinary(Elementwise::Multiply<Real>(),*mask_);
314 dv.axpy(-one,*mask_);
315 // dv[i] = 0 if \chi[i] = 1 else -3
316
317 mask_->set(g);
318 mask_->applyUnary(Elementwise::Sign<Real>());
319 dv.plus(*mask_); // dv[i] = sgn(g[i]) if \chi[i] = 1 else dv[i] <= -2
320
321 // Set the dv elements that = 0 to sgn(u + l - 2x).
322 mask_->set(*upper_);
323 mask_->plus(*lower_);
324 mask_->axpy(-two,x);
325 mask_->applyUnary(Elementwise::Sign<Real>());
326 dv.applyBinary(setZeroEntry_,*mask_);
327
328 // Set the dv elements that = 0 to 1.
329 mask_->setScalar(one);
330 dv.applyBinary(setZeroEntry_,*mask_);
331
332 // Set the dv elements that are <= -2 to 0.
333 mask_->set(dv);
334 dv.applyBinary(Active(-two),*mask_);
335
336 dv.applyBinary(Elementwise::Multiply<Real>(), g);
337 dv.applyBinary(Elementwise::Multiply<Real>(), v);
338}
339
340} // namespace ROL
341
342#endif
Objective_SerialSimOpt(const Ptr< Obj > &obj, const V &ui) z0 zero)()
Provides the interface to apply upper and lower bound constraints.
Ptr< Vector< Real > > upper_
void activateLower(void)
Turn on lower bound.
Ptr< Vector< Real > > lower_
void activateUpper(void)
Turn on upper bound.
Elementwise::ReductionMin< Real > minimum_
void buildScalingFunction(Vector< Real > &d, const Vector< Real > &x, const Vector< Real > &g) const
void applyInverseScalingFunction(Vector< Real > &dv, const Vector< Real > &v, const Vector< Real > &x, const Vector< Real > &g) const override
Apply inverse scaling function.
void project(Vector< Real > &x) override
Project optimization variables onto the bounds.
void pruneUpperActive(Vector< Real > &v, const Vector< Real > &x, Real eps=Real(0)) override
Set variables to zero if they correspond to the upper -active set.
void pruneLowerActive(Vector< Real > &v, const Vector< Real > &x, Real eps=Real(0)) override
Set variables to zero if they correspond to the lower -active set.
Ptr< Vector< Real > > mask_
void projectInterior(Vector< Real > &x) override
Project optimization variables into the interior of the feasible set.
void applyScalingFunctionJacobian(Vector< Real > &dv, const Vector< Real > &v, const Vector< Real > &x, const Vector< Real > &g) const override
Apply scaling function Jacobian.
Bounds(const Vector< Real > &x, bool isLower=true, Real scale=1, Real feasTol=std::sqrt(ROL_EPSILON< Real >()))
bool isFeasible(const Vector< Real > &v) override
Check if the vector, v, is feasible.
Defines the linear algebra or vector space interface.
virtual void set(const Vector &x)
Set where .
virtual void setScalar(const Real C)
Set where .
virtual void applyBinary(const Elementwise::BinaryFunction< Real > &f, const Vector &x)
virtual void plus(const Vector &x)=0
Compute , where .
virtual ROL::Ptr< Vector > clone() const =0
Clone to make a new (uninitialized) vector.
virtual void axpy(const Real alpha, const Vector &x)
Compute where .
ROL::BlockOperator2Diagonal BlockOperator2 apply(V &Hv, const V &v, Real &tol) const
Real ROL_INF(void)