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