ROL
ROL_Gaussian.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_GAUSSIAN_HPP
11 #define ROL_GAUSSIAN_HPP
12 
13 #include "ROL_Distribution.hpp"
14 #include "ROL_ParameterList.hpp"
15 
16 namespace ROL {
17 
18 template<class Real>
19 class Gaussian : public Distribution<Real> {
20 private:
21  Real mean_;
22  Real variance_;
23 
24  std::vector<Real> a_;
25  std::vector<Real> b_;
26  std::vector<Real> c_;
27  std::vector<Real> d_;
28 
29  Real erfi(const Real p) const {
30  const Real zero(0), half(0.5), one(1), two(2), pi(ROL::ScalarTraits<Real>::pi());
31  Real val(0), z(0);
32  if ( std::abs(p) > static_cast<Real>(0.7) ) {
33  Real sgn = (p < zero) ? -one : one;
34  z = std::sqrt(-std::log((one-sgn*p)*half));
35  val = sgn*(((c_[3]*z+c_[2])*z+c_[1])*z + c_[0])/((d_[1]*z+d_[0])*z + one);
36  }
37  else {
38  z = p*p;
39  val = p*(((a_[3]*z+a_[2])*z+a_[1])*z + a_[0])/((((b_[3]*z+b_[2])*z+b_[1])*z+b_[0])*z+one);
40  }
41  val -= (erf(val)-p)/(two/std::sqrt(pi) * std::exp(-val*val));
42  val -= (erf(val)-p)/(two/std::sqrt(pi) * std::exp(-val*val));
43  return val;
44  }
45 
46 public:
47 
48  Gaussian(const Real mean = 0., const Real variance = 1.)
49  : mean_(mean), variance_((variance>0.) ? variance : 1.) {
50  a_.clear(); a_.resize(4,0.); b_.clear(); b_.resize(4,0.);
51  c_.clear(); c_.resize(4,0.); d_.clear(); d_.resize(2,0.);
52  a_[0] = 0.886226899; a_[1] = -1.645349621; a_[2] = 0.914624893; a_[3] = -0.140543331;
53  b_[0] = -2.118377725; b_[1] = 1.442710462; b_[2] = -0.329097515; b_[3] = 0.012229801;
54  c_[0] = -1.970840454; c_[1] = -1.624906493; c_[2] = 3.429567803; c_[3] = 1.641345311;
55  d_[0] = 3.543889200; d_[1] = 1.637067800;
56  }
57 
58  Gaussian(ROL::ParameterList &parlist) {
59  mean_ = parlist.sublist("SOL").sublist("Distribution").sublist("Gaussian").get("Mean",0.);
60  variance_ = parlist.sublist("SOL").sublist("Distribution").sublist("Gaussian").get("Variance",1.);
61  variance_ = (variance_ > 0.) ? variance_ : 1.;
62  a_.clear(); a_.resize(4,0.); b_.clear(); b_.resize(4,0.);
63  c_.clear(); c_.resize(4,0.); d_.clear(); d_.resize(2,0.);
64  a_[0] = 0.886226899; a_[1] = -1.645349621; a_[2] = 0.914624893; a_[3] = -0.140543331;
65  b_[0] = -2.118377725; b_[1] = 1.442710462; b_[2] = -0.329097515; b_[3] = 0.012229801;
66  c_[0] = -1.970840454; c_[1] = -1.624906493; c_[2] = 3.429567803; c_[3] = 1.641345311;
67  d_[0] = 3.543889200; d_[1] = 1.637067800;
68  }
69 
70  Real evaluatePDF(const Real input) const {
71  return std::exp(-std::pow(input-mean_,2)/(2.*variance_))/(std::sqrt(2.*ROL::ScalarTraits<Real>::pi()*variance_));
72  }
73 
74  Real evaluateCDF(const Real input) const {
75  const Real half(0.5), one(1), two(2);
76  return half*(one+erf((input-mean_)/std::sqrt(two*variance_)));
77  }
78 
79  Real integrateCDF(const Real input) const {
80  ROL_TEST_FOR_EXCEPTION( true, std::invalid_argument,
81  ">>> ERROR (ROL::Gaussian): Gaussian integrateCDF not implemented!");
82  }
83 
84  Real invertCDF(const Real input) const {
85  //return std::sqrt(2.*variance_)*erfi(2.*input-1.) + mean_;
86  const Real zero(0), half(0.5), one(1), eight(8);
87  const Real dev(std::sqrt(variance_)), eps(1.24419211485e-15);
88  // Set lower and upper bounds to the mean plus/minus 8 standard
89  // -- deviations this ensures that 1-eps probability mass is
90  // -- covered by the interval.
91  const Real lVal = mean_ - eight*dev;
92  const Real uVal = mean_ + eight*dev;
93  // If the input is outside of the interval (half*eps,1-half*eps)
94  // -- then set the return value to be either the lower or
95  // -- upper bound. This case can occur with probability eps.
96  if ( input <= half*eps ) { return lVal; }
97  if ( input >= one-half*eps ) { return uVal; }
98  // Determine maximum number of iterations.
99  // -- maxit is set to the number of iterations required to
100  // -- ensure that |b-a| < eps after maxit iterations.
101  size_t maxit = static_cast<size_t>(one-std::log2(eps/(eight*dev)));
102  maxit = (maxit < 1 ? 100 : maxit);
103  // Run bisection to solve CDF(x) = input.
104  Real a = (input < half ? lVal : mean_);
105  Real b = (input < half ? mean_ : uVal );
106  Real c = half*(a+b);
107  Real fa = evaluateCDF(a) - input;
108  Real fc = evaluateCDF(c) - input;
109  Real sa = ((fa < zero) ? -one : ((fa > zero) ? one : zero));
110  Real sc = ((fc < zero) ? -one : ((fc > zero) ? one : zero));
111  for (size_t i = 0; i < maxit; ++i) {
112  if ( std::abs(fc) < eps || (b-a)*half < eps ) {
113  break;
114  }
115  if ( sc == sa ) { a = c; fa = fc; sa = sc; }
116  else { b = c; }
117  // Compute interval midpoint.
118  c = (a+b)*half;
119  fc = evaluateCDF(c) - input;
120  sc = ((fc < zero) ? -one : ((fc > zero) ? one : zero));
121  }
122  return c;
123  }
124 
125  Real moment(const size_t m) const {
126  Real val = 0.;
127  switch(m) {
128  case 1: val = mean_; break;
129  case 2: val = std::pow(mean_,2) + variance_; break;
130  case 3: val = std::pow(mean_,3)
131  + 3.*mean_*variance_; break;
132  case 4: val = std::pow(mean_,4)
133  + 6.*std::pow(mean_,2)*variance_
134  + 3.*std::pow(variance_,2); break;
135  case 5: val = std::pow(mean_,5)
136  + 10.*std::pow(mean_,3)*variance_
137  + 15.*mean_*std::pow(variance_,2); break;
138  case 6: val = std::pow(mean_,6)
139  + 15.*std::pow(mean_,4)*variance_
140  + 45.*std::pow(mean_*variance_,2)
141  + 15.*std::pow(variance_,3); break;
142  case 7: val = std::pow(mean_,7)
143  + 21.*std::pow(mean_,5)*variance_
144  + 105.*std::pow(mean_,3)*std::pow(variance_,2)
145  + 105.*mean_*std::pow(variance_,3); break;
146  case 8: val = std::pow(mean_,8)
147  + 28.*std::pow(mean_,6)*variance_
148  + 210.*std::pow(mean_,4)*std::pow(variance_,2)
149  + 420.*std::pow(mean_,2)*std::pow(variance_,3)
150  + 105.*std::pow(variance_,4); break;
151  default:
152  ROL_TEST_FOR_EXCEPTION( true, std::invalid_argument,
153  ">>> ERROR (ROL::Distribution): Gaussian moment not implemented for m > 8!");
154  }
155  return val;
156  }
157 
158  Real lowerBound(void) const {
159  return ROL_NINF<Real>();
160  }
161 
162  Real upperBound(void) const {
163  return ROL_INF<Real>();
164  }
165 
166  void test(std::ostream &outStream = std::cout ) const {
167  size_t size = 1;
168  std::vector<Real> X(size,4.*(Real)rand()/(Real)RAND_MAX - 2.);
169  std::vector<int> T(size,0);
170  Distribution<Real>::test(X,T,outStream);
171  }
172 };
173 
174 }
175 
176 #endif
std::vector< Real > b_
Objective_SerialSimOpt(const Ptr< Obj > &obj, const V &ui) z0_ zero()
Gaussian(ROL::ParameterList &parlist)
Gaussian(const Real mean=0., const Real variance=1.)
Real lowerBound(void) const
std::vector< Real > a_
Real integrateCDF(const Real input) const
Real evaluateCDF(const Real input) const
virtual void test(std::ostream &outStream=std::cout) const
std::vector< Real > d_
Real moment(const size_t m) const
std::vector< Real > c_
Real invertCDF(const Real input) const
Real upperBound(void) const
Real erfi(const Real p) const
Real evaluatePDF(const Real input) const
void test(std::ostream &outStream=std::cout) const