ROL
ROL_Bounds_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_BOUNDS_DEF_H
11 #define ROL_BOUNDS_DEF_H
12 
13 #include "ROL_BoundConstraint.hpp"
14 
22 namespace ROL {
23 
24 template<typename Real>
25 Bounds<Real>::Bounds(const Vector<Real> &x, bool isLower, Real scale, Real feasTol)
26  : scale_(scale), feasTol_(feasTol), mask_(x.clone()), min_diff_(ROL_INF<Real>()) {
27  lower_ = x.clone();
28  upper_ = x.clone();
29  if (isLower) {
30  lower_->set(x);
31  upper_->applyUnary(Elementwise::Fill<Real>( BoundConstraint<Real>::computeInf(x)));
33  }
34  else {
35  lower_->applyUnary(Elementwise::Fill<Real>(-BoundConstraint<Real>::computeInf(x)));
36  upper_->set(x);
38  }
39 }
40 
41 template<typename Real>
42 Bounds<Real>::Bounds(const Ptr<Vector<Real>> &x_lo, const Ptr<Vector<Real>> &x_up,
43  const Real scale, const Real feasTol)
44  : scale_(scale), feasTol_(feasTol), mask_(x_lo->clone()) {
45  lower_ = x_lo;
46  upper_ = x_up;
47  const Real half(0.5), one(1);
48  // Compute difference between upper and lower bounds
49  mask_->set(*upper_);
50  mask_->axpy(-one,*lower_);
51  // Compute minimum difference
52  min_diff_ = mask_->reduce(minimum_);
53  min_diff_ *= half;
54 }
55 
56 template<typename Real>
58  struct Lesser : public Elementwise::BinaryFunction<Real> {
59  Real apply(const Real &xc, const Real &yc) const { return xc<yc ? xc : yc; }
60  } lesser;
61 
62  struct Greater : public Elementwise::BinaryFunction<Real> {
63  Real apply(const Real &xc, const Real &yc) const { return xc>yc ? xc : yc; }
64  } greater;
65 
67  x.applyBinary(lesser, *upper_); // Set x to the elementwise minimum of x and upper_
68  }
70  x.applyBinary(greater,*lower_); // Set x to the elementwise maximum of x and lower_
71  }
72 }
73 
74 template<typename Real>
76  // Make vector strictly feasible
77  // Lower feasibility
79  class LowerFeasible : public Elementwise::BinaryFunction<Real> {
80  private:
81  const Real eps_;
82  const Real diff_;
83  public:
84  LowerFeasible(const Real eps, const Real diff)
85  : eps_(eps), diff_(diff) {}
86  Real apply( const Real &xc, const Real &yc ) const {
87  const Real tol = static_cast<Real>(100)*ROL_EPSILON<Real>();
88  const Real one(1);
89  Real val = ((yc <-tol) ? yc*(one-eps_)
90  : ((yc > tol) ? yc*(one+eps_)
91  : yc+eps_));
92  val = std::min(yc+eps_*diff_, val);
93  return xc < val ? val : xc;
94  }
95  };
96  x.applyBinary(LowerFeasible(feasTol_,min_diff_), *lower_);
97  }
98  // Upper feasibility
100  class UpperFeasible : public Elementwise::BinaryFunction<Real> {
101  private:
102  const Real eps_;
103  const Real diff_;
104  public:
105  UpperFeasible(const Real eps, const Real diff)
106  : eps_(eps), diff_(diff) {}
107  Real apply( const Real &xc, const Real &yc ) const {
108  const Real tol = static_cast<Real>(100)*ROL_EPSILON<Real>();
109  const Real one(1);
110  Real val = ((yc <-tol) ? yc*(one+eps_)
111  : ((yc > tol) ? yc*(one-eps_)
112  : yc-eps_));
113  val = std::max(yc-eps_*diff_, val);
114  return xc > val ? val : xc;
115  }
116  };
117  x.applyBinary(UpperFeasible(feasTol_,min_diff_), *upper_);
118  }
119 }
120 
121 template<typename Real>
124  Real one(1), epsn(std::min(scale_*eps,static_cast<Real>(0.1)*min_diff_));
125 
126  mask_->set(*upper_);
127  mask_->axpy(-one,x);
128 
129  Active op(epsn);
130  v.applyBinary(op,*mask_);
131  }
132 }
133 
134 template<typename Real>
135 void Bounds<Real>::pruneUpperActive( Vector<Real> &v, const Vector<Real> &g, const Vector<Real> &x, Real xeps, Real geps ) {
137  Real one(1), epsn(std::min(scale_*xeps,static_cast<Real>(0.1)*min_diff_));
138 
139  mask_->set(*upper_);
140  mask_->axpy(-one,x);
141 
142  UpperBinding op(epsn,geps);
143  mask_->applyBinary(op,g);
144 
145  v.applyBinary(prune_,*mask_);
146  }
147 }
148 
149 template<typename Real>
152  Real one(1), epsn(std::min(scale_*eps,static_cast<Real>(0.1)*min_diff_));
153 
154  mask_->set(x);
155  mask_->axpy(-one,*lower_);
156 
157  Active op(epsn);
158  v.applyBinary(op,*mask_);
159  }
160 }
161 
162 template<typename Real>
163 void Bounds<Real>::pruneLowerActive( Vector<Real> &v, const Vector<Real> &g, const Vector<Real> &x, Real xeps, Real geps ) {
165  Real one(1), epsn(std::min(scale_*xeps,static_cast<Real>(0.1)*min_diff_));
166 
167  mask_->set(x);
168  mask_->axpy(-one,*lower_);
169 
170  LowerBinding op(epsn,geps);
171  mask_->applyBinary(op,g);
172 
173  v.applyBinary(prune_,*mask_);
174  }
175 }
176 
177 template<typename Real>
179  const Real half(0.5);
180  bool flagU = false, flagL = false;
182  mask_->set(v);
183  mask_->applyBinary(isGreater_,*upper_);
184  flagU = mask_->reduce(maximum_) > half ? true : false;
185  }
187  mask_->set(*lower_);
188  mask_->applyBinary(isGreater_,v);
189  flagL = mask_->reduce(maximum_) > half ? true : false;
190  }
191  return ((flagU || flagL) ? false : true);
192 }
193 
194 template<typename Real>
196  // TODO: Cache values?
197 
198  const Real zero(0), one(1);
199 
200  // This implementation handles -l and/or u = infinity properly.
201 
202  /*
203  The first if statement defining d_II (Equation 4.4 of [Ulbrich et al. 1999])
204  is equivalent to x - l <= min(u - x, -g) = - max(g, x - u).
205  * Our x, l, u represent u, a, b in the paper.
206  * Indeed, if x - l = -g < u - x, then min{|g|,c} = min{x - l, u - x, c}.
207  */
208 
209  d.set(x);
210  d.axpy(-one,*upper_);
211  d.applyBinary(Elementwise::Max<Real>(),g);
212  d.plus(x);
213  d.axpy(-one,*lower_); // = x - l + max(g, x - u)
214 
215  mask_->set(x);
216  mask_->axpy(-one,*lower_);
217  mask_->applyBinary(Elementwise::Min<Real>(),g);
218  mask_->plus(x);
219  mask_->axpy(-one,*upper_);
220  mask_->scale(-one); // = u - x - min(g, x - l)
221 
222  mask_->applyBinary(Elementwise::Min<Real>(),d);
223 
224  d.setScalar(-one);
225  d.applyBinary(Active(zero),*mask_);
226  mask_->setScalar(one);
227  d.plus(*mask_);
228  // d[i] = 1 when one of the if conditions in (4.4) are true else 0
229 
230  mask_->set(g);
231  mask_->applyUnary(Elementwise::AbsoluteValue<Real>());
232  d.applyBinary(Elementwise::Multiply<Real>(),*mask_);
233  // d[i] = |g[i]| when one of the if conditions in (4.4) are true else 0
234 
235  // Compute min(x - l, u - x).
236  // * We use the identity min(p, q) = min(p + r, q + r) - r to handle the case
237  // where -l or u = infinity.
238  mask_->set(x);
239  mask_->axpy(-one,*lower_);
240  mask_->plus(x);
241  mask_->applyBinary(Elementwise::Min<Real>(),*upper_);
242  mask_->axpy(-one,x); // = min(x - l, u - x)
243 
244  // When one of the if conditions in (4.4) are true |g|[i] >= (*mask_)[i].
245  // Thus by taking
246  d.applyBinary(Elementwise::Max<Real>(),*mask_);
247  // d_II follows as min(d, c), i.e.,
248  mask_->set(*upper_);
249  mask_->axpy(-one,*lower_);
250  mask_->applyUnary(buildC_);
251  d.applyBinary(Elementwise::Min<Real>(),*mask_);
252 }
253 
254 template<typename Real>
257  dv.applyBinary(Elementwise::DivideAndInvert<Real>(),v);
258 }
259 
260 template<typename Real>
262  const Real one(1), two(2), three(3);
263 
264  // This implementation builds Equation (5.1) of [Ulbrich et al. 1999].
265 
266  Bounds<Real>::buildScalingFunction(dv, x, g); // dv = d_II
267 
268  mask_->set(*upper_);
269  mask_->axpy(-one,*lower_);
270  mask_->applyUnary(buildC_); // = c
271  mask_->axpy(-one,dv); // = c - d_II
272  dv.setScalar(one);
273  dv.applyBinary(Active(0),*mask_); // = \chi
274 
275  mask_->setScalar(three);
276  dv.applyBinary(Elementwise::Multiply<Real>(),*mask_);
277  dv.axpy(-one,*mask_);
278  // dv[i] = 0 if \chi[i] = 1 else -3
279 
280  mask_->set(g);
281  mask_->applyUnary(Elementwise::Sign<Real>());
282  dv.plus(*mask_); // dv[i] = sgn(g[i]) if \chi[i] = 1 else dv[i] <= -2
283 
284  // Set the dv elements that = 0 to sgn(u + l - 2x).
285  mask_->set(*upper_);
286  mask_->plus(*lower_);
287  mask_->axpy(-two,x);
288  mask_->applyUnary(Elementwise::Sign<Real>());
289  dv.applyBinary(setZeroEntry_,*mask_);
290 
291  // Set the dv elements that = 0 to 1.
292  mask_->setScalar(one);
293  dv.applyBinary(setZeroEntry_,*mask_);
294 
295  // Set the dv elements that are <= -2 to 0.
296  mask_->set(dv);
297  dv.applyBinary(Active(-two),*mask_);
298 
299  dv.applyBinary(Elementwise::Multiply<Real>(), g);
300  dv.applyBinary(Elementwise::Multiply<Real>(), v);
301 }
302 
303 } // namespace ROL
304 
305 #endif
Ptr< Vector< Real > > upper_
virtual ROL::Ptr< Vector > clone() const =0
Clone to make a new (uninitialized) vector.
virtual void plus(const Vector &x)=0
Compute , where .
virtual void axpy(const Real alpha, const Vector &x)
Compute where .
Definition: ROL_Vector.hpp:119
void activateLower(void)
Turn on lower bound.
virtual void applyBinary(const Elementwise::BinaryFunction< Real > &f, const Vector &x)
Definition: ROL_Vector.hpp:214
Real ROL_INF(void)
Definition: ROL_Types.hpp:71
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.
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:46
Objective_SerialSimOpt(const Ptr< Obj > &obj, const V &ui) z0_ zero()
void project(Vector< Real > &x) override
Project optimization variables onto the bounds.
Elementwise::ReductionMin< Real > minimum_
Definition: ROL_Bounds.hpp:37
void activateUpper(void)
Turn on upper bound.
Ptr< Vector< Real > > lower_
virtual void setScalar(const Real C)
Set where .
Definition: ROL_Vector.hpp:241
Bounds(const Vector< Real > &x, bool isLower=true, Real scale=1, Real feasTol=std::sqrt(ROL_EPSILON< Real >()))
void buildScalingFunction(Vector< Real > &d, const Vector< Real > &x, const Vector< Real > &g) const
bool isFeasible(const Vector< Real > &v) override
Check if the vector, v, is feasible.
Provides the interface to apply upper and lower bound constraints.
ROL::DiagonalOperator apply
virtual void set(const Vector &x)
Set where .
Definition: ROL_Vector.hpp:175
Real min_diff_
Definition: ROL_Bounds.hpp:35
void applyInverseScalingFunction(Vector< Real > &dv, const Vector< Real > &v, const Vector< Real > &x, const Vector< Real > &g) const override
Apply inverse scaling function.
Ptr< Vector< Real > > mask_
Definition: ROL_Bounds.hpp:33
void projectInterior(Vector< Real > &x) override
Project optimization variables into the interior of the feasible 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.
void applyScalingFunctionJacobian(Vector< Real > &dv, const Vector< Real > &v, const Vector< Real > &x, const Vector< Real > &g) const override
Apply scaling function Jacobian.