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
50namespace ROL {
51
52template<class Real>
53class Gaussian : public Distribution<Real> {
54private:
55 Real mean_;
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
80public:
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 }
117
118 Real invertCDF(const Real input) const {
119 //return std::sqrt(2.*variance_)*erfi(2.*input-1.) + mean_;
120 const Real zero(0), half(0.5), one(1), eight(8);
121 const Real dev(std::sqrt(variance_)), eps(1.24419211485e-15);
122 // Set lower and upper bounds to the mean plus/minus 8 standard
123 // -- deviations this ensures that 1-eps probability mass is
124 // -- covered by the interval.
125 const Real lVal = mean_ - eight*dev;
126 const Real uVal = mean_ + eight*dev;
127 // If the input is outside of the interval (half*eps,1-half*eps)
128 // -- then set the return value to be either the lower or
129 // -- upper bound. This case can occur with probability eps.
130 if ( input <= half*eps ) { return lVal; }
131 if ( input >= one-half*eps ) { return uVal; }
132 // Determine maximum number of iterations.
133 // -- maxit is set to the number of iterations required to
134 // -- ensure that |b-a| < eps after maxit iterations.
135 size_t maxit = static_cast<size_t>(one-std::log2(eps/(eight*dev)));
136 maxit = (maxit < 1 ? 100 : maxit);
137 // Run bisection to solve CDF(x) = input.
138 Real a = (input < half ? lVal : mean_);
139 Real b = (input < half ? mean_ : uVal );
140 Real c = half*(a+b);
141 Real fa = evaluateCDF(a) - input;
142 Real fc = evaluateCDF(c) - input;
143 Real sa = ((fa < zero) ? -one : ((fa > zero) ? one : zero));
144 Real sc = ((fc < zero) ? -one : ((fc > zero) ? one : zero));
145 for (size_t i = 0; i < maxit; ++i) {
146 if ( std::abs(fc) < eps || (b-a)*half < eps ) {
147 break;
148 }
149 if ( sc == sa ) { a = c; fa = fc; sa = sc; }
150 else { b = c; }
151 // Compute interval midpoint.
152 c = (a+b)*half;
153 fc = evaluateCDF(c) - input;
154 sc = ((fc < zero) ? -one : ((fc > zero) ? one : zero));
155 }
156 return c;
157 }
158
159 Real moment(const size_t m) const {
160 Real val = 0.;
161 switch(m) {
162 case 1: val = mean_; break;
163 case 2: val = std::pow(mean_,2) + variance_; break;
164 case 3: val = std::pow(mean_,3)
165 + 3.*mean_*variance_; break;
166 case 4: val = std::pow(mean_,4)
167 + 6.*std::pow(mean_,2)*variance_
168 + 3.*std::pow(variance_,2); break;
169 case 5: val = std::pow(mean_,5)
170 + 10.*std::pow(mean_,3)*variance_
171 + 15.*mean_*std::pow(variance_,2); break;
172 case 6: val = std::pow(mean_,6)
173 + 15.*std::pow(mean_,4)*variance_
174 + 45.*std::pow(mean_*variance_,2)
175 + 15.*std::pow(variance_,3); break;
176 case 7: val = std::pow(mean_,7)
177 + 21.*std::pow(mean_,5)*variance_
178 + 105.*std::pow(mean_,3)*std::pow(variance_,2)
179 + 105.*mean_*std::pow(variance_,3); break;
180 case 8: val = std::pow(mean_,8)
181 + 28.*std::pow(mean_,6)*variance_
182 + 210.*std::pow(mean_,4)*std::pow(variance_,2)
183 + 420.*std::pow(mean_,2)*std::pow(variance_,3)
184 + 105.*std::pow(variance_,4); break;
185 default:
186 ROL_TEST_FOR_EXCEPTION( true, std::invalid_argument,
187 ">>> ERROR (ROL::Distribution): Gaussian moment not implemented for m > 8!");
188 }
189 return val;
190 }
191
192 Real lowerBound(void) const {
193 return ROL_NINF<Real>();
194 }
195
196 Real upperBound(void) const {
197 return ROL_INF<Real>();
198 }
199
200 void test(std::ostream &outStream = std::cout ) const {
201 size_t size = 1;
202 std::vector<Real> X(size,4.*(Real)rand()/(Real)RAND_MAX - 2.);
203 std::vector<int> T(size,0);
204 Distribution<Real>::test(X,T,outStream);
205 }
206};
207
208}
209
210#endif
Objective_SerialSimOpt(const Ptr< Obj > &obj, const V &ui) z0_ zero()
virtual void test(std::ostream &outStream=std::cout) const
Real invertCDF(const Real input) const
Real lowerBound(void) const
Real erfi(const Real p) const
Real evaluateCDF(const Real input) const
std::vector< Real > a_
Gaussian(const Real mean=0., const Real variance=1.)
Gaussian(ROL::ParameterList &parlist)
Real integrateCDF(const Real input) const
Real upperBound(void) const
std::vector< Real > c_
std::vector< Real > b_
void test(std::ostream &outStream=std::cout) const
Real evaluatePDF(const Real input) const
std::vector< Real > d_
Real moment(const size_t m) const