10 #ifndef ROL_GAUSSIAN_HPP
11 #define ROL_GAUSSIAN_HPP
14 #include "ROL_ParameterList.hpp"
29 Real
erfi(
const Real p)
const {
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);
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);
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));
48 Gaussian(
const Real mean = 0.,
const Real 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;
59 mean_ = parlist.sublist(
"SOL").sublist(
"Distribution").sublist(
"Gaussian").get(
"Mean",0.);
60 variance_ = parlist.sublist(
"SOL").sublist(
"Distribution").sublist(
"Gaussian").get(
"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;
75 const Real half(0.5), one(1), two(2);
80 ROL_TEST_FOR_EXCEPTION(
true, std::invalid_argument,
81 ">>> ERROR (ROL::Gaussian): Gaussian integrateCDF not implemented!");
86 const Real
zero(0), half(0.5), one(1), eight(8);
87 const Real dev(std::sqrt(
variance_)), eps(1.24419211485e-15);
91 const Real lVal =
mean_ - eight*dev;
92 const Real uVal =
mean_ + eight*dev;
96 if ( input <= half*eps ) {
return lVal; }
97 if ( input >= one-half*eps ) {
return uVal; }
101 size_t maxit =
static_cast<size_t>(one-std::log2(eps/(eight*dev)));
102 maxit = (maxit < 1 ? 100 : maxit);
104 Real a = (input < half ? lVal :
mean_);
105 Real b = (input < half ?
mean_ : uVal );
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 ) {
115 if ( sc == sa ) { a = c; fa = fc; sa = sc; }
128 case 1: val =
mean_;
break;
130 case 3: val = std::pow(
mean_,3)
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;
152 ROL_TEST_FOR_EXCEPTION(
true, std::invalid_argument,
153 ">>> ERROR (ROL::Distribution): Gaussian moment not implemented for m > 8!");
159 return ROL_NINF<Real>();
163 return ROL_INF<Real>();
166 void test(std::ostream &outStream = std::cout )
const {
168 std::vector<Real> X(size,4.*(Real)rand()/(Real)RAND_MAX - 2.);
169 std::vector<int> T(size,0);
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
Real integrateCDF(const Real input) const
Real evaluateCDF(const Real input) const
virtual void test(std::ostream &outStream=std::cout) const
Real moment(const size_t m) const
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