44 #ifndef ROL_GAUSSIAN_HPP
45 #define ROL_GAUSSIAN_HPP
48 #include "ROL_ParameterList.hpp"
63 Real
erfi(
const Real p)
const {
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);
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);
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));
82 Gaussian(
const Real mean = 0.,
const Real 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;
93 mean_ = parlist.sublist(
"SOL").sublist(
"Distribution").sublist(
"Gaussian").get(
"Mean",0.);
94 variance_ = parlist.sublist(
"SOL").sublist(
"Distribution").sublist(
"Gaussian").get(
"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;
109 const Real half(0.5), one(1), two(2);
114 ROL_TEST_FOR_EXCEPTION(
true, std::invalid_argument,
115 ">>> ERROR (ROL::Gaussian): Gaussian integrateCDF not implemented!");
116 return ((input <
mean_) ? 0.0 : input);
121 const Real
zero(0), half(0.5), one(1), eight(8);
122 const Real dev(std::sqrt(
variance_)), eps(1.24419211485e-15);
126 const Real lVal =
mean_ - eight*dev;
127 const Real uVal =
mean_ + eight*dev;
131 if ( input <= half*eps ) {
return lVal; }
132 if ( input >= one-half*eps ) {
return uVal; }
136 size_t maxit =
static_cast<size_t>(one-std::log2(eps/(eight*dev)));
137 maxit = (maxit < 1 ? 100 : maxit);
139 Real a = (input < half ? lVal :
mean_);
140 Real b = (input < half ?
mean_ : uVal );
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 ) {
150 if ( sc == sa ) { a = c; fa = fc; sa = sc; }
163 case 1: val =
mean_;
break;
165 case 3: val = std::pow(
mean_,3)
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;
187 ROL_TEST_FOR_EXCEPTION(
true, std::invalid_argument,
188 ">>> ERROR (ROL::Distribution): Gaussian moment not implemented for m > 8!");
194 return ROL_NINF<Real>();
198 return ROL_INF<Real>();
201 void test(std::ostream &outStream = std::cout )
const {
203 std::vector<Real> X(size,4.*(Real)rand()/(Real)RAND_MAX - 2.);
204 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